Accelerated parametric chamfer alignment using a parallel, pipelined GPU realization

Abstract

Parametric chamfer alignment (PChA) is commonly employed for aligning an observed set of points with a corresponding set of reference points. PChA estimates optimal geometric transformation parameters that minimize an objective function formulated as the sum of the squared distances from each transformed observed point to its closest reference point. A distance transform enables efficient computation of the (squared) distances, and the objective function minimization is commonly performed via the Levenberg–Marquardt (LM) nonlinear least squares iterative optimization algorithm. The point-wise computations of the objective function, gradient, and Hessian approximation required for the LM iterations make PChA computationally demanding for large-scale datasets. We propose an acceleration of the PChA via a parallelized and pipelined realization that is particularly well suited for large-scale datasets and for modern GPU architectures. Specifically, we partition the observed points among the GPU blocks and decompose the expensive LM calculations in correspondence with the GPU’s single-instruction multiple-thread architecture to significantly speed up this bottleneck step for PChA on large-scale datasets. Additionally, by reordering computations, we propose a novel pipelining of the LM algorithm that offers further speedup by exploiting the low arithmetic latency of the GPU compared with its high global memory access latency. Results obtained on two different platforms for both 2D and 3D large-scale point datasets from our ongoing research demonstrate that the proposed PChA GPU implementation provides a significant speedup over its single CPU counterpart.

Appendix: Compositional approach for projective transformation estimation

The Jacobian matrix elements $\{J^c_{j,l}\}$ associated with the projective transformation in Table 1 require a division operation per-element, which is computationally expensive. To simplify the Jacobian calculations, we adopt a compositional approach [31] that eliminates the division operations and also enables further simplifications. The compositional approach for our 2D point set alignment by projective transformation proceeds as follows:

• The projective transformation defined by the current estimate ${\varvec{\alpha }}_t$ of the parameters is applied to each OS point ${\mathbf {p}}_j$ to obtain a corresponding warped point ${\mathbf {p'}}_j \,=\, {\mathcal {T}}_{{\varvec{\alpha }}_t}\left( {\mathbf {p}}_j\right)$.

• Each LM iteration, then estimates the incremental parameter update ${\varvec{\delta }}$ that minimizes

  $$\begin{aligned} f({\varvec{\delta }}) \,=\, \sum \limits _{j=1}^{N_p} \left|| {\mathbf {r}}\left( {\mathcal {T}}_{\left( {\varvec{\alpha }}^I+{\varvec{\delta }}\right) }\left( {\mathbf {p'}}_j\right) \right) \right||^2 , \end{aligned}$$

  (12)

  where ${\varvec{\alpha }}^I\,=\,\left[ 1,0,0,0,1,0,0,0 \right]$ is the parameter vector that corresponds to the identity transformation, i.e., ${\mathcal {T}}_{{\varvec{\alpha }}^I} \left( {\mathbf {p'}}_j \right) =\,{\mathbf {p'}}_j$.

• The updated projective transformation is obtained as

  $$\begin{aligned} {\mathcal {T}}_{{\varvec{\alpha }}_{t+1}} \,=\, {\mathcal {T}}_{{\varvec{\alpha }}_t} \circ {\mathcal {T}}_{{\varvec{\delta }}}, \end{aligned}$$

  (13)

  where $\circ$ denotes composition, or equivalently multiplication of the corresponding matrix representations.

Considerable simplification of the Jacobian matrix calculation is obtained because the calculation is performed at ${\varvec{\alpha }}^I$, where the term $w_j$ in Table 1 becomes unity, eliminating the need for division operations. Specifically, the Jacobian matrix ${\mathbf {J}}_j$ at the transformed point ${\mathbf {p'}}_j\equiv {\mathcal {T}}_{{\varvec{\alpha }}^I}({\mathbf {p'}}_j)$, is computed as

$$\begin{aligned} {\mathbf {J}}_j \,=\, \frac{\partial {\mathcal {T}}_{{\varvec{\alpha }}}({\mathbf {p'}}_j)}{\partial {\varvec{\alpha }}}\bigg |_{{\varvec{\alpha }}\,=\,{\varvec{\alpha }}^I} \,=\, \begin{pmatrix} x'&{}y'&{}1&{}0&{}0&{}0&{}-x'^2&{}-x'y'\\ 0&{}0&{}0&{}x'&{}y'&{}1&{}-x'y'&{}-y'^2 \end{pmatrix}. \end{aligned}$$

(14)

The Hessian matrix approximation elements are shown in Table 2, where additional simplifications are also incorporated.