An Unsupervised, Iterative N-Dimensional Point-Set Registration Algorithm

TL;DR

This paper introduces an unsupervised, iterative algorithm for aligning unlabeled N-dimensional point clouds, leveraging linear least squares to consider all point pairings and refine registration without prior correspondence knowledge.

Contribution

It presents a novel unsupervised iterative registration method that handles unlabeled N-dimensional point clouds using linear least squares and pairwise consideration.

Findings

Effective alignment of unlabeled point clouds in multiple dimensions

Iterative process converges without prior point correspondence

Applicable to high-dimensional data sets

Abstract

An unsupervised, iterative point-set registration algorithm for an unlabeled (i.e. correspondence between points is unknown) N-dimensional Euclidean point-cloud is proposed. It is based on linear least squares, and considers all possible point pairings and iteratively aligns the two sets until the number of point pairs does not exceed the maximum number of allowable one-to-one pairings.

Equations

$${\mathbf{u}}_i\in {\mathbb{R}}^N\in \mathbf{U}i=1,\dots ,N_U$$

$${\mathbf{v}}_k\in {\mathbb{R}}^N\in \mathbf{V}k=1,\dots ,N_V$$

$${\mathbf{u}}_i^{{}^{\prime }}=\mathbf{L}{\mathbf{u}}_i+\mathbf{t},$$

$$e(\mathbf{U},\mathbf{V};\mathbf{L},\mathbf{t})=\frac{{\sum }_{i=1}^{N_U}{\sum }_{k=1}^{N_V}{m}_{ik}({\mathbf{u}}_i^{{}^{\prime }}-{\mathbf{v}}_k{)}^T({\mathbf{u}}_i^{{}^{\prime }}-{\mathbf{v}}_k)}{{\sum }_{i=1}^{N_U}{\sum }_{k=1}^{N_V}{m}_{ik}}$$

$$\begin{gathered} e(\mathbf{U},\mathbf{V};\mathbf{L},\mathbf{t})=\sum _{i=1}^{N_U}\sum _{k=1}^{N_V}{m}_{ik}[{\mathbf{v}}_k^T{\mathbf{v}}_k-2\text{Tr}\left(\mathbf{L}{\mathbf{u}}_i{\mathbf{v}}_k^T\right)+2\text{Tr}\left(\mathbf{L}{\mathbf{u}}_i{\mathbf{t}}^T\right) \\ +{\mathbf{u}}_i^T{\mathbf{u}}_i-2{\mathbf{t}}^T{\mathbf{v}}_k+{\mathbf{t}}^T\mathbf{t}], \end{gathered}$$

$$\underset{\mathbf{L},\mathbf{t}}{\mathrm{minimize}}$$

$$\mathbf{L}{\mathbf{L}}^T={\mathbf{I}}_{N\times N},$$

$$\text{det}\mathbf{L}=1.$$

$$\mathcal{L}(\mathbf{L},\mathbf{t},\mathbf{\alpha },\lambda )=e(\mathbf{U},\mathbf{V};\mathbf{L},\mathbf{t})+\text{Tr}\left(\mathbf{\alpha }(\mathbf{L}{\mathbf{L}}^T-{\mathbf{I}}_{N\times N})\right)+\lambda \left(\text{det}\mathbf{L}-1\right),$$

$$\overline{\mathbf{u}}=\sum _{i=1}^{N_U}\sum _{k=1}^{N_V}{m}_{ik}{\mathbf{u}}_i$$

$$\overline{\mathbf{v}}=\sum _{i=1}^{N_U}\sum _{k=1}^{N_V}{m}_{ik}{\mathbf{v}}_k$$

$$\hat{\mathbf{t}}=\overline{\mathbf{v}}-\mathbf{L}\overline{\mathbf{u}}$$

$$\frac{\partial \mathcal{L}}{\partial \mathbf{L}}=2\left(\overline{\mathbf{v}}{\overline{\mathbf{u}}}^T-\sum _{i=1}^{N_U}\sum _{k=1}^{N_V}{m}_{ik}{\mathbf{v}}_k{\mathbf{u}}_i^T+\mathbf{\alpha L}\right)+\lambda \mathbf{L}=0_{N\times N},$$

