ROBUST RESECTION MODEL FOR ALIGNING THE MOBILE MAPPING SYSTEMS TRAJECTORIES AT DEGRADED AND DENIED URBAN ENVIRONMENTS

: Mobile mapping systems MMS equipped with cameras and laser scanners are widely used nowadays for different geospatial applications with centimetric accuracy either in project wise or national wise scales. The achieved positioning accuracy is very much related to the navigation unit, namely the GNSS and IMU onboard. Accordingly, in GNSS denied and degraded environments, the absolute positioning accuracy is worsened to few meters in some cases. Frequently, ground control points GCPs of a high positioning accuracy are used to align the MMS trajectories and to improve the accuracy when needed. The best way to integrate the MMS trajectories to the GCPs is by measuring them on the MMS images where the positioning accuracy is dropped. MMS images are mostly spherical panoramic (equirectangular) images and sometimes perspective and, in both types, it is required to precisely determine the images orientation in what is called as space resection or camera pose determination. For perspective images, the pose is conventionally determined by collinearity equations or by using projection and fundamental matrices. Whereas for equirectangular panoramic images it is based on resecting vertical and horizontal angles. However, there is still a challenge in the state–of–the–art of image pose determination because of the model nonlinearity and the sensitivity to proper initialization and spatial distribution of the points. In this research, a generic method is presented to solve the pose resection problem for the perspective and equirectangular images using oblique angles. The oblique angles are derived from the measured image coordinates and based on spherical trigonometry rules and vector geometry. The developed algorithm has proven to be highly stable and steadily converge to the global minimum. This is related to the robust geometric constraint offered by the oblique angles that are enclosed between the object points and the camera. As a result, the MMS trajectories are realigned accurately to the GCPs and the absolute accuracy is highly refined. Four experimental tests are presented where the results show the efficiency of the proposed angular based model in different cases of simulated and real data with different image types.


INTRODUCTION
Mobile Mapping Systems MMS are widely used nowadays for different applications in urban and rural environments like for highways and roads engineering projects, façades dimensional surveying, traffic signs detection, as built surveying applications, infrastructure monitoring and many other applications (Karel Sukup and Jan Sukup, 2010). Mostly, those mentioned applications require a high level of positioning accuracy which might be a challenge in urban canyons, tunnelling roads and dense forest areas of denied GNSS environments. Accordingly, ground control points GCPs either surveyed in the field or acquired from aerial images are used to realign the MMS trajectories to centimetre or decimetre accuracy levels. Normally the GCPs are measured on the MMS images and a space resection model is then applied. The output of space resection or camera pose computations are the exterior orientation parameters namely: the camera position and its rotation angles at the instant of capture. It should be noted that the MMS can be equipped either with a frame camera or with a 360 panoramic camera like the Ladybug camera which is widely used in the MMS nowadays (FLIR 2019) Hence, in the case of spherical panoramic imaging, an angular based model is used for solving the resection problem. In principle, a spherical projection is used to project the whole sphere on a flat surface and result in an equirectangular (equidistant projection) image, with a 2:1 ratio where the width is exactly twice the height. This is logical, since it covers 360 horizontally and 180° vertically and allowing the viewer to look in every direction, including the zenith and nadir (Figure 1). Therefore, every point in the equirectangular panorama has a latitude (vertical angle) and longitude (horizontal angle) coordinates. The horizontal and vertical angles are used to solve the pose of the equirectangular panorama by space resection. However, the resection model is nonlinear and sensitive to initial values that should not be far by more than 1' to have a proper convergent solution (Shepherd 1982, Schofield 2007, Alsadik 2016). On the other hand, when a frame camera is used, then the collinearity resection model is used. Collinearity is considered as the standard mathematical model in photogrammetry and machine vision to compute the images orientation within the bundle adjustment (Lourakis andArgyros 2004, McGlone, Mikhail et al. 2004). The collinearity concept is based on assuming the object point, its image coordinates, and the camera lens are collinear, and all located on a straight line ( Figure 2).

