A Novel Absolute Orientation Method Using Local Similarities Representation

Absolute orientation is an important method in the field of photogrammetry. The technique is used to transform points between a local coordinate reference system and a global (geodetic) reference system. The classical transformation method uses a single set of similarity transformation parameters. However, the root mean square error (RMSE) of the classical method is large, especially for large-scale aerial photogrammetry analyses in which the points used are triangulated through free-net bundle adjustment. To improve the transformation accuracy, this study proposes a novel absolute orientation method in which the transformation uses various sets of local similarities. A Triangular Irregular Network (TIN) model is applied to divide the Ground Control Points (GCPs) into numerous triangles. Local similarities can then be computed using the three vertices of each triangle. These local similarities are combined to formulate the new transformation based on a weighting function. Both simulated and real data sets were used to assess the accuracy of the proposed method. The proposed method yields significantly improved plane and z-direction transformed point accuracies compared with the classical method. On a real data set with a mapping scale of 1:30,000 for a 53 km × 35 km study area, the plane and z RMSEs can be reduced from 1.2 m and 12.4 m to 0.4 m and 3.2 m, respectively.


Introduction
Absolute orientation is a fundamental procedure in the field of photogrammetry and plays a crucial role in the transformation of three-dimensional (3D) points from a local to a global (geodetic) coordinate reference system.The technique has been widely used in numerous applications, such as the construction of 3D city models using multi-view stereo mobile mapping systems [1][2][3], the generation of 3D maps [4,5], the reconstruction of buildings [6][7][8][9] and 3D human models, and the manufacture of 3D physical models of large statues [10,11].
The classical absolute orientation method accurately represents the geometric relationships between 3D points in the local and global coordinate reference systems by means of a similarity transformation, which consists of one scaling parameter, three rotation angles and three translation vector parameters [12].The problem to be solved can be defined as a non-linear least-squares problem based on these seven unknown parameters for three or more pairs of 3D points in two coordinate reference systems.The solution can be obtained via a non-linear optimization method, in which the optimal values in the unknown variable space can be acquired through sequential iteration from some initial conditions.Blostein et al. [13] studied the global minima of this non-linear optimization problem using an iterative solution method to yield precise estimations.Although the iterative method can produce accurate estimations, the transformation parameters require initialization.In addition, divergence may occur if the Hessian matrix is singular.Many studies have focused on obtaining a closed-form solution without initialization.Horn et al. [14] proposed a closed-form absolute orientation solution using orthonormal matrices.However, this algorithm often produces an incorrect rotation matrix.Aurn et al. [15] solved this problem by utilizing another closed-form method based on singular value decomposition (SVD) to stably determine the transformation parameters.
In these previous methods, a single similarity transformation with a single set of seven parameters was used to represent the relationships between objects in two coordinate reference systems, as applied to aerial photogrammetry.However, if the study area is sufficiently large, these methods produce inaccurate point transformations.These inaccuracies arise primarily because point transformations with large deformation in a large-scale study case cannot satisfy the requirements of rigid-body transformation, as represented by the single set of seven similarity transformation parameters used in the classical method.Therefore, it is infeasible to accurately represent the transformation of two such point sets using a single similarity transformation of this type, and thus, a larger root mean square error (RMSE) is incurred.This study proposes a new transformation method that combines various local similarities rather than relying on a single similarity transformation.A local similarity refers to a set of similarity transformations computed using three pairs of 3D points in a local region.A Triangular Irregular Network (TIN) model is used to divide the given pairs of 3D points into various triangles.Local similarities can then be calculated using the three vertices of each triangle.The new transformation is expressed in terms of the local similarities calculated from all triangles.These similarities are combined by means of a weighted power function.
This paper is organized as follows.Section 2 briefly introduces the classical absolute orientation method based on a single similarity transformation.Section 3 introduces the proposed novel absolute orientation method based on local similarities.The performance of the proposed method is demonstrated using six simulated data sets and a real data set in Section 4. Section 5 presents the conclusions.