$$\mathbf{Z}=\sum _{i=1}^{N_U}\sum _{k=1}^{N_V}{m}_{ik}{\mathbf{v}}_k{\mathbf{u}}_i^T-\overline{\mathbf{v}}{\overline{\mathbf{u}}}^T,$$

$$\mathbf{Z}=\left(\begin{array}{ccccc}\text{cov}(u^1,v^1)& \text{cov}(u^1,v^2)& \text{cov}(u^1,v^3)& \cdots & \text{cov}(u^1,v^N)\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ \text{cov}(u^N,v^1)& \text{cov}(u^N,v^2)& \text{cov}(u^N,v^3)& \cdots & \text{cov}(u^N,v^N)\end{array}\right),$$

$$\text{cov}(u^j,v^{j^{\prime }})=\sum _{i=1}^{N_U}\sum _{k=1}^{N_V}{m}_{ik}{u}_i^j{v}_k^{j^{\prime }}-\left(\sum _{i=1}^{N_U}\sum _{k=1}^{N_V}{m}_{ik}{u}_i^j\right)\left(\sum _{i=1}^{N_U}\sum _{k=1}^{N_V}{m}_{ik}{v}_k^{j^{\prime }}\right)$$

$$\hat{\mathbf{L}}={\left(\mathbf{\alpha }+\frac{\lambda }{2}{\mathbf{I}}_{N\times N}\right)}^{-1}\mathbf{Z}$$

$$\mathbf{L}{\mathbf{L}}^T={\left(\mathbf{\alpha }+\frac{\lambda }{2}{\mathbf{I}}_{N\times N}\right)}^{-1}\mathbf{Z}{\mathbf{Z}}^T{\left(\mathbf{\alpha }+\frac{\lambda }{2}{\mathbf{I}}_{N\times N}\right)}^{-1}={\mathbf{I}}_{N\times N},$$

$$\hat{\mathbf{L}}=\pm {\left(\sqrt{\mathbf{Z}{\mathbf{Z}}^T}\right)}^{-1}\mathbf{Z}$$

$$\mathbf{Z}{\mathbf{Z}}^T=\mathbf{P}\mathbf{D}{\mathbf{P}}^T,$$

$${\left(\mathbf{P}{\mathbf{D}}^{\frac{1}{2}}{\mathbf{P}}^T\right)}^2=\mathbf{P}{\mathbf{D}}^{\frac{1}{2}}{\mathbf{P}}^T\mathbf{P}{\mathbf{D}}^{\frac{1}{2}}{\mathbf{P}}^T=\mathbf{P}\mathbf{D}{\mathbf{P}}^T=\mathbf{Z}{\mathbf{Z}}^T$$

$$\hat{\mathbf{L}}=\mathbf{P}{\left(\sqrt{\mathbf{D}}\right)}_{\text{PD}}^{-1}{\mathbf{P}}^T\mathbf{Z}$$

$$\hat{\mathbf{L}}={\left(\sqrt{\mathbf{Z}{\mathbf{Z}}^T}\right)}_{\text{PD}}^{-1}\mathbf{Z}$$

$$\hat{\mathbf{t}}=\overline{\mathbf{v}}-{\left(\sqrt{\mathbf{Z}{\mathbf{Z}}^T}\right)}_{\text{PD}}^{-1}\mathbf{Z}\overline{\mathbf{u}}$$

$$\hat{\mathbf{L}}=\mathbf{P}{\left(\sqrt{\mathbf{D}}\right)}_{\text{PD}}^{-1}{\mathbf{P}}^T\mathbf{Z}$$

$$\hat{\mathbf{t}}=\overline{\mathbf{v}}-\mathbf{P}{\left(\sqrt{\mathbf{D}}\right)}_{\text{PD}}^{-1}{\mathbf{P}}^T\mathbf{Z}\overline{\mathbf{u}}$$

$$\begin{gathered} e_{\text{min}}=\left(\sum _{i=1}^{N_U}\sum _{k=1}^{N_V}{m}_{ik}{\mathbf{u}}_i^T{\mathbf{u}}_i-{\overline{\mathbf{u}}}^T\overline{\mathbf{u}}\right)+\left(\sum _{i=1}^{N_U}\sum _{k=1}^{N_V}{m}_{ik}{\mathbf{v}}_i^T{\mathbf{v}}_i-{\overline{\mathbf{v}}}^T\overline{\mathbf{v}}\right) \\ -2\text{Tr}\left({\left(\sqrt{\mathbf{Z}{\mathbf{Z}}^T}\right)}_{\text{PD}}\right) \end{gathered}$$