Figure 2. Collinearity equations principle.
A minimum of three reference points or GCPs are required to solve the resection problem (Luhman, Robson et al. 2014). Basically, collinearity equations are nonlinear and normally solved by starting from a proper set of initial pose values which are crucial in the solution. Different methods are introduced in the field of photogrammetry and machine vision to solve the nonlinearity in the collinearity model and to efficiently solve with a minimum number of control points. From a machine vision perspective, the camera pose is commonly solved as the perspective n-point (PnP) problem and this direct solution approach is mainly developed with three points (P3P) or four points (P4P) methods (Horaud, Conio et al. 1989, Grafarend and Shan 1997, Quan and Lan 1999, Gao and Chen 2001, Kneip, Scaramuzza et al. 2011. However, the PnP direct methods result in more than one solution and the correct solution should be determined. Furthermore, these methods are limited to a specific number of reference points (n=3, 4, or 5) and do not consider the redundancy in observations which is supposed to strengthen the solution from a statistical viewpoint (Alsadik 2016). Another common approach is the well-known direct linear transformation (DLT) (Marzan and Karara 1975) and it is based on the projective relations between the image and the object space by estimating the so called projection matrix. A minimum of five reference points is necessary to solve the system of DLT equations in a direct linear solution (Luhman, Robson et al. 2014). In summary, the mentioned camera orientation approaches are either nonlinear solutions and dependent on the proper initial values or linear, non-rigorous solutions which cannot handle statistically the redundant observations. Moreover, these methods are sensitive to the poor spatial distribution of the points and when the interior orientation of the camera is also required. Specifically, bad distribution is present when the points are placed on or close to a straight line. Similarly, DLT can result in a weak conditioned system of equations if all reference points are positioned on a common plane (Luhman, Robson et al. 2014). Accordingly, the aim of this research is to develop a generic mathematical model for perspective and equirectangular panoramic image orientation determination that: 1) consider the redundancy in observations and 2) robust to improper initial values. This means to introduce a solution with a geometric stability, reliability and converge to global minimum even with improper initial values. The developed model is based on using angular conditions represented by oblique angles instead of using horizontal and vertical angles. Spherical trigonometry laws will be used to derive the oblique angles while vectors geometry will be applied to develop the mathematical model that relates: the observed (derived) angles, the camera unknown orientation parameters, and the coordinates of the object points. The MMS absolute accuracy improvement, through the realignment of the trajectories, will highly benefit from such proposed resection model.
Following in section 2, the method of deriving the oblique angles and their use in the mathematical model will be presented. In section 3, different experimental tests will be presented for panoramic images and perspective images. Finally, in section 4, the discussion of results and conclusions will be presented.

METHODOLOGY
The proposed method of the image resection determination is based on deriving oblique angles from the observed image coordinates in the viewing cameras. A spherical trigonometry law is used to apply the oblique angle derivation as will be discussed in the following section 2.1. Afterward, the development of the proposed mathematical model of the image pose determination will be shown in section 2.2.

Oblique angle derivation
A spherical triangle ABC with two vertical angles ( ), one horizontal angle ( ), and one oblique angle is shown in Figure 3. The cosine rule in equation 1 can be used to solve the spherical triangle ABC and to compute the oblique angle (Murray 1908).
The described solution of the spherical triangle and the derived oblique angle can be used to solve the problem of pose determination and the triangulation problem as will be shown in section 2.2 and section 2.3 respectively. A related point to consider is the minimum required number of GCPs. Three points are necessary to define the oblique angle in the camera pose problem. Accordingly, if every oblique angle formulates one observation equation then a minimum of three oblique angles is necessary to define the object space coordinates of the camera and the same applies for the image triangulation problem (Alsadik 2016).

