Optimization problem formulation

If too many transformations are left as free variables for the optimizer to find, the system will be under-determined, and the optimization routine will fail: complaining about a "not positive definite" (singular in this case) Hessian.

Example: we have 1 stationary camera observing 10 chessboards. We want to be able to uniquely represent the transformation between each chessboard and the camera, for a total of 10 transformations. If we optimize a camera pose \(T_\mathrm{cr}\) camera and 10 separate chessboard poses \(T_\mathrm{rf}\) for each chessboard, we will have 11 transformations in the optimization vector. Since 11 > 10, we have more variables in the optimization vector than are needed to uniquely describe the world. So the system is under-determined, and the optimization will fail.

In a vanilla calibration problem such as this, we would address this by fixing the reference coordinate system to one of the camera or chessboard frames. The mrcal convention is to fix the reference coordinate system to camera 0, setting \(T_\mathrm{cr}\) to the identity transformation. In the above example, this would reduce the number of transformations being optimized from 11 to 10, which resolves the issue.

Any other method of making the optimization variables unique is valid also. For instance, the chessboard poses might be known. In that case we don't need to optimize any \(T_\mathrm{rf}\), and solving for the \(T_\mathrm{cr}\) is enough.

The reference coordinate system is a single coordinate system common to the whole optimization problem that all the objects in the world can use to localize themselves. It does not have any physical meaning beyond that. In particular, the reference coordinate system is not attached to any fixed object in the world. Thus noise in the chessboard observations would shift the reference coordinate system, just as it would shift the camera and chessboard coordinate systems. The projection uncertainty documentation talks about this in depth.

Detections of these corners serve as the input features to mrcal. mrcal is a purely geometrical toolkit, and this vision problem must be handled by another library. A number of tools are available to detect chessboard corners. These did not work well for my use cases, so I recommend mrgingham for all corner detections. mrgingham is fast and is able to find chessboard corners subject to very un-pinhole-like projections. At this time it has two limitations that will be lifted eventually:

• It more or less assumes a grid of 10x10 corners (i.e. 11x11 squares)
• It requires all the corners to be observed in order to report the detections from an image. Incomplete chessboard observations aren't supported

If these are unacceptable, any other detector may be used instead.

When given an image of a chessboard, the detector is directly observing the feature we actually care about: the corner. Another common calibration board style is a grid of circles, where the feature of interest is the center of each circle. When given an image of such a grid of circles, the detector either

• detects the contour at the edge of each circle
• finds the pixel blob comprising each circle observation

and from either of these, the detector infers the circle center. This can work when looking at head-on images, but when given tilted images subjected to non-pinhole lens behaviors, getting accurate circle centers from outer contours or blobs is hard. The resulting inaccuracies in the detections of circle centers will introduce biases into the solve that aren't modeled by the projection uncertainty routine, so chessboards are strongly recommended in favor of circle grids.

mrcal assumes independent noise on each point observation, so correlated sources of point observations (such as corners of an apriltag) are also not appropriate sources of data. Apriltag centers would work, however.

The calibration object is assumed to be nominally planar. However, large calibration boards used for calibration of wide lenses are never flat: temperature and humidity effects deform the board strongly-enough to affect the calibration. mrcal currently supports a simple 2-parameter deformation model. This model uses two axis-aligned parabolic factors. Let the chessboard grid span \([-1,1]\) along the \(x\) and \(y\) axes. Then I define the non-planar deformation as \(z \equiv k_x (1 - x^2) + k_y (1 - y^2)\) with \(k_x\) and \(k_y\) being the two deformation factors being optimized by the solver. If the board were flat, \(k_x\) and \(k_y\) would be 0, and thus we would have \(z=0\) everywhere. The deflection at the edges is 0, and is strongest at the center.

Empirically, this appears to work well: I get better-fitting solves, and less systematic error. And the optimal deformation factors \(k_x\), \(k_y\) are consistent between different calibrations.

Clearly, this does not work for especially strong or asymmetric deflections. There's a richer 5-parameter deformation model in a not-yet-released branch that appears to work even for asymmetric deflections. This needs more testing, and has not yet been released. Talk to Dima if you want to play with it.

When analyzing the behavior of the optimizer it is often useful to pick out particular elements of the full optimization vector \(\vec b\). mrcal provides a number of functions to report the index and size of the block of \(\vec b\) that contains specific data. In C:

