Week 9: Automatic differentiation but not automatic comprehension

When I saw that the function taking the longest was from Ceres Solver, I was very surprised to see the familiar name and quite happy that it turned out to be so important after all the hours I was forced to labor for it. The operation that I looked into this week was bundle adjust. To start off, let’s revisit how reconstruction works. First, a pair of images is added to the reconstruction, kicking off the grow_reconstruction function. Using the matching features, the RANSAC random sampling algorithm is used to create a triangulation, a representation of the geometric model using triangles. To increase the accuracy of the triangulation, it is bundle adjusted to the prior reconstruction by minimizing the discrepancy in their projections. If there’s GPS or GCP (ground control points) data as well, those are considered first to be absolutely accurate. 

Ceres Solver

Ceres Solver, the library I was struggling with several weeks ago to setup, is a framework for solving optimization problems. To use Ceres Solver, you need to setup the problem with a cost function and a parameter block. A cost function is the customizable measure for difference between expected and actual values. A parameter block represents either a single parameter or a group of scalars that define a parameter, like the 3 components of a translation vector and the 4 components of a quaternion for camera poses. Typically the least squares method is used on the cost function to calculate the residual block. Then, the parameter block is iteratively tweaked until the residual block is minimized. Since we are working with experimental data, which might have outliers, we can introduce loss functions to reduce the impact of the outliers on the residuals. 

To calculate the change in the residual block after tweaking the parameter block, we need to know the Jacobian. The Jacobian is a matrix representation of the partial derivatives of a function with respect to each of its parameters, for each of the functions defined with respect to those parameters. Rather than analytically calculating and storing the Jacobian, there is a method called automatic differentiation. 

Automatic differentiation is based on the concept of dual numbers. Similar to imaginary numbers, dual numbers introduce another unit e which has the property that e^2 = 0 but e is an infinitesimally small number that approaches but does not equal 0. We can make use of this unit to calculate the derivative of a function by substituting the variable which we are taking the derivative with respect to, with the sum of that variable and e. By solving for the function, the derivative is the coefficient of e. For example, if we had the function

f(x) = x(x + 1)

We know that the derivative of f(x) must be

f’(x) = 2x + 1

If we substitute x + e for x,

f(x + e) = (x + e)(x + e + 1)

            = x^2 + (2x + 1)e + x + e^2

In Ceres Solver, these dual numbers are implemented with the class Jet. 

Bundle Adjust

The primary function of bundle adjust is to setup the minimization problem. The linear solver used for bundle adjust is usually SPARSE_SCHUR because the solution tends to be sparse due to the lack of correlation between the parameters for different points and cameras. Here is a sample of how much data was in each of the parameter blocks: 

There are cost functions for prior errors, relative motion/rotation errors, common position errors, up vector errors, pan/tilt/roll angle errors, and linear motion errors. After adding all the cost functions, parameter blocks, and loss functions, the ceres::Solve() function is called to initialize the process. Each complete run took about 10 seconds on my 60 images. 

I didn’t get to do that much this week because I was bedridden with the flu for most of it (never been sick so often in the period of a few months.) For the upcoming week, I’ll be looking at the function in compute_depthmaps, PatchMatchUpdatePixel().