Single Similarity Transformation
The classical absolute orientation method searches for the optimal rigid-body transformation to perform a 3D point transformation from a local to a global coordinate reference system.This process is commonly expressed in terms of a similarity transformation with only seven degrees of freedom.A single set of similarity transformation parameters is estimated to transform points from the local to the global coordinate system.The classical method is also referred to as the single similarity transformation in this study.The single similarity transformation can be written as follows: where P L and P G refer to corresponding points in the local and global coordinate reference systems, respectively, and s (scaling parameter), R (rotation matrix) and t (translation vectors) constitute a single set of seven unknown similarity parameters.If an object can be continually rotated about the y-, xand z-axes [16], as illustrated in Figure 1, then the detailed rotation matrix can be written as follows: The classical method of transforming local 3D points (blue) to global 3D points (red), in which a single similarity transformation (s, R, t) is used to fully represent the transformation between the two coordinate systems, is illustrated in Figure 2.
Three equations can be formulated for each pair of 3D points.Therefore, it is necessary to estimate the seven parameters using at least two and one-third pairs of 3D points.Redundant observations exist for a single set of similarity transformation parameters when three or more pairs are available.The unknown parameters can be optimized using an iterative optimization method.The optimization problem can also be solved by means of a closed-form solution that makes full use of SVD.Detailed SVD solutions have been provided in previous studies [14,17].
using the classic absolute orientation method.

Local Similarities Representation
The measurements for the classical non-linear least-squares problem involve all 3D points in the global coordinate reference system.These points are used to optimize the seven similarity transformation parameters for the absolute orientation problem.Here, these 3D points are considered to play the same role in the point transformations, producing similar constraints.However, the distribution of these points can be expected to impact the accuracy with which the absolute orientation problem can be solved.
This section discusses the newly proposed transformation method, which uses the local similarities representation related to the geometric topography of the Ground Control Points (GCPs).A TIN model is used to construct the GCP topography, and various local similarities are computed based on each triangle in the TIN model.The new transformation is then estimated based on these local similarities.Unlike in the classical absolute orientation method, the number of sets of estimated transformation parameters is not one but rather equal to the number of local points, which should be transformed into the global coordinate reference.These sets of estimated transformation parameters are used to directly transform points in the local coordinate reference system into the global coordinate reference system.A weighting function is proposed to combine the various local similarities based on the point distribution.

Generation of the Triangular Irregular Network Model
TINs are digital data structures that are widely used to represent surfaces in geographic information system (GIS) and photogrammetry studies [18,19].The first TIN model was proposed in the early 1970s as a simple method for building a surface using a set of irregularly spaced points.TIN models can be broadly applied because of the manner in which they describe and preserve spatial relationships [20,21].A TIN model is used in the proposed transformation method to divide a region into numerous adjacent triangles using the given GCPs, which are located on the vertices of these triangles in the local coordinate reference system.Seven similarity transformation parameters can then be retrieved using the three pairs of GCPs from each triangle.The similarity transformations thus retrieved from each triangle are called the local similarities, which best represent the relationships between two corresponding triangles in the two coordinate reference systems.The TIN model is constructed from the GCPs such that the Delaunay triangle criterion is satisfied at the local scale, meaning that a circle that passes through the three GCPs of a given triangle will not contain any other points.The minimum interior angles of all triangles are maximized when all such circles passing through three GCPs satisfy the Delaunay criterion.
The triangulation growth algorithm [22] illustrated in Figure 3 is used to generate the TIN model from the GCPs.A simple example of TIN construction using 15 GCPs is explained as follows: (1) initial Delaunay triangle generation.First, locate a GCP "A" that is close to the geometric center of all GCPs.Then, search for the initial base line AB by finding another point "B" that is closest to "A".The initial Delaunay triangle is generated by finding a third point "C" that satisfies the Delaunay triangle criterion; (2) triangular network growth.The three closest points "D", "E" and "F" that satisfy the Delaunay criterion for each edge of the initial triangle are identified; (3) TIN model construction.The triangular network is developed by repeating (2) until all of the GCPs have been incorporated as triangle vertices.

