IMPROVED CONJUGATE GRADIENT BUNDLE ADJUSTMENT OF DUNHUANG WALL PAINTING IMAGES

Bundle adjustment with additional parameters is identified as a critical step for precise orthoimage generation and 3D reconstruction of Dunhuang wall paintings. Due to the introduction of self-calibration parameters and quasi-planar constraints, the structure of coefficient matrix of the reduced normal equation is banded-bordered, making the solving process of bundle adjustment complex. In this paper, Conjugate Gradient Bundle Adjustment (CGBA) method is deduced by calculus of variations. A preconditioning method based on improved incomplete Cholesky factorization is adopt to reduce the condition number of coefficient matrix, as well as to accelerate the iteration rate of CGBA. Both theoretical analysis and experimental results comparison with conventional method indicate that, the proposed method can effectively conquer the ill-conditioned problem of normal equation and improve the calculation efficiency of bundle adjustment with additional parameters considerably, while maintaining the actual accuracy.


INTRODUCTION
Bundle adjustment is the problem of refining a visual reconstruction to produce jointly optimal 3D structure and viewing parameter (camera pose and/or calibration) estimates (Triggs et al., 2000).In close-range photogrammetry, since the photographic structure is diverse and complicated, bundle adjustment is applied in case of redundancy due to the combination of different types of geometric constraints in object space, which include topology constraints (such as object point, line and plane constraints, etc.) and object constraints (such as parallelism, perpendicularity and symmetry, etc.) (van den Heuvel, 1998).In the "Digital Dunhuang" project, the paintings usually lie on a near-planar wall surface and the forward overlap between adjacent images is smaller than 60%.Near-planarity constraints are combined with bundle adjustment to control the model's error propagation (Zhang et al., 2011;Zhang, Hu and Huang 2012).Hence, the coefficient matrix of reduced normal equation is no longer regular sparse-banded matrix.It's difficult to design fast solution method of bundle adjustment with additional parameters.

Review of Related Work
Taking various kinds of geometric constraints into account, McGlone (1995) constructed a unified least square adjustment model with additional parameters and solved this model using generalized inverse method.Wang and Clarke (2001) proposed a separate adjustment model based on separate least squares estimation of the camera's exterior orientation parameters and object point coordinates, which can reduce both the operation time and the memory consumption with guaranteed precision.Bartoli (2002) constructed a bundle adjustment model based on quasi-linear optimization which preserved the original cost function.Different projection models, a large number of views and points, as well as missing data are handled in a unified framework.Ni, Steedly and Dellaert (2007) presented an out-ofcore bundle adjustment model by partitioning the whole scene into several smaller scenes which possess their own local coordinate systems.This model is straightforward to parallelize and can take good advantage of the physical memory.All the methods mentioned above are focus on efficiency improvement of bundle adjustment.
Conjugate gradient method is evolved on the basis of steepest descent algorithm and is characterized by its simplicity, easy realization and low computational cost (Hestenes and Stiefel, 1952).Whereas, its application was limited seriously due to the sensitivity to rounding error.In recent years, this problem is solved effectively by combining conjugate gradient method with a variety of preconditioning methods.Thus, conjugate gradient method is widely used for solving large-scale sparse linear equations and it has been introduced into the field of bundle adjustment.Due to the property of invariance to orthogonal transformation, Byröd and Åström (2009) first combined its multi-scale representation with standard preconditioners to speed up the convergence rate of conjugate gradient bundle adjustment (CGBA).Then they further optimized CGBA by Levenberg-Marquardt method (Byröd and Åström, 2010).QR factorization based block preconditioner was presented in accordance with the sparse structure of normal equation.Agarwal et al. (2010) designed an inexact Newton type method for solving large-scale bundle adjustment, in which conjugate gradient method is used for calculating Newton step.Jian et al. (2011) adapt the idea of subgraph preconditioned conjugate gradient to solve large scale bundle adjustment problems, and proposed a generalized subgraph preconditioning technique to avoid the over-estimation of uncertainty.Jeong et al. (2012) used block-based preconditioned CGBA on the reduced camera system.Experimental results demonstrated that it is a more efficient way to perform the procedures of the bundler as well as a faster way to converge.

