An unsupervised, iterative $N$-dimensional point-set registration algorithm
UDC 517.9 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 is proposed. The algorithm considers all possible point pairings and iteratively aligns the two sets until the n...
Gespeichert in:
| Datum: | 2022 |
|---|---|
| Hauptverfasser: | , , |
| Format: | Artikel |
| Sprache: | Englisch |
| Veröffentlicht: |
Institute of Mathematics, NAS of Ukraine
2022
|
| Online Zugang: | https://umj.imath.kiev.ua/index.php/umj/article/view/6969 |
| Tags: |
Tag hinzufügen
Keine Tags, Fügen Sie den ersten Tag hinzu!
|
| Назва журналу: | Ukrains’kyi Matematychnyi Zhurnal |
| Завантажити файл: |  |
Institution
Ukrains’kyi Matematychnyi Zhurnal| _version_ | 1860512571054882816 |
|---|---|
| author | Hosseinbor, P. Zhdanov, R. Ushveridze, A. Hosseinbor, P. Zhdanov, R. Ushveridze, A. |
| author_facet | Hosseinbor, P. Zhdanov, R. Ushveridze, A. Hosseinbor, P. Zhdanov, R. Ushveridze, A. |
| author_sort | Hosseinbor, P. |
| baseUrl_str | https://umj.imath.kiev.ua/index.php/umj/oai |
| collection | OJS |
| datestamp_date | 2025-03-31T08:44:52Z |
| description | UDC 517.9
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 is 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.
  |
