Application of simulated annealing algorithm for 3D coordinate transformation problem solution

1 Introduction

A three-dimensional coordinate transformation is commonly used in photogrammetry, engineering geodesy, geographical information systems or computer visualisations of spatial objects. It consists of converting the spatial coordinates from the original coordinate system to the target one. Similarity transformation is the most frequently used mathematical model, which determines seven parameters, of which three determine translation, three rotation and one change in scale. This method is also known as Helmert transformation. Although it was proposed in 1893 [1] (as a 2D transformation) it is still being analysed and improved [2,3]. The Helmert transformation of point i coordinates is described by equation (1):

where a_i = [x_i, y_i, z_i]^T is the position vector of a given point in the original system, b_i = [X_i, Y_i, Z_i] is the position vector of the point in the target system, t = [ΔX, ΔY, ΔZ]^T is the translation vector and λ defines scale parameter. R is a 3 × 3 rotation matrix. Its elements are the functions of three angles of rotation.

Usually, to determine the values of transformation parameters a group of common points is used, for which coordinates are known in both systems. Helmert transformation, which is the subject of this study, requires at least three such points. In some cases it is possible to unambiguously determine the seven transformation parameters on the basis of XYZ coordinates of two common points and one coordinate (e.g. Z) of a third point. In no event can the points be collinear.

There are many algorithms used to determine transformation parameters based on the coordinates of common points. They can be divided into two groups: analytical procedures and iterative numerical procedures. The first group covers the algorithms like Procrustes algorithm [4,5] proposed by Grafarend and Awange, or algorithm based on distribution of eigenvalues and eigenvectors [6,7]. Analytical algorithms are considered to be complex, therefore they are not customarily applied in practice. They do, however, have the advantage of being able to determine the definitive values of transformation parameters in a direct way. Moreover, they do not require the knowledge of approximated values of those parameters.

Iterative algorithms, being easier to compute, are more commonly used. They usually utilise linearisation of the dependencies binding the transformed coordinates with the transformation parameters. This method, however, requires the knowledge of approximated values of those parameters. While determining approximate values of the parameters describing translation or scale change is relatively easy, determination of rotation requires more advanced methods. This is particularly true for large angles of rotation. This category should include the algorithm utilising modified Gauss–Newton method [8], Pareto algorithm [9], algorithms analysed in reference [10] or algorithm using total least squares (TLS) method [11].

2 Parameterisation of rotation

Definition of parameters describing the rotation components is of particular importance in case of iterative algorithms. Euler angles (φ, ψ, θ) corresponding to each subsequent rotation about the axis of the coordinate system are customarily used. In photogrammetry (especially aerial photogrammetry) ω, φ, κ angles are used. The problematic aspect of using angles of rotation lies in the fact that their value is discontinued after exceeding 360°. Another common characteristic of the abovementioned solutions is the need to calculate trigonometric functions in order to identify elements corresponding to these angles in the rotation matrix. Admittedly, operational efficiency of currently produced computers allows for a slightly less restrictive approach to computational optimisation, nevertheless, in the case of iterative algorithms, by nature requiring multiple iterations of computational tasks, the question of cost/time of the numerical operations should be taken into consideration.

The solution utilising quaternions to describe rotation of the coordinate system is devoid of the mentioned shortcomings. William Rowan Hamilton, an Irish astronomer living in the eighteenth century, is the author of quaternion algebra [12]. Some contribution to the development of the idea of using quaternions to describe transformations in 3D space was made by Olinde Rodrigues, a French banker living in the same time period [13]. Interesting details about this discovery can be found in reference [14].

The details on quaternion theory as well as its application are commonly known. An inquisitive reader may look into a book by Hanson [15]. In the later part of this paper, the author will narrow the scope of the study to the key aspects of applying quaternions to describe transformations of 3D space.

A quaternion is a four-element numerical structure, which when simplified can be expressed as the following equation:

Its elements fulfil the following equation:

A rotation matrix corresponding to the quaternion (Rodrigues’ matrix) can be defined as follows:

Each of the four quaternion elements defining Rodrigues’ matrix may take values in the range of 〈−1;1〉 (meaningful when quaternion elements are acquired with non-algebraic method). Important characteristic of the rotation matrix (4) in relation to quaternion is the symmetry with respect to zero. It can be easily checked that quaternions Q and −Q return the same rotation matrix.

It has to be noted that in contrast to the rotation matrix created based on angles, in order to calculate all nine elements of rotation matrix it is not necessary to carry out numerically cost-heavy trigonometric calculations. It is a gain in spite of the fact that we have to deal with one more rotation parameter and respect constraint (3). Another important advantage of Rodrigues’ matrix is linked to the fact that the functions describing the variability of individual quaternion elements in rotation function are continuous (there is no problem of exceeding 360° as it was the case for angles).

The above-mentioned advantages resulted in Rodrigues’ idea becoming widely used in various, often distant scientific and technological disciplines like computer vision [16], astronomy [17,18] or chemistry [19,20].

3 Simulated annealing

Dissemination of computer calculation methods and the related decrease in the cost of computational operations caused a shift in the way we look at estimation tasks and the algorithms used for those purposes. This new approach resulted in the appearance of new methods from the Monte Carlo or Metaheuristics Family, where numerous repetitions of a computation sequence are required to obtain the results.

Simulated annealing (SA), also known as Monte Carlo annealing, probabilistic hill climbing or stochastic relaxation, is an algorithm which belongs to metaheuristic methods. The idea behind SA comes from thermodynamics and reflects the process of solidification of a liquid into a crystalline solid. As the cooling process progresses molecule mobility decreases. If the cooling rate is sufficiently slow, the molecules can achieve the state of mutual order corresponding to the lowest energy state (e.g. create crystal lattice).

SA algorithm was first published in the study by Metropolis et al. [21]. Since then it has been refined by other scientists, like Kirkpatrick et al. [22]. More information on SA algorithm can be found in reference [23]. The most interesting characteristic of annealing process reflected in the algorithm structure is the ability to avoid local minima corresponding to pseudocrystalline state with energy level higher than minimal. However, this will only take place if the decrease in temperature is sufficiently slow. The idea of solution seeking for one-dimensional objective function is illustrated in Figure 1.

Figure 1

Idea of solution seeking for one-dimensional objective function.

Figure 2

Flowchart of the SA method.

The properties of SA algorithm, especially its ability to solve complex computational problems with local minima, led to a situation, where it is used in fields sometimes very distant from thermodynamics. The famous travelling salesman problem is frequently used as an example. Metaheuristic methods have been used to solve some optimisation problems in geodesy and surveying techniques e.g. references [24,25,26] but uses of SA in this field are limited. They include research on geodetic network design [27,28] or terrestrial laser scanning survey design [29]. Although this publication concerns an issue (coordinate transformation) of a much broader spectrum of applications, it proposes a new use for SA algorithm in the field of surveying techniques.

SA is an iterative algorithm, in which constant temperature change of a cooling liquid is replaced by incremented changes realised in subsequent iterations. Its use for solving estimation tasks requires a definition of several essential elements:

• Objective function corresponding to the molecular energy level during annealing process that will be minimalised

• Cooling scheme comprised assumed initial temperature and dependency defining temperature drop after each iteration

• Molecule movement model corresponding to the actual temperature

• Termination criteria for the iterative process. The condition can be formulated on the basis of:

  • final temperature (minimal),

  • acceptable value of the objective function,

  • range of molecule movement that usually corresponds to estimated model parameter changes defining the objective function.

By defining: x_i – vector of the current model parameters in the ith iteration of the estimation task, f(x_i) – objective function defined for the current parameters, and Δx_i – change in parameter vector in the ith iteration, the SA algorithm can be presented using the following scheme:

4 Concept of a hybrid algorithm

