Bundle Block Adjustment with 3D Natural Cubic Splines

Point-based methods undertaken by experienced human operators are very effective for traditional photogrammetric activities, but they are not appropriate in the autonomous environment of digital photogrammetry. To develop more reliable and accurate techniques, higher level objects with linear features accommodating elements other than points are alternatively adopted for aerial triangulation. Even though recent advanced algorithms provide accurate and reliable linear feature extraction, the use of such features that can consist of complex curve forms is more difficult than extracting a discrete set of points. Control points that are the initial input data, and break points that are end points of segmented curves, are readily obtained. Employment of high level features increases the feasibility of using geometric information and provides access to appropriate analytical solutions for advanced computer technology.

where M is the local structure matrix and α is a parameter so that 25 .
with I as an image. The local structure matrix M is: where 2 2 , and . y x g g C  A corner is detected when: where H thr is the threshold parameter on corner strength. The Harris operator searches points where variations in two orthogonal directions are large using the local autocorrelation function and provides good repeatability under varying rotation, scale, and illumination. The Förstner corner detector is also based on the covariance matrix for the gradient at a target point. Marr [9] proposes the zero-crossing edge point detector utilizing second order rather than first order directional derivatives. The maximum of first order derivatives indicates the location of an edge whereas it is the zero of second order derivatives that indicates an edge. Physical boundaries of objects of straight linear features and the candidate feature can be used in higher applications. A reason for developing curve features is that they will be prior to, and a fundamental aspect of, the next highest features such as surfaces, areas, and 3D volumes that consist of free-form linear features. Line-based photogrammetry is most suitable in the development of robust and accurate techniques for automation. If linear features are employed as control features, they provide advantages over points in the automation of aerial triangulation. Photogrammetry based on the manual measurement and identification of conjugate points is less reliable than line-based photogrammetry because it has the limitations of occlusion (visibility), ambiguity (repetitive patterns), and semantic information when considering the need for reliable and effective automation. The manual identification of corresponding entities within two images is crucial in the automation of point based photogrammetric tasks. No knowledge of the point-to-point correspondence is required in line-based photogrammetry. In addition, point features do not carry information about the scene whereas linear features contain the semantic information related to real object features. Additional information concerning linear features can increase the redundancy of the point system.

Literature Review
A review of related works begins with those using methods of pose estimation in imagery based on linear features that appear in most man-made objects such as buildings, roads, and parking lots. Over the years, a number of researchers in photogrammetry and computer vision have used line instead of point features; for example, Masry [25], Heikkila [26], Kubik [27], Petsa and Patias [28], Gülch [29], Wiman and Axelsson [30], Chen and Shibasaki [31], Habib [32], Heuvel [33], Tommaselli [34], Vosselman and Veldhuis [35], Förstner [36], Smith and Park [37], Schenk [38], Tangelder et al. [39], and Parian and Gruen [40]. Mulawa and Mikhail [41] originally proved the feasibility of linear features for close-range photogrammetric applications such as space intersection and resection, and relative and absolute orientation. This was the first step in employing linear feature-based methods in close-range photogrammetric applications. Mulawa [42] later developed linear feature-based methods for different sensors.
Whereas straight linear features and conic sections can be represented as unique mathematical expressions, free-form lines in nature cannot be described by algebraic equations. Hence, Mikhail and Weerawong [43] used splines and polylines to represent free-form lines as analytical expressions. Tommaselli and Tozzi [44] proposed that the accuracy of the straight line parameter be a subpixel with the representation of four degrees of freedom in an infinite line. Many researchers in photogrammetry have described straight lines as infinite lines using minimal representation to reduce unknown parameters. The main consideration in straight line expression is in the singularities. Habib et al. [45] made a bundle block adjustment using a 3D point set lying on control linear features instead of traditional control points. EOPs were reconstructed hierarchically employing automatic single photo resection (SPR).
Habib et al. [46] summarized linear features extracted from a mobile mapping system, a GIS database, and maps for various photogrammetric applications such as SPR, triangulation, digital camera calibration, image matching, 3D reconstruction, image to image registration, and surface to surface registration. In their work, matched linear feature primitives were utilized in space intersection for the reconstruction of object space features, and linear features in the object space were used as control features in triangulation and digital camera calibration.
Mikhail [47] and Habib et al. [48] accomplished the geometrical modeling and the perspective transformation of linear features within a triangulation process. Linear features were used to recover relative orientation parameters. Habib et al. proposed a free-form line in object space by a sequence of 3D points along the object space line.
Schenk [49] extended the concept of aerial triangulation from point features to linear features. The line equation of six dependent parameters replaced the point-based collinearity equation: where a real variable is t, the start point (X A ,Y A ,Z A ), and the direction vector (a,b,c).
The traditional point-based collinearity equation was extended to line features: 33 with x p ,y p as photo coordinates, f the focal length, X C ,Y C ,Z C the camera perspective center, and r ij the elements of the 3D orthogonal rotation matrix. The extended collinearity equation with six parameters was derived as the line expression of four parameters (ф,θ,x 0 ,y 0 ) because a 3D straight line has only four independent parameters. Two constraints are required to solve a common form of the 3D straight equations using six parameters determined by two vectors: where z is a real variable. The advantage of the 3D straight line using four independent parameters is that it reduces the computation and time complexity in processes such as bundle block adjustment. The collinearity equation as the straight line function of four parameters was developed: 33 where X, Y, and Z were defined in equation (7).
Zalmanson [50] updated EOPs using the correspondence between the parametric control free-form line in object space and the projected 2D free-form line in image space. The hierarchical approach, the modified iteratively close point (ICP) method, was developed to estimate curve parameters. The ray lies on the free-form line whose parametric equation represented by one parameter follows. Besl and McKay [51] employed the ICP algorithm to solve a matching problem of point sets, free-form curves, surfaces, and terrain models in 2D and 3D space. In their work, an ICP algorithm was executed without prior knowledge of the correspondence between points. The ICP method affected the Zalmanson's dissertation on the development of the recovery of EOPs using 3D free-form lines in photogrammetry. Euclidean 3D transformation was then employed in a search for the closest entity in the geometric data set. Rabbani et al. [52] utilized the ICP method in the registration of Lidar point clouds to divide them into four categories (spheres, planes, cylinders, and tori) by direct and indirect methods. The parametric curve Γ(t) = [X(t) Y(t) Z(t)] T was obtained by minimizing the Euclidian distance between two parametric curves: with two independent variables l and t as in (11).
Akav et al. [53] employed planar free-form curves for aerial triangulation with the ICP method. Because the effect of the Z parameter as compared with that of X or Y was large in a normal plane equation aX + bY + cZ = 1, a different plane representation was developed to avoid numerical problems: with  the angle from the XY plane,  the angle around the Z axis, n the unit vector of plane normal, and D the distance between the plane and the origin. Five relative orientation parameters and three planar parameters were obtained by using the homography mapping system, which searched for the conjugate point in an image corresponding to a point in the other image.
Lin [54] proposed the method of autonomous recovery of exterior orientation parameters by an extension of the traditional point-based modified iterated Hough transform (MIHT) to the 3D free-form linear feature-based MIHT. Straight polylines were generalized for matching primitives in the pose estimation because the mathematical representation of straight lines is much clearer than the algebraic expression of conic sections and splines.
Gruen and Akca [55] matched 3D curves whose forms were defined by a cubic spline using matching least squares. Subpixels were localized by the matching, and the quality of the localization was decided by the geometry of image patches. Two free-form lines were defined as: where   3  2  1  0  3  2  1  0   ,  ,  ,  ,  ,  ,  ,  ], The Taylor expansion was employed to adopt the Gauss-Markov model:

3D Natural Cubic Splines
The choice of the right feature model is important in the development of a feature-based approach because an ambiguous feature representation leads to unstable adjustment. A spline is a piecewise polynomial function in the n of vector graphics. Splines are widely used for data fitting in computer science because of the resultant simplicity in curve reconstruction. Complex figures are well approximated through curve fitting and a spline lends strength to the accuracy evaluation, data interpolation, and curve smoothing. One of the important properties of a spline is that it can easily be morphed. A spline represents a 2D or 3D continuous line within a sequence of pixels and segmentation. The relationship between pixels and lines is applied to a bundle block adjustment or a functional representation. A spline of degree 0 is the simplest spline, a linear spline has degree 1, a quadratic spline has degree 2, and a common natural cubic spline has degree 3 with continuity C 2 . The geometrical meaning of continuity C 2 is that the first and second derivatives are proportional at joint points and the parametric importance of continuity C 2 is that the first and second derivatives are equal at connected points.
The number of break points that are the determination of a set of piecewise cubic functions varies depending upon the spline parameters. A natural cubic degree guarantees a second-order continuity, which means that the first and second order derivatives of two consecutive natural cubic splines are continuous at the break point. The intervals for a natural cubic spline do not need to be the same as the distance of every two consecutive data points. The best intervals are chosen by a least squares method. In general, the total number of break points is less than that of original input points. The algorithm of a natural cubic spline is as below.
Generate the break point (control point) set for the spline of the original input data. Calculate the maximum distance between the approximated spline and the original input data while the maximum distance > the threshold of the maximum distance. Add the break point to the break point set at the location of the maximum distance. Compare the maximum distance with the threshold.
A larger threshold makes for more break points with a more accurate spline to the original input data. N piecewise cubic polynomial functions between two adjacent break points are defined from the N + 1 break points. There is a separate cubic polynomial for each segment with its own coefficients.
The strength of this approach is that segmented lines represent a free-form line with analytical parameters. The number of break points is reduced and the input error should be absorbed by a mathematical model, especially in the expression of points on a straight line. A natural cubic spline is a data-independent curve fitting. The disadvantage is that the whole curve shape depends on all of the passing points, and changing any one of these alters the entire curve.
The correspondence between the 3D curve in the object space coordinate system and its projected 2D curve in the image coordinate system is implemented using an accommodating natural cubic spline curve feature because of its boundary conditions that retain zero second derivatives at the end points. A natural cubic spline is composed of a sequence of cubic polynomial segments as in Figure 1 with x 0 ,x 1 ,…,x n as the n + 1 control points and X 0 ,X 1 ,…,X n-1 as the ground coordinates of n segments.