$$\sum _{i=1}^{N_U}\sum _{k=1}^{N_V}{m}_{ik}{\mathbf{u}}_i^T{\mathbf{u}}_i-{\overline{\mathbf{u}}}^T\overline{\mathbf{u}}=\sum _{j=1}^N{\sigma }_{u^j}^2,$$

$${\sigma }_{u^j}^2=\sum _{i=1}^{N_U}\sum _{k=1}^{N_V}{m}_{ik}{\left(u_i^j\right)}^2-{\left(\sum _{i=1}^{N_U}\sum _{k=1}^{N_V}{m}_{ik}{u}_i^j\right)}^2$$

$$\begin{array}{rl}{e}_{\text{min}}& =\sum _{j=1}^N\left({\sigma }_{u^j}^2+{\sigma }_{v^j}^2\right)-2\text{Tr}\left({\left(\sqrt{\mathbf{Z}{\mathbf{Z}}^T}\right)}_{\text{PD}}\right)\\ & =\sum _{j=1}^N\left({\sigma }_{u^j}^2+{\sigma }_{v^j}^2\right)-2\text{Tr}\left(\mathbf{P}{\left(\sqrt{\mathbf{D}}\right)}_{\text{PD}}{\mathbf{P}}^T\right)\end{array}$$

$$e^{\prime }=\sum _{j=1}^N\left({\sigma }_{u^j}^2+{\sigma }_{v^j}^2\right)+2\text{Tr}\left({\left(\sqrt{\mathbf{Z}{\mathbf{Z}}^T}\right)}_{\text{PD}}\right)$$

$${\mathbf{u}}_i^{{}^{\prime }}=s\mathbf{L}{\mathbf{u}}_i+\mathbf{t}$$

$$\begin{gathered} e(\mathbf{U},\mathbf{V};\mathbf{L},\mathbf{t},s)=\sum _{i=1}^{N_U}\sum _{k=1}^{N_V}{m}_{ik}[{\mathbf{v}}_k^T{\mathbf{v}}_k-2s\text{Tr}\left(\mathbf{L}{\mathbf{u}}_i{\mathbf{v}}_k^T\right)+2s\text{Tr}\left(\mathbf{L}{\mathbf{u}}_i{\mathbf{t}}^T\right) \\ +s^2{\mathbf{u}}_i^T{\mathbf{u}}_i-2{\mathbf{t}}^T{\mathbf{v}}_k+{\mathbf{t}}^T\mathbf{t}] \end{gathered}$$

Full text

0.1 Introduction

Point-set registration has been extensively studied in computer vision. For labeled data, an important class of solutions is linear least squared techniques [1, 2, 3, 4]. Unlike the aforementioned work of [1, 2, 3], the paper in [4] treated the more general case of two point-sets of unequal size and do not assume correspondence. They first established correspondence by numerically determining the matching pairs support between the two points sets (i.e. finding an optimal subset of pairings between the two sets), and then derived (analytical) least squared solutions to the transformation parameters that optimally align the (labeled) optimal subset of pairings. However, their approach is computationally expensive, having quartic polynomial (average-case) complexity.

For unlabeled data, the alignment of two point patterns is a two-part problem; both the correspondence and the optimal affine transform that minimizes some dissimilarity metric between the two point-sets need to be determined. Various point pattern matching algorithms for unlabeled data have been proposed in [5, 6, 7, 8, 9], but the necessary optimization for each is numerically based. The works of [5, 6, 7] determined the optimal affine transform and correspondence simultaneously by numerically solving a constrained least squares problem. The method in [9] models each point-set as a Gaussian mixture, and determines the appropriate alignment by (non-linear) optimization of the L2 distance between the two Gaussian mixtures.