Estimation of the New Transformation Based on Local Similarities
A set of seven similarity parameters can be effectively estimated from three GCPs, as stated in Section 2. Various sets of seven similarity transformation parameters can thus be obtained, each from a different set of three triangle vertices.Such a set of similarity transformation parameters based on three triangle vertices is called a local similarity, and these local similarities are used to estimate the new transformation proposed in this study.
Suppose that n triangles exist for a given point set containing N pairs of GCPs and their corresponding local points.Various sets of local transformation parameters (s i , R i , t i ) can be retrieved from each triangle i(1 ≤ i ≤ n) using the method described in Section 2. Each triangle exerts a different influence when transforming a local coordinate reference Check Point (CP) P L to the global coordinate CP P G .Thus, the local coordinate P L can be transformed to the global coordinate P G , the proposed transformation, which is a function of the local similarities, as follows: where w i (P j L ) is the weight of the i th triangle, which is determined by a weighting function that will be discussed in the next section.This function defines the contributions of the various local similarities to the point transformation of P j L .The eight GCPs and the CP (A) illustrated in Figure 4 provide a simple example.Eight triangles are constructed based on the TIN model, in which each vertex corresponds to a GCP.The three vertices of each triangle are used to obtain eight local similarities, which can be written as {(s i , R i , t i ), i = 1, 2, ..., 8}.Meanwhile, the weights of the eight triangles are denoted by w ∆1,2,5 , w ∆2,5,6 , w ∆5,6,8 , w ∆5,7,8 , w ∆4,5,7 , w ∆3,4,7 , w ∆1,3,4 and w ∆1,4,5 .Point A, represented by P A L in the local coordinate reference system, can be transformed into point P A G in the global coordinate reference via: The weighting values should satisfy the following constraint:

Weighting Function for Combining Local Similarities
The weighting function introduced in Equation ( 3) is used to combine the local similarities calculated from all triangles in the TIN model to estimate the new transformation.This weighting function can be expressed as a function of the point distribution, which can be measured based on the distances between a CP and the three vertices of each triangle in the local coordinate reference system.The weights calculated using the weighting function should be small for triangles far from the CP and large for triangles closer to the CP.
The red dashed lines in Figure 3 represent the distances between a CP and the three vertices of a nearby triangle, whereas the blue dashed lines represent the distances between a CP and the three vertices of a distant triangle.The weight I ∆1,2,5 of triangle ∆ 1,2,5 with respect to a given CP can be expressed as a power function: where , and d 5A represent the distances between the three vertices of the triangle and the CP.The power-law index q is set to 60 for optimal accuracy in this study.The influence of the power-law index on the accuracy is analyzed in detail in Section 4.3.
The constraint in Equation ( 5) can then be satisfied as follows.For n triangles, the normalized weight of the i th triangle with respect to the j th CP P i L is given by

Experiments and Discussion
To analyze the performance of the absolute orientation method, two sets of 3D points in two different coordinate reference systems should be provided.For the analysis presented in this paper, the bundle adjustment (BA) method, a type of aerial triangulation, was utilized to generate such 3D points using two-dimensional (2D) image points as input to serve as the measurements for the BA non-linear optimization problem.Note that the proposed method is also valid for other applications besides aerial triangulation.For example, 3D LiDAR points can be transformed between two coordinate systems using the proposed absolute orientation method.
This section reports the use of six simulated data sets and one real data set representing an aerial photogrammetric terrain scenario to assess the accuracies of the classical absolute orientation method and the proposed method.For convenience, 3D points in the global coordinate reference system are called global points, whereas points in the local coordinate reference system are called local points.

