A MODIFIED METHOD FOR IMAGE TRIANGULATION USING INCLINED ANGLES

The ongoing technical improvements in photogrammetry, Geomatics, computer vision (CV), and robotics offer new possibilities for many applications requiring efficient acquisition of three-dimensional data. Image orientation is one of these important techniques in many applications like mapping, precise measurements, 3D modeling and navigation. Image orientation comprises three main techniques of resection, intersection (triangulation) and relative orientation, which are conventionally solved by collinearity equations or by using projection and fundamental matrices. However, different problems still exist in the state – of –the –art of image orientation because of the nonlinearity and the sensitivity to proper initialization and spatial distribution of the points. In this research, a modified method is presented to solve the triangulation problem using inclined angles derived from the measured image coordinates and based on spherical trigonometry rules and vector geometry. The developed procedure shows promising results compared to collinearity approach and to converge to the global minimum even when starting from far approximations. This is based on the strong geometric constraint offered by the inclined angles that are enclosed between the object points and the camera stations. Numerical evaluations with perspective and panoramic images are presented and compared with the conventional solution of collinearity equations. The results show the efficiency of the developed model and the convergence of the solution to global minimum even with improper starting values.


INTRODUCTION
Collinearity equations are considered as the standard model used in photogrammetry and computer vision to compute the image orientation within bundle adjustment (Lourakis and Argyros 2004;McGlone et al. 2004).The concept is based on the central projection with an ideal situation that the object point, its image coordinates, and the camera perspective center are collinear and lying on a straight line.The exterior and interior orientation parameters of the camera are efficiently represented within these collinearity equations.The calculation of the exterior orientation parameters is based on observing a minimum of three reference points by resection.Whereas intersection or triangulation is to determine the space coordinates of target points by knowing the exterior orientation of at least two viewing cameras (Luhman et al. 2014).Basically, collinearity equations are nonlinear and normally solved by starting from a proper set of approximate values.
Many researchers in the field of photogrammetry and computer vision have tried to avoid the indirect nonlinear solution of collinearity and to solve with a minimum number of reference points.They introduced during the last three decades different direct closed form solutions for resection, triangulation and relative orientation problems.From a computer vision perspective, the space resection problem is solved as the perspective n-point (PnP) problem and this approach is mainly developed with three points P3P or four points P4P methods (Gao and Chen 2001;Grafarend and Shan 1997;Horaud et al. 1989;Quan and Lan 1999).However, these methods have resulted in more than one solution and a decision should be made to find the unique solution.Moreover, these methods are not considering the redundancy in observations that is supposed to strengthen the solution from a statistical viewpoint.
Another approach is called direct linear transformation DLT (Marzan and Karara 1975) and is based on the projective relations between the image and the object space by estimating the so called projection matrix P. A minimum of five reference points is necessary to solve the system of DLT equations in a direct linear solution (Luhman et al. 2014).For image triangulation which is the subject of this paper, Hartley and Sturm (1997) introduced an extensive study of the available triangulation methods and evaluated the best methods by studying their stability against observation noise.For perspective image triangulation problem, a rigorous direct linear solution can be derived with collinearity (Alsadik 2013).
Generally, the challenge of image orientation with the nonlinear collinearity methods is to initialize the solution with the proper approximate values of orientation parameters.Starting from improper approximate values can mostly lead to a wrong divergent solution.On the other hand, collinearity is not always applicable in its standard form like in the case of panoramic imaging.
The aim of this research is to develop a mathematical model of image triangulation that is appropriate to handle for poor approximations and solution instabilities.Especially when collinearity cannot be applied as in the case of spherical or cylindrical equirectangular panoramic imagery.This means to try introducing a solution with a geometric stability, reliability and converge to global minimum even with improper initial values or when the points are badly distributed.The developed model is based on using angular conditions represented by inclined angles instead of collinearity conditions.Based on land surveying principles, intersecting horizontal and vertical angles are sensitive to improper approximate values.A shifted approximate point coordinates of more than 1 ′ in directions is probably failing to converge (Shepherd 1982).Therefore, inclined angles are used instead of horizontal and vertical angles for image orientation because they offer a higher geometric constraint on the solution and converge steadily to the global minimum.
The following section 2 introduces the method of deriving the inclined angles and how to use the observed image coordinates in the derivation and image triangulation.Section 3 illustrates the experimental results and section 4 presents the discussion of results and conclusions.