Unlike the numerical optimization schemes discussed above, a much more computationally efficient approach would be an analytical optimization scheme for unlabeled data, which is proposed in this paper. Our alignment approach is similar to the labeled techniques of [1, 2, 3, 4] in that we derive analytical solutions for the optimal affine transform via linear least squares. But unlike them, we do not assume or establish correspondence prior to registration, but use registration to establish correspondence. To the best of our knowledge, our derivation of the (closed-form) linear least square solutions for the registration of two unlabeled point-sets, though remarkably simple, is absent from the available literature. The presented derivation is generalized to any $N$-dimensional point-set (i.e. each point in a point-set resides in $\mathbb{R}^{N}$), and the obtained solutions are then used to create an unsupervised, iterative $N$-dimensional point-set registration algorithm. The $N=2$ case was shown in [10], and utilized in the context of fingerprint matching.

The paper is organized as follows. In Section II, we lay the theoretical foundations of the proposed algorithm; in Section III, we describe its numerical implementation; and in Section IV, we conclude with a discussion of the algorithm and its potential applications.

0.2 Theory

Consider two $N$-dimensional point sets ${\bf U}$ and ${\bf V}$ comprising $N_{U}$ and $N_{V}$ singular points, respectively. The Cartesian coordinates of the singular points will be expressed as an $N$-dimensional vector:

[TABLE]

The elements of ${\bf u}_{i}$ and ${\bf v}_{k}$ are denoted as $u_{i}^{j}$ and $v_{k}^{j}$, respectively, where $j=1,...,N$. ${\bf U}$ and ${\bf V}$ can be interpreted as matrices whose columns are formed by the vectors ${\bf u}_{i}$ and ${\bf v}_{k}$, respectively, i.e. ${\bf U}\in\mathbb{R}^{N_{U}\times N}$ and ${\bf V}\in\mathbb{R}^{N_{V}\times N}$; and $u^{j}$ and $v^{j}$ can be interpreted as the features of ${\bf U}$ and ${\bf V}$, respectively.

We want to register point-set ${\bf U}$ to ${\bf V}$. There are $N_{U}N_{V}$ possible matching (cross) pairs and at most $\min\{N_{U},N_{V}\}$ one-to-one matching pairs. Let $m_{ik}$ denote the weight of a matching pair; the weight can be interpreted as a probability that the points ${\bf u}_{i}$ and ${\bf v}_{k}$ match locally. We assume that the (initial) $m_{ik}$ for each (cross) pair has been determined prior to alignment. One way of computing $m_{ik}$ is discussed in [10] within the context of fingerprint matching.

0.2.1 Case I: No Scale

We apply a global rotation and translation to point set ${\bf U}$ such that

[TABLE]

where $L$ is the $N\times N$ rotation matrix and ${\bf t}$ is the $N$-d vector of translation parameters. Since ${\bf L}$ is a rotation matrix, it is orthogonal, i.e. ${\bf L}{\bf L}^{T}={\bf L}^{T}{\bf L}={\bf I}_{N\times N}$, and has a determinant of 1.

The measure of closeness of the transformed point set ${\bf U}$ and the template set ${\bf V}$ may be defined as the weighted sum of the squared distances between their points (i.e. Euclidean distance metric):

[TABLE]

Interpreting $m_{ik}$ as the probability of a match between ${\bf u}_{i}$ and ${\bf v}_{k}$ leads to $\sum_{i=1}^{N_{U}}\sum_{k=1}^{N_{V}}m_{ik}=1$. Expanding out the product yields

[TABLE]

where ’Tr’ denotes the trace operator. We seek the optimal values of parameters ${\bf L}$ and ${\bf t}$ that minimize $e({\bf U},{\bf V};{\bf L},{\bf t})$, so the constrained optimization problem to be solved is

[TABLE]

The Lagrangian for the above (constrained) optimization problem is

[TABLE]

where $\boldsymbol{\alpha}$ is an $N\times N$ symmetric matrix of Lagrangian multipliers (it is symmetric because ${\bf L}{\bf L}^{T}$ is symmetric, and consequently contains $N(N+1)/2$ Lagrangian multipliers), and $\lambda$ is a scalar Lagrangian multiplier.

Before proceeding further, let us define the weighted average coordinate vectors as

[TABLE]

Optimizing Eq. (4) with respect to ${\bf t}$ yields

[TABLE]

Substituting Eq.(5) back into our Lagrangian and then partial differentiating it with respect to ${\bf L}$ yields

[TABLE]