Simulated Data Sets
For an aerial triangulation application, the absolute orientation method followed by free-net BA can be used to obtain an initialization for the BA solution with control points (BA-CP).The absolute orientation method can be used to transform 3D tie points and camera poses in the local coordinate reference system to ensure that the initial unknown parameters for BA-CP are located in the global coordinate reference system.
The main workflow for these simulations is explained as follows.When image points of 3D tie points and GCPs are available, these image points can be input into a free-net BA model to simultaneously optimize the exterior orientation (EO) of the camera and the 3D points by minimizing the total difference between the predicted and measured image points.Then, the various 3D points corresponding to the input image points are triangulated.For free-net BA, the origin of the coordinate reference system is often defined as the location of the first camera; thus, these triangulated 3D points should be located in the local coordinate reference system.Now, the local and global coordinates of the GCPs are known, and thus, the transformation parameter estimation is complete.Finally, the global CPs are used to assess the point transformation accuracy of the absolute orientation method.
To obtain 2D image points, it is necessary to generate a synthetic terrain with 3D tie points, GCPs and CPs in the global coordinate reference system.The terrains of the first four data sets in Table 1 correspond to village topologies, ranging from approximately 0 m to 50 m in elevation, as illustrated in Figure 5. Additionally, mountain topologies were also simulated in the test data, corresponding to Test 5 and Test 6, as shown in Figure 6.The maximum elevations of the two mountain topologies are approximately 300 m and 500 m, respectively.The numbers of GCPs, CPs and tie points are listed in the eighth, ninth and tenth columns, respectively.In Figures 5 and 6, GCPs and CPs are indicated by red triangles and green circles, respectively.Both the GCPs and CPs are evenly distributed throughout the study area, which encompasses a 200 km × 100 km area.In addition, a Gaussian noise σ = 5 mm were imposed into the planimetric and elevation coordinates of GCPs and CPs in all simulated datasets.The second column refers to the flight height, and the mapping scales are calculated based on a fixed focal length of 0.12 m for the camera used.The areas considered in the six cases are equal, as listed in the fourth column of Table 1.To generate 3D local points when 3D tie points are available in the global coordinate reference system, 2D image points are first obtained and then input into the BA model.The use of a camera with a maximum image resolution of 7680 × 13,824 pixels that is free of geometric distortion is assumed.The camera is placed at the specified flight height to generate each sequence of aerial images.The forward and side overlaps are as large as 60% and 30%, respectively.The flight trajectories are pre-defined, i.e., the EOs are known.Then, the 3D tie points, GCPs and CPs can be projected onto image films to form 2D image points based on the collinearity equation.These projected image points have zero errors.To make the simulated data more realistic, some noise is added to the projected image points; here, a Gaussian noise of σ = 0.3 pixels is randomly added.These image points with Gaussian noise are input into a free-net BA model to triangulate the local 3D points.Here, ParallaxBA [23] is used as the free-net BA model because of its good performance, including accurate estimation and a rapid convergence rate.The performance of the proposed method was assessed and validated via comparison with the results of the classical method.The CP accuracies achieved on the six simulated data sets are presented in Table 2, where CAO and PAO refer to the classical absolute orientation method and the proposed absolute orientation method, respectively.Columns three through six represent the RMSEs in the x-, yand z-directions.The proposed method yielded improved accuracies for all data sets compared with the classical method.For example, the plane and z RMSE values for Test 1, corresponding to a flight height of 2500 m, were reduced from 0.433 m and 0.866 m to 0.021 m and 0.077 m, respectively, corresponding to increases in the plane and z accuracies by factors of 21 and 11, respectively.This accuracy increase can be attributed to the consideration of the deformation between the two coordinate reference systems, which is not addressed in the classical method.Noise is present in the extracted image points, which introduces error that is propagated into the triangulated points.Thus, reconstructed terrain deformation occurs between the local and global coordinate reference systems.The classical transformation between the two coordinate systems uses a rigid-body transformation and assumes that a single similarity transformation with seven degrees of freedom can sufficiently represent the transformation of deformed objects.However, this assumption cannot simultaneously satisfy the transformation requirements for two deformed objects.Therefore, the proposed method utilizes a representation based on local similarities to express the relationship between the two coordinate systems instead of only a single similarity transformation, thereby decreasing the deformation induced by the noise in the image points.In addition, another method, BA-CP, was tested for comparison with the absolute orientation methods, as shown in Table 2.It is evident that BA-CP showed the best accuracy among the three methods, although that of the proposed method is only slightly poorer.Clearly, for aerial triangulation, the BA-CP method cannot be entirely replaced by the absolute orientation method.However, the primary purpose of this paper is to improve the accuracy of the point transformation that can be achieved using the classical absolute orientation method.
Finally, the CP errors in the plane and z-directions are presented to illustrate the error distribution.Only the plane and z errors from Test 1 are shown in Figures 7 and 8.The green triangles and purple curves in Figure 7 represent the CPs and contour lines, respectively.The blue and red arrows represent the error directions and magnitudes.The horizontal and vertical axes are scaled by factors of 10,000 and 100,000, respectively, for clarity.The z-direction error distributions of each CP were assessed for both methods, and the results are plotted in Figure 8.The blue and red points represent the residual errors of the classical and proposed methods, respectively.Figures 7 and 8 demonstrate that the plane and z errors of the proposed method are generally smaller than those of the classical method.
The time consumed by the classical and proposed absolute orientation methods are presented in Table 3.The efficiency of the proposed method is found to be slightly poorer than that of the classical one.The primary reason for this finding is that only one set of similarity transformation parameters must be calculated in the classical method, whereas it is necessary to compute numerous sets of transformation parameters using Equation (3) in the proposed method.

