Iterative parametric point (IPP) algorithm

Input: A set of profile maps ${\bf j}_k$, $k=1,\ldots ,L$, acquired from different viewpoints. Initial estimates for the unknown parameters.

Output: Refined estimates for ${\bf a}$, or ${\bf c}$, or ${\bf a}$ and ${\bf c}$.

Iterate until convergence,

1.

Transform the maps ${\bf j}_k$, $k=1,\ldots ,L-1$, onto all the other maps ${\bf j}_l$ for which l>k using the current estimates of the unknown parameters.

2.

For each pair of overlapping maps ${\bf j}_k$ and ${\bf j}_l$, k<l, construct an interpolated image ${\bf\tilde j}_l$ on S_k from the pixel values of ${\bf j}_l$ at the intermediate locations ${\bf s}_{kl}$ hit by the transformed map ${\bf j}_{kl}$. For each pixel $({\bf s}_k,j_k({\bf s}_k))$ of ${\bf j}_k$, the corresponding pixel on ${\bf j}_l$ is chosen to be $({\bf s}_{k},\tilde j_l({\bf s}_{k}))$.

3.

Find compatible matches and compute the weighting images ${\bf w}_{kl}$ on S_k. These include

• an adaptive weighting based on the statistical distribution of the difference in the direction of the surface normal at the corresponding points; ${\bf w}^{(\hbox{n})}_{kl}=\chi _{<0}\left ({\bf n}_{\eta _{kl}}-\eta _{kl,\hbox{max}}^{(r)}\right )$, where $\chi$ is a characteristic function, ${\bf n}_{\eta _{kl}}$ is an image of the angles $\eta _{kl}$ between the surface normals, and $\eta _{kl,\max}^{(r)}$ is a threshold set for $\eta _{kl}$ at the rth iteration. Denote the weighted sample mean and deviation of the values of ${\bf n}_{\eta _{kl}}$ at the rth iteration by $\nu _{kl}$ and $\varrho _{kl}$, respectively. The threshold is then updated to
  $$\begin{align*}\eta _{kl,\max}^{(r+1)}=&(\nu _{kl}+3\varrho _{kl})\cdot\chi _{<\e... ...n _2)}(\nu _{kl})+\nu _{kl}\cdot\chi _{\geq\epsilon _2}(\nu _{kl}), \end{align*}$$
  where $0<\epsilon _0<\epsilon _1<\epsilon _2$ are constants that are chosen according to an expected accuracy of the surface normal compatibility.
• an adaptive weighting based on the statistical distribution of the distance between the corresponding points; ${\bf w}^{(\hbox{d})}_{kl}$ [Zhang, 1994].
• a weighting according to a decay function (or zero weights) at locations near edges;
  
  $${\bf w}^{(\hbox{e})}_{kl}=\big ({\bf 1}_k\wedge\vert\vert{\bf... ...ge\vert\vert{\bf\tilde e}_l\vert +1-\epsilon\vert ^{-1}\big ),$$
  
  
  where $\wedge$ stands for the binary minimum operation, ${\bf e}_k$ and ${\bf\tilde e}_l$ are Laplacian images, and $\epsilon >0$.
• a weighting image taking into account the precision of the distance between the corresponding points; $w_{kl}^{(\hbox{v})}({\bf s})=s_\sigma (\hbox{Var}(d_{kl}({\bf s})))^{-1/2}$, where the variance is estimated using the rules of first order error propagation and the scale $s_\sigma$ is selected according to the noise level of the data so that the weights are around one.

4.

Update simultaneously all the parameters by taking one step according to the Levenberg-Marquardt algorithm towards the minimum of the mean of the squares of weighted distances between the corresponding points given by

$$f_1({\bf j}_1,\ldots ,{\bf j}_L,{\bf a},{\bf c})=\sum _{l=2}^L{\sum _{k=1}^{l-1}{{\bf w}_{kl}^2\bullet{\bf d}_{kl}^2/K}},$$ (3)

where the distance image ${\bf d}_{kl}={\bf j}_{kl}-{\bf\tilde j}_l$, $K=\sum _{l=2}^L{\sum _{k=1}^{l-1}{{\bf w}_{kl}^2\bullet{\bf 1}_k}}$, and ${\bf 1}_k$ is a unit image on S_k having ones at all locations.

Remarks:

• For high accuracy, the interpolation errors should be as equal as possible within the overlapping areas. Edge areas are removed or weighted according to a decay function. A bicubic interpolation is used in smooth areas and a bilinear one elsewhere.
• Other distance images include ${\bf d}_{kl}'=[h(({\bf r}_{kl}-{\bf r}_l^\ast )^2)]^{1/2}$, where ${\bf r}=({\bf x},{\bf y},{\bf z})$, the square of a multi valued image is a unary homogeneous operation, and the function h sums up the components of a multi valued image, e.g., $h(({\bf x},{\bf y},{\bf z}))={\bf x}+{\bf y}+{\bf z}$. The corresponding point for $({\bf s}_k,{\bf r}_k({\bf s}_k))$ on S_k is chosen to be the one $({\bf s}_l,{\bf r}_l({\bf s}_l))$ on S_l which is closest as measured by the Euclidean distance between ${\bf r}_{kl}({\bf s}_k)$ and ${\bf r}_l({\bf s}_l)$. For each ${\bf s}_k\in S_k$, the pixel value of ${\bf r}_l^\ast$ equals the one of ${\bf r}_l$ which has been defined as the corresponding point by the closest point operation. This distance image is used in the standard iterative closest point (ICP) algorithm [Besl and McKay, 1992].
• In sequential registration, the merit function to be minimized is given by
  

  $$f_2({\bf j}_1,\ldots ,{\bf j}_l,{\bf c}_1,\ldots ,{\bf c}_l,{... ...) =\sum _{k=1}^{l-1}{{\bf w}_{kl}^2\bullet{\bf d}_{kl}^2/K_l},$$ (4)

  
  
  where $K_l=\sum _{k=1}^{l-1}{{\bf w}_{kl}^2\bullet{\bf 1}_k}$. The unknown parameters ${\bf a}_{1l}$ can be updated using a closed form solution based on the unit quaternions.