Genetic Algorithm in the Computation of the Camera External Orientation Rudolf Urban, Martin Štroner Department of Special Geodesy Faculty of Civil Engineering, Czech Technical University in Prague, Thákurova 7, 166 00 Prague 6, Czech Republic rudolf.urban@fsv.cvut.cz Abstract The article addresses the solution of the external orientation of the camera by means of a generic algorithm which replaces complicated calculation models using the matrix inverse. The computation requires the knowledge of four control points in the spatial coordinate system and the image coordinate system. The computation procedure fits very well computer-based solutions thanks to it being very simple. Keywords: Photogrammetry, Adjustment, Algorithms, Camera, Observations, Automation 1. Introduction Genetic algorithms simulating nature’s natural development as a subset of global optimization algorithms allow simple computations of optimum values based on a selected evaluation stan- dard. The price for this simplicity, on the other hand, is very high computational demands. In more complex computations using the Least Square Method the convergent iteration compu- tation of a set of non-linear equations, for example, requires the determination of sufficiently accurate approximate values. In numerous geodetic problems, these may be determined with ease, e.g. in the area of geodetic network adjustment, but in other applications such as pho- togrammetry or laser scanning, or in geometric primitives fitting in particular, this is a rather complicated and frequently also unreliable option. Together with their computational force supported by the latest computer technology, genetic algorithms theoretically offer algorith- mically very simple methods allowing finding initial values for subsequent optimizations using the Least Square Method. Based on the knowledge of its mathematical nature, however, the whole problem must be simplified so that the solution time does not pose limitations to practi- cal usability. The computation of the elements of the camera external orientation is a problem frequently solved in photogrammetry, there can be noted e.g. [1], [2], [3]. This problem with the availability of known approximate values and redundant numbers of measurements, it may be iteratively solved using the Least Square Method. The computation of approximate values has recently been solved several times (see [4], [5]), here a very simple computational procedure using the genetic algorithm is described, which will apply in the computation of the first part of the determination of the external orientation elements, namely in the com- putation of the projection centre coordinates. The computational algorithm is complemented with a simple computation of the camera rotation angles. 2. Genetic algorithms and their potential for use in geodetic computations Genetic algorithms belong to global optimization techniques; they are based on a heuristic procedure applying the principles of evolution biology to solving complex problems. Evolution Geoinformatics FCE CTU 9, 2012 5 Urban, R., Štroner, M.: Genetic Algorith in Camera External Orientation processes like inheritance, mutation, crossover and natural selection are simulated there. 2.1. Principle of genetic algorithms The principle of how a genetic algorithm works is that it gradually generates generations of different solutions to a given problem. As these solutions (populations) undergo a develop- ment, the solutions improve. The solution is traditionally represented by binary numbers, but other expressions are also used, like fields, matrices etc. In the first generation, the popu- lation is composed of more solutions created usually by random generation. In the transition into a new generation, an evaluation is computed for each solution of the population using the evaluation (so-called fitness) function, which expresses the quality of the given solution and represents natural selection. According to this evaluation, suitable solutions are selected to be reproduced (copied) into the next generation and further modified by means of muta- tion (random change in part of the solution) and crossover (exchange of parts of the solution between individual solutions). Thus, a new generation is generated and the solution is itera- tively repeated while the population gradually improves, i.e. the quality (optimality) of the solution improves. The direction of selection and optimization are given by the evaluation function. The principles and applications of genetic algorithms and evolution procedures in general may also be found e.g. in [6], [7] or [8]. The advantage of using these algorithms is the fact that the generation of the following solutions is a very simple random process using the generator of pseudorandom numbers; there is no need for solving optimization calculations, you just need to define the evaluation function. 2.2. General flowchart of genetic algorithms The general process of the genetic algorithm’s functioning may be described in the following six points: 1. Initiation: the zero population is generated, i.e. various random solutions are generated in the whole population numbers. 2. Start of iteration • Selection of high-quality individuals with the best evaluation. • By means of crossover, mutation and reproduction a new population is generated from selected individuals (Generation Procedure, besides, another procedure is used where only one new solution is generated to replace the worst solution in the previous generation). • The evaluation function for the new population is computed. 3. End of cycle: unless the end condition is fulfilled, the procedure again continues with point 2. 4. End of optimization: the solution with the best evaluation is the sought solution. The end condition may be specified by the computation time, the evaluation quality of the best solution (provided it may be defined), the total number of cycles or the total number of generations. Geoinformatics FCE CTU 9, 2012 6 Urban, R., Štroner, M.: Genetic Algorith in Camera External Orientation 2.3. Advantages and disadvantages of genetic algorithms Among the principal advantages of genetic algorithms there is the fact that no mathematical solution of the problem is necessary – what suffices is the knowledge of the evaluation function, if correctly set up the algorithms are resistant to finding only the local optimum and are very universal. The disadvantage is the problem of finding the “accurate” optimum solution; another limiting factor may be a large number of computations of the evaluation function (if it sets high demands for computing). Yet another disadvantage is its excessive universality and resulting extreme computational demands; for this reason hybrid algorithms are frequently used instead where the knowledge of the solved problem is applied. 2.4. Application to solutions of geodetic optimization problems Optimization problems in the field of geodesy and cartography are usually solved using the Least Square Method, i.e. the optimization problem [pvv] = vT · P · v = min, where vi are corrections of measured variables (i = 1...n) and pi is the measurement weight, v is the correction vector and P is the weight matrix. For non-linear systems of equations, the iteration solution (for more detail see [9]) is known in the form dx = ( AT · P · A )−1 AT · P · l ′ , (1) where A is the Jacobian matrix (plan matrix) and l′ is the vector of reduced measurements. The corrections are determined from the formula v = A · dx + l′. (2) Part of the solution is the matrix inverse of N = AT · P·A where numerical problems occur during the solution if large numbers of measurements are solved (as discussed e.g. in [10] and applied in software in [11]). Using the generic algorithm no inversion is necessary, only corrections are computed and, therefore, numerical problems are practically eliminated in this case. Besides, the performance of the convergent iterative computation requires sufficiently accurate approximate values. In the case of not so extensive problems in common calculations in coordinate systems in geodesy or in smaller purpose-built networks this may be ensured by standard calculations of coordinates without any adjustment; in extensive solved systems like photogrammetric computations (see the solved problem) or in the area of 3D scanning where e.g. fitting of geometric primitives onto hundreds of thousands of points are solved, on the contrary, problems arise. Provided they are appropriately applied, genetic algorithms may be one of the ways of how to obtain these approximate values and successively perform the final optimization using standard methods. Another area of use is optimization according to a criterion which cannot be solved so simply as the Least Square Method. An example may be L1-minimization used in geodesy mainly for detecting large measurement errors, i.e. the solution of a set of equations under the condition [|v|] = ∑ |v| = min., or minimization according to another standard selected by a mere change in the evaluation function. A fitting application of genetic algorithms that cannot be solved by any other means is for example the optimization of the geodetic network configuration (see [12]). It must, however, be pointed out that if a given problem lends itself to a satisfactory mathematical solution, in this case Geoinformatics FCE CTU 9, 2012 7 Urban, R., Štroner, M.: Genetic Algorith in Camera External Orientation Figure 1: Illustrative schema of used coordinate systems (In photogrammetry is general used the coordinate system, which z axis is oriented to zenith and y is oriented allong the normal of the image plane). the application of genetic algorithms is generally not beneficial as it poses extreme demands for computing. 3. Calculation procedure of the camera projection centre position The first part of the calculation of the external orientation elements is the determination of the position of the camera projection centre. The underlying idea of the calculation is a tetrahedron (Fig. 1), which is delimited by the entrance pupil P and three control points A, B, C. The symbol H refers to the principal point of the image. Then, respective direction vectors of position vectors (a, b, c), the lengths of position vectors from the entrance pupil (Ra,Rb,Rc) and the apex angles of position vectors (α,β,γ) in the calculated triangle are displayed. The distance of the entrance pupil from the principal point of the image is the camera constant (denoted in further calculations by symbol f). After the third dimension is added into the image coordinate system, the direction vectors of position vectors become known and are defined by the formula a = (xA yA f)T , (3) b = (xB yB f)T , (4) c = (xC yC f)T , (5) where xi, yi(i = A,B,C) are point coordinates in the image. Before any calculations are made, point coordinates in the image must be reduced removing the coordinates of the principal point of the image and further corrected for the effect of lens distortion. The direction vectors of position vectors allow calculating the values of apex Geoinformatics FCE CTU 9, 2012 8 Urban, R., Štroner, M.: Genetic Algorith in Camera External Orientation angles of the position vectors using the formulae: cos α = (aT · b) (|a| · |b|) , (6) cos β = (bT ·c) (|b| · |c|) , (7) cos γ = (cT ·a) (|c|·|a|) . (8) 3.1. Calculation of one solution The basic idea is the solution of a tetrahedron which in general leads to maximally two different solutions of the entrance pupil position in the respective half-plane. To be able to quite accurately determine which solution is correct, another point must be added into the calculation [13] thus three more equations must be solved. Therefore, six apex angles which may be calculated among four control points enter the calculation. 3.2. Application of a genetic algorithm In terms of classification, a hybrid genetic algorithm is used here. The position of the entrance pupil is sought quite randomly. A random vector of a selected length is generated from the previous position. The length of the vector is decreased based on the evaluation criterion so that the optimum speed and optimum computation accuracy may be achieved. In this case, the zero generation point may conveniently be determined by graphic means - located into the correct half-plane and situated at the previously estimated camera distance from the object so that the computation will not take too long. The initial point and the random vector serve for the generation of a new point (new generation), which is checked against the evaluation criterion. If the new point meets the criterion of higher quality than the previous generation, it is adopted as the current solution. The procedure continues by the generation of another random vector and all is repeated again. If the criterion fails to be of higher quality, the algorithm continues from the previous position. The length of the random vector is decreased, always after a selected number of generations which fail to meet the evaluation criterion. It is, therefore, appropriate to select a relatively large step at the start of the computation so that the computation may converge faster. In this way, the algorithm will gradually determine the correct solution of the entrance pupil position. So that the algorithm may function correctly, it is necessary to appropriately select the number of unsuitable generations for the modification of the vector’s length. The greater the number, the higher-quality the computation is, being however, also longer. The above-described algorithm may be ranked in the hybrid group. Due to the absence of a greater population and crossover the algorithm is highly sensitive to local minima; in this case, however, it may be used as there is only one minimum within the solved half-plane being thus of a global type as can be seen in Fig. 2 (the chart of contour lines was created for the algorithm testing data presented in Tab. 1). The procedure is also called “Hill Climbing” and its description including further alternatives (e.g. “Hill Climbing“ with random restarts) less sensitive to the local minimum may be found in [8]. Although the evaluation criterion checks in a simple way the approximation quality, nev- ertheless it cannot be interpreted in terms of the quality of the entrance pupil positioning. Geoinformatics FCE CTU 9, 2012 9 Urban, R., Štroner, M.: Genetic Algorith in Camera External Orientation Figure 2: Global minimum of the evaluation criterion in the plane of the correct solution Therefore, the computation is terminated after the reduction of the random vector length below a certain level which thus indirectly determines the interval within which the correct solution is found. 3.3. Evaluation criterion The principal step of the genetic algorithm consists in the specification of the criterion which will reveal whether the result is better or worse. For the computation of the external orienta- tion, the criterion will depend on the magnitude of apex angles of position vectors. To be able to compare deviations from the correct solution, the angles calculated from coordinates must be subtracted from the correct angles from the image coordinates and the camera constant. These deviations have different signs and there are six of them in all. So that the evaluation criterion may be assessed with the best result, it is expedient to compare only one argument or one digit, doing so the resulting deviation may be characterized by the sum of individual squared deviations in the angle (which corresponds to the solution using the Least Square Method). A criterion set up in this way will decide whether the new position of the entrance pupil is of higher quality than the preceding one. The algorithm will, therefore, always retain in its memory only two positions of the entrance pupil and two evaluation criteria, which it successively evaluates. 4. Calculation of angles of rotation Only three identical points are necessary for the calculation of angles of rotation. This calculation procedure is based on the publication by B. K. P. Horn [15] and [16] and allows simple determination without iteration. Geoinformatics FCE CTU 9, 2012 10 Urban, R., Štroner, M.: Genetic Algorith in Camera External Orientation Figure 3: Schematic representation of three points in two coordinate systems and (According to [15]) 4.1. Theoretical basis A schematic graphic representation of the problem is in Fig. 3. The coordinates of three points in two coordinate systems (X1,X2,X3,x1,x2,x3) are known and the task is to find three angles of rotation (or the orthonormal matrix of rotation R). Identical triples of points lend themselves to setting up two mutually equivalent triples of orthonormal vectors (unit, normal to each other) so that a pair of points are used for the calculation of a vector which is normalized, then by means of a third point and Gram-Smidt orthogonalization (described e.g. in [14]) a vector normal to the first one and via a vector product a third vector are created. The mutual relationship shared by these two triples of orthonormal vectors is only rotation in space, which may be simply determined. The calculation procedure of the triple of orthonormal vectors Vx for x is (for VXthe procedure is identical, just with use of) v1 = x2 − x1, v2 = x3 − x1, (9) e1 = v1/||v1||, where ||v1|| = √ (vT1 � v1), (10) e2 = o2/||o2||, where o2 = v2 − (vT2 � e1) � e1, (11) e3 = e1 × e2, (12) Vx = (e1 e2 e3). (13) Since the following holds true R · Vx = VX, (14) then it holds true for the calculation of R that: R · Vx · V−1x = VX ·V −1 x , (15) hence R · E = VX · V−1x , (16) Geoinformatics FCE CTU 9, 2012 11 Urban, R., Štroner, M.: Genetic Algorith in Camera External Orientation Figure 4: Example of correct orientation of coordinate systems corresponding to formulas and hence R = VX · V−1x . (17) Since all matrices in formula (18) are orthogonal, it also holds true for V−1x that V−1x = VTx hence R = VX · VTx . (18) 4.2. Identical triples of points While calculating the angles of rotation the correct turning of the coordinate systems is of key importance (Fig. 3). Further, the coordinates of three control points in the image coordinate system must be recalculated into the space behind the entrance pupil. This recalculation may be performed using the entrance pupil P, the position vectors of control points (Ra,Rb,Rc) and direction vectors of the position vectors (a, b, c), all modified so that they correspond to the correct turning of the coordinate systems. Modified direction vectors according to Fig. 4 a = (xA −x0, −(yA −y0), −f)T , (19) b = (xB −x0, −(yB −y0), −f)T , (20) c = (xC −x0, −(yC −y0), −f)T (21) where xi, yi are image coordinates of points and x0, y0 are image coordinates of the principal point in the image. The differences in the coordinates of control points behind the pupil are obtained by multi- plying the unit direction vector (the vector marked with index 0 on the bottom right side) by the length of the position vector and by subsequent adding of the result to the coordinates of the entrance pupil P. Thus, the matrix definition of point A looks as follows  X P A Y PA ZPA   =   XPYP ZP   + Ra · a. (22) Geoinformatics FCE CTU 9, 2012 12 Urban, R., Štroner, M.: Genetic Algorith in Camera External Orientation Hence the input matrices of two mutually rotated systems before orthogonalization are X =   XA XB XCYA YB YC ZA ZB ZC   , x =   X P A X P B X P C Y PA Y P B Y P C ZPA Z P B Z P C   . (23) 5. Testing of the computational algorithm The algorithm was practically tested on more different cases; the most easily explainable one is presented here. Four control points were selected on a calibration field with known spatial coordinates for the experiment (Tab. 1). The used apparatus was the Lumenera Lu125C camera (camera constant f = 2445.8997, principal point in the image x0 = 677.1816,y0 = 504.3293). Prior to the application of the computational algorithm the effects of radial and tangential distortion were removed from the set of image coordinates. After the computation of the complete external orientation the results were subjected to the iteration computation of the projection transformation to confirm whether the accuracy of the results would be sufficient for the iterative solution of the external orientation in bundle adjustment. The computations were performed on a desktop computer fitted with the Intel Pentium D processor (3.0 GHz) with 2GB RAM in the Scilab 5.0.3 programme. Point No. X [m] Y [m] Z [m] x [pixel] y [pixel] 1 5001.22710 98.67664 997.50504 551.11 895.69 2 5001.55797 99.00517 997.49921 1129.16 371.59 3 5001.08773 99.15847 997.48235 338.45 74.27 4 5001.55615 99.10514 996.99086 980.61 179.56 Table 1: Coordinates of control points The testing of the computation was performed with the initial setup of the vector’s length as 0.5 m, the number of incorrect solutions of 15 for the reduction of the vector’s length into one half and with a tolerance for the iteration cycle termination with the reduction of the vector’s length to a value of 1E-15. Different starts in the respective half-space were selected for testing. Tab. 2 presents the deviations of the start from the correct solution in individual coordinate axes (dX, dY, dZ) and, further, for illustration, also the spatial distance from the correct solution, the number of iterations and the total time necessary for the computation of the entrance pupil position. The resulting position of the entrance pupil from the iteration solution was identical for all tested points (X = 5001.198 m, Y = 99.139 m, Z = 998.924 m). The resulting position of the entrance pupil from the projection transformation adjustment was (X = 5001.199 m, Y = 99.138 m, Z = 998.925 m). The rotation matrix was nearly identical due to minimum differences in the position of the entrance pupil. Hence the resulting rotation matrix is R =   0.9973281 −0.0332701 −0.06503720.0429059 0.9873119 0.1528864 0.0591255 −0.1552684 0.9861014   Tab. 2 clearly shows that the computation speed depends on the distance from the correct solution (and also on the number of iterations). However, even in the case of a considerably Geoinformatics FCE CTU 9, 2012 13 Urban, R., Štroner, M.: Genetic Algorith in Camera External Orientation Time [s] dX [m] dY [m] dZ [m] Distance [m] No. of iterations 0.238 -1.481 -2.491 0.631 2.966 489 0.141 -1.526 -2.410 2.629 3.879 288 0.189 -1.571 -2.329 4.626 5.413 381 0.153 -1.906 -0.539 0.542 2.054 306 0.126 -1.952 -0.458 2.540 3.236 255 0.135 -1.997 -0.377 4.538 4.972 279 0.191 -2.332 1.413 0.453 2.764 393 0.122 -2.377 1.494 2.451 3.727 223 0.162 -2.423 1.575 4.449 5.305 315 0.165 0.473 -2.064 0.658 2.217 315 0.199 0.428 -1.983 2.655 3.342 412 0.135 0.383 -1.902 4.653 5.041 280 0.134 0.047 -0.112 0.569 0.582 272 0.155 0.002 -0.031 2.567 2.567 321 0.461 -0.043 0.050 4.565 4.565 956 0.158 -0.379 1.840 0.480 1.939 324 0.133 -0.424 1.922 2.478 3.164 273 0.172 -0.469 2.003 4.476 4.926 356 0.425 2.427 -1.636 0.685 3.006 876 0.178 2.381 -1.555 2.682 3.910 367 0.131 2.337 -1.474 4.680 5.435 272 0.123 2.001 0.316 0.596 2.111 255 0.157 1.956 0.397 2.594 3.272 324 0.323 1.910 0.478 4.591 4.996 661 0.135 1.575 2.268 0.507 2.807 278 0.237 1.530 2.349 2.505 3.759 492 0.159 1.485 2.430 4.503 5.328 325 0.321 18.802 10.861 41.075 46.461 648 0.245 3.802 20.861 21.075 29.896 501 0.292 -1.199 -4.140 51.075 51.257 609 0.555 -21.198 -14.139 101.075 104.238 1135 Table 2: Results of testing random starting points Geoinformatics FCE CTU 9, 2012 14 Urban, R., Štroner, M.: Genetic Algorith in Camera External Orientation outlying start the computational algorithm reaches time well below 1 second, which may be seen in the last tested point in the table. 6. Conclusion The article presents a non-traditional method of the computation of approximate values of elements of the camera external orientation from four control points which is very simple and, at the same time, highly efficient if computer technology is used. It applies the genetic method where there is no need for the calculation of the matrix inverse or any other complicated matrix calculations. It is solely based on a repetitive generation of a random vector with a length depending on the number of previous unsuccessful attempts. The method was tested on practical examples; the results are reliable and the time of their determination depends on the accuracy of the estimated initial position of the camera entrance pupil. The procedure converges from a practically randomly set initial position. Acknowledgements The article was written with support from the grant of Czech Technical University in Prague No. SGS12/051/OHK1/1T/11 “Optimization of 3D data acquisition and processing for the needs of engineering geodesy”. References [1] Grunert, J. A.: Das Pothenotische Problem in erweiterter Gestalt nebst Über seine Anwendungen in der Geodäsie. Grunerts Archiv für Mathematik und Physik, Band 1, 1841, pp. 238-248,(German). [2] Haralick, R. - Lee, C. - Ottenberg, K. - Nolle, M.: Review and Analysis of Solutions of the Three Point Perspective Pose Estimation Problem. Intational Journal of Computer Vision, 13, 3, 331-356, 1994. [3] Lepetit, V. - Moreno-Noguer, F. – Fua, P.: EPnP: An Accurate O(n) Solution to the PnP Problem. International Journal of Computer Vision 81(2): 155-166, 2009. [4] Urban, R.: Solution of the Camera External Orientation from four Control Points. Pro- ceedings of the Juniorstav 2011 Conference. Brno: Vysoké učení technické v Brně, Fakulta stavební, part 1, 376 p. 2011. ISBN 978-80-214-4232-0. (in Czech) [5] Koska, B. - Pospíšil, J. - Obr, V.: Eliminations of Some Defects of the Digital Cameras Used in the Laser Scanning Systems. In: INGEO 2008 – Bratislava, 2008. ISBN 978-80- 227-2971-0. [6] Holland, J. H.: Adaptation in Natural and Artificial Systems, University of Michigan Press, Ann Arbor, 1975. [7] Mitchell, M.: An Introduction to Genetic Algorithms, MIT Press, Cambridge, MA, 1996. [8] Weise, T.: Global Optimization Algorithms - Theory and Application. Electronic mono- graph, on-line available at http://www.it-weise.de/projects/book.pdf, 20.2.2012. Geoinformatics FCE CTU 9, 2012 15 Urban, R., Štroner, M.: Genetic Algorith in Camera External Orientation [9] Böhm, J. - Radouch, V. - Hampacher, M.: Theory of Errors and Adjustment Calculus. Geodetický a kartografický podnik Praha, 2nd edition, Prague, 1990. ISBN 80-7011-056- 2. (in Czech) [10] Čepek, A. - Pytel, J.: A Note on Numerical Solutions of Least Squares Adjustment in GNU Project Gama In: Interfacing Geostatistics and GIS. Berlin: Springer-Verlag, 2009, pp. 179-193. ISBN 978-3-540-33235-0. [11] Čepek, A.: Program GNU Gama. http://www.gnu.org/software/gama/gama.cs.html. 14.4.2012. [12] Berné, J. L. - Baselga, S.: First-Order Design of Geodetic Networks Using the Simulated Annealing Method. Journal of Geodesy, 78, Springer-Verlag, 2004. [13] Kraus, K.: Photogrammetry Volume 2 - Advanced Methods and Applications. Dümmler, Bonn, Germany, 4th edition, 1997. ISBN 3-427-78694-3. [14] Štroner, M. - Pospíšil, J.: Terrestrial Scanning Systems. 1st edition. Praha: Česká tech- nika - nakladatelství ČVUT, 2008. 187 p. ISBN 978-80-01-04141-3. (in Czech) [15] Horn, B. K. P.: Closed-Form Solution of Absolute Orientation Using Unit Quaternions. Journal of the Optical Society A, 4, 629–642, 1987. [16] Horn, B. K. P. - Hilden, H. M. - Negahdaripour, S.: Closed-Form Solution of Absolute Orientation Using Orthonormal Matrices. Journal of the Optical Society of America A. Vol. 5 Issue 7, pp.1127-1135, 1988. Geoinformatics FCE CTU 9, 2012 16