• mrcal_state_index_intrinsics(): Return the index in the optimization vector of the intrinsics of camera i
• mrcal_state_index_extrinsics(): Return the index in the optimization vector of the extrinsics of camera i
• mrcal_state_index_frames(): Return the index in the optimization vector of the pose of frame i
• mrcal_state_index_points(): Return the index in the optimization vector of the position of point i
• mrcal_state_index_calobject_warp(): Return the index in the optimization vector of the calibration object warp
• mrcal_num_states_intrinsics(): Get the number of intrinsics parameters in the optimization vector
• mrcal_num_states_extrinsics(): Get the number of extrinsics parameters in the optimization vector
• mrcal_num_states_frames(): Get the number of calibration object pose parameters in the optimization vector
• mrcal_num_states_points(): Get the number of point-position parameters in the optimization vector
• mrcal_num_states_calobject_warp(): Get the number of parameters in the optimization vector for the board warp
• mrcal_num_states(): Get the full length of the optimization vector

And in Python:

• mrcal.state_index_intrinsics(): Return the index in the optimization vector of the intrinsics of camera i
• mrcal.state_index_extrinsics(): Return the index in the optimization vector of the extrinsics of camera i
• mrcal.state_index_frames(): Return the index in the optimization vector of the pose of frame i
• mrcal.state_index_points(): Return the index in the optimization vector of the position of point i
• mrcal.state_index_calobject_warp(): Return the index in the optimization vector of the calibration object warp
• mrcal.num_states_intrinsics(): Get the number of intrinsics parameters in the optimization vector
• mrcal.num_states_extrinsics(): Get the number of extrinsics parameters in the optimization vector
• mrcal.num_states_frames(): Get the number of calibration object pose parameters in the optimization vector
• mrcal.num_states_points(): Get the number of point-position parameters in the optimization vector
• mrcal.num_states_calobject_warp(): Get the number of parameters in the optimization vector for the board warp
• mrcal.num_states(): Get the full length of the optimization vector

If plotting a whole vector of state (or a vector of measurements), it is really helpful to annotate the plot to make it clear which variables correspond to each block of state (or measurements). mrcal provides helper functions to help with this:

• mrcal.plotoptions_state_boundaries(): Return the set plot options for gnuplotlib to show the state boundaries
• mrcal.plotoptions_measurement_boundaries(): Return the set plot options for gnuplotlib to show the measurement boundaries

The nonlinear least squares-solving library used by mrcal is libdogleg, which implements Powell's dogleg method. This is a trust-region algorithm that represents the trust region as a ball in state space. I.e. the radius of this trust region is the same in every direction. And that means that the optimization will work best when each state variable in \(\vec b\) affects the cost function \(E\) evenly. Example of what we don't want: camera positions measured in km, while the chessboard positions are measured in mm, with both sets of these very different numbers stored in \(\vec b\).

Clearly getting identical behavior from each variable is impossible, but we can scale the elements of \(\vec b\) to keep things more or less even. mrcal applies this scaling, and the libdogleg optimization library never sees the full state vector \(\vec b\), but the scaled vector \(\vec b_\mathrm{packed}\). Similarly, it never sees the full jacobian \(J \equiv \frac{\partial \vec x}{\partial \vec b}\), but rather \(J_\mathrm{packed} \equiv \frac{\partial \vec x}{\partial \vec b_\mathrm{packed}}\). This means that the optimization callback functions report packed state. These are

• mrcal_optimizer_callback() routine in the C API
• mrcal.optimizer_callback() routine in the Python API

To pack or unpack an array of state, mrcal provides some routines. In C:

• mrcal_pack_solver_state_vector(): Scales a state vector to the packed, unitless form used by the optimizer
• mrcal_unpack_solver_state_vector(): Scales a state vector from the packed, unitless form used by the optimizer

And in Python:

• mrcal.pack_state(): Scales a state vector to the packed, unitless form used by the optimizer
• mrcal.unpack_state(): Scales a state vector from the packed, unitless form used by the optimizer

In addition to \(\vec x_\mathrm{observations}\), the measurement vector contains regularization terms. These are mostly-insignificant terms that are meant to improve the convergence of the solver. They are also aphysical, and cause a bias in the solution, so mrcal is careful to keep these small-enough to not break anything noticeably. The behavior of these terms is likely to change in the future, so I don't document these in detail; please consult the sources. Currently the logic is at the end of the optimizer_callback() function in mrcal.c.

It is possible to control whether a solve does/does not include regularization terms with the do_apply_regularization bit in mrcal_problem_selections_t or the do_apply_regularization key in the call to mrcal.optimize().

When analyzing the behavior of the optimizer it is often useful to pick out particular elements of the full measurement vector \(\vec x\). mrcal provides a number of functions to report the index and size of the block of \(\vec x\) that contains specific data. In C:

• mrcal_measurement_index_boards(): Return the measurement index of the start of a given board observation
• mrcal_measurement_index_points(): Return the measurement index of the start of a given point observation
• mrcal_measurement_index_regularization(): Return the index of the start of the regularization measurements
• mrcal_num_measurements_boards(): Return how many measurements we have from calibration object observations
• mrcal_num_measurements_points(): Return how many measurements we have from point observations
• mrcal_num_measurements_regularization(): Return how many measurements we have from regularization
• mrcal_measurements(): Return how many measurements we have in the full optimization problem