| doi_str_mv | 10.37863/umzh.v74i3.6969 |
| first_indexed | 2026-03-24T03:30:54Z |
| format | Article |
| fulltext |
DOI: 10.37863/umzh.v74i3.6969
UDC 517.9
P. Hosseinbor (BIO-key International, Eagan, MN, USA),
R. Zhdanov (CyberOptics, Minneapolis, MN, USA),
A. Ushveridze (Algostream, Plymouth, MN, USA)
AN UNSUPERVISED, ITERATIVE \bfitN -DIMENSIONAL POINT-SET
REGISTRATION ALGORITHM
НЕКОНТРОЛЬОВАНИЙ IТЕРАЦIЙНИЙ \bfitN -ВИМIРНИЙ АЛГОРИТМ
РЕЄСТРАЦIЇ ТОЧОК
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 is 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.
Запропоновано неконтрольований iтерацiйний N -вимiрний алгоритм реєстрацiї точок для непомiчених даних (тобто
спiввiдношення мiж точками невiдоме) на основi лiнiйного методу найменших квадратiв. Алгоритм використовує
всi можливi парування точок та iтерацiйно вирiвнює цi двi множини поки кiлькiсть точок не досягне максимуму
однозначних парувань.
1. Introduction. Point-set registration seeks to align an N-dimensional query point-set (or point
cloud) to some N-dimensional reference reference point-set via some mathematical transform, and
has been extensively studied in computer vision. For labeled data, an important class of solutions
is linear least squared techniques [1, 2, 6, 11]. Unlike the aforementioned work of [1, 6, 11], the
paper in [2] 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 corre-
spondence and the optimal affine transform that minimizes some dissimilarity metric between the
two point-sets need to be determined. One of the most important class of solutions for the unlabeled
case is iterative closest point (ICP).
Various point pattern matching algorithms for unlabeled data have been proposed in [3, 5, 7 – 10],
but the necessary optimization for each is numerically based. The works of [3, 8, 9] determined the
optimal affine transform and correspondence simultaneously by numerically solving a constrained
least squares problem. The method in [5] models each point-set as a Gaussian mixture, and deter-
mines 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 ef-
ficient 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, 6, 11] 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
c\bigcirc P. HOSSEINBOR, R. ZHDANOV, A. USHVERIDZE, 2022
ISSN 1027-3190. Укр. мат. журн., 2022, т. 74, № 3 427
428 P. HOSSEINBOR, R. ZHDANOV, A. USHVERIDZE
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 \BbbR N ), and the obtained solutions are then used to create an un-
supervised, iterative N-dimensional point-set registration algorithm. The N = 2 case was shown in
[4], and utilized in the context of fingerprint matching.
The iterative closest point algorithm (ICP) is probably the most closely related method in the lit-
erature to our’s because both are iterative, pertinent to unlabeled data, and employ linear least squares
to estimate the optimal alignment parameters. However, there is are major differences between the
two: ICP inherently assumes that correspondence between points can be inferred from the distances
between points.
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.
2. Theory. Consider two N-dimensional point sets \bfU and \bfV comprising NU and NV singular
points, respectively. We refer to \bfU and \bfV as the query and reference point-sets, respectively. The
Cartesian coordinates of the singular points will be expressed as an N-dimensional vector:
\bfu i \in \BbbR N \in \bfU , i = 1, . . . , NU ,
\bfv k \in \BbbR N \in \bfV , k = 1, . . . , NV .
The elements of \bfu i and \bfv k are denoted as uji and vjk, respectively, where j = 1, . . . , N. \bfU and \bfV
can be interpreted as matrices whose columns are formed by the vectors \bfu i and \bfv k, respectively,
i.e., \bfU \in \BbbR NU\times N and \bfV \in \BbbR NV \times N ; and uj and vj can be interpreted as the features of \bfU and \bfV ,
respectively.
We want to register point-set \bfU to \bfV . There are NUNV possible matching (cross) pairs and at
most \mathrm{m}\mathrm{i}\mathrm{n}\{ NU , NV \} one-to-one matching pairs. Let mik denote the weight of a matching pair; the
weight can be interpreted as a probability that the points \bfu i and \bfv k match locally. We assume that
the (initial) mik for each (cross) pair has been determined prior to alignment. One way of computing
mik is discussed in [4] within the context of fingerprint matching.
2.1. Case I: No Scale. We apply a global rotation and translation to point set \bfU such that
\bfu \prime
i = \bfL \bfu i + \bft , (1)
where L is the N \times N rotation matrix and \bft is the N-d vector of translation parameters. Since \bfL is
a rotation matrix, it is orthogonal, i.e., \bfL \bfL T = \bfL T\bfL = \bfI N\times N , and has a determinant of 1.
The measure of closeness of the transformed point set \bfU and the template set \bfV may be defined
as the weighted sum of the squared distances between their points (i.e., Euclidean distance metric):
e(\bfU ,\bfV ;\bfL , \bft ) =
\sum NU
i=1
\sum NV
k=1
mik(\bfu
\prime
i - \bfv k)
T (\bfu \prime
i - \bfv k)\sum NU
i=1
\sum NV
k=1
mik
. (2)
Interpreting mik as the probability of a match between \bfu i and \bfv k leads to
\sum NU
i=1
\sum NV
k=1
mik = 1.
Expanding out the product yields
e(\bfU ,\bfV ;\bfL , \bft ) =
NU\sum
i=1
NV\sum
k=1
mik
\Bigl[
\bfv T
k \bfv k - 2\mathrm{T}\mathrm{r}
\bigl(
\bfL \bfu i\bfv
T
k
\bigr)
+ 2\mathrm{T}\mathrm{r}
\bigl(
\bfL \bfu i\bft
T
\bigr)
+ \bfu T
i \bfu i - 2\bft T\bfv k + \bft T \bft
\Bigr]
, (3)
ISSN 1027-3190. Укр. мат. журн., 2022, т. 74, № 3
AN UNSUPERVISED, ITERATIVE N-DIMENSIONAL POINT-SET REGISTRATION ALGORITHM 429
where “Tr” denotes the trace operator. We seek the optimal values of parameters \bfL and \bft that
minimize e(\bfU ,\bfV ;\bfL , \bft ), so the constrained optimization problem to be solved is
minimize
\bfL ,\bft
e(\bfU ,\bfV ;\bfL , \bft )
subject to \bfL \bfL T = \bfI N\times N ,
\mathrm{d}\mathrm{e}\mathrm{t}\bfL = 1.
The Lagrangian for the above (constrained) optimization problem is
\scrL (\bfL , \bft ,\bfitalpha , \lambda ) = e(\bfU ,\bfV ;\bfL , \bft ) + \mathrm{T}\mathrm{r}
\bigl(
\bfitalpha (\bfL \bfL T - \bfI N\times N )
\bigr)
+ \lambda (\mathrm{d}\mathrm{e}\mathrm{t}\bfL - 1) , (4)
where \bfitalpha is an N \times N symmetric matrix of Lagrangian multipliers (it is symmetric because \bfL \bfL 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
\bfitu =
NU\sum
i=1
NV\sum
k=1
mik\bfu i,
\bfitv =
NU\sum
i=1
NV\sum
k=1
mik\bfv k.
Optimizing Eq. (4) with respect to \bft yields
\^\bfitt = \bfitv - \bfitL \bfitu . (5)
Substituting Eq. (5) back into our Lagrangian and then partial differentiating it with respect to \bfL
yields
\partial \scrL
\partial \bfL
= 2
\Biggl(
\bfitv \bfitu T -
NU\sum
i=1
NV\sum
k=1
mik\bfv k\bfu
T
i +\bfitalpha \bfitL
\Biggr)
+ \lambda \bfL = \bfzero N\times N , (6)
where \bfzero N\times N denotes an N \times N matrix of zeros. Let
\bfZ =
NU\sum
i=1
NV\sum
k=1
mik\bfv k\bfu
T
i - \bfitv \bfitu T ,
which is a (N \times N) cross-covariance matrix that we assume to be invertible. Specifically,
\bfZ =
\left( \mathrm{c}\mathrm{o}\mathrm{v}(u1, v1) \mathrm{c}\mathrm{o}\mathrm{v}(u1, v2) \mathrm{c}\mathrm{o}\mathrm{v}(u1, v3) \cdot \cdot \cdot \mathrm{c}\mathrm{o}\mathrm{v}(u1, vN )
...
...
...
...
...
\mathrm{c}\mathrm{o}\mathrm{v}(uN , v1) \mathrm{c}\mathrm{o}\mathrm{v}(uN , v2) \mathrm{c}\mathrm{o}\mathrm{v}(uN , v3) \cdot \cdot \cdot \mathrm{c}\mathrm{o}\mathrm{v}(uN , vN )
\right) ,
where \mathrm{c}\mathrm{o}\mathrm{v}(uj , vj
\prime
) denotes the covariance between features uj and vj
\prime
:
\mathrm{c}\mathrm{o}\mathrm{v}(uj , vj
\prime
) =
NU\sum
i=1
NV\sum
k=1
miku
j
iv
j\prime
k -
\Biggl(
NU\sum
i=1
NV\sum
k=1
miku
j
i
\Biggr) \Biggl(
NU\sum
i=1
NV\sum
k=1
mikv
j\prime
k
\Biggr)
. (7)
Solving Eq. (6) for \bfL yields
\^\bfitL =
\biggl(
\bfitalpha +
\lambda
2
\bfI N\times N
\biggr) - 1
\bfZ .
Since \bfL is an orthogonal matrix, we have
ISSN 1027-3190. Укр. мат. журн., 2022, т. 74, № 3
430 P. HOSSEINBOR, R. ZHDANOV, A. USHVERIDZE
\bfL \bfL T =
\biggl(
\bfitalpha +
\lambda
2
\bfI N\times N
\biggr) - 1
\bfitZ \bfitZ T
\biggl(
\bfitalpha +
\lambda
2
\bfI N\times N
\biggr) - 1
= \bfI N\times N ,
which upon isolating the cross-covariance, \bfZ , yields
\biggl(
\bfitalpha +
\lambda
2
\bfI N\times N
\biggr) 2
= \bfitZ \bfitZ T , so
\^\bfitL = \pm
\Bigl( \sqrt{}
\bfitZ \bfitZ T
\Bigr) - 1
\bfitZ . (8)
The minimum of the error function corresponds to
\Bigl( \surd
\bfitZ \bfitZ T
\Bigr) - 1
\bfitZ , so we discard the solution
-
\Bigl( \surd
\bfitZ \bfitZ T
\Bigr) - 1
\bfitZ .
The matrix \bfitZ \bfitZ T exhibits several interesting properties:
1. It is real symmetric. It is real because the elements of \bfitZ are covariances, which are always
real by definition of covariance. And it is symmetric because \bfitZ \bfitZ T =
\bigl(
\bfitZ \bfitZ T
\bigr) T
.
2. It is diagonalizable as a consequence of being real symmetric. We will exploit this fact later.
3. It is positive-definite, which we now prove. First,
\bfx T\bfitZ \bfitZ T\bfx = (\bfitZ T\bfx )T\bfitZ T\bfx = \| \bfitZ T\bfx \| 2 \geq 0 \forall \bfx \in \BbbR N \setminus \bfzero ,
which establishes that \bfitZ \bfitZ T is inherently positive semidefinite. Now recall that we further assume
that \bfitZ is invertible, which means \bfitZ \bfx = \bfzero if and only if \bfx = \bfzero . Thus, \| \bfitZ T\bfx \| 2 > 0 \forall \bfx \in \BbbR N \setminus \bfzero .
4. It’s eigenvalues are all positive as a consequence of positive-definiteness.
Eq. (8) requires taking the square root of the matrix \bfitZ \bfitZ T . According to the spectral theorem,
any real symmetric matrix can be diagonalized by an orthogonal matrix:
\bfitZ \bfitZ T = \bfP \bfD \bfP T , (9)
where \bfP is an N \times N orthogonal matrix and \bfD is the N \times N diagonal matrix; \bfP and \bfD are
formed by the eigenvectors and eigenvalues, respectively, of \bfitZ \bfitZ T . The square root of \bfitZ \bfitZ T is then
\bfP \bfD
1
2\bfP T because \Bigl(
\bfP \bfD
1
2\bfP T
\Bigr) 2
= \bfP \bfD
1
2\bfP T\bfP \bfD
1
2\bfP T = \bfP \bfD \bfP T = \bfitZ \bfitZ T .
Now, any N \times N matrix with N distinct eigenvalues has 2N 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 \bfitZ \bfitZ T to extract only its positive-definite square-root, so \bfD
1
2 is formed by the positive square roots
of the eigenvalues of \bfitZ \bfitZ T . To this end, we employ the notation
\bigl( \surd
arg
\bigr)
PD
to refer to the positive-
definite square root of any square matrix. Taking all this into account, Eq. (8) can be expressed as
\^\bfitL = \bfP
\bigl( \surd
\bfD
\bigr) - 1
PD
\bfP T\bfitZ . (10)
In summary, the optimal alignment parameters are
\^\bfitL =
\Bigl( \surd
\bfitZ \bfitZ T
\Bigr) - 1
PD
\bfitZ ,
\^\bfitt = \bfitv -
\Bigl( \surd
\bfitZ \bfitZ T
\Bigr) - 1
PD
\bfitZ \bfitu .
ISSN 1027-3190. Укр. мат. журн., 2022, т. 74, № 3
AN UNSUPERVISED, ITERATIVE N-DIMENSIONAL POINT-SET REGISTRATION ALGORITHM 431
Or equivalently,
\^\bfitL = \bfP
\Bigl( \surd
\bfD
\Bigr) - 1
PD
\bfP T\bfitZ ,
\^\bfitt = \bfitv - \bfP
\Bigl( \surd
\bfD
\Bigr) - 1
PD
\bfP T\bfitZ \bfitu .
The minimum squared error is found by substituting the optimal alignment parameters back into
Eq. (3), which yields
emin =
\Biggl(
NU\sum
i=1
NV\sum
k=1
mik\bfu
T
i \bfu i - \bfitu T \bfitu
\Biggr)
+
+
\Biggl(
NU\sum
i=1
NV\sum
k=1
mik\bfv
T
i \bfv i - \bfitv T \bfitv
\Biggr)
- 2\mathrm{T}\mathrm{r}
\Bigl( \Bigl( \sqrt{}
\bfitZ \bfitZ T
\Bigr)
PD
\Bigr)
.
Note that
NU\sum
i=1
NV\sum
k=1
mik\bfu
T
i \bfu i - \bfitu T \bfitu =
N\sum
j=1
\sigma 2
uj , (11)
where
\sigma 2
uj =
NU\sum
i=1
NV\sum
k=1
mik
\bigl(
uji
\bigr) 2 - \Biggl( NU\sum
i=1
NV\sum
k=1
miku
j
i
\Biggr) 2
.
In other words, \sigma 2
uj is the variance of feature uj , and likewise \sigma 2
vj
is the variance of feature vj . So
we can rewrite the minimum squared error as
emin =
N\sum
j=1
\bigl(
\sigma 2
uj + \sigma 2
vj
\bigr)
- 2\mathrm{T}\mathrm{r}
\Bigl( \Bigl( \sqrt{}
\bfitZ \bfitZ T
\Bigr)
PD
\Bigr)
=
=
N\sum
j=1
\bigl(
\sigma 2
uj + \sigma 2
vj
\bigr)
- 2\mathrm{T}\mathrm{r}
\Bigl( \bigl( \surd
\bfD
\bigr)
PD
\Bigr)
, (12)
where we exploited the fact that \mathrm{T}\mathrm{r}
\Bigl(
\bfP
\bigl( \surd
\bfD
\bigr)
PD
\bfP T
\Bigr)
= \mathrm{T}\mathrm{r}
\Bigl( \bigl( \surd
\bfD
\bigr)
PD
\Bigr)
.
Note that had we used \^\bfitL = -
\Bigl( \surd
\bfitZ \bfitZ T
\Bigr) - 1
PD
\bfitZ , the generated error would be
e\prime =
N\sum
j=1
\bigl(
\sigma 2
uj + \sigma 2
vj
\bigr)
+ 2\mathrm{T}\mathrm{r}
\Bigl( \Bigl( \sqrt{}
\bfitZ \bfitZ T
\Bigr)
PD
\Bigr)
.
The eigenvalues of a positive-definite matrix all always positive, which means \mathrm{T}\mathrm{r}
\Bigl( \Bigl( \surd
\bfitZ \bfitZ T
\Bigr)
PD
\Bigr)
>
> 0, so e\prime > emin.
2.2. Case II: Uniform Scale. If we include a uniform scale, s, into our affine transform, we
then have
\bfu \prime
i = s\bfL \bfu i + \bft . (13)
So our cost function becomes
ISSN 1027-3190. Укр. мат. журн., 2022, т. 74, № 3
432 P. HOSSEINBOR, R. ZHDANOV, A. USHVERIDZE
e(\bfU ,\bfV ;\bfL , \bft , s) =
NU\sum
i=1
NV\sum
k=1
mik
\Bigl[
\bfv T
k \bfv k - 2s \mathrm{T}\mathrm{r}
\bigl(
\bfL \bfu i\bfv
T
k
\bigr)
+
+2s \mathrm{T}\mathrm{r}
\bigl(
\bfL \bfu i\bft
T
\bigr)
+ s2\bfu T
i \bfu i - 2\bft T\bfv k + \bft T \bft
\Bigr]
(14)
and our resulting optimization problem is
minimize
\bfL ,\bft ,s
e(\bfU ,\bfV ;\bfL , \bft , s)
subject to \bfL \bfL T = \bfI N\times N ,
\mathrm{d}\mathrm{e}\mathrm{t}\bfL = 1.
The Lagrangian for the above (constrained) optimization problem is
\scrL (\bfL , \bft , s,\bfitalpha , \lambda ) = e(\bfU ,\bfV ;\bfL , \bft , s) + \mathrm{T}\mathrm{r}
\bigl(
\bfitalpha (\bfL \bfL T - \bfI N\times N )
\bigr)
+ \lambda (\mathrm{d}\mathrm{e}\mathrm{t}\bfL - 1) . (15)
Optimizing Eq. (15) with respect to \bft yields
\^\bfitt = \bfitv - s\bfitL \bfitu . (16)
Substituting Eq. (16) back into our Lagrangian and then partial differentiating it with respect to \bfL
yields
\partial \scrL
\partial \bfL
= 2 (\bfitalpha \bfitL - s\bfZ ) + \lambda \bfL = \bfzero N\times N ,
which solving for \bfL yields
\^\bfitL = s
\biggl(
\bfitalpha +
\lambda
2
\bfI N\times N
\biggr) - 1
\bfZ .
Exploiting the fact that \bfL is an orthogonal matrix, we obtain
\^\bfitL =
\Bigl( \sqrt{}
\bfitZ \bfitZ T
\Bigr) - 1
PD
\bfitZ , (17)
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) into the Lagrangian, Eq. (15), and optimizing it with respect
to s yields
\partial \scrL
\partial s
= 2
\Biggl(
s
\Biggl(
NU\sum
i=1
NV\sum
k=1
mik\bfu
T
i \bfu i - \bfitu T \bfitu
\Biggr)
- \mathrm{T}\mathrm{r}
\Bigl( \Bigl( \sqrt{}
\bfitZ \bfitZ T
\Bigr)
PD
\Bigr) \Biggr)
= 0. (18)
Note that
NU\sum
i=1
NV\sum
k=1
mik\bfu
T
i \bfu i - \bfitu T \bfitu =
N\sum
j=1
\sigma 2
uj , (19)
where
\sigma 2
uj =
NU\sum
i=1
NV\sum
k=1
mik
\bigl(
uji
\bigr) 2 - \Biggl( NU\sum
i=1
NV\sum
k=1
miku
j
i
\Biggr) 2
.
ISSN 1027-3190. Укр. мат. журн., 2022, т. 74, № 3
AN UNSUPERVISED, ITERATIVE N-DIMENSIONAL POINT-SET REGISTRATION ALGORITHM 433
In other words, \sigma 2
uj is the variance of feature uj . Solving Eq. (18) for s yields
\^s =
\mathrm{T}\mathrm{r}
\Bigl( \Bigl( \surd
\bfitZ \bfitZ T
\Bigr)
PD
\Bigr)
\sum N
j=1
\sigma 2
uj
. (20)
Since both the numerator and denominator of Eq. (20) are always positive, it is always guaranteed
that \^s > 0; this is reasonable because a negative scale defies physical interpretation.
So the optimal alignment parameters are
\^\bfitL =
\Bigl( \surd
\bfitZ \bfitZ T
\Bigr) - 1
PD
\bfitZ ,
\^s =
\mathrm{T}\mathrm{r}
\Bigl( \Bigl( \surd
\bfitZ \bfitZ T
\Bigr)
PD
\Bigr)
\sum N
j=1
\sigma 2
uj
,
\^\bfitt = \bfitv -
\mathrm{T}\mathrm{r}
\Bigl( \Bigl( \surd
\bfitZ \bfitZ T
\Bigr)
PD
\Bigr) \Bigl( \surd
\bfitZ \bfitZ T
\Bigr) - 1
PD
\bfitZ \bfitu \sum N
j=1
\sigma 2
uj
.
Or equivalently,
\^\bfitL = \bfP
\bigl( \surd
\bfD
\bigr) - 1
PD
\bfP T\bfitZ ,
\^s =
\mathrm{T}\mathrm{r}
\Bigl( \bigl( \surd
\bfD
\bigr)
PD
\Bigr)
\sum N
j=1
\sigma 2
uj
,
\^\bfitt = \bfitv -
\mathrm{T}\mathrm{r}
\Bigl( \bigl( \surd
\bfD
\bigr)
PD
\Bigr)
\bfP
\bigl( \surd
\bfD
\bigr) - 1
PD
\bfP T\bfitZ \bfitu \sum N
j=1
\sigma 2
uj
.
The minimum squared error is found by substituting the optimal alignment parameters back into
Eq. (14), which yields
emin =
N\sum
j=1
\sigma 2
vj -
\Bigl[
\mathrm{T}\mathrm{r}
\Bigl( \Bigl( \surd
\bfitZ \bfitZ T
\Bigr)
PD
\Bigr) \Bigr] 2
\sum N
j=1
\sigma 2
uj
=
=
N\sum
j=1
\sigma 2
vj -
\Bigl[
\mathrm{T}\mathrm{r}
\Bigl( \bigl( \surd
\bfD
\bigr)
PD
\Bigr) \Bigr] 2
\sum N
j=1
\sigma 2
uj
. (21)
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 mik is separable, i.e.,
mik = \sigma i\gamma k, which means that
ISSN 1027-3190. Укр. мат. журн., 2022, т. 74, № 3
434 P. HOSSEINBOR, R. ZHDANOV, A. USHVERIDZE
\bfitu =
NU\sum
i=1
\sigma i\bfu i,
\bfitv =
NV\sum
k=1
\gamma k\bfv k.
Consequently, the cross-covariance matrix becomes
\bfZ =
NV\sum
k=1
\gamma k\bfv k
NU\sum
i=1
\sigma i\bfu
T
i - \bfitv \bfitu T = \bfzero N\times N .
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., mik = c for
all i and k, where c is a real constant.
2.4. Relation to Labeled Scenario. In the case where the correspondence is known a priori, we
have mik = wi\gamma ik, where
\gamma ik =
\left\{ 0, points i and k do not correspond,
1, points i and k do correspond.
So our cost function reduces to
D(\bfU ,\bfV ; a, b, sx, sy, \theta ) =
\sum N
n=1
wn(\bfu
\prime
n - \bfv n)
T (\bfu \prime
n - \bfv k)\sum N
n=1
wi
, (22)
where N \leq \mathrm{m}\mathrm{i}\mathrm{n}\{ NU , NV \} is the number of matching pairs between the two sets and wn 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, 6, 11].
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 \bfM as the data structure (say an array) storing all NUNV potential pairings. After
alignment, the Euclidean distance between any two (cross) points forming a pair is
\Delta ik =
\sqrt{} \bigl(
\^s\^\bfitL \bfu i +\^\bfitt - \bfv k
\bigr) T \bigl(
\^s\^\bfitL \bfu i +\^\bfitt - \bfv k
\bigr)
.
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,
\Delta ik > T,
then the pair is an outlier and we remove it from the array \bfM . Otherwise, we recompute its weight as
ISSN 1027-3190. Укр. мат. журн., 2022, т. 74, № 3
AN UNSUPERVISED, ITERATIVE N-DIMENSIONAL POINT-SET REGISTRATION ALGORITHM 435
mik = 1 - \Delta ik
T
. (23)
Upon iterating across every pair in \bfM , 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., length(\bfM ) < \mathrm{m}\mathrm{i}\mathrm{n} (NU , NV ). If this happens to be the case, then we are
done; the remaining pairs in \bfM form the optimal matching pairs. Otherwise, we repeat Stage II
using the updated \bfM 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 > \mathrm{m}\mathrm{i}\mathrm{n}(NU , NV ) and T > 0 do
2: Align query object to template
3: for each pair Mj in array M do
4: compute weighted sum, \Delta j , of radial displacements
5: if \Delta j > T then
6: remove pair Mj from array M
7: else
8: compute new weight of Mj
9: end if
10: end for
11: if no pairs are removed then
12: T := T - \epsilon
13: end if
14: end while
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-parameters. The computational complexity of the proposed
algorithm is then \scrO (INUNV ). 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 algorith-
mic (i.e., computational complexity) perspective. In such a scenario, the total number of attempted
alignments is I(T, \epsilon ) =
\biggl\lceil
T
\epsilon
\biggr\rceil
, so the complexity is \scrO
\biggl( \biggl\lceil
T
\epsilon
\biggr\rceil
NUNV
\biggr)
. Note, a worst computational
complexity than this could be achieved; for example, having more than
\biggl\lceil
T
\epsilon
\biggr\rceil
iterations until con-
vergence is reached is completely plausible because any given threshold may be used more than
once.
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
ISSN 1027-3190. Укр. мат. журн., 2022, т. 74, № 3
436 P. HOSSEINBOR, R. ZHDANOV, A. USHVERIDZE
emin =
N\sum
j=1
\sigma 2
vj -
\Bigl[
\mathrm{T}\mathrm{r}
\Bigl(
\bfP
\bigl( \surd
\bfD
\bigr)
PD
\bfP T
\Bigr) \Bigr] 2
\sum N
j=1
\sigma 2
uj
.
It would ideally be small for genuine matches and large for dissimilar point-sets.
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 of the algorithm, i.e., the optimal matching point-pairs, may not necessarily
be injective. In fact, a point in \bfU may match to more than one point in \bfV , or more than one point
in \bfU may match to the same point in \bfV . 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, from 3D com-
puter vision and graphics to face and fingerprint recognition. For example, assuming \bfU and \bfV
describe two different sets of individuals, we may be interested in identifying those pairs of indi-
viduals that are most likely to become friends. Another example is where \bfU describing a set of
individuals and \bfV 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.
References
1. K. S. Arun, T. S. Huang, S. D. Blostein, Least-squares fitting of two 3-D point sets, IEEE Trans. Pattern Anal.
Machine Intell., 9, 698 – 700 (1987).
2. S. Chang, F. Cheng, W. Hsu, G. Wu, Fast algorithm for point pattern matching: Invariant to translations, rotations
and scale changes, Pattern Recognition, 30, 311 – 320 (1997).
3. S. Gold, A. Rangarajan, C. P. Lu, S. Pappu, E. Mjolsness, New algorithms for 2d and 3d point matching: pose
estimation and correspondence, Pattern Recognition, 31, 1019 – 1031 (1998).
4. A. P. Hosseinbor, R. Zhdanov, A. Ushveridze, An unsupervised 2d point-set registration algorithm for unlabeled
feature points: Application to fingerprint matching, Pattern Recognition Letters, 100, 137 – 143 (2017).
5. B. Jian, B. C. Vemuri, Robust point set registration using gaussian mixture models, IEEE PAMI, 33, 1633 – 1645
(2010).
6. W. Kabsch, A solution for the best rotation to relate two sets of vectors, Acta Cryst., 32, 922 – 923 (1976).
7. S. Lan, Z. Guo, J. You, A non-rigid registration method with application to distorted fingerprint matching, Pattern
Recognition, 95, 48 – 57 (2019).
8. A. Rangarajan, H. Chui, F. L. Bookstein, The softassign procrustes matching algorithm, IPMI, 29 – 42 (1997).
9. A. Rangarajan, H. Chui, E. Mjolsness, S. Pappu, L. Davachi, P. S. Goldman-Rakic, J. S. Duncan, A robust point
matching algorithm for autoradiograph alignment, Med. Image Anal., 1, 379 – 398 (1997).
10. Y. Tsin, T. Kanade, A correlation-based approach to robust point set registration, ECCV, 558 – 569 (2004).
11. S. Umeyama, Least-squares estimation of transformation parameters between two point sets, IEEE Trans. Pattern
Anal. Machine Intell., 13, 376 – 380 (1991).
Received 28.10.21
ISSN 1027-3190. Укр. мат. журн., 2022, т. 74, № 3
|
| id | umjimathkievua-article-6969 |
| institution | Ukrains’kyi Matematychnyi Zhurnal |
| keywords_txt_mv | keywords |
| language | English |
| last_indexed | 2026-03-24T03:30:54Z |
| publishDate | 2022 |
| publisher | Institute of Mathematics, NAS of Ukraine |
| record_format | ojs |
| resource_txt_mv | umjimathkievua/3d/b049a0b2ef80687c292c9c317d3ac43d.pdf |
| spelling | umjimathkievua-article-69692025-03-31T08:44:52Z An unsupervised, iterative $N$-dimensional point-set registration algorithm An unsupervised, iterative $N$-dimensional point-set registration algorithm Hosseinbor, P. Zhdanov, R. Ushveridze, A. Hosseinbor, P. Zhdanov, R. Ushveridze, A. Unsupervised method point-registration algorithm least square optimization UDC 517.9 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 is 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. &nbsp; УДК 517.9 Неконтрольований ітераційний $n$-вимірний алгоритм реєстрації точок Запропоновано неконтрольований ітераційний $N$-вимірний алгоритм реєстрації точок для непомічених даних (тобто співвідношення між точками невідоме) на основі лінійного методу найменших квадратів. Алгоритм використовує всі можливі парування точок та ітераційно вирівнює ці дві множини поки кількість точок не досягне максимуму однозначних парувань. Institute of Mathematics, NAS of Ukraine 2022-04-26 Article Article application/pdf https://umj.imath.kiev.ua/index.php/umj/article/view/6969 10.37863/umzh.v74i3.6969 Ukrains’kyi Matematychnyi Zhurnal; Vol. 74 No. 3 (2022); 427-436 Український математичний журнал; Том 74 № 3 (2022); 427-436 1027-3190 en https://umj.imath.kiev.ua/index.php/umj/article/view/6969/9210 Copyright (c) 2022 Renat Zhdanov |
| spellingShingle | Hosseinbor, P. Zhdanov, R. Ushveridze, A. Hosseinbor, P. Zhdanov, R. Ushveridze, A. An unsupervised, iterative $N$-dimensional point-set registration algorithm |
| title | An unsupervised, iterative $N$-dimensional point-set registration algorithm |
| title_alt | An unsupervised, iterative $N$-dimensional point-set registration algorithm |
| title_full | An unsupervised, iterative $N$-dimensional point-set registration algorithm |
| title_fullStr | An unsupervised, iterative $N$-dimensional point-set registration algorithm |
| title_full_unstemmed | An unsupervised, iterative $N$-dimensional point-set registration algorithm |
| title_short | An unsupervised, iterative $N$-dimensional point-set registration algorithm |
| title_sort | unsupervised, iterative $n$-dimensional point-set registration algorithm |
| topic_facet | Unsupervised method point-registration algorithm least square optimization |
| url | https://umj.imath.kiev.ua/index.php/umj/article/view/6969 |
| work_keys_str_mv | AT hosseinborp anunsupervisediterativendimensionalpointsetregistrationalgorithm AT zhdanovr anunsupervisediterativendimensionalpointsetregistrationalgorithm AT ushveridzea anunsupervisediterativendimensionalpointsetregistrationalgorithm AT hosseinborp anunsupervisediterativendimensionalpointsetregistrationalgorithm AT zhdanovr anunsupervisediterativendimensionalpointsetregistrationalgorithm AT ushveridzea anunsupervisediterativendimensionalpointsetregistrationalgorithm AT hosseinborp unsupervisediterativendimensionalpointsetregistrationalgorithm AT zhdanovr unsupervisediterativendimensionalpointsetregistrationalgorithm AT ushveridzea unsupervisediterativendimensionalpointsetregistrationalgorithm AT hosseinborp unsupervisediterativendimensionalpointsetregistrationalgorithm AT zhdanovr unsupervisediterativendimensionalpointsetregistrationalgorithm AT ushveridzea unsupervisediterativendimensionalpointsetregistrationalgorithm |