where $\boldsymbol{0}_{N\times N}$ denotes an $N\times N$ matrix of zeros. Let

[TABLE]

which is a ($N\times N$) cross-covariance matrix. Specifically,

[TABLE]

where $\text{cov}(u^{j},v^{j^{\prime}})$ denotes the covariance between features $u^{j}$ and $v^{j^{\prime}}$:

[TABLE]

Solving Eq. (6) for ${\bf L}$ yields

[TABLE]

Since ${\bf L}$ is an orthogonal matrix, we have

[TABLE]

which upon isolating the cross-covariance, ${\bf Z}$, yields $\left(\boldsymbol{\alpha}+\frac{\lambda}{2}{\bf I}_{N\times N}\right)^{2}=\boldsymbol{Z}\boldsymbol{Z}^{T}$, so

[TABLE]

The minimum of the error function corresponds to $\left(\sqrt{\boldsymbol{Z}\boldsymbol{Z}^{T}}\right)^{-1}\boldsymbol{Z}$, so we discard the solution $-\left(\sqrt{\boldsymbol{Z}\boldsymbol{Z}^{T}}\right)^{-1}\boldsymbol{Z}$.

The matrix $\boldsymbol{Z}\boldsymbol{Z}^{T}$ exhibits several interesting properties:

It is real symmetric. It is real because the elements of $\boldsymbol{Z}$ are covariances, which are always real by definition of covariance. And it is symmetric because $\boldsymbol{Z}\boldsymbol{Z}^{T}=\left(\boldsymbol{Z}\boldsymbol{Z}^{T}\right)^{T}$. 2. 2.

It is diagonalizable as a consequence of being real symmetric. We will exploit this fact later. 3. 3.

It is positive-definite because ${\bf x}^{T}\boldsymbol{Z}\boldsymbol{Z}^{T}{\bf x}=(\boldsymbol{Z}^{T}{\bf x})^{T}\boldsymbol{Z}^{T}{\bf x}=||\boldsymbol{Z}^{T}{\bf x}||^{2}>0$ for any ${\bf x}\in\mathbb{R}^{N}$. 4. 4.

It’s eigenvalues are all positive as a consequence of positive-definiteness.

Eq. (8) requires taking the square root of the matrix $\boldsymbol{Z}\boldsymbol{Z}^{T}$. According to the spectral theorem, any real symmetric matrix can be diagonalized by an orthogonal matrix:

[TABLE]

where ${\bf P}$ is an $N\times N$ orthogonal matrix and ${\bf D}$ is the $N\times N$ diagonal matrix; ${\bf P}$ and ${\bf D}$ are formed by the eigenvectors and eigenvalues, respectively, of $\boldsymbol{Z}\boldsymbol{Z}^{T}$.The square root of $\boldsymbol{Z}\boldsymbol{Z}^{T}$ is then ${\bf P}{\bf D}^{\frac{1}{2}}{\bf P}^{T}$ because

[TABLE]

Now, any $N\times N$ matrix with $N$ distinct eigenvalues has $2^{N}$ square roots because the square root of each eigenvalue can be either positive or negative. If such a matrix is further positive-definite, then it has precisely one positive-definite square root; the positive-definite square root corresponds to the case where only the positive square root of each eigenvalue is taken. We exploit these properties of $\boldsymbol{Z}\boldsymbol{Z}^{T}$ to extract only its positive-definite square-root, so ${\bf D}^{\frac{1}{2}}$ is formed by the positive square roots of the eigenvalues of $\boldsymbol{Z}\boldsymbol{Z}^{T}$. To this end, we employ the notation $\left(\sqrt{\text{arg}}\right)_{\textit{PD}}$ to refer to the positive-definite square root of any square matrix. Taking all this into account, Eq. (8) can be expressed as

[TABLE]

In summary, the optimal alignment parameters are

[TABLE]

Or equivalently,

[TABLE]

The minimum squared error is found by substituting the optimal alignment parameters back into Eq. (3), which yields

[TABLE]

Note that

[TABLE]

where

[TABLE]

In other words, $\sigma_{u^{j}}^{2}$ is the variance of feature $u^{j}$, and likewise $\sigma_{v^{j}}^{2}$ is the variance of feature $v^{j}$. So we can rewrite the minimum squared error as

[TABLE]