The most important advantage of the SA algorithm is the ability to find a solution (global minimum) for a task where local minima are present. Slow cooling rate and non-zero probability of accepting worse solution are crucial in such a case. It is not without significance that the starting point of the procedure has no significant impact on the speed and possibilities of finding a solution. This means that in practice there is no need to know the approximate solution and that the starting point of the process can be assumed freely. These advantages lose their significance in the last stage of the process of finding a solution, when the current solution is the close vicinity of the sought minimum for which the objective function changes in a regular manner. It results from the fact that, in contrast to gradient methods, SA procedure uses only the value of the objective function ignoring completely the character of its distribution. For the same reason and due to the stochastic character of the solution search process itself, when using SA method such solution will always be an approximation. By proper control of termination criteria it is only possible to cause the difference between the obtained result and the global minimum to be negligible for the given task.

For the above reasons, a concept of a hybrid algorithm was born. It is comprised of two extremely distinct algorithms executed sequentially depending on the phase of solution seeking process. First the SA algorithm is used. Its ultimate goal is to find an approximation for the quaternion parameter values that will make it possible to apply the chosen gradient method in the next stage. In this case, it is the widely known least squares (LS) method using linear dependency between the coordinates of the common points and the seven parameters of 3D transformation.

Another advantage of LS lies in the consistent evaluation or result accuracy estimation – both in transformation parameters as well as its results (transformed coordinates). As the LS does not include mechanisms of dealing with gross errors, in practical applications it should be supported by the quality control routine. It can be based on LS itself [30]. It is also possible to apply another robust method such as the previously mentioned TLS [11].

5 Numerical realisation of the hybrid algorithm

Spatial coordinate transformation is described by equation (1). Assuming that the a_i vector is a position vector (related to the origin of the coordinate system) the t translation vector describes the translation of the origin point of the original coordinate system. Therefore, the t vector corresponds to the position of that point in the target system. This situation is not favourable, because it is not possible to determine the approximate t components based on the coordinates of the common points – especially if the points lie within a considerable distance from the origin of the coordinate system. Consequently, in practical realisations the concept of a base point for a group of common points is utilised. In this case equation (1) takes the form (11)

where a₀ and b₀ vectors accordingly represent the position of the base point in the original and target system. The u = [dx, dy, dz] vector is a residual vector of the translation. The centre of gravity of the common points or a selected point belonging to this group can be used as the base point. Such approach makes it possible to easily determine all transformation parameters excluding the rotation components.

6 Numerical test

The way the proposed algorithm operates will be illustrated on two numerical examples.

The first test object consists of eight points. Four of them (1–4) are common (reference) points and four others (5–8) are control points. Control points are not used to calculate transformation parameters. Figure 3 shows distribution of both groups of points on XOY plane.

Figure 3

Test object 1 points placement.

Spatial coordinates of all points are presented in Table 1.

Solution seeking process for the SA method is best illustrated using the so-called Markov chain. It consists of a sequence of states (in this case, sets of q₂, q₃, q₄ values) corresponding to subsequent, at the time best, objective function values. Due to the stochastic character of the solution seeking process every time the procedure is stared it takes different way to approach the final solution. For that reason it is not possible to obtain two identical Markov chains by two independent starts of the SA procedure. Analysing Markov chains of SA is essential while optimising the key parameters of the procedure in the particular task.

To better illustrate how the SA algorithm is working out approximate rotation parameters the procedure was repeated three times. The Markov chains describing course of the solution seeking process for the test 1 were juxtaposed in Table 2. In the case of SA method, the total number of iterations was constant, as it is dependent on the initial value of the standard deviation value for quaternion elements – σ₀, standard deviation decreasing rate – β and termination condition for the iterative process. For the values assumed here (previous chapter) it was equal to 230. The it. nr column presents the iteration numbers for which the SA procedure found the better solution.

The last (final) set of quaternion parameters set is marked in bold. As it can be seen, each of the three trials delivered quite similar parameters values. Differences in q₂, q₃ and q₄ don’t exceed 10% of their values.

The unknowns vector norm of the LS first iteration (also marked in bold) can be treated as an assessment of the SA result. The final solution being the result of application of LS was always identical. Residuals on common (bold) and control points are presented in Table 3.

Due to proper number and distribution of the first test object common points, finding transformation parameters is a well-defined and relatively easy task. To test the effectiveness of the SA algorithm in a more demanding case another test object was invented. The second test object consists of three common points and three control points. The placement of common points was selected in a very disadvantageous way. They form an elongated triangle. The layout of all the points is presented in Figure 4 and their coordinates in Table 4. All the measures in Figure 4 are in meters.