Real Data Set
The local and global 3D points of the real data set used to assess the performance of the absolute orientation methods were generated in a manner similar to that described in Section 4.1.However, unlike in the case of the simulated data sets, the 2D image points could be extracted using point feature extraction technology.
The real data set considered here includes 737 Digital Mapping Camera System (DMC) images captured at a height of approximately 3500 m and encompassing an area of 53 km × 35 km.The mapping scale for this data set is 1:30,000.The images are distributed among 12 tracks, of which the forward and side overlaps are as large as 60% and 30%, respectively.Fifty-two 3D points are available in the global Beijing-54 coordinate reference system [24].In addition, image points for these 3D points are also available.Thirty-two evenly distributed 3D points were chosen as GCPs, and the remaining points were treated as CPs.The 3D coordinates of GCPs and CPs were supplied by a surveying and mapping corporation.Their planimetric accuracies were better than ±10 mm, while elevation accuracies measured by traditional leveling survey were up to ±5 mm.The GCP and CP distributions are shown in Figure 9.A feature extraction and matching algorithm was used to extract 2D image point associations to link adjacent images and provide constraints for the non-linear least-squares optimization problem addressed by the absolute orientation method.A Scale Invariant Feature Transformation (L 2 -SIFT) algorithm [25] was used because it can effectively process large-scale aerial photogrammetry data.The local 3D points could be triangulated in the local coordinate reference system based on the extracted image points and 3D image points used to solve the BA least-squares problem.The coordinate origin is commonly defined based on the first camera.These triangulated 3D local points corresponded to the input image point associations.Thus, the local points corresponding to the GCPs and CPs could also be obtained.
Based on the criteria of the TIN generation illustrated in Figure 3, we divided all GCPs into the subsets of numerous irregular TIN triangles.For each irregular TIN triangle, the local similarity could be optimized using three pair of the vertices of the triangle.Thus, the number of estimated local similarities were equal to that of the triangles.To achieve the transformation for each CP, the global similarity of the CP can be estimated by combining all local similarities based on the previously described power function.The RMSEs in all four directions were evaluated after the transformation of the local points into the Beijing-54 coordinate reference system using the proposed method, as shown in Table 4.The second row of Table 4 presents the results for the classical method.The RMSEs of the proposed method are presented in the third row.The x-, yand z-direction RMSEs of the proposed method were approximately two, three and three times smaller, respectively, than those of the classical method.The plane-and z-direction error distributions for 22 CPs are shown in Figures 10 and 11, where the errors of the proposed method are again less than those of the classical method.In Figure 10, the xand y-axes are scaled by a factor of 5000.The amounts of time consumed by the classical and proposed absolute orientation methods were 0.002 and 0.026 s, respectively.Similar to the experimental summary presented in Table 2, the results of the two absolute orientation methods are compared with those of the BA-CP method in Table 4.It is clear that the accuracy of the proposed method is very close to that of BA-CP.

Analysis of the Influence of the Power-Law Index on the Accuracy of the Proposed Method
The power-law index in Equation ( 6) plays a very important role in the weighting function, which assigns different weights to different triangles depending on the distance between the CP and the triangle.In this section, the influence of this power-law index on the accuracy of the proposed method is analyzed based on the results of the experiments.Here, the first four data sets in Table 1 and the real data set in Table 2 are used.
The plane and height RMSEs achieved using the proposed method when the power law-index is set to 0, 3, 5, 7, 10, 20, 30, 40, 50, 60, 70, 80, 90 and 100 are listed in Table 5.In addition, the plane and height RMSEs are plotted as a function of the power-law index in Figure 12.
The RMSEs are very large when the index is set to zero and the weights calculated using Equation ( 6) are always one.This means that the weighting is unaffected by the distance between the CP and the triangle.This situation disobeys the underlying assumption of the proposed method, namely, that a closer triangle (at a smaller distance) should have a greater contribution to the transformation of a CP.As the power-law index increases, the accuracies continuously improve, primarily because larger and smaller contributions are assigned to closer and farther triangles, respectively, when this parameter is large.However, once the index becomes very large (e.g., 70 for the simulated data sets), the transformation fails.The primary reason for this failure is that the denominator in Equation ( 6) tends toward infinity, causing the weight of each triangle to become zero when the power-law index is large, and consequently, the normalization of the weights in Equation ( 7) fails because the denominator also becomes equal to zero.In summary, when the power-law index is set to 60, the best accuracy can be obtained for data presented in Table 1.Thus, the power-law index 60 is used for the experiments in Sections 4.1 and 4.2 of this study.