Note that had we used $\hat{\boldsymbol{L}}=-\left(\sqrt{\boldsymbol{Z}\boldsymbol{Z}^{T}}\right)_{\textit{PD}}^{-1}\boldsymbol{Z}$, the generated error would be

[TABLE]

The eigenvalues of a positive-definite matrix all always positive, which means $\text{Tr}\left(\left(\sqrt{\boldsymbol{Z}\boldsymbol{Z}^{T}}\right)_{\textit{PD}}\;\right)>0$, so $e^{\prime}>e_{\text{min}}$.

0.2.2 Case II: Uniform Scale

If we include a uniform scale, $s$, into our affine transform, we then have

[TABLE]

So our cost function becomes

[TABLE]

and our resulting optimization problem is

[TABLE]

The Lagrangian for the above (constrained) optimization problem is

[TABLE]

Optimizing Eq. (15) with respect to ${\bf t}$ yields

[TABLE]

Substituting Eq.(16) back into our Lagrangian and then partial differentiating it with respect to ${\bf L}$ yields

[TABLE]

which solving for ${\bf L}$ yields

[TABLE]

Exploiting the fact that ${\bf L}$ is an orthogonal matrix, we obtain

[TABLE]

which is the same rotation matrix we obtained in the non-scalar case. Thus, the inclusion of a fixed scale into our affine transform does not affect the rotation matrix.

Substituting both Eqs. (16) and (17) in Eq. (15) results in $\mathcal{L}({\bf L},{\bf t},s,\boldsymbol{\alpha},\lambda)=\mathcal{L}(s,\boldsymbol{\alpha},\lambda)$. Optimizing the Lagrangian with respect to $s$ yields

[TABLE]

Note that

[TABLE]

where

[TABLE]

In other words, $\sigma_{u^{j}}^{2}$ is the variance of feature $u^{j}$. Solving Eq. (18) for $s$ yields

[TABLE]

Since both the numerator and denominator of Eq. (20) are always positive, it is always guaranteed that $\hat{s}>0$; this is reasonable because a negative scale defies physical interpretation.

So the optimal alignment parameters are

[TABLE]

Or equivalently,

[TABLE]

The minimum squared error is found by substituting the optimal alignment parameters back into Eq. (14), which yields

[TABLE]

0.2.3 Uncoupled Weights: An Ill-Posed Problem

We now examine a specific case that makes the minimization problem ill-posed, i.e. no solution exists. Assume that $m_{ik}$ is separable, i.e. $m_{ik}=\sigma_{i}\gamma_{k}$, which means that

[TABLE]

Consequently, the cross-covariance matrix becomes

[TABLE]

As a result, an infinite number of rotation matrices are admissible. Thus, a unique solution for the minimization problem is guaranteed if and only if the weight term is coupled between the two point-sets. A special sub-case is when the weight term is fixed for all possible pairings, i.e. $m_{ik}=c$ for all $i$ and $k$, where $c$ is a real constant.

0.2.4 Relation to Labeled Scenario

In the case where the correspondence is known a priori, we have $m_{ik}=w_{i}\gamma_{ik}$, where

[TABLE]

So our cost function reduces to

[TABLE]

where $N\leq\min\{N_{U},N_{V}\}$ is the number of matching pairs between the two sets and $w_{n}$ is the ”strength” of the match between the two points forming pair $n$. The solution to Eq. (22) is the same as that derived in unlabeled scenario, except that now the (weighted) averages, variances, and covariances are no longer coupled. Eq. (22) is the cost function minimized in [1, 2, 3, 4].

0.3 Numerical Implementation

We have derived the closed-form solutions for the optimal alignment parameters for two different cases: the absence of scale and the presence of a uniform scale. Without loss of generality, the present discussion will focus on the case of uniform scale.

Denote ${\bf M}$ as the data structure (say an array) storing all $N_{U}N_{V}$ potential pairings. After alignment, the Euclidean distance between any two (cross) points forming a pair is

[TABLE]

If a pair constitutes a genuine match, then ideally $\Delta_{ik}$ will be small, and we will want to keep it. And if the pair is spurious, then $\Delta_{ik}$ will be large, and we will want to discard it. To do that, we need to compare each pair’s $\Delta_{ik}$ to some threshold, $T$. If for a given pair,

