Материал: [2.1] 3D Imaging, Analysis and Applications-Springer-Verlag London (2012)

where λ is the inverse step length. Gradient descent has the nice property that, unless a local minimum has been reached, one can always decrease the error by making the step length small enough. On the other hand, gradient descent is known to be slow and rather inefficient.

The Levenberg-Marquardt algorithm combines both Gauss-Newton and gradient descent updates in a relatively simple way:

A large value of λ yields a small, safe, gradient-descent step while a small value of λ favor large and more accurate steps of Gauss-Newton that make convergence to a local minimum faster. The art of a Levenberg-Marquardt algorithm implementation is in tuning λ after each iteration to ensure rapid progress even where Gauss-Newton fails. The now standard implementation is to multiply λ by 10 if the error increases and to divide it by 10 if the error decreases (with an upper bound at 108 and a lower bound at 10−4 for instance). In order to make the method robust to outliers, one may attenuate the influence of points with a large error by replacing the square error function by an M-estimator ε and an Iterative Reweighted Least Squared (IRLS)- like reweighting procedure. For instance, the following robust functions can be used:

6.6.2 Computing the Derivatives

An important issue in how Levenberg-Marquardt is applied to ICP is the one of computing the derivatives of the error function. The simplest approach is based on using finite differencing, assuming that the error function is smooth. However, this leads to a cost of p extra function evaluations per inner loop. In [32] a more effective solution was proposed based on the distance transform which also drastically improves the computational efficiency. The distance transform is defined as:

where x X and X is a discrete grid representing the volume which encloses the model-view M. Indeed, each data-point di can be easily associated to grid-points by obtaining the residual error ei = X(di ) in one shot.5 In other words, LM-ICP merges

5Note that the volume is discretized into integer values, therefore the data-point di should be rounded to recover X(di ).

the two main steps of ICP, namely closest point computation and transformation estimation, in a single step. Note further that when the mapping x → ε2( x ) is monotonic, we obtain that Dε (x) = ε2( D(x) ), so existing algorithms to compute D may be used to compute Dε , without requiring knowledge of the form of ε.

By combining Eq. (6.12) with Eq. (6.20) the new formulation of the registration problem becomes:

i=1

This formulation makes it much easier to compute the derivatives of E. In fact, since the distance transform is computed in a discrete form, it is possible to compute finite differences derivatives. More specifically,

In practice, xDε remains constant through the minimization, and we get:

Nd

aE(a) = xDε T (a, di ) aTT (a, di ). (6.22)

i=1

Note that the computation of aTT (a, di ) depends on the rigid transformation parametrization being used. In [32], the author proposed to model rotations by unitary quaternions for which the derivatives can be easily computed analytically. Finally, in order to compute the derivatives using matrix operators the Jacobian matrix is defined as Ji,j = ( xDε (T (a, di )) · aTj T (a, di )), where aj T (a, di ) =

6.6.3 The Case of Quaternions

Let the quaternion be defined by q = [s, v] where s and v are the scalar and vectorial components respectively [97]. Let d be the point on which the rotation must be

applied. To this aim such a point must be represented in quaternion space, leading to r = [0, d]. Therefore, the rotated point is obtained by:

r = qrq−1.

By multiplying in quaternion space6 we obtain:

r = 0, s2d + (d · v) · v + 2s(v × d) + v × (v × d) .

We represent this rotated point as:

r = [0, Tx , Ty , Tz],

where:

Tx = s2dx + (dx vx + dy vy + dzvz)vx + 2s(vy dz − vzdy ) + vy (vx dy − vy dx )

−vz(vzdx − vx dz)

=s2dx + vx2dx + vx vy dy + vx vzdz + 2svy dz − 2svzdy + vx vy dy − vy2dx

−vz2dx + vx vzdz

=s2 + vx2 − vy2 − vz2 dx + 2(vx vy − svz)dy + 2(vx vz + svy )dz

Ty = s2dy + (dx vx + dy vy + dzvz)vy + 2s(vzdx − vx dz) + vz(vy dz − vzdy )

−vx (vx dy − vy dx )

=s2dy + vx vy dx + vy2dy + vy vzdz + 2svzdx − 2svx dz + vy vzdz − vz2dy

−vx2dy + vx vy dx

=2(vx vy + svz)dx + s2 − vx2 + vy2 − vz2 dy + 2(vy vz − svx )dz

Tz = s2dz + (dx vx + dy vy + dzvz)vz + 2s(vx dy − vy dx ) + vx (vzdx − vx dz)

−vy (vy dz − vzdy )

=s2dz + vx vzdx + vy vzdy + vz2dz + 2svx dy − 2svy dx + vx vzdx − vx2dz

−vy2dz + vy vzdy

=2(vx vz − svy )dx + 2(vy vz + svx )dy + s2 − vx2 − vy2 + vz2 dz.

Now we introduce the translation component [tx , ty , tz] and normalize the quaternion by obtaining:

6A multiplication between two quaternions q and q is defined as [ss − v · v , v × v + sv + s v].