Conclusions
A new absolute orientation method is proposed to accurately transform points between local and global (geodetic) coordinate reference systems.The transformation uses multiple local similarities instead of a single similarity transformation.The GCPs are divided into numerous subsets based on the generation of a TIN.The local similarities can be individually estimated by minimizing the differences between the fitted points and the three vertices of a triangle.The new transformation is obtained by combining the various local similarities.The results suggest that the x-, yand z-direction RMSEs can be improved by factors as large as 22, 20 and 21, respectively, based on a simulated data set acquired at a height of 2500 m.Additionally, for a real data set, these RMSEs can be reduced by factors of approximately two, three and three, respectively, compared with the classical method.In addition, the accuracies of both the classical and proposed absolute orientation methods are compared with that of another method, BA-CP.Clearly, the BA-CP method cannot be replaced by absolute orientation methods as the method that offers the best accuracy for aerial triangulation.However, for other applications in which image points are not available, the BA-CP method will fail, and the proposed method can provide a feasible solution for point transformation.
At present, the automatic determination of the optimal power-law index is a complex task.Future work will primarily focus on how to robustly determine this parameter.

Figure 1 .
Figure 1.Three rotation angles and the corresponding axes.

Figure 3 .
Figure 3.The TIN model generation procedure using GCPs for the absolute orientation problem.

Figure 4 .
Figure 4.A simple example of the proposed transformation based on local similarities.

Figure 5 .
Figure 5.The GCP and CP distributions for the first to fourth simulated data sets.(a) Test 1; (b) Test 2; (c) Test 3; (d) Test 4. The red triangles and green circles represent GCPs and CPs, respectively.Additionally, Tie Points (TPs) are drawn with blue quadrangles.

Figure 6 .
Figure 6.The GCP and CP distributions for the fifth and sixth simulated data sets.(a) Test 5; (b) Test 6.The red triangles and green circles represent GCPs and CPs, respectively.Additionally, TPs are drawn with blue quadrangles.

Figure 7 .
Figure 7.The 2D plane error distributions for each CP from the simulated data set for Test 1 assessed for the classical absolute orientation method (CAO) and the proposed absolute orientation method (PAO).

Figure 8 .
Figure 8.The z-direction errors for each CP from the simulated Test 1 data set assessed for the classical absolute orientation method (CAO) and the proposed absolute orientation method (PAO).

Figure 9 .
Figure 9.The GCP and CP distributions for the real data set.(a) Side View.(b) Top View.The red triangles and green circles represent GCPs and CPs, respectively.Additionally, TPs are drawn with blue quadrangles.

Figure 10 .
Figure 10.The 2D plane error distributions of 22 CPs in the real data set assessed for the CAO method and the PAO method.

Figure 11 .
Figure 11.The height error distributions of each CP in the real data set assessed for the classical absolute orientation method (CAO) and the proposed absolute orientation method (PAO).

Figure 12 .
Figure 12.The plane and height RMSEs as functions of the power-law index in Equation (6).(a) Plane RMSEs; (b) Z RMSEs.

Table 1 .
The overall parameters of the four simulation data sets.

Table 2 .
For simulation data sets, accuracies of CPs are resulted by the classical and proposed methods as well as bundle adjustment with additional GCPs parameters (BA-CP).

Table 3 .
Time consumed by the classical and proposed method using the simulated data sets.

Table 4 .
RMSEs achieved for the real data set using the classical and proposed methods as well as bundle adjustment with additional GCPs parameters (BA-CP).