Contribution of This Paper
In this paper, a novel method of conjugate gradient bundle adjustment with preconditioning based on improved incomplete Cholesky factorization (improved CGBA) is developed for precise orthoimage generation and 3D reconstruction of the Dunhuang wall paintings, which processes the following advantages: 1.It is suitable for solving large-scale normal equation with additional parameters introduced by self-calibration and geometric constraints; 2. It can effectively reduce the condition number of coefficient matrix and improve the structure of normal equation in bundle adjustment; 3. It can accelerate the iteration rate and improve the calculation efficiency of bundle adjustment considerably, while maintaining the actual accuracy.This paper first illustrates the modelling and solving methods of the improved CGBA, then makes theoretical analysis and experimental comparison with conventional method.The remaining sections are organized as follows: In Section2, CGBA method is deduced based on calculus of variations.In Section 3, improved incomplete Cholesky factorization method is proposed for preconditioning, and the block diagonal matrix of the coefficient matrix is used as preconditioning matrix.The detailed solving process of improved CGBA method with preconditioning is put forward with five steps.In Section 4, the proposed method is compared with Gauss elimination method on space complexity and time complexity.Experiments are carried out by four datasets of Dunhuang wall painting images to validate the effectiveness and efficiency of the proposed method.Finally, conclusions are drawn in Section 5.

Basic Principles of CGBA
Conjugate gradient method is an iterative approximation method for solving large sparse systems with nonzero entries occurring in predictable patterns, which frequently arise in the solution of boundary-value problems.Good results are obtained when an effective preconditioning method is adopted for matrix computation (Burden and Faires 2000).
In bundle adjustment, conjugate gradient method is combined with the steepest descend method to construct a sequence of conjugated directions by gradients at known points, and further to search the minimum value of objective function along these directions.The normal equation of bundle adjustment can be written as follows: where X = unknown matrix N = coefficient matrix, it is a n-order positive definite symmetric matrix W = constant matrix In n-dimensional vector space, there certainly exists a linear independent vector group 1 2 , , n R R R  , which satisfies: is called a conjugated vector group about N , i.e. i R and j R are conjugated about N .They can be constructed by Gram-Schmidt process (Burden and Faires, 2000).