Extended Collinearity Equation Model for Splines
Collinearity equations are the commonly used condition equations to determine relative orientation. The space intersection calculates a point location in object space using the projection ray intersection from two or more images, and the space resection determines the coordinates of a point on an image and the EOPs with respect to the object space coordinate system. The space intersection and the space resection are the fundamental operations in photogrammetry for further processes such as triangulation. The basic concept of the collinearity equation is that all points on the image, the perspective center, and the corresponding point in the object space are all on a straight line. The relationship between the image and object coordinate systems is expressed by three position parameters and three orientation parameters. Collinearity equations play an important role in photogrammetry because each control point in object space produces two collinearity equations for every photograph in which the point appears. If m points appear in n images, then 2mn collinearity equations can be employed in the bundle block adjustment. The extended collinearity equations relating a natural cubic spline in object space with ground coordinates (X i (t),Y i (t),Z i (t)) with image space having photo coordinates (x pi ,y pi ) are seen as (16). A natural cubic spline allows the utilization of the collinearity model for expressing orientation parameters and curve parameters as below: 33 with (x pi ,y pi ) as the photo coordinate, f the focal length, X C ,Y C ,Z C the camera perspective center, and r ij the elements of the 3D orthogonal rotation matrix T R by the angular elements (,,) of EOPs. The extended collinearity equations can be written as follows: To recover the 3D natural cubic spline parameters and the exterior orientation parameters in a bundle block adjustment, a nonlinear mathematical model of the extended collinearity equation is differentiated. The models of exterior orientation recovery are classified into linear and nonlinear Sensors 2009, 9 9639 methods. Whereas linear methods decrease the computation load, the accuracy and reliability of linear algorithms are not excellent. Nonlinear methods are more accurate and predictable. However, nonlinear methods require initial estimates and they increase the computational complexity. The relationship between a point in image space and a corresponding point in object space is established by the extended collinearity equation. Prior knowledge of the correspondences between individual points in the 3D object space and their projected features in the 2D image space is not required in extended collinearity equations with 3D natural splines. One point on a cubic spline has 19 parameters (X C ,Y C ,Z C ,,,,a 0 ,a 1 ,a 2 ,a 3 ,b 0 ,b 1 ,b 2 ,b 3 ,c 0 ,c 1 ,c 2 ,c 3 ,t). The differentials of (17) are derived by (18): with differentials of du, dv, and dw (19).

Arc-Length Parameterization of 3D Natural Cubic Splines
The assumption made in bundle block adjustment by the Gauss-Markov model is that all the estimated parameters are uncorrelated. Hence, the design matrix of the adjustment must be full rank, nonsingular, and normal. However, because the spline parameters are not independent of their location parameters, additional observations are required to obtain parameter estimations. In a point-based approach, the point location relationship between image and object space is established for the pose estimation to include the fundamental camera position and orientation, the remote sensing, and the computer vision. The coordinates of a point are necessary for the space intersection and resection. To remove any rank deficiency caused by datum defects in point-based photogrammetry, some constraints are adopted to estimate the unknown parameters. The most common constraints are coplanarity, symmetry, perpendicularity, and parallelism. The minimum number of constraints is equal to the rank deficiency of the system. Inner constraints are often used in a photogrammetric network, which can be applied to both the object features and the camera orientation parameters. Angle or distance condition equations provide information on the relativity between observations in object space and points in image space. Absolute information can be obtained from the fixed control points.
In this research, an arc-length parameterization is applied as an additional condition equation to solve the rank deficient problem in extended collinearity equations using 3D natural cubic splines. The concept of differentiable parameterization is that the arc length of a curve can be divided into minute pieces and these can be summed such that each piece will be approximately linear. The sum of the squares of derivatives is the same with a velocity vector because a parametric curve can be considered as a point trajectory. A velocity vector describes the path of a curve and movement characteristics. If a particle on a curve moves at a constant rate, the curve is parameterized by the arc length. While the extended collinearity equation provides the only information, curves have additional geometric constraints such as arc length, tangent of location, and curvature. These support space resection under the assumption of properly accounting for additional independent observations in both the image and object space. Arc length in object space is determined by a geometric integration using a construction from the differentiable parameterization of a spline. Arc length in image space is calculated by a geometric integration of a construction from the differentiable parameterization of the photo coordinates derived from a spline in the object space: where f is the focal length and: Because the problem of the arc-length parameterization of splines has no analytical solution, several numerical approximations of reparameterization techniques for splines or other curve representations have been developed. While most curves are not parameterized for arc length, the arc length of a B-spline can be reparameterized by adjusting its knots. Wang et al. [56] approximated the parameterized arc length of spline curves by generating a new curve that accurately approximated the original spline curve to reduce the computation complexity of the arc-length parameterization. They showed that the approximation of the arc-length parameterization works well in a variety of real-time applications including driving simulations.
Guenter and Parent [57] employed a hierarchical approach algorithm to develop a linear search arc-length subdivision table for parameterized curves to reduce the arc-length computation time. A table of the correspondence between parameter t and the arc length can be established to accelerate the arc-length computation. After dividing the parameter range into intervals, the arc length of each interval is computed for mapping parameters to the arc length. A reference table for various intervals of arc length can be developed. Another method of arc-length approximation is to use explicit functions such as the Bé zier curve, which has advantages in fast function evaluations. Adaptive Gaussian integrations employ a recursive method that starts from a few samples and adds on more as necessary. Adaptive Gaussian integration also uses a table that maps curves or spline parameter values according to the arc-length values.
Nasri et al. [58] proposed an arc-length approximation method of circles and piecewise circular splines generated by control polygons or points using a recursive subdivision algorithm. While B-splines have various tangents over the curve depending upon the arc-length parameterization, circular splines have constant tangents whose vectors are useful in arc-length computation.
Simpson's rule is the numerical approximation of definite integrals. The geometric integration of the arc length in the image space can be calculated by this rule as follows: 2  20  1  19  3  18  2  17  1  16   0  15  3  14  2  13  1  12  0  11  3 being the approximate parameter using the following:  0  0  0  0  0   ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  , , , and with a e being the stochastic error of the arc length between two locations with zero expectation. A 1 ,…A 20 denote the partial derivatives of the arc-length parameterization of a 3D natural cubic spline.