[TABLE]

then the pair is an outlier and we remove it from the array ${\bf M}$. Otherwise, we recompute its weight as

[TABLE]

Upon iterating across every pair in ${\bf M}$, we count the number of pairs left in the array. If no pairs have been removed, this is indicative of the threshold being too large, so we repeat Step 2 using the threshold $T:=T-\epsilon$, where $0<\epsilon<T$.

If pairs have been removed, then we need to check if the convergence criterion has been met. We define convergence as when the number of pairs does not exceed the maximum number of allowable one-to-one matches, i.e. $\text{length}({\bf M})<\min{\left(N_{U},N_{V}\right)}$. If this happens to be the case, then we are done; the remaining pairs in ${\bf M}$ form the optimal matching pairs. Otherwise, we repeat Stage II using the updated ${\bf M}$ and new weights.

Note that the algorithm consists of two hyper-parameters: the max threshold, $T$, and the threshold spacing, $\epsilon$. They have to be tuned appropriately with respect to the data.

We now summarize the proposed algorithm as a pseudocode:

1:  while number of pairs $>\min(N_{U},N_{V})$ and $T>0$ do

2:     Align query object to template

3:     for each pair $M_{j}$ in array M do

4:        compute weighted sum, $\Delta_{j}$, of radial displacements

5:        if $\Delta_{j}>T$ then

6:            remove pair $M_{j}$ from array M

7:        else

8:           compute new weight of $M_{j}$

9:        end if

10:     end for

11:     if no pairs are removed then

12:        $T:=T-\epsilon$

13:     end if

14:  end while

0.3.1 Computational Complexity

Denote $I(T,\epsilon)$ as the number of iterations until convergence is reached; it is a function of the two hyper-paramters. The computational complexity of the proposed algorithm is then $\mathcal{O}(IN_{U}N_{V})$. In the best case scenario, the algorithm will achieve convergence after a single iteration, i.e. $I=1$.

Let us look at the scenario where no pairs will be removed during each iteration of alignment. Such a scenario is the worst-case from a practical standpoint, but not necessarily from an algorithmic (i.e. computational complexity) perspective. In such a scenario, the total number of attempted alignments is $I(T,\epsilon)=\left\lceil{\frac{T}{\epsilon}}\right\rceil$, so the complexity is $\mathcal{O}(\left\lceil{\frac{T}{\epsilon}}\right\rceil N_{U}N_{V})$. Note, a worst computational complexity than this could be achieved; for example, having more than $\left\lceil{\frac{T}{\epsilon}}\right\rceil$ iterations until convergence is reached is completely plausible because any given threshold may be used more than once.

0.3.2 Similarity Score

Based on the set of (optimal) matching point pairs outputted by the proposed algorithm, a similarity score measuring the strength of the match between the two point-sets can be computed from the minimum squared error generated by the optimal alignment parameters, which recall for the case of uniform scale is

[TABLE]

It would ideally be small for genuine matches and large for dissimilar point-sets.

0.4 Discussion

An unsupervised, iterative $N$-dimensional point-set registration algorithm for unlabeled data (i.e. correspondence between points is unknown) and based on linear least squares has been proposed. The algorithm considers all possible point pairings and iteratively aligns the two sets until the number of point pairs does not exceed the maximum number of allowable one-to-one pairings.

Note that the output the of algorithm, i.e. the optimal matching point-pairs, may not necessarily be injective. In fact, a point in ${\bf U}$ may match to more than one point in ${\bf V}$, or more than one point in ${\bf U}$ may match to the same point in ${\bf V}$. Such a situation frequently arises in fingerprint matching where due to image noise and/or faulty image processing, a minutia in the query image may genuinely correspond to more than one minutia in the reference image.

The proposed algorithm may be utilized in a wide variety of matching problems, beyond just fingerprint matching. For example, assuming ${\bf U}$ and ${\bf V}$ describe two different sets of individuals, we may be interested in identifying those pairs of individuals that are most likely to become friends. Another example is where ${\bf U}$ describing a set of individuals and ${\bf V}$ a set of companies, and we seek to find which company a given individual best matches to. In these two examples, injectivity may not be a realistic assumption; there may be more than one company that a given individual would fit well in, or there may be one company where more than one individual would be a great fit.