CGBA Derivation by Calculus of Variations
If the initial value of the unknown matrix X in Equation ( 1) is set to 0 X , the normal equation is transformed into: are still expressed by X and W , respectively.Then the expression of the normal equation remains the same as Equation ( 1).The problem of solving Equation ( 1) can be transformed into searching the minimum value of the following function by calculus of variations (Burden and Faires 2000).
The solution to Equation ( 1 , which is equivalent to X when ( ) S X reaches the minimum value.
Assuming that we have searched the solution of Equation ( 4) . Then the search direction in step The problem is converted to solve the minimum value of In Step 2, as proved by mathematical induction (Burden and Faires 2000), . Hence, 1 k   in Equation ( 7) can be further simplified as: , the expression of k  in Equation ( 6) can be further simplified as: There are n linearly independent vectors at most in ndimensional vector space.At least one of 0 1 , ,... n r r r is zero.Hence in theory, at most n times iterations are required in searching the solutions of normal equation by CGBA.

Solution method of CGBA
The orthogonality of residual vectors and the conjugacy of search directions cannot be guaranteed in solving large-scale sparse equations, because there exists some rounding error in floating-point numbers arithmetic operation.It is meaningful to solve the normal equation by CGBA with iteration.The solution method of CGBA is represented as the following steps: . Then return to Step 3 and judge again.

Improved Incomplete Cholesky Factorization
During the computation process of CGBA of large and complex image blocks, such as Dunhuang wall painting images, an effective preconditioning method is needed for improving the structure of normal equation and reducing the condition number of coefficient matrix.
Based on the principle of numerical analysis (Burden, 2000), the model precision of CGBA in step k can be estimated as: where As can be seen from Equation ( 10), the larger the condition number K is, the more the iteration times are, and the lower the convergent rate is.
As for the precise orthoimage generation and 3D reconstruction of Dunhuang wall paintings, the test paintings lies on a nearplanar wall surface.A pre-calibrated digital camera is fixed on an orbital platform that is near-parallel to the painting.The images are captured along regular strips, with the forward overlap and sidelap between adjacent images maintained at 50%.The principal optical axes of the images are parallel with each other and perpendicular to the painting; and the origin of the object space coordinate system is set at the lower left corner of the painting, with the X-axis parallel to the track, the Y-axis pointing to the zenith, and the Z-axis perpendicular to the painting (Zhang, Hu and Huang 2012).The photographic structure of the test paintings is shown in Figure 1.Since N is a n-order positive definite symmetric matrix, it is meaningful to combine incomplete Cholesky factorization with CGBA to reduce the condition number of the normal equation.This preconditioning method is described as follows: First of all, assuming that there is a positive definite symmetric matrix  , which can be decomposed by incomplete Cholesky factorization as follows: Where L is a lower triangular matrix.
Then, the normal equation ( 1) is equivalent to Equation (12) as: . F is also a positive definite symmetric matrix.Equation ( 11) can be simplified as: The original problem is transformed into solving Equation ( 13).
Since the condition number of F is much lower than that of N , the convergence rate of Equation ( 13) can be accelerated obviously.Therefore, the equation (1) needs to optimize for solving Y as follows: 1 In order to calculate unknowns X in Equation (1) directly, set and a variable . By introducing them into iterative computation, the solving method of CGBA based on improved incomplete Cholesky factorization is determined.Next, k Z can be calculated by solving The derivation process of improved preconditioning method indicates that the preconditioning matrix M is required as close as possible to N while possessing simple structure.

Algorithm Design of Preconditioning Matrix
It's reasonable for the block diagonal matrix of the coefficient matrix N in Equation ( 1) to be used as the preconditioning matrix, because one simple way of solving bundle adjustment is by iterative computation between resection and intersection (Mikhail, Bethel and McGlone, 2001).To illustrate this, the structure of Equation (1) can be described as: t in turn.Based on the properties of continuous function (Burden and Faires 2000), if the initial values of the parameters are given close to the truth values, the process in each iteration can be expressed by solving the following equation approximately.
The coefficient matrix of Equation ( 15) is an approximation of N .It is also an n-order positive definite symmetric matrix while possesses simple structure.Hence, it can be treated as the preconditioning matrix of CGBA.

Solution Method of Improved CGBA
Above all, the solution method of CGBA with preconditioning based on improved incomplete Cholesky factorization can be summarized as the following steps: 1. Initialize the unknown matrix , and . Calculate the value of k Z by solving . Then return to Step 3 and judge again.

Comparisons on Space and Time Complexity
To evaluate the performance of improved CGBA of Dunhuang wall painting images, the commonly used Gaussian elimination bundle adjustment (GEBA) (Triggs et al. 2000) is used for algorithm comparison and analysis.
Suppose that the corrections of exterior orientation parameters and additional parameters ( , ) T t E in Equation ( 15) is represented by Q .Then, Equation ( 14) can be written as: Since the number of unknowns of Q are usually far less than those of object point coordinates D , Q is usually solved first by eliminating D from Equation ( 16).The reduced normal equation of bundle adjustment is: where is called a schur complement matrix.
Due to the introduction of self-calibration parameters and quasiplanar constraints in bundle adjustment of Dunhuang wall painting images, H does not meet with regular banded structure.Therefore in GEBA, Equation ( 17) is solved by calculating the inverse matrix of H in the following calculation.
In the solving process of improved CGBA and GEBA, both coefficient matrix and constant matrix are required to store in memory.The space and time complexity of improved CGBA and GEBA are compared and analysed as follows: Suppose there are p images and q object points in bundle adjustment model, the number of unknowns of the exterior orientation parameters and the object point coordinates are 6 p and 3q , respectively.Besides, suppose that each object point is imaged on k images, l additional parameters are introduced into this model.The coefficient matrix of Equation ( 17) is a banded-bordered matrix with the border width l .The storage amount of improved CGBA is: Compared with improved CGBA, GEBA requires additional storage for A and B in Equation ( 17) as well during the calculation.Because A is formed by 3*3 matrix on the principal diagonal, we do not need to store the whole matrix A and B. On assumption of sparse storage the storage amount of GEBA is: As can be seen from Equation ( 18) and Equation ( 19), the space complexity of improved CGBA, which represented by storage amount, is little lower than that of GEBA.
The time complexity of improved CGBA and GEBA are measured by calculation amount, which are mainly determined by multiplication times.The difference between them is mainly caused by solving the normal equation and the reduced normal equation.In solving the normal equation, the calculation amount of CGBA at each iteration includes one multiplication of Equation ( 16) and vectors of the unknowns, as well as the calculation of preconditioning equation, which in total is: As illustrated in Section 3.2, n-time iterations at most are required in searching the solutions of normal equation by CGBA.The maximum iteration times is 6 3 p q l   , and the maximum calculation amount of CGBA is: In practice, the iteration times of improved CGBA with preconditioning based on incomplete Cholesky factorization are far less than n.The actual calculation amount of improved CGBA is about: By contrast, the calculation amount of GEBA include not only the construction and solution of Equation ( 17), but also the calculation of object point coordinates.Hence, the amount of GEBA is: 3 ( ) (162 27 27) (18 3 9) [( 6) ] By comparison between Equation ( 22) and Equation ( 23), the time complexity of improved CGBA, which is represented by calculation amount, is much lower than that of GEBA.

Experimental Data and Result Analysis
Four datasets of Dunhuang wall painting images are used for experiments of algorithm comparison between improved CGBA and GEBA.In Dataset I and Dataset II, eight self-calibration parameters of Brown model (Mikhail, Bethel and McGlone, 2001) are introduced into bundle adjustment, which include the corrections of three interior orientation parameters, three radial distortion parameters and two tangential distortion parameters.
In Dataset III and Dataset IV, four quasi-planar constraint parameters (Zhang, Hu and Huang, 2012) are introduced into bundle adjustment.
To verify the average actual accuracy of bundle adjustment, high-precision control and check points are evenly distributed and measured by electronic total station.Specifications of the four datasets are illustrated in details in Experimental results in Table 2 show that the RMS errors and average actual accuracy of the four datasets by improved CGBA and GEBA are relatively consistent.CGBA is proved to be correct and accurate.
In terms of operational efficiency, when the data amount is relatively small, the iteration times of improved CGBA is more than those of GEBA, and their operation times are nearly the same.With the data amount increase greatly, particularly when there are plenty of object points, the operation time of improved CGBA is far less than that of GEBA, which illustrates that improved CGBA is more suitable for solving large-scale bundle adjustment equations.
To verify the necessity of preconditioning method based on improved incomplete Cholesky factorization, the four datasets are solved respectively by improved CGBA and CGBA, i.e.CGBA with or without preconditioning.Their differences are mainly described by the condition number of matrix of the normal equation and the average iteration times.Experimental results comparison between improved CGBA and CGBA are shown in Table 3.
Since the order of equation is too large, it is difficult to calculate the condition number of normal equation directly.In this paper, the maximum and minimum eigenvalues of the coefficient matrix are firstly solved by the power method and the inverse power method, respectively.Then, the condition number is estimated by the two eigenvalues according to the positive definite symmetric property of coefficient matrix.The sign "--" represents the situations that the iteration of the normal equation is non-convergence or the solution to the normal equation fluctuates continually during iterations.

Methods
As shown in Table 3, the condition numbers of coefficient matrices of the four datasets by CGBA reach up to 10 22 at least, which is caused by the small overlap and quasi-planar imaging relationship of Dunhuang wall painting images.Under this condition, the rounding error of floating-point numbers arithmetic operation indicates higher sensitivity to the solution, the iteration convergence rate of bundle adjustment is slow, or even we cannot get a convergence solution.
Table 3 also shows that the preconditioning method based on improved incomplete Cholesky factorization can significantly decrease the condition number and the numerical sensitivity of normal equation, thus improving not only the equation structure, but also the computing performance of improved CGBA.The iteration times of improved CGBA are proved far less than n as theoretically illustrated in Section 4.1.

CONCLUSIONS
In this work, a novel conjugate gradient bundle adjustment with preconditioning based on improved incomplete Cholesky factorization has been proposed and applied for precise orthoimage generation and 3D reconstruction of Dunhuang wall paintings.Both theoretical analysis and experimental comparison demonstrate that this method can dramatically improve the structure of the normal equation and accelerate the convergence of bundle adjustment.This method is proved significant in practice and can be further modified and applied for other different 3D reconstruction applications.
a small value  .If k r   , then continue to doStep 4; Else, stop the iteration and output k X

Figure 1 .
Figure 1.Photographic Structure of Dunhuang Wall Paintings (Zhang et al. 2011) . The initial value of Y is replaced by 14) where D = the corrections of object point coordinates t = the corrections of exterior orientation parameters E = the corrections of additional parametersThis iterative computation method is to solve three a small value  , If k r   , then continue to do Step 4; Else stop the iteration and output k

Table 3 .
Experimental results of improved CGBA and CGBA Notes: a In the solving process of CGBA, the normal equation needs to calculate at each iteration.Average iteration times refers to the average iterative steps during iterations. b