Mathematical model
The computations start by deriving the observed oblique angles enclosed between the camera and two control points and . For perspective images, this is possible by using the observed image coordinates referenced to the principle point of both control points as follows: where observed image coordinates of the points and respectively. camera focal length.
While for the equirectangular panoramic images, one horizontal angle and two vertical angles are easily computed by applying equations 7 and 8 assuming the images are aligned on the horizontal plane at the panorama centre. (7) , , where = the observed columns ( -coordinates) in pixels for the control points and respectively. , the observed rows (y-coordinates) in pixels for the control points and .
the pixel size in degrees.
Then, for both types of images, the oblique angle can be computed using the spherical trigonometry rules as mentioned in Equation 1. It should be noted that a high redundancy of the observed angles can be attained based on the number of control points . As shown in Figure 4a, three GCPs produce three observed oblique angles while the addition of a fourth GCP will result in six observed oblique angles as shown in Figure 4b and so on. Therefore, the possible number of oblique angles can be calculated based on the number of measured GCPs using the following equation 9: a) b) Figure 4. The composed oblique angles related to three and four observed control points. a) three observed oblique angles. b) six observed oblique angles.
The high redundancy of the observed oblique angles is expected to strengthen the stability of the solution and the convergence to an optimal minimum as will be shown in the experimental tests. So far, the oblique angles are derived from the measured imag coordinates. The same oblique angles can also be computed using the approximate camera 3D coordinates applying the vectors dot product as shown in equation 10.
where the difference in coordinates between the viewing camera and points . = spatial distances between the camera and the observed control points respectively.
Equation 10 can be reformulated as shown in equation 11 to further compute the partial derivatives ( ) to the unknown coordinates of the camera station ( ) as follows: ( 1 1 )   (12) The partial derivatives of equation 12 are necessary to solve the nonlinearity and redundancy of the observation equations system. The solution will be applied by least squares adjustment using either a Gauss-Newton technique or Levenberg Marquardt (Marquardt 1963, Ghilani and where = vector of corrections to the camera coordinates . matrix of partial derivatives to the camera coordinates , .
cofactor weight matrix of the observations. covariance matrix of the camera coordinates.

Proposed method workflow
The workflow diagram of Figure 5 is shown to summarize the implementation of the developed algorithm and a MATLAB code of the proposed resection model is listed in the appendix as well.
The algorithm is applicable for both: perspective and panoramic images where the GCPs are measured on the images in pixels. The image coordinates are transformed into horizontal and vertical angles which are finally converted into derived oblique angles using spherical trigonometry.
Since the sky is open with no surrounding objects exist nearby aeroplanes, aerial images absolute accuracy is not degraded by GNSS multipath and line of sight errors. Accordingly, well identified objects on both: aerial images and the MMS images can serve for GCPs of centimetric accuracy. Another approach is to have GCPs by the traditional field surveying which is of a higher accuracy, but it requires labour force, time and cost. Figure 5. Oblique angular-based resection workflow.

Pose estimation of the panoramic images
In this section, we introduce two experimental tests to use the proposed method of panoramic image pose determination by oblique angles. The first test is applied on simulated data while the second test is applied on real panoramic image of a mobile mapping system. As mentioned, the images are assumed aligned on the horizontal plane at the panorama centre.

First test -simulated data
To test the developed method for panoramic images, a simulated reality scaled 3D environment is created using Blender tool (Community 2018) where the exact (true) locations are known. The true camera position P is shown as a red arrow in Figure 6a while a panoramic equirectangular image is created with dimension of 4800 2400 pixels as shown in Figure 6b. Reference points A, B, C, and D on well identified corners are selected and measured on the panoramic image to compute its pose by the developed method of oblique angles (Table 1). The camera P coordinates are simulated to be at m m, and m. The efficiency of the developed method is validated through checking whether the computed camera P position is the same as designed in the simulated environment or not.   The computed adjusted coordinates are: m m m. The slight differences with respect to the true values are related to the imperfection of the manual measurements. However, in practice, automated point measurements can be applied.

Second test -real MMS trajectory data
The second test is applied on an MMS data consist of 129 recordings at a degraded GNSS urban environment shown in Figure 8a. The MMS is equipped with a panoramic camera producing equirectangular images of a 4800 2400 pixels resolution (CycloMedia 2019). The suggested resection model is applied to determine the position of every image recording along the trajectory by the realignment to five predefined GCPs measured from Leica CityMapper camera with an accuracy of 10 cm on average (Figure 8b). For illustration, one detailed image resection is given in the following description. Given: five control points are measured from aerial oblique images ( Figure 8) and then their measured pixel coordinates on the MMS panoramic image are measured as shown in Table 3. The horizontal, vertical and then the oblique angles are computed between the panoramic image P and the control points based on the observed pixel coordinates as illustrated in Table 4  To test the robustness of the suggested algorithm, challenging starting values of (0,0,0) for the camera coordinates are selected. The solution converged to the optimal values after 22 iterations as shown in Figure 9 to 92255.78 m, 437597.07 m, and 2.65 m.

Pose estimation of the perspective images
Two tests are applied to verify the proposed camera pose estimation using oblique angles-based method. The first test is applied in a simulated scene where the true values are known. Since real MMS perspective image data was not available, the second test is applied to an aerial image as published in (Easa 2010).

First test -simulated data
To evaluate the computed image orientation, a ground truth data is created in a simulation environment. The simulated image is rendered where four reference targets are mounted on a building façade (Figure 10a). The coordinates of the targets are fixed as listed in Table 5 and the rendered image is created assuming a Canon 500D camera with a focal length of 18mm (Figure 10b). The camera is assumed free of the lens distortion.
ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume V-1-2020, 2020 XXIV ISPRS Congress (2020 edition) a) b) Figure 10. a) Simulated model with four control targets. b) Synthetic Image captured by Canon 500D camera.
The camera pose computations are applied using the proposed oblique angular based method. The approximate values of the camera space coordinate ( ) and the rotation angles ( , , ) are selected by starting from far improper values to verify the method robustness. The GCPs XYZ coordinates and their manually measured image coordinates are given as follows in Table 5 -10 95 7 -10 -10 90 Table 6. Initial camera orientation parameters.
The computation of the P camera coordinates is applied using the proposed angular based method. The computed horizontal angles , vertical angles ( , and the derived oblique angles are shown in Table 7 using the observed image coordinates of the GCPs given in Table 5.  19798 APC 37.50925 29.50303 9.89955 -13.54558 APD 24.48941 5.76123 9.89955 -13.91349 BPC 13.16006 4.59969 -1.19798 -13.54558 BPD 32.88362 30.66457 -1.19798 -13.91349 CPD 34.22699 35.26427 -13.54558 -13.91349 Table 7. Derived angles of the simulated test.
Then, the least squares adjustment is iterated until the corrections are converged to zero. Figure 11 illustrates the convergence of the solution to optimal minimum. Further, collinearity equations are reduced to only three unknown rotation angles while fixing the determined camera coordinates.

Second test -Aerial image data
The second test is applied to an aerial perspective image which is listed in (Easa 2010). Four GCPs and their measured image positions are given in Firstly, the camera position is computed using the suggested method with least squares adjustment and the solution is rapidly converged to optimal minimum in three iterations as shown in Figure 12. Then the rotation angles are computed using collinearity condition where the camera position is fixed. It's worth to mention that even when selecting challenging approximate values for the camera exterior parameters of zero values as: , we still get a robust solution which converges to optimal minimum after 12 iterations as illustrated in Figure 13. Then, a reduced collinearity model can be used for computing the angular rotations.

DISCUSSION AND CONCLUSIONS
In this paper, a developed method for the equirectangular and perspective image pose determination is presented. This is aimed ,among many potential applications, for improving the mobile mapping trajectories using predefined control points. The suggested camera pose estimation is based on using derived oblique angles and showed promising results. In the four presented tests, the method was robust to converge to the optimal solution despite the far initial values given to the nonlinear model. For demonstration, Figure 14 illustrate the black dots which represent the initial values given to the panoramic image position within an area extent of 20 square kilometres and where the optimal value is located in the middle (red dot). Using the proposed algorithm and starting from any initial point coordinates represented by the black dots will converge to the optimal coordinates as shown by the blue arrows. On the other hand, if we applied the conventional solution using horizontal and vertical angles instead of the oblique angles, then the solution couldn't converge far by 10 meters away from the optimal coordinates. Figure 14. Robust convergence of the proposed method toward optimal value at the centre.
For perspective images, the proposed method showed a successful integration with the collinearity equations to solve the six pose parameters as shown in section 3.2 of the simulated and real experimental tests. Future work can investigate the efficiency of the proposed mathematical model to include the modelling of the interior camera calibration parameters. A comparison of performance with the other methods of pose determination for perspective images can also be investigated.

ACKNOWLEDGMENT
The author would like to introduce his gratitude to CycloMedia company for supporting with some test data.