Model Integration
The objective of bundle block adjustment is twofold, namely to calculate the exterior orientation parameters of a block of images and also the coordinates of the ground features in object space. In the determination of orientation parameters, additional interior conditions such as lens distortion, atmospheric refraction, and principal point offset can be obtained by self-calibration. In general, orientation parameters are determined by bundle block adjustment using a large number of control points. This establishment of control points, however, means expensive fieldwork, so an economical and accurate adjustment method is required. Linear features have several advantages to complement points in that they are useful for higher level tasks and they are easily extracted in man-made Sensors 2009, 9 9643 environments. The line photogrammetric bundle adjustment in this research aims at the estimation of exterior orientation parameters and 3D natural cubic spline parameters using the correspondence between splines in object space and spline observations of multiple images in image space. Nonlinear functions of orientation parameters, spline parameters, and spline location parameters are represented by extended collinearity and arc-length parameterization equations. Five observation equations are produced by each two points, and these are four extended collinearity equations (21) and one arc-length parameterization equation (24). An integrated model provides not only for the recovery of the image orientation parameters but also enables surface reconstruction using 3D curves. Of course, as the equation system of the integrated model has seven datum defects, control information about the coordinate system is required to obtain parameters. This is a step toward higher level vision tasks such as object recognition and surface reconstruction. In the case of straight lines and conic sections, tangents are additional observations in the integrated model. Conic sections, like points, provide good mathematical constraints because such sections provide nonsingular second degree equations. Such equations provide information for reconstruction and transformation and conic sections are divided by the eccentricity e. Because such sections can adopt more constraints than points and straight line features, they are useful for close range photogrammetric applications. In addition, conic sections have strength in correspondence establishment between 3D sections in object space and their counterpart features in 2D projected image space.
Ji et al. [59] employed conic sections for the recovery of EOPs, and Heikkila used them for camera calibration. A Hough transformation reduces the time complexity of conic section extraction using five parameter spaces for a SPR, camera calibration, and triangulation.
Parameters are linearized in the previous sections and the Gauss-Markov model is employed for the unknown parameter estimation. The equation system of the integrated model is described as:  dc  dc  dc  dc  db  db  db  db  da  da  da  da   3  2  1  0  3  2  1  0  3 , m as the number of images, n the number of points on a spline segment, k the kth image, and i the ith spline segment. Because the equation system of the integrated model has seven datum defects, the control information for the coordinate system is required to obtain seven transformation parameters. In a general photogrammetric network, the rank deficiency referred to as datum defects is seven. Estimates of the unknown parameters are obtained by the least squares solution, which minimizes the sum of squared deviations. A nonlinear least squares system is required in a conventional nonlinear photogrammetric solution to obtain orientation parameters. are considered as stochastic constraints, the reduction of the normal equation matrix can be applied. Control information is implemented as stochastic constraints in a bundle block adjustment. The distribution and quality of control features depend on the number and the density of control features, the number of tie features, and the degree of overlap of the tie features. If adding stochastic constraints removes the rank deficiency of the Gauss-Markov model, bundle adjustment can be implemented employing only the extended collinearity equations for the 3D natural cubic splines. Fixed exterior orientation parameters, control splines, or control spline location parameters can be stochastic constraints.

Evaluation of Bundle Block Adjustment
Bundle block adjustment must be followed by an evaluation postadjustment analysis to check the suitability of project specifications and requirements. Iteratively reweighed least squares and least median of squares are the appropriate implementation of a statistical evaluation that removes poor observations. The important element affecting bundle block adjustment is the geometry of aerial images. Generally, the previous flight plan is adopted to obtain suitable results. A simulation bundle block adjustment is implemented before employing a flight plan within the new project design because such a simulation can reduce the effect of error measurements.
A qualitative evaluation that allows the operator to recognize the adjustment characteristics is often used after bundle block adjustment. The sizes of the residuals in images are drawn for the evaluation. The image residuals can be points or long lines and if all image residuals have the same orientation, then the image has a systematic error such as atmospheric refraction or an orientation parameter error. In addition, a lack of flatness in the focal plane may cause systematic errors in the image space, which affects the accuracy of a bundle block adjustment. Distortions are different from one location to another in the entire image space. The topographic measurement of the focal plane can correct the lack of focal plane flatness. Image coordinate errors are correlated in the case of systematic image errors. A poor measurement can result in an indicated opposite residual direction or an exaggerated residual.
The three main elements in the statistical evaluation of bundle block adjustments are precision, accuracy, and reliability. Precision is calculated employing parameter variances and covariances, because a small variance indicates that the estimated values have a small range and a large variance means that the estimated values are not calculated properly. The range of the parameter variance is from zero, in the case of error free parameters, to infinity, in the case of completely unknown parameters. A dispersion matrix may contain diagonal elements that are parameter variances. These and any off-diagonal elements are covariances between two parameters. Accuracy can be verified using check points that are not contained in bundle block adjustment like control points. Reliability can be confirmed from other redundant observations. The extended collinearity equations are a mathematical model for bundle block adjustment. The mathematical model consists of both functional and stochastic models. The functional one represents the geometrical properties and the stochastic one describes the statistical properties. Repeated measurements at the same location in the image space are represented with respect to the functional model and the redundant observations of image locations in the image space are expressed with respect to the stochastic model. While the Gauss-Markov model uses indirect observations, condition equations such as coordinate transformations and the coplanarity condition can be employed in the adjustment.
The Gauss-Markov model and the condition equation can be combined into the Gauss-Helmert model. In addition, functional constraints such as points having the same height or straight railroad segments can be added into the block adjustment.
The difference between condition and constraint equations is that condition equations consist of observations and parameters, and constraint equations consist of only parameters. With the advance of technology, the photogrammetrical input data has increased so adequate formulation of adjustment is required. All the variables are involved in the mathematical equations and the weight matrix of the variables changes from zero to infinity depending upon the variances. Variables with near to zero weight are considered as unknown parameters and variables with close to infinite weight are considered as constants. Most actual observations exist between the two boundary cases. Assessment by postadjustment analysis is important in photogrammetry to evaluate the results. One of the assessment methods is to compare the estimated variance with the two-tailed confidence interval based on the normal distribution. The two-tailed confidence interval is computed by a reference variance 2 0  with 2  distribution as: where r is degrees of freedom and  is a confidence coefficient (or a confidence level). If 2 0  has a value outside of the interval, we can assume that the mathematical model of adjustment is incorrect through the wrong formulation or linearization, blunders, or systematic errors.

