Bundle adjustment

The observed points ${\tt m}_{ki}$ being fixed, a specific residual $r_{ki}=D({\tt m}_{ki},{\bf\hat{P}}_k {\tt\hat{M}}_i)^2$ is only dependent on the point $i$-th point and the $k$-th projection matrix. This results in a very sparse matrix for the Jacobian. This is illustrated in figure 5.6 for 3 views and 4 points.

Figure 5.6: Sparse structure of Jacobian for bundle adjustment.

Because of the block structure of the Jacobian solving the normal equations ${\bf J}^\top{\bf Jx=b}$ have a structure as seen in figure 5.7.

Figure 5.7: Block structure of normal equations.

It is possible to write down explicit formulas for each block. Let us first introduce the following notation:

$${\bf U}_k = \sum_i (\frac{\partial \hat{\tt m}_{ki}}{\partial \hat{\bf P}_k})^\top \frac{\partial \hat{\tt m}_{ki}}{\partial \hat{\bf P}_k}$$ (E16)

$${\bf V}_i = \sum_k (\frac{\partial \hat{\tt m}_{ki}}{\partial \hat{\tt M}_i})^\top \frac{\partial \hat{\tt m}_{ki}}{\partial \hat{\tt M}_i}$$ (E17)

$${\bf W}_{ki} = (\frac{\partial \hat{\tt m}_{ki}}{\partial \hat{\bf P}_k})^\top \frac{\partial \hat{\tt m}_{ki}}{\partial \hat{\tt M}_i}$$ (E18)

$$\epsilon({\bf P}_k) = \sum_i (\frac{\partial \hat{\tt m}_{ki}}{\partial \hat{\bf P}_k})^\top \epsilon_{ki}$$ (E19)

$$\epsilon({\tt M}_i) = \sum_i (\frac{\partial \hat{\tt m}_{ki}}{\partial \hat{\tt M}_i})^\top \epsilon_{ki}$$ (E20)

with $\frac{\partial \hat{\tt m}_{ki}}{\partial \hat{\bf P}_k}$ and $\frac{\partial \hat{\tt m}_{ki}}{\partial \hat{\tt M}_i}$ matrices containing the partial derivatives from the coordinates of $\hat{\tt m}_{ki}$ to the parameters of $\hat{\bf P}_k$ and $\hat{\tt M}_i$ respectively. In this case the normal equations can be rewritten as

$$\left[\begin{array}{cc} {\bf U} & {\bf W} \\ {\bf W}^\top &... ...{c} \epsilon({\bf P}) \\ \epsilon({\tt M}) \end{array} \right]$$ (E21)

where the matrices ${\bf U,V,W}, \Delta({\bf P}), \Delta({\tt M}), \epsilon({\bf P})$ and $\epsilon({\tt M})$ are composed of the blocks defined previously. Assuming ${\bf V}$ is invertible both sides of equation 5.21 multiplied on the left with

$$\left[ \begin{array}{cc} {\bf I} & {\bf -WV}^{-1} \\ {\bf0}& {\bf I} \end{array}\right]$$

to obtain

$$\left[ \begin{array}{cc} {\bf U-WV}^{-1}{\bf W}^\top & {\bf0... ...^{-1}\epsilon({\tt M}) \\ \epsilon({\tt M}) \end{array}\right]$$ (E22)

This can be separated in two groups of equations. The first one is

$$\left( {\bf U-WV}^{-1}{\bf W}^\top \right) \Delta({\bf P}) = \epsilon({\bf P}) - {\bf WV}^{-1}\epsilon({\tt M})$$ (E23)

and can be used to solve for $\Delta({\bf P})$. The solution can be substituted in the other group of equations:

$$\Delta({\tt M})= {\bf V}^{-1}\left( \epsilon({\bf M}) - {\bf W}^\top \Delta({\bf P})\right)$$ (E24)

Note that due to the sparse block structure of ${\bf V}$ its inverse can be computed very efficiently.

The only computationally expensive step consist of solving equation 5.23. This is however much smaller than the original problem. For 20 views and 2000 points the problem is reduced from solving 6000 unknowns concurrently to more or less 200 unknowns. To simplify the notations the normal equations were used in this presentation. It is however simple to extend this to the augmented normal equations.