Bundle Block Adjustment (BBA) is the problem of estimating unknown elements associated with a network of images that are connected to each other with tie points. The initial geometrical structure of the network must be known prior to performing the bundle adjustment calculations.
As the name of the method suggests, the main foundation of BBA is related to the equations of collinearity in which we imagined a straight ray of light that established a connection between image space and model (or ground).
During BBA calculations, all possible equations are gathered under non-linear least-square optimization umbrella to obtain an accurate estimate of the unknowns. Here we solve a convex optimization problem of root finding type. Newton’s optimization method based on the sum of squares of image residual errors (or it is better to say residuals of all observational equations) is usually used to address the problem in hand. In this process, all the unknown parameters associated with a network of images are extracted simultaneously, so some inconsistent parameters such as lens distortions, along with other parameters such as exterior orientation elements and coordinates of ground points, are simultaneously solved. This step is usually considered as the last step (or one of the last steps) in the problem of extracting 3D geometric information from 2D images. In BBA problem, the main observational equations are
1- Equations of the collinearity condition: we have two equations for each model point in each image (see relative orientation page).
2- Equations for ground control points: control points in the BBA are divided into three categories: 1- height control points, 2- planar control points, and 3- complete control points. For each element of a control point, an equation is added to the cost function. In the calibration problem, we saw the control points in the model space as the corners of a raster pattern or coded target.
3- Equations related to scale bars (known lengths in the model space). Scale bars are objects whose length of one side is known to us with acceptable accuracy. Scale bars are used to obtain the unknown scale of the photogrammetric model. For each scale bar, an observation equation is added to BBA equations.
4- Equations of Global Navigation Satellite Systems (GNSS) readings on the measuring platform: Usually, the measuring platform used in imaging is equipped with a locator sensor. In this case, the observations of locations at the moment of image capture are added as new equations to the collection of equations. An important point in this regard is the necessity of synchronizing the location sensor observations with the camera shutter (see more information in time synchronization page). If there is a time shift, this shift is also considered as an unknown parameter.
5- Equations of Inertial Measurement Unit [IMU] readings: This device is responsible for measuring the direction of the platform. The measurements of the IMU device are also considered as additional observational equations. Involving IMU observations enable us to perform direct geo-referencing.
Unknowns
The unknowns of BBA problem are
1- Elements of camera calibration: interior orientation parameters are potentially considered among the unknowns of BBA. Geometric strength of a network determines how many of these elements could be potentially estimated. In a strong network, all IOPs are included in the calculations, while in a weak network, all or part of the parameters are excluded from the calculation process and their values are assumed to be constant.
2- Elements of the internal geometry of the camera system: In a calibration process in which more than one camera involves, the internal geometry of the camera is entered as internal constraints in the calculation process.
3- Elements of exterior orientations of images: the positions and orientations of the images are considered unknown in the BBA.
4- 3D coordinates of tie points and GCPs.
5- Elements of the relative orientation of the optical center of the camera in relation to the phase center of the GNSS antenna, if observations from this sensor are available.
6- The elements of the relative rotations of the camera compared to the flatness measuring device. If there is an IMU in the measurement platform, three elements of relative rotation are added to the list of unknowns.
7- Model coordinates of nodal points
8- Unknown elements of control points
Observational Equations
Now that we described observational equations and unknowns, we state a mathematical formulation for those equations.
1- The main equations in solving the bundle problem are collinearity equations. It is possible that the parameters of several different sensors are involved in the problem at the same time. In this case, we must be careful to consider the index of unknowns correctly.
1- For GCPs we consider the equations
2- Sometimes there are one or more scale-bars with known details in the model or ground coordinate system. In this case, the observation equations for the known lengths are written as
In the equation above, the symbol SB is used to show the length of a scale bar.
3- To transfer the coordinates read by the positioning sensor, it is necessary to consider the relative parameters of the optical center with respect to the center of the sensor’s antenna. In this case, the following equations can be used
In the above equation, π0_πΊππ_π‘ is the coordinates read by the position sensor at time (t), π_πΆ_π‘ is the position of the optical center of the camera at time (t), π _πΆ_π‘ is the rotation matrix of the camera at time (t), π _πΆ_πΌ is the relative rotations of the camera with respect to the IMU measuring device, and π_πΆ_πΌ is the relative position of the camera with respect to the center of the positioning sensor antenna.
4- In order to transfer the rotations read by the IMU device, the relative rotations of the camera with respect to the device should be considered as follows
After forming the system of observational equations, we form the relevant matrices. It should be noted that placement of unknowns will affect the structure of matrices. In an optimal state, instead of directly forming the Jacobian matrix of coefficients (A), the elements of the matrix A^T.A are formed instead. This technique is important because it can significantly save computing resources and system storage. We will examine this method in detail later in this page.
Initial estimation of network structure
The initial estimation of the network structure is problem of estimating the network with sufficient accuracy in such a way that the exterior orientation elements of the network images will become known with sufficient accuracy. We need a suitable tool and strategy for the initial estimation of the network. A suitable tool, for example, is relative orientation methos. A suitable strategy for forming the network is an optimal way by which we can estimate the network with a good initial accuracy with our limited computing resources. Choosing an appropriate strategy in forming the network is very important, because an appropriate strategy in a logical timing provides us with a solid structure, while an inappropriate strategy can simply lead to divergence of the optimization operation.
To better understand the problem, consider the following example. In this example, 21 images were taken in one direction by a sensor installed on the drone.
A simple way to build a network is to find tie points between successive images and connect successive stereo pairs. This type of connection is the simplest ways to form the above network structure, but at the same time, it causes the resulting network to have very weak geometric strength.
The second strategy that we can express for this problem is connecting each image to two or three consecutive images that lie after it. By using this strategy, we usually expect to achieve a more robust structure in estimating the network structure, if there are problems in this case as well. One of the problems that may happen here is the problem of accumulating errors after connecting several consecutive images.
In the figure above, we can see the three-dimensional rotational ellipses showing the uncertainty about the position of the images with the above strategy with a visual magnification factor. In this figure, we can see that the farther the network images are from the first image, the largest the error ellipses become. One way to improve the second strategy that was mentioned is to identify closed loops in the network structure. In this way, if we guess that at some point the images have returned to the first place, we can establish a connection between the last image and the first image. The detection of rings can potentially play an effective role in reducing the size of error rotational ellipses. The second most effective suggestion for improving the network structure is to detect adjacent bands with overlap. In this case, we first identify the bands and calculate the horizontal overlaps, then connect the images in the two detected bands that have sufficient horizontal overlap.
Internal constraints and coordinate system definition
In BBA, three different modes can be considered to define a coordinate system
1- Minimum coordinate system (Minimum Constraint Block Bundle Adjustment)
In this case, we are satisfied with the minimum information about the coordinate system and we make the grid. A minimal definition of the coordinate system can be, for example, similar to the relative orientation in the form of defining the coordinate center corresponding to the first image in the network, and considering the largest optical center element of the farthest image from the reference image as a constant. Defining the minimal coordinate system is possible in many situations. It should be noted that the definition of the coordinate system will directly affect the shape and size of the three-dimensional rotational ellipses. Therefore, it seems necessary to be careful in choosing a proper coordinate system.
2- coordinate system corresponding to the center of gravity of optical centers (Inner Constraint Block Bundle Adjustment)
In this case, the center of gravity of the image centers is considered as the coordinate center, and sufficient restrictions are defined to prevent the rotation of the light centers in this system. In many cases, this coordinate system definition leads to acceptable results.
3- coordinate system corresponding to ground control points (Over constraint Block Bundle Adjustment)
In this case, the coordinate system corresponding to the ground control points is considered, and all rotational and translational parameters of the images are considered as free parameters.
4- coordinate system corresponding to optical centers (Over constraint Block Bundle Adjustment)
In this case, the coordinate system corresponding to the readings of the positions of the optical centers of the images is considered, and all the rotational and translations of the images are considered as free parameters. This mode is similar to the previous mode and it is also possible to combine this definition with the previous mode.
Large Decimal Numbers in a Global Coordinate System
Sometimes entering large decimal numbers in the bundle calculation is problematic. To solve this issue, a three-dimensional conformal transformation can be used to reduce numbers using a local auxiliary coordinate system. For this purpose, we consider one of the points as the origin of local coordinates to access smaller numbers. We can also consider the largest length in order to calculate a suitable scale.
Calculation of degrees of freedom
Example: In the following example, we calculate the degrees of freedom in a bundle adjustment. In this example, six images, four tie points, two complete ground control points, one planar control point, and one height control point are considered. We can see the visibility status of each point in the table below.
Now, using the definition of observation equations, we first count the number of possible equations:
The number of observation equations in the
1- first image: 5×2=10
2- second image: 4×2=8
3- third image: 5×2=10
4- fourth image: 3×2=6
5- fifth image: 4×2=8
6- sixth image: 6×2=12
7- The number of observation equations for complete control points is 2×3=6
So we have 60 equations in total. Now we count the number of unknowns in the problem:
EOPs: 6×6=36
Ground coordinates of tie points: 4×3=12
z component of planar control points: 1×1=1
x and y components of height control points: 1×2=2
Therefore, we have a total of 43 unknowns. As a result, there are 17 degrees of freedom in this problem.
Solving Bundle Adjustment using Non-Linear Least Squares and Sub-Matrix Sparse Method
According to the system of observational equations that were mentioned in this page, a system of equations is formed to divide the bundle, which in a simplified linear form is as
We use the least-square method with the following cost function
The first term is a weighted sum of squares of image residuals, and the second term is Lagrange coefficients. The partial derivatives of this function should be equal to zero at the solution so
Now we can write a set of constrained equations to solve them
Finally, one step of optimization would be in form of
The arrangement of the unknowns in the matrix X Μ is very important in the structure of the matrix A^T.(B.P^+.B^T)^+.A. In an optimal state, first the unknowns of the calibration parameters are considered and then the unknowns of the coordinates of the tie points are placed. In this case, the matrix A^T.(B.P^+.B^T)^+.A takes a structure like below
In the above structure, if the number of tie points is relatively large, the sub-matrix H_3 will have a much larger dimension than H_1. As we can see, the matrix H_3 has a manageable sparse structure, so it is clear that the inverse calculation of this matrix is ββrelatively simple considering this internal structure. For the reasons mentioned, we use the sub-matrices to invert the (B.P^+.B^T)^+.A matrix as follows
An effective solution for the above equations is obtained if we first find the initial part of the inverse of the matrix as
Now we put the result of the above equation in (1) as
Therefore, we simply get to the desired results as
The second part of the inverse matrix is calculated by the following equation as
In the next step, we calculate the final two parts of the inverse matrix, although it is possible to skip their calculation and obtain the unknown coordinates of the tie points by multi-image intersection. In case there are computational limitations or the dimensions of the problem are large, the mentioned technique can be used, otherwise, in order to calculate the last part effectively, we put the results of equation (7) in equation (4) as follows
The last element is calculated as