Figure 4

Test object 2 points placement.

The coordinates in system 2 were obtained by rotation and translation of the established points in system 1 using equation (1), with rotation matrix R created on the basis of commonly known equations used in photogrammetry presented i.a. in reference [32] (equation (2.3)) for rotation angles ω = −167°, φ = 195°, κ = 45° respectively. A translation vector t = [500.00,500.00,200.00] was set. Coordinates in system 2 were distorted with errors generated based on normal distribution with standard deviation σ_xyz = ±0.02 m.

As in the previous test the procedure was repeated three times. The Markov chains describing course of the solution seeking process for test 2 were juxtaposed in Table 5.

Residuals on common and control points are presented in Table 6.

Despite the significant disadvantage resulting from the placement of common points the SA procedure was able to deliver angular transformation parameters that were good enough, so that the subsequently introduced LS procedure was able to provide the exact solution after a few iterations. The Markov chains (columns q₂, q₃ and q₄) presented in Table 5 point to a high diversity in the way the solution was reached. Relatively big differences present themselves in the final quaternion component values. It is especially true of q₃ component. The previous test and other experiments carried out by the author that are not presented in this paper indicate that these discrepancies are closely linked to point set geometry and can be reduced to small extent by increasing the precision of the SA method (e.g. by slowing down the standard deviation decreasing rate or delaying the termination of solution seeking process by decreasing value of σ_end).

Disadvantageous geometry of the common points is also responsible for large residuals on the control points shown in Table 6.

It is worth remembering that even though all tests presented in Table 5 were successful, due to the random nature of the SA procedure itself one cannot eliminate the possibility of a failure, which will be understood as obtaining a solution distant from the final one. The probability of a failure increases the more the SA parameters are optimised for processing speed. In such case an attempt to continue the process using LS will lead to a situation, where two subsequent iterations i and i + 1 will produce an inequality |δx_i+1| ≥ |δx_i|, which will be interpreted as an error and will terminate the LS iteration process. If SA parameters have been selected correctly, a repeated attempt to solve the task is the best countermeasure.

7 Conclusions

The performed computational tests confirmed the predicted effectiveness of the SA procedure in solving a fundamental 3D transformation problem, namely to determine approximate values of parameters defining the change in orientation. While the first numerical test used to check and demonstrate the algorithm using four regularly distributed common points can be considered as easy, the second one utilised an extremely disadvantageous placement of three common points. Poor geometry of common points caused the distribution of the objective function in the space of the sought parameters to be less clear. In result, quaternion elements being rotation parameters obtained as a result of SA procedure in the second test were much more differentiated among subsequent trials.

It should be noted that the length of the Markov chains presented in Tables 2 and 5 did not exceed 10 nodes for 230 iterations. Such nature of the process results from the used parameters which caused the SA procedure to progress relatively quickly (β = 0.99) and finish relatively early (σ_end = 0.1). The selection of parameters proved sufficient even in the case of unfavourable geometry. The unfavourable geometry of the common points in test 2 also caused that, on average, one more iteration of LS were needed to obtain the final solution.

The author is aware of the shortcomings of the LS method, of which the most often referenced one is insufficient robustness in case of presence of gross errors. That is why its application in the described algorithm is of demonstrative value only. Nothing prevents software used for real 3D tasks from utilising a robust estimation method e.g. TLS.

Due to the high probable diversity of cases as well as high dependency of the course of the SA procedure and the precision of the delivered results on its control parameters, one should provide an advanced user the possibility to modify these parameters in order to optimise the course of calculations for specific cases.

In the case of a numerous common points their number can be limited in order to accelerate the SA algorithm operation process (only approximated values are needed).

Parameters of SA procedure were optimised for potentially least optimal case and were assumed as constant. It seems that in case of more advantageous distribution of common points and larger number of them it is possible to accelerate the solution seeking process. Formulating rules for customisation of SA parameters with regard to the task’s characteristic needs requires further research.