METHODOLOGY
The proposed method of image triangulation is based on deriving inclined angles from the observed image coordinates in the viewing cameras.A spherical trigonometry law is used to apply the derivation as discussed in the following section 2.1.Afterward, the derivation of the proposed mathematical model of image triangulation and the least square adjustment will be shown in section 2.2.

Inclined angle derivation
Given a spherical triangle ABC (Fig. 1) with two known vertical angles (β 1 , β 1 ) and one known horizontal angle (θ ), the inclined angle γ may be computed from the cosine rule (Murray 1908).The cosine rule in equation 1 can be used to solve the spherical triangle ABC and to compute the inclined angle γ. cos( γ) = cos(θ) cos(β 1 ) cos(β 2 ) + sin ( β 1 )sin (β 2 ) (1) The described solution of the spherical triangle of inclined angle is useful to solve the problem of image triangulation as will be shown in section 2.2.From the aforementioned discussion of the spherical triangle, three points are necessary to define the inclined angle in the orientation problem, the viewing camera center C1, the object point A and the adjacent camera C2 as shown in Fig. 2. Accordingly, if every inclined angle formulate an observation equation then a minimum of three inclined angles is necessary to define the object space coordinates XYZ of a point.
As stated, inclined angles are defined by knowing the horizontal and vertical angles to the target point.These angles are indirectly observed by knowing the image coordinates.The horizontal angle θ is defined as the difference between direction angles α 1 and α 2 .In image space, these angles (α, β) are computed by using the image coordinates x, y of the target point and the camera focal length f as shown in equations 2, 3, 4 and 5. Fig. 3a illustrates the derivation of these angles for a single camera based on the measured image coordinates with respect to the image center (principal point p.p) system.However, cameras in reality are oriented in an angular rotation, which is not parallel to the space coordinate system (XYZ) as shown in Fig. 3b.The angular rotations can be represented as Euler's angles ω, φ, k around the X-Y-Z axis respectively (Wolf and DeWitt 2000).
The exterior orientation angles ω, φ, k should be considered in the computations of the correct angles (α, β) to finally estimate the inclined angles.It is worth to mention that the derivation is based on a right-handed system.(2) where M ω , M φ , M k are the rotation matrices with respect to orientation angles ω, φ, k respectively and xa, ya refer to the image coordinates and f is the focal length.
The second pair of angular directions (α 2 , β 2 ) is computed between the viewing camera C1 and the adjacent camera C2.The computation is simply as follows: It must be noted that the aforementioned derivation is only applicable to the perspective images.On the other hand, the angles are easily computed in panoramic images as will be shown in section three.Finally, the inclined angle γ can be derived as shown in equation 1 after computing the horizontal angle θ which is the absolute difference between direction angles α 1 and α 2 as mentioned earlier.

The mathematical model of image triangulation with inclined angles
To apply the image orientation with inclined angles, a mathematical model should be developed to relate the inclined angles with camera orientations and the object points.This is done by using two vectors dot product as shown in Fig. 4 and equation 11.
Where This high redundancy in the observed angles is expected to strengthen the stability of the solution and convergence to an optimal minimum as will be shown in the numerical examples.
To solve the non-linearity and redundancy of the observation equations 12, least square adjustment using a Gauss-Markov model is applied (Ghilani and Wolf 2006) as shown in equation 14.This nonlinear system of equations is normally linearized by using Taylor's expansion and neglecting higher order terms.
where ∆: Vector of corrections δX P , δY P , δZ P .B: Matrix of partial derivatives to unknowns X P , Y P , Z P .F: Vector of function value of discrepancy between the approximate and correct values of unknowns.Q e : Cofactor matrix of observations.N: Normal equation matrix.
With respect to equation 14 and Fig. 5, the coefficient matrix B of the partial derivatives to the unknown parameters of image orientation will be formed as: The partial derivatives can be formulated as:

EXPERIMENTAL TESTS
To evaluate the developed method, three numerical tests were implemented.The first test was applied to three intersected street view panoramas to a manhole control point.For perspective images, two tests are applied.A simulated experiment of five intersecting cameras and a real test of the intersection of twelve cameras triangulated at a control point.The first two tests were applied with a focal length camera of 18mm and image format 22.3×14.8mm 2 .
Panoramic image intersection: Three intersected spherical panoramas in a street mapping system (CycloMedia.B.V. Technology) are shown in Fig. 6 where the vehicle equipped with a GNSS system and IMU.A ground control point GCP is measured on each panorama (red dot) and listed in Table 1 to determine finally its position by intersection.The image coordinates are transformed into the image center before computing the angular observations as: x i = c i − 2400 − 0.5 y i = 1200 − r i − 0.5 Obviously, angular value per pixel is 0.075 degrees and these results in computing the horizontal and vertical HV angles as follows ( As mentioned before, the conventional intersection of horizontal and vertical angles (HV method) is not successful and will diverge with this improper starting value.Therefore, non-linear least square adjustment of inclined angles is applied to finally result in the adjusted position.The solution iterations are shown in Fig. 7. Table 3. Simulated data of a triangulation problem.The maximum number of 20 angles is observed as shown in equation 13.An approximate starting value of the nonlinear least square adjustment was selected to be far from the true coordinates [1000; 1500; 500] to investigate the stability of the solution based on the measured image coordinates of Table 3.To validate the efficiency of the developed method, the solution was compared with the nonlinear collinearity image triangulation.A normally distributed random noise between 1 and 10 pixels were generated and added to the image coordinates.This is to investigate the deviation magnitude away from the true value of the object point coordinates as shown in Fig. 10.The shift in the object point coordinates was computed with respect to a noise increment of one pixel.The graph shows a close result between the proposed method of inclined angles and the collinearity equation method.
Figure 10.Comparison between the proposed method and collinearity by estimating shift in the object space (y-axis) with respect to the added normally distributed image noise (x-axis).

Perspective images -Second test:
The second test is applied in a real imaging case of a target point marked on a building façade and measured with a total station for checking and found to be (975.524, 20044.271, 302.718m) in , , and  respectively.Twelve images are triangulated to the object point as shown in Fig. 11 with a Canon 500D camera.The interior orientation parameters of the camera are: f=18.1mm,frame size= 4752×3168 pixels, pixel size=4.7 microns while the lens distortion parameters are neglected.As shown in equation 13, a maximum of 132 inclined angles is derived from the twelve images and processed as the observed values to the target object point.
Figure 11.Triangulation -second test  The initial values to run the HV method is based on the output of the inclined angle method to ensure the solution convergence.
Figure 12.The convergence of corrections to zero values of the second test.
Similarly to the first test, a validation of the developed method is applied by adding noise (1-10 pixels) to the image coordinates.The solution was compared with the collinearity image triangulation.The shift in the object point coordinates was computed with respect to a noise increment of one pixel.
The graph (Fig. 13) shows a close result between the proposed method of inclined angles and the collinearity equation method.
The graph shows a maximum shift of 6 cm resulted in the point position when adding a noise of 10 pixels.
Figure 13.Comparison between the shifts resulted after adding a random noise to the image coordinates of the second test.

DISCUSSION AND CONCLUSIONS
In this paper, we presented a modified algorithm for image triangulation using inclined angles.These inclined angles were derived from the measured image coordinates using spherical trigonometry rules.Dot product of two vectors was used to model the inclined angles enclosed between the points and their viewing cameras as illustrated in equation 12.The derived inclined angels offered a higher redundancy than collinearity equations with the same number of points and cameras.Therefore, least square adjustment (equation 14) was used in the computations of triangulation.The partial derivatives of the object points were listed in section 2.2.The minimum number of three cameras was needed to solve the triangulation.In addition to orient frame-perspective images, the method was promising in the application of panoramic image intersection of Fig. 6.
The inclined angles were efficient and stable even when noise added to the observations (Fig. 10).However, the linear solution of triangulation with collinearity is still preferable since it is a direct method without the need for starting values and can be solved with two stereo images.However, the proposed method can also be used to solve panoramic image intersection and geodetic intersection problems.
The developed method of triangulation as shown in the results is stable and converges efficiently to their true values.In practice, starting from close approximations is more efficient to solve the image triangulation problem.However, similarly to collinearity, the developed method was also in the situations where the cameras or reference points were located on a straight line as mentioned previously.Either this may result in zero valued inclined angles or rank defected system of equations and the method fails.Many issues can be investigated for future work like the extension of this approach to solve the resection problem and relative orientation.Further, to combine collinearity with inclined angles to make the solution more robust or to develop direct linear solutions.Camera calibration is also to be considered with the developed model.The use of other forms of rotation instead of Euler angles is to be explored like the use of quaternion for more computational simplicity.

Figure 4 .
Figure 4. Vectors dot product.v1 ⃑⃑⃑⃑ .v2 ⃑⃑⃑⃑ = |v1||v2| cos γ (11) where |v1|and |v2| represent the vectors length and γ is the inclined angle enclosed between them.On this basis, and with reference to Fig.5, the mathematical model can developed as shown in equation 12 for inclined angle γ enclosed by two vectors connecting the viewing camera C1 with camera C2 and object point P.

Figure 5 .
Figure 5. Inclined angle derivation between two cameras and the object point.Accordingly, every derived inclined angle is formulated in one observation equation.Hence, three intersected images at least are needed in image triangulation to compute the object point coordinates X P , Y P , Z P and three reference points in resection for b X P = ∆X − q (X P o − Xc) b Y P = ∆Y − q (Y P o − Yc) (16) b Z P = ∆Z − q (Z P o − Zc) Where P o = approximate coordinates of point P. q = cos γ * l /l C 1 P o

Figure 6 .
Figure 6.Three panoramic images triangulation.Orientation angles ω and φ are set to zero in all images, while the format size of each image is 4800×2400 pixels.A starting value of GCP point coordinates is: [1 1 1].The image coordinates are transformed into the image center before computing the angular observations as:

Figure 7 .
Figure 7. Corrections to GCP point by inclined angles triangulation from the second iteration.

Figure 8 .
Figure 8. Triangulation test of five simulated cameras.

Fig. 9
shows 20 iterations for the visualization of the solution convergence.It shows a fast convergence after four iterations to zero values and to reach the global optimum of the coordinates of P with least squares.It must be noted that the condition number of the normal equation matrix was stable and near unity value in the iterations.

Figure 9 .
Figure 9.The log plot of convergence of iterations to zero values in the first test.
As shown in the first test, a far approximate coordinates of [1, 1, 1] are selected to test the stability and convergence of the triangulation with inclined angles.The correction converges to zero after eight iterations and the solution reached the global minimum as shown in Fig.12.The second iteration shows a large correction value because we started from a far approximation as mentioned earlier.The standard deviations after the triangulation is computed in meters in three methods as follows: with Inclined angles:  =0.012,  =0.005,  =0.139 with H-V angles:  =0.034,  =0.016,  =0.014 with Collinearity:  =0.034,  =0.019,  =0.015 Accordingly, improper starting values like[-1000; -1000; 500]   in the simulated test of Fig.8produced a wrong triangulation despite the convergence of the solution to a local minimum.This problem is caused because the intersecting angles can be satisfied in the opposite direction as shown in Fig.14along the dotted red line.

Figure 14 .
Figure 14.The possible solutions of triangulated angles along the dotted red circle.

Table 2 )
: Image Horizontal angles Vertical angles

Table 2 .
horizontal and vertical angles computed from each panorama to GCP point.

Table 4
lists the predetermined camera exterior orientations with the measured pixel coordinates of the object point in each image.

Table 4 .
Orientations and measurements of the second test.