Pose Estimation with an ICP Algorithm
In the previous spline segment case, the correspondence between spline segments in the image and the object space was assumed. In the present consideration, it is not known which image points belong to which spline segment. The ICP algorithm can be utilized for the recovery of EOPs because the initial estimated parameters of the relative pose can be obtained from the orientation data for general photogrammetric tasks. The original ICP algorithm steps are as follows. The closest point operators search the associate point using the nearest neighboring algorithm and then the transformation parameters are estimated using a mean square cost function. The point is transformed by the estimated parameters and this step is iteratively established towards convergence into a local minimum of the mean square distance. The transformation, which includes translation and rotation between two clouds of points, is estimated iteratively towards convergence into a global minimum. In other words, the iterative calculation of the mean square errors is terminated when a local minimum falls below a predefined threshold. A small global minimum or a fluctuated curve requires more memory-intensive and time-consuming computation. In every iteration step, a local minimum is calculated with varying transformation parameters, but convergence into a global minimum with the correct transformation parameters is not always the result.
By the definition of a natural cubic spline, each parametric equation of a spline segment )) ( ( t S i can be expressed as: as the object space coordinates and , as the coefficients of the ith spline segment. The ray from the perspective center (X C ,Y C ,Z C ) to the image point (x p ,y p ,-f) is: where: EOPs at the kth iteration.
A point on the ray searches the closest to a natural cubic spline by minimizing the following target function for every spline segment. Transformation parameters related to an image point and its closest spline segment can be established using the least squares method: The global minimum of ) Substituting (28) and (29) into (31) and taking the derivatives with respect to l and t leads to: Convergence into a global minimum does not exist because equation (32) is not a linear system in l and t. The relationship between an image space point and its corresponding spline segment cannot be established with the minimization method.

Experiments and Results
This section demonstrates the feasibility and the performance of the proposed model for the acquisition of spline parameters, spline location parameters, and image orientation parameters based on control and tie splines in the object space within the simulated and real data sets. In general photogrammetric tasks, the correspondence between image edge features must be established either automatically or manually, but in this study correspondence between image edge features is not required. In a series of six experiments with the synthetic data set, the first test recovers spline parameters and spline location parameters in an error free EOPs case. The second test recovers the partial spline parameters related to the spline shape. The third procedure estimates the spline location parameters with error free EOPs. The fourth step calculates EOPs and spline location parameters, followed by the fifth step that estimates EOPs with full controlled splines in which the parametric curves used as control features are assumed to be error free. In the last experiment, EOPs and tie spline parameters are obtained using the control spline.
Object space knowledge concerning splines, their relationships, and the orientation information of images can be considered as control information. Spline parameters in a partial control spline or orientation parameters can be considered as stochastic constraints in the integrated adjustment model. The starting point of a spline is considered to be a known parameter in the partial control spline in which a 0 ,b 0 , and c 0 of the X, Y, and Z coordinates of a spline are known. The number of unknowns is displayed in Table 1 and Figure 3, where n is the number of points in the object space, t shows the number of spline location parameters, and m represents the number of overlapped images in the target area.

Unknown EOP
Partial control spline 6m + 9(n-1) + t Full control spline 6m + t Four points on a spline segment in one image are the only independent observations so additional points on the same segment do not provide nonredundant information to reduce the overall deficiency of the EOP and spline parameter recovery. To verify the information content of an image spline, we demonstrate that any five points on a spline segment generate a dependent set of extended collinearity equations. Any combination of four points yielding eight collinearity equations are independent observations, but five points bearing 10 collinearity equations produce a dependent set of observations Sensors 2009, 9 9649 related to the correspondence between a natural cubic spline in the image and the object space. More than four point observations on an image spline segment increase the redundancy related to the accuracy but do not decrease the overall rank deficiency of the proposed adjustment system. In the same fashion, the case using a polynomial of degree 2 can be implemented. Three points on a quadratic polynomial curve in one image are the only independent sets, so additional points on the same curve segment are a dependent observation. More than the independent point observations on a polynomial increase the redundancy related to the accuracy, but they do not provide nonredundant information.
The amount of information carried by a natural cubic spline can be calculated with the redundancy budget. Every spline segment has 12 parameters and every point measured on a spline segment contributes one additional parameter. Let n be the number of points measured on one spline segment in the image space and m be the number of images that contain a tie spline. 2nm collinearity equations and m (n − 1), the arc-length parameterizations, are equations and 12 (the number of one spline segment parameters) + nm (the number of spline location parameters) are unknowns. The redundancy is 2nm − m − 12 for one spline segment, so that if two images (m = 2) are used for bundle block adjustment, the redundancy is 4n − 14. Four points are required to determine spline and spline location parameters, in which case one spline segment and one degree of freedom to the overall redundancy budget is solved by each point measurement with the extended collinearity equation. Arc-length parameterization also contributes one degree of freedom to the overall redundancy budget. The fifth point does not provide additional information to reduce the overall deficiency but only strengthens the spline parameters. This means it increases the overall precision of the estimated parameters.
This fact shows the advantage of adopting splines in which the number of degrees of freedom is four because in straight tie lines only two points per line are independent. Independent information, the number of degrees of freedom of a straight line, is two from two points or a point with its tangent direction. A redundancy is r = 2m − 4 with a line expression of four parameters because there are 2 nm collinearity equations and the unknowns are 4 + nm [49]. Only two points (n = 2) are available to determine four line parameters with two images (m = 2) so at least three images must contain a tie line. The information content of t tie lines on m images is t (2m − 4). One straight line adds two degrees of freedom to the redundancy budget and at least three lines are required in the space resection. An additional point on a straight line does not provide additional information to reduce the rank deficiency of the recovery of EOPs but only contributes image line coefficients. If spline location parameters or spline parameters enter the integrated adjustment model through stochastic constraints, employing extended collinearity equations is enough to solve the system without the arc-length parameterization.
The redundancy budget of a tie point is r = 2m − 3 so tie points provide one more independent equation than the tie lines. However, using tie points requires a semiautomatic matching procedure to identify the tie points on all the images, and using linear features provides a more reliable and accurate basis for object recognition, pose determination, and other higher photogrammetric activities than using point features.

Synthetic Data Description
To evaluate the new bundle block adjustment model using natural cubic splines, an analysis of the sensitivity and robustness of the model is required. The model suitability can be verified by using the Sensors 2009, 9 9650 estimated parameters with a dispersion matrix that includes standard deviations and correlations. The accuracy of bundle block adjustment is determined by the geometry of a complete block of images and the quality of the position and attitude information of a camera. A novel approach is a simulation of the bundle block adjustment. This is required prior to an actual experiment with real data in order to evaluate the performance of the proposed algorithms. Such a simulation can control the measurement errors to minimize random noise affecting the overall geometry of a block. Individual observations are generated based on the general situation of bundle block adjustment in order to estimate the properties of the proposed algorithms. A simulation allows adjustment for geometric problems or conditions with various experiments. A spline is derived via three ground control points (3232, 4261, 18), (3335, 4343, 52), and (3373,4387,34). Several factors that affect the estimates of exterior orientation parameters, spline parameters, and spline location parameters are discerned using the proposed bundle block adjustment model together with both the simulated image and the real image blocks.

Experiments with Error Free EOPs
Spline parameters and spline location parameters are dependent upon various controls, and the unknowns can be obtained by a combined model of extended collinearity equations and the arc-length parameterization equations of splines. Splines in the object space are considered as tie lines in the same fashion as tie points in a conventional bundle block adjustment. Data on the exterior orientation parameters is considered as control information in this experiment. A well-known fact in employing the least squares system is that good initial estimates of true values make the system swiftly convergent towards the correct solution.
Normally distributed random noise is added to points in the image space coordinate system in all the experiments. This has a zero mean and  = 5m standard deviation. Generally, the larger the noise level the more accurate are the approximations required to achieve the ultimate convergence of the results. A worst case scenario for estimation is that the large noise level causes the proposed model not to converge towards the specific estimates because the convergence radius is then proportional to the noise level. The parameter estimation is sensitive to the noise of the image measurement. Error propagation related to the noise in image space observation is one of the most important elements in the estimation theory. The proposed bundle block adjustment can be evaluated statistically using the variances and the covariances of parameters because a small variance indicates that the estimated values have a small range and a large variance means that the estimates are not properly calculated. The range of parameter variance is from zero in the case of error free parameters to infinity with completely unknown parameters. The result of one spline segment is expressed in Table 3 with  0 as the initial values and ˆ as the estimates. The estimated spline and spline location parameters along with their standard deviations are established without the knowledge of the point-to-point correspondence. Table 3. Spline parameter and spline location parameter recovery. Image 1  Image 2  Image 3   If no random noise is added to image points, the estimates converge to the true values. The quality of initial estimates is important in the least squares system because it determines the iteration number of the system and the accuracy of the convergence. The assumption is that two points on one spline segment are measured in each image so the total number of equations is 2  6 (the number of images)  2 (the number of points) + 6 (the number of the arc length), and the total number of unknowns is 12 (the number of spline parameters) + 12 (the number of spline location parameters). The redundancy (=the number of equations − the number of parameters), that is, the degrees of freedom, is six. While some of the geometric constraints such as slope and distance observations are dependent on the extended collinearity equations using splines, other constraints such as slope and arc length increase the nonredundant information in the adjustment to reduce the overall rank deficiency of the system.

Spline location parameters
The coplanarity approach is another mathematical model of the perspective relationship between the image and the object space features. The projection plane defined by the perspective center in the image space and the plane including the straight line in the object space are identical. Because the coplanarity condition is only for straight lines, the coplanarity approach cannot be extended to curves. Object space knowledge about the starting point of a spline can be employed in bundle block adjustment. Because the control information about a starting point is available for only three parameters of a total of 12 unknown parameters to a spline, a spline with control information about a starting point is called a partial control spline. Three spline parameters related to the starting point of a spline are set to stochastic constraints and the result is seen in Table 4. The total number of equations is 2  6 (the number of images)  2 (the number of points) + 6 (the number of the arc length) = 30, and the total number of unknowns is 9 (the number of partial spline parameters) + 12 (the number of spline location parameters) = 21 so the redundancy is nine. A convergence of partial spline and spline location parameters has been archived with a partial control spline. Table 4. Partial spline parameter and spline location parameter recovery.
Spline location parameters Image 1 Image 2 Image 3  In the next experiment, spline location parameters are estimated with known EOPs and a full control spline. Because spline parameters and spline location parameters are dependent upon other parameters, the unknowns can be obtained from the model of an observation equation with stochastic constraints. In this experiment, spline parameters are set to stochastic constraints and the result is seen in Table 5. The total number of equations is 2  6 (the number of images)  3 (the number of points) = 36, and the total number of unknowns is 18 (the number of spline location parameters) so the redundancy is 18. Because spline location parameters are independent of each other, the arc-length parameterization is not required. The result indicates that a convergence of spline location parameters has been achieved with fixed spline parameters considered as stochastic constraints. The proposed model is robust with respect to the initial approximations of spline parameters. The uncertain information related to the representation of a natural cubic spline is described in the dispersion matrix.

Recovery of EOPs and Spline Parameters
The object space knowledge of splines is available to recover the exterior orientation parameters in a bundle block adjustment. Control spline and partial control spline approaches are applied to verify the feasibility of using control information with splines. In both cases, equations of the arc-length parameterization are not necessary if we have enough equations to solve the system because spline parameters are independent of each other. In the experiment for a full control spline, the total number of equations is 2  6 (the number of images)  4 (the number of points) + 3 (the number of arc lengths)  6 (the number of images) = 66, and the total number of unknowns is 36 (the number of EOPs) + 24 (the number of spline location parameters) = 60. The redundancy is six. In the case of the partial control spline with one spline segment, the total number of equations is 2  6 (the number of images)  4 (the number of points) + 3 (the number of arc lengths)  6 (the number of images) = 66, and the total number of unknowns is 36 (the number of EOPs) + 9 (the number of partial spline parameters) + 24 (the number of spline location parameters) = 69. Thus, one more segment is required to solve the underdetermined system. The total number of equations using two spline segments is 2  6 (the number of images)  4 (the number of points)  2 (the number of spline segments) + 3 (the number of arc lengths)  6 (the number of images)  2 (the number of spline segments) = 132, and the total number of unknowns is 36 (the number of EOPs) + 9 (the number of partial spline parameters)  2 (the number of spline segments) + 24 (the number of spline location parameters)  2 (the number of spline segments) = 102. The redundancy is 30. A convergence of the EOPs of an image block and the spline parameters has been achieved in both experiments. Table 6 expresses the convergence achievement of EOPs and spline location parameters. The correlation coefficient between parameter X C and  is high (  1) in the dispersion matrix, that is, two parameters are highly correlated among the EOPs. The correlation coefficient between parameters Y C and  is approximately 0.85. In general, the correlation coefficient between parameters X C and  is higher than between parameters Y C and .
Because a control spline provides the object space information about the coordinate system having datum defects of seven, tie spline parameters and EOPs can be recovered simultaneously. In the experiment of combined splines, the total number of equations is 2  6 (the number of images)  3 (the number of points)  2 (the number of splines) + 12 (the number of arc lengths)  2 (the number of splines) = 96, and the total number of unknowns is 36 (the number of EOPs) + 12 (the number of tie spline parameters) + 18 (the number of tie spline location parameters) + 18 (the number of control spline location parameters) = 84.  Knowledge of object space information about a spline referred to as a full control spline is available prior to aerial triangulation. A control spline is considered to be a stochastic constraint in the proposed Sensors 2009, 9 9656 adjustment model and the representation of a control spline is the same as that of a tie spline. The result for combined splines that demonstrates the feasibility of using tie splines and control splines for bundle block adjustment is illustrated in Table 7.
Iteration with an incorrect spline segment in which a spline in the image space does not lie on the projection of a 3D spline in the object space results in a divergence of the system. A control spline is taken to be error free, but in reality this assumption is not correct. The accuracy of control splines is propagated into the proposed bundle block adjustment algorithm, but initial data such as a GIS database, maps, or orthophotos cannot be without error.

Tests with Real Data
In this section, actual experiments with real data are undertaken to verify the feasibility of the proposed bundle block adjustment algorithm using splines for the recovery of EOPs and spline parameters. Medium scale aerial images covering the area of Jakobshavn Isbrae in West Greenland are employed for this study. The aerial photographs were obtained by Kort and Matrikelstyrelsen (KMS: Danish National Survey and Cadastre) in 1985. KMS established aerial triangulation using GPS ground control points with a  1 pixel root mean square error under favorable circumstances and images were oriented to the WGS84 reference frame. Technical information on the aerial images is described in Table 8. The diapositive films were scanned with a RasterMaster photogrammetric precision scanner, which has a maximum image resolution of 12 μm and a scan dimension of 23 cm  23 cm to obtain digital images for a softcopy workstation as seen in Figure 6. The first experiment is the recovery of spline parameters with known EOPs obtained by manual operation using a softcopy workstation. A spline consists of four parts and the second segment parameters are recovered. The total number of equations is 2  3 (the number of images)  3 (the number of points) + 2 (the number of arc lengths)  3 (the number of images) = 24, and the total number of unknowns is 12 (the number of spline parameters) + 9 (the number of spline location parameters) = 21 so the redundancy is three. Table 9 shows the convergence achievement of spline and spline location parameters. Estimation of spline parameters including their location parameters is established by the relationship between splines in the object space and their projection in the image space without the knowledge of the point-to-point correspondence. Because bundle block adjustment using splines does not require conjugate points generated by point-to-point correspondence knowledge, a more robust and flexible matching algorithm can be adopted. Table 10 shows the available object space information without knowledge of the point-to-point correspondence (the full control spline). All locations are assumed as lying on the second spline segment and the second spline segment as calculated from the softcopy workstation is used as control information.  Even though edge detectors are often used in digital photogrammetry and remote sensing software, the control points are extracted manually because edge detection is not our main goal. Among the three segments, the second spline segment is used for the EOP recovery. The information of the control spline is obtained by a manual operation using the softcopy workstation with an estimated accuracy of ±1 pixel. The convergence radius of the proposed iterative algorithm is proportional to the estimated accuracy level. The image coordinate system is converted into the photo coordinate system using the interior orientation parameters from KMS. The association between a point on a 3D spline segment and a point on a 2D image is not established in this study. Of course, 3D spline measurement in the stereo model using the softcopy workstation cannot be without error so the accuracy of the control spline is propagated into the recovery of EOPs. The result is illustrated in Table 11. The spline control information is utilized as stochastic constraints in the adjustment model. Because adding these constraints removes the rank deficiency of the Gauss-Markov model corresponding to spline parameters that are dependent upon spline location parameters, a bundle block adjustment can be made using only the extended collinearity equations for natural cubic splines.

Conclusions
In this paper, traditional least squares of a bundle block adjustment process have been augmented by support splines instead of conventional point features. Estimation of EOPs and spline parameters including location parameters is established by the relationship between splines in the object space and their projection into the image space without any knowledge of the point-to-point correspondence. Because bundle block adjustment using splines does not require conjugate points generated by the point-to-point correspondence knowledge, a more reliable and flexible matching algorithm can be adopted. Point-based aerial triangulation with experienced human operators is effective for traditional photogrammetric activities but is not appropriate within the autonomous environment of digital photogrammetry. Feature-based aerial triangulation is suitable for the development of reliable and accurate automation techniques. If linear features are employed as control features, they provide advantages over point features in aerial triangulation automation. Point-based aerial triangulation based on manual measurement and the identification of conjugate points is less reliable than feature-based aerial triangulation because it has the limitations of visibility (occlusion), ambiguity (repetitive patterns), and semantic information in the light of robust and appropriate automation. Automation of aerial triangulation and pose estimation is obstructed by the correspondence problem, but the employment of splines is one way to overcome occlusion and ambiguity issues. The manual identification of corresponding entities in two images is crucial in the automation of photogrammetric tasks. A further problem of point-based approaches is their weak geometric constraints as compared with feature-based methods, so accurate initial values for the unknown parameters are required. Feature-based aerial triangulation can be implemented without conjugate points because the measured points in each image are not the conjugate points in this proposed adjustment model. Thus, tie splines that do not appear in all the overlapped images together can be employed in feature-based aerial triangulation. Another advantage of employing splines is that the adoption of high level features increases the feasibility of geometric information and provides an appropriate analytical solution that emphasizes the redundancy of aerial triangulation. 3D linear features expressed by 3D natural cubic splines are employed as the mathematical model of linear features in the object space and its counterpart in the projected image space for bundle block adjustment. To solve overparameterization of 3D natural cubic splines, arc-length parameterization using Simpson's rule is developed, and in the case of straight lines and conic sections, spline tangents can be additional equations to the overparameterized system. Photogrammetric triangulation by the proposed model, including the extended collinearity and arc-length parameterization equations, is developed to show the feasibility of tie and control splines for the estimation of the exterior orientation of multiple images, splines, and spline location parameters. A useful stochastic constraint for a spline segment is examined for its utility to become a full or partial control spline such as known EOPs with a tie, partial control, and full control spline, and unknown EOPs with a partial and full control spline. In addition, the information content of an image spline is calculated and the feasibility of a tie spline and a control spline for a block adjustment is described. A simulation bundle block adjustment is implemented prior to the actual experiment with real data in order to evaluate the performance of the proposed algorithms. A simulation can control the measurement errors so that random noises minimally affect the overall geometry of a block. The individual observations are generated based on the general situation of bundle block adjustment to estimate the properties of the proposed algorithms. A simulation allows adjustment for geometric problems or varying conditions within individual experiments.