# Two-View 3D reconstruction using the sparse Levenberg–Marquardt algorithm

I have trouble implementing the Levenberg-Marguardt algorithm as described in the book Multi View Geometry in Computer Vision. To be more specific I have trouble calculating the partial derivatives that are needed in order to pursue with the main algorithm implementation.

Consider having the projection vector in two views $\mathbf{X}=[\mathbf{x},\mathbf{y},\mathbf{x}',\mathbf{y}']$ where the $x,y$ are the projected coordinates in the first view and $\mathbf{x}′,\mathbf{y}′$ are the coordinates in the second view respectively. Also lets say we are given an initial estimated $\mathbf{A}$ Matrix with 12 parameters describing the rotation and the translation of the second view with the respect to the first view and an initial estimation for the point in the 3D coordinates $\mathbf{b}=[\mathbf{X},\mathbf{Y},\mathbf{Z}]$.

How can the matrices $\partial\mathbf{X}/\partial{A}$ and $\partial\mathbf{X}/\partial{b}$ can be calculated ?

PS: I tried searching for numerical differentiation methods but I had hard time understanding them i.e finite difference.

• Your question is more suited for Computational Science (SCICOMP.SE). Feb 3 '17 at 10:21
• Also, you should try not to cut-and-paste equations into here. They don't work... and make them unreadable. Closing for clarity issues.
– Peter K.
Feb 3 '17 at 14:35
• @Gilles thank you for pointing my mistake out. I have been asking the very same question to both stackoverflow and mathexchange thats why i had to copy and paste the body of the question. Im really sorry Feb 3 '17 at 15:45

Along with many other derivative based descent algorithms, the Levenberg-Marquardt algorithm relies on the partial derivative matrix, a.k.a. Jacobian Matrix, which is the matrix of all first-order partial derivatives of a vector-valued function:

\begin{align} J_i &= \frac{\partial f(x_i, \beta)}{\beta} \\ \mathbf{J} &= \begin{bmatrix} \dfrac{\partial f_1}{\partial x_1} & \cdots & \dfrac{\partial f_1}{\partial x_n}\\ \vdots & \ddots & \vdots\\ \dfrac{\partial f_m}{\partial x_1} & \cdots & \dfrac{\partial f_m}{\partial x_n} \end{bmatrix} \end{align}

Now how you compute Jacobian depends on what you actually want from the minimization. You could for example do a full bundle adjustment where Jacobian might have an involved expression. But, let me give an example of the simplest scheme I could think of. Triangulation from multiple views is generally written in the form of minimization of the reprojection error. Here is the cost function for the two view problem:

$$C(X) = d(x,\hat{x})^2 + d(x',\hat{x}')^2$$ where $x$ and $x'$ are points in 2D views and $\hat{x}$ and $\hat{x}'$ are the projections of 3D points. In the matrix notation, the error function reads:

$$\begin{pmatrix} \epsilon_1 \\ \vdots \\ \epsilon_m \\ \end{pmatrix} = \begin{pmatrix} \|x_1^a-P_a X^1\|+ \|x_1^b-P_b X^1\|\\ \vdots \\ \|x_m^a-P_a X^m\|+\|x_m^b-P_b X^m\| \\ \end{pmatrix}$$ where $P_a, P_b$ are the projection matrices of cameras $a$ and $b$. Since we are searching for the 3D points only $\{\mathbf{X}_i\}$, the Jacobian is $Mx3$:

\begin{align} \mathbf{J} &= \begin{bmatrix} \dfrac{\partial \epsilon_1}{\partial x} & \dfrac{\partial \epsilon_1}{\partial y} & \dfrac{\partial \epsilon_1}{\partial z}\\ \vdots & \ddots & \vdots\\ \dfrac{\partial \epsilon_m}{\partial x} & \dfrac{\partial \epsilon_m}{\partial y} & \dfrac{\partial \epsilon_m}{\partial z} \end{bmatrix} \end{align}

Of course this problem requires a good initial solution, which could be obtained from a linear closed form solution. The quality of the non-linear minimizer is highly dependent on the initialization and the quality of 2D points.

For the rest, and different parameterizations, check out the links.

Here are some remarks:

1. If you are not looking for re-inventing the wheel, then I would recommend Ceres Solver, which includes Automatic Differentiation. This saves the headache of Jacobians (and is much more accurate than numerical derivatives - numerical approximations are highly discouraged). Ceres also includes samples for bundle adjustment.

2. If you are also optimizing for the pose (camera orientation), then I do not recommend operating on the 12 parameter vector of rotation and translations. Rather choose a representation such as quaternions or the exponential map (twist coordinates). As this preserves the manifold-ness, you will get much higher quality solutions, without ambiguities. Some examples are given in the links I provided.