And in Python:

• mrcal.measurement_index_boards(): Return the measurement index of the start of a given board observation
• mrcal.measurement_index_points(): Return the measurement index of the start of a given point observation
• mrcal.measurement_index_regularization(): Return the index of the start of the regularization measurements
• mrcal.num_measurements_boards(): Return how many measurements we have from calibration object observations
• mrcal.num_measurements_points(): Return how many measurements we have from point observations
• mrcal.num_measurements_regularization(): Return how many measurements we have from regularization
• mrcal.num_measurements(): Return how many measurements we have in the full optimization problem

I solve the calibration problem using Ordinary Least Squares, minimizing the discrepancies between pixel observations and their predictions. The pixel observations \(\vec q_\mathrm{ref}\) are noisy, and I assume that they are zero-mean, independent and normally-distributed. In particular, I treat the 2 values in each observation (\(x\) and \(y\)) as two independent measurements. I have no prior proof that the noise truly meets all those criteria, but empirical evidence suggests that these are all reasonable assumptions. And they simplify lots of analyses that we want to do. In order to propagate the input noise, we need to quantify it: for the \(i\) -th observed point, what is \(\mathrm{Var}\left(\vec q_{\mathrm{ref}_i}\right)\)?

Chessboard corner detectors often make it easy to infer the relative accuracy levels between the different corners, as opposed to an absolute noise level for each one. Thus the implementation splits the observed noise into two parts:

• The baseline standard deviation of the noise \(\sigma\). This is one value that applies to all the observations
• The scale \(s_i\) applied to that baseline. These are different for each observed corner

The mrgingham corner detector, in particular, reports the resolution used in detecting each corner as a decimation level: level-0 is "full-resolution", level-1 is "half-resolution" and so on. From that decimation level we get the relative scale

\[ s_i \equiv 2^{\mathrm{level}} \]

and we can define the 2x2 variance for each observed corner

\[ \mathrm{Var}\left( \vec q_{\mathrm{ref}_i} \right) = s_i^2 \sigma^2 I \]

and the variance for all the pixel observations

\[\mathrm{Var}\left(\vec q_\mathrm{ref}\right) = \mathrm{diag}\left(s_i^2\right) \sigma^2 \]

The remaining piece is to compute \(\sigma\), but this is hard to measure directly. There's an attempt in mrgingham, but it doesn't obviously work well. Thus the current method is to estimate \(\sigma\) from the solve residuals:

\[\sigma = \sqrt{\mathrm{Var}\left( \vec x^*_\mathrm{observations} \right)}\]

where \(\vec x^*\) is the measurement vector at the optimum. If all the assumptions are satisfied and we have enough data, then the input noise is the only source of error, and the only thing at affects the optimal \(\vec x\), and this way of evaluating \(\sigma\) is reasonable. If we have too little data, we're going to be overfitting, and the \(\sigma\) computed using the above method will be too low. The current thought is that this will be hard to do with chessboards and even a half-hearted calibration (we'll always have enough data), so mrcal currently doesn't try to catch this case. This might be added in the future.

We know where each point was observed in reality, and we know where the state vector \(\vec b\) predicts each one would have been observed. So we can construct a vector of errors \(\vec q_\mathrm{err} \equiv \vec q_\mathrm{predicted}\left( \vec b \right) - \vec q_\mathrm{ref}\).

For the purposes of optimization we want to weight the errors of uncertain observations less than confident ones, and to do that we can use the same \(s_i\) scale factor we computed earlier. For point \(i\) I define the weight

\[w_i \equiv \frac{1}{s_i} \]

Let's construct a diagonal matrix of all these weights: \(W \equiv \mathrm{diag}\left( \vec w \right)\). Then the measurement vector is

\[ \vec x_\mathrm{observations} \equiv W q_\mathrm{err} = W \left( \vec q_\mathrm{predicted}\left( \vec b \right) - \vec q_\mathrm{ref} \right) \]

If we assume that the model fits the data and that we have enough data to not overfit (both reasonable assumptions), then

\[\mathrm{Var}\left( \vec x_\mathrm{observations} \right) = W \mathrm{Var}\left( \vec q_\mathrm{ref} \right) W^T = \sigma^2 I \]

where \(\sigma\) is the input noise we're propagating. Furthermore, \(\vec x_\mathrm{observations}\) is homoscedastic: each element as the same variance. I make two more (reasonable) assumptions:

• The rest of the measurement vector \(\vec x\) (regularization) is insignificant
• I consider the linear problem at the local linearization of my nonlinear system

And then I can make a larger statement: the optimal parameter vector we compute from the least-squares optimization is the maximum-likelihood estimate of the true solution.