Procrustes Analysis with Deformations: A Closed-Form Solution by Eigenvalue Decomposition

Generalized Procrustes Analysis (GPA) is the problem of bringing multiple shapes into a common reference by estimating transformations. GPA has been extensively studied for the Euclidean and affine transformations. We introduce GPA with deformable transformations, which forms a much wider and difficult problem. We specifically study a class of transformations called the Linear Basis Warps, which contains the affine transformation and most of the usual deformation models, such as the Thin-Plate Spline (TPS). GPA with deformations is a nonconvex underconstrained problem. We resolve the fundamental ambiguities of deformable GPA using two shape constraints requiring the eigenvalues of the shape covariance. These eigenvalues can be computed independently as a prior or posterior. We give a closed-form and optimal solution to deformable GPA based on an eigenvalue decomposition. This solution handles regularization, favoring smooth deformation fields. It requires the transformation model to satisfy a fundamental property of free-translations, which asserts that the model can implement any translation. We show that this property fortunately holds true for most common transformation models, including the affine and TPS models. For the other models, we give another closed-form solution to GPA, which agrees exactly with the first solution for models with free-translation. We give pseudo-code for computing our solution, leading to the proposed DefGPA method, which is fast, globally optimal and widely applicable. We validate our method and compare it to previous work on six diverse 2D and 3D datasets, with special care taken to choose the hyperparameters from cross-validation.

This work was supported by ANR via the TOPACS project (ANR-19-CE45-0015). We thank the authors of the public datasets which we could use in our experiments. We appreciate the valuable comments of the anonymous reviewers that have helped improving the quality of the manuscript.

Authors and Affiliations

1. ENCOV, TGI, Institut Pascal, UMR6602 CNRS, Université Clermont Auvergne, Clermont-Ferrand, France

   Fang Bai

2. Department of Clinical Research and Innovation, CHU de Clermont-Ferrand, ENCOV, TGI, Institut Pascal, UMR6602 CNRS, Université Clermont Auvergne, Clermont-Ferrand, France

   Adrien Bartoli

Communicated by Xavier Pennec.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

The Thin-Plate Spline

1.1 General Mapping Definition

The Thin-Plate Spline (TPS) is a nonlinear mapping driven by a set of chosen control centers. Let \({\varvec{c}}_1, {\varvec{c}}_2, \dots , {\varvec{c}}_l\), with \({\varvec{c}}_i \in {\mathbb {R}}^{d}\) \((d=2,3)\), be the l control centers of TPS. We present the TPS model in 3D; the 2D case can be analogously derived by removing the z coordinate and letting the corresponding TPS kernel function be \(\phi \left( r\right) = r^2 \log \left( r^2\right) \) (Bookstein 1997). Strictly speaking, the TPS only occurs in 2D, but is part of a general family of functions called the harmonic splines. For simplicity, we call the function TPS independently of the dimension.

We denote a point in 3D by \({\varvec{p}}^{\intercal } = \left[ x, y,z\right] \). The TPS function \(\tau ({\varvec{p}})\) is a \({\mathbb {R}}^3 \rightarrowtail {\mathbb {R}}\) mapping:

where \(\phi \left( \cdot \right) \) is the TPS kernel function. In 3D, we have \(\phi \left( r\right) = - \vert r \vert \).

We collect the parameters in a single vector \(\varvec{\eta }^{\intercal } = \left[ {\varvec{w}}^{\intercal }, {\varvec{a}}^{\intercal }\right] \), with the coefficient in the nonlinear part as \({\varvec{w}}^{\intercal } = \left[ w_1, w_2, \cdots , w_l\right] \), and that in the linear part as \({\varvec{a}}^{\intercal } = \left[ a_1, a_2, a_3, a_4\right] \). Let \(\tilde{\varvec{p}}^{\intercal } = \left[ x, y, z, 1\right] \). Collect the nonlinear components in a vector:

Then the TPS function compacts in the vector form as:

We define \(\tilde{\varvec{c}}_i^{\intercal } = \left[ {\varvec{c}}_i^{\intercal }, 1\right] \) \((i \in \left[ 1:l \right] )\), and \(\tilde{\varvec{C}} = \left[ \tilde{\varvec{c}}_1, \tilde{\varvec{c}}_2, \dots , \tilde{\varvec{c}}_l\right] \). The transformation parameters in the nonlinear part of the TPS model are constrained with respect to the control centers, i.e., \( \tilde{\varvec{C}} {\varvec{w}} = {\varvec{0}} \). To solve the parameters of TPS, the l control centers \({\varvec{c}}_1, {\varvec{c}}_2, \dots , {\varvec{c}}_l\), are mapped to l scalar outputs \(y_1, y_2, \dots , y_l\), by \(\tau ({\varvec{c}}_k) = y_k\) \((k \in \left[ 1:l \right] )\). We collect the outputs in a vector \({\varvec{y}}^{\intercal } = \left[ y_1, y_2, \dots , y_l\right] \).

1.2 Standard Form as a Linear Regression Model

While there are \(l+4\) parameters in \(\varvec{\eta }\), we will show in the following that the 4 constraints in \( \tilde{\varvec{C}} {\varvec{w}} = {\varvec{0}} \) can be absorbed in a reparameterization of the TPS model, by using l parameters only.

To proceed, we define a matrix:

where we use the shorthand \( \phi _{ij} = \phi \left( \Vert {\varvec{c}}_i - {\varvec{c}}_j \Vert \right) \). The original diagonal elements of \({\varvec{K}}_{\lambda }\) are all zeros, i.e., \(\phi _{11} = \phi _{22} = \cdots = \phi _{ll} = 0\). \(\lambda \) is an adjustable parameter acting as an internal smoothing parameter to improve the conditioning of the matrix.

The l TPS mappings of control centers \(\tau ({\varvec{c}}_k) = y_k\) \((k \in \left[ 1:l \right] )\), and the 4 parameter constraints \( \tilde{\varvec{C}} {\varvec{w}} = {\varvec{0}} \) can be collected in the following matrix form:

which admits a closed-form solution:

with the matrix \(\varvec{{\mathcal {E}}}_{\lambda } \in {\mathbb {R}}^{\left( l+4\right) \times l}\) defined as:

Therefore the TPS function writes:

Let \(\varvec{\beta } ({\varvec{p}}) = \varvec{{\mathcal {E}}}_{\lambda }^{\intercal } \begin{bmatrix} \varvec{\phi }_{{\varvec{p}}} \\ \tilde{\varvec{p}} \end{bmatrix}\). Then \( \tau ({\varvec{p}}) = {\varvec{y}}^{\intercal } \varvec{\beta } ({\varvec{p}}) \).

The 3D TPS warp model \({\mathcal {W}}\left( {\varvec{p}}\right) \,:\, {\mathbb {R}}^3 \rightarrowtail {\mathbb {R}}^3\) is obtained by stacking 3 TPS functions together, one for each dimension:

This is the standard form of the general warp model defined in Eq. (16), with \({\varvec{W}} = \left[ {\varvec{y}}_{x}, {\varvec{y}}_{y}, {\varvec{y}}_{z} \right] \in {\mathbb {R}}^{l \times 3}\) being a set of new transformation parameters that are free of constraints.

1.3 Regularization by Bending Energy Matrix

Let \(\varvec{\bar{{\mathcal {E}}}}_{\lambda } = {\varvec{K}}_{\lambda }^{-1} - {\varvec{K}}_{\lambda }^{-1} \tilde{\varvec{C}}^{\intercal } \left( \tilde{\varvec{C}} {\varvec{K}}_{\lambda }^{-1} \tilde{\varvec{C}}^{\intercal }\right) ^{-1} \tilde{\varvec{C}} {\varvec{K}}_{\lambda }^{-1}\). The \(l \times l\) matrix \(\varvec{\bar{{\mathcal {E}}}}_{\lambda }\) is the bending energy matrix of the TPS warp (Bookstein 1989), which satisfies \( \varvec{\bar{{\mathcal {E}}}}_{\lambda } {\varvec{K}}_{\lambda } \varvec{\bar{{\mathcal {E}}}}_{\lambda } = \varvec{\bar{{\mathcal {E}}}}_{\lambda } \). The bending energy matrix \(\varvec{\bar{{\mathcal {E}}}}_{\lambda }\) is positive semidefinite, with rank \(l - 4\). The bending energy of a single TPS function is proportional to \({\varvec{w}}^{\intercal } {\varvec{K}}_{\lambda } {\varvec{w}}\), with \({\varvec{w}} = \varvec{\bar{{\mathcal {E}}}}_{\lambda } {\varvec{y}}\), which is given by:

For the 3D TPS warp model, the overall bending energy is:

where \(\sqrt{ \varvec{\bar{{\mathcal {E}}}}_{\lambda } }\) is the matrix square root of \(\varvec{\bar{{\mathcal {E}}}}_{\lambda }\). Then the bending energy of the TPS warp writes into the standard form \( {\mathcal {R}} ({\varvec{W}}) = \Vert \sqrt{ \varvec{\bar{{\mathcal {E}}}}_{\lambda } } {\varvec{W}} \Vert _F^2 \), given in Eq. (18).

1.4 The Thin-Plate Spline Satisfies \(\varvec{{\mathcal {B}}}\left( \cdot \right) ^{\intercal } {\varvec{x}} = {\varvec{1}}\) and \({\varvec{Z}} {\varvec{x}} = {\varvec{0}}\)

Let an arbitrary point cloud \({\varvec{D}}\) be \( {\varvec{D}} = \left[ {\varvec{p}}_1,\, {\varvec{p}}_2, \dots , {\varvec{p}}_m \right] \). Then by Eq. (33) and the point-cloud operation defined in Sect. 5.1.2, we write:

Hence \(\varvec{{\mathcal {B}}}({\varvec{D}})^{\intercal } \) has structure:

1.4.1 \(\varvec{{\mathcal {B}}}(\cdot )^{\intercal } {\varvec{x}} = {\varvec{1}}\)

It suffices to show such an \({\varvec{x}}\) satisfies \( \varvec{{\mathcal {E}}}_{\lambda } {\varvec{x}} = \begin{bmatrix} {\varvec{0}} \\ 1 \end{bmatrix} \). Because:

such an \({\varvec{x}}\) indeed exists, which is \({\varvec{x}} = \tilde{\varvec{C}}^{\intercal } \begin{bmatrix} {\varvec{0}} \\ 1 \end{bmatrix} \).

1.4.2 \({\varvec{Z}} {\varvec{x}} = {\varvec{0}}\)

We need to show:

The proof is immediate by noting that \({\varvec{Z}}^{\intercal } {\varvec{Z}} = \varvec{\bar{{\mathcal {E}}}}_{\lambda } \), and \( \varvec{\bar{{\mathcal {E}}}}_{\lambda } {\varvec{x}} = \varvec{\bar{{\mathcal {E}}}}_{\lambda } \tilde{\varvec{C}}^{\intercal } \begin{bmatrix} {\varvec{0}} \\ 1 \end{bmatrix} = {\varvec{O}} \begin{bmatrix} {\varvec{0}} \\ 1 \end{bmatrix} = {\varvec{0}} \).

1.5 The Thin-Plate Spline is Invariant to the Coordinate Transformation

Given a datum shape \({\varvec{D}}\), we apply the same rigid transformation \(({\varvec{R}},\, {\varvec{t}})\) to each datum point \({\varvec{p}}_i\), and denote the transformed datum point as \({\varvec{p}}_i' = {\varvec{R}} {\varvec{p}}_i + {\varvec{t}}\). Each control point \({\varvec{c}}_i\) is transformed to \( {\varvec{c}}_i' = {\varvec{R}} {\varvec{c}}_i + {\varvec{t}} \). The matrix \({\varvec{K}}_{\lambda }\) is invariant by replacing \({\varvec{c}}_i\) with \({\varvec{c}}_i'\) as its elements are radial basis functions. By replacing \(\tilde{\varvec{C}}\) by \(\tilde{\varvec{C}}' = \left[ \tilde{\varvec{c}}_1', \tilde{\varvec{c}}_2', \dots , \tilde{\varvec{c}}_l'\right] \), the matrix \(\varvec{{\mathcal {E}}}_{\lambda }\) becomes:

\(\varvec{\phi }_{{\varvec{p}}_i}\) is also invariant by simultaneously replacing \({\varvec{c}}_i\) with \({\varvec{c}}_i'\) and \({\varvec{p}}_i\) with \({\varvec{p}}_i'\). Notice that \(\tilde{\varvec{p}}_i' = {\varvec{T}} \tilde{\varvec{p}}_i\), thus:

This proves \(\varvec{{\mathcal {B}}}\left( {\varvec{D}}'\right) = \varvec{{\mathcal {B}}}\left( {\varvec{D}}\right) \) where \({\varvec{D}}' = {\varvec{R}} {\varvec{D}} + {\varvec{t}} {\varvec{1}}^{\intercal }\). Denote \(\varvec{\bar{{\mathcal {E}}}}_{\lambda }'\) to be the top \(l \times l\) sub-block of \(\varvec{{\mathcal {E}}}_{\lambda }'\). It is obvious that \(\varvec{\bar{{\mathcal {E}}}}_{\lambda }' = \varvec{\bar{{\mathcal {E}}}}_{\lambda }\) thus the regularization matrix \({\varvec{Z}}' = \sqrt{ \varvec{\bar{{\mathcal {E}}}}_{\lambda }' } = \sqrt{ \varvec{\bar{{\mathcal {E}}}}_{\lambda } } = {\varvec{Z}} \) remains unchanged.

Lemmas and Proof of Theorem 2

1.1 Lemmas on Positive Semidefinite Matrices

Lemma 7

For any \({\varvec{M}} = \sum _{i=1}^{n} {\varvec{M}}_i\), if \({\varvec{M}}_i \succeq {\varvec{O}}\) for each \(i \in \left[ 1:n \right] \), then \({\varvec{M}} \succeq {\varvec{O}}\).

Lemma 8

For any \({\varvec{M}} \succeq {\varvec{O}}\), we have \({\varvec{x}}^{\intercal } {\varvec{M}} {\varvec{x}} = 0 \iff {\varvec{M}} {\varvec{x}} = {\varvec{0}}\).

Proof

Given \({\varvec{M}} \succeq {\varvec{O}}\), there exist a matrix \(\sqrt{{\varvec{M}}}\), such that \({\varvec{M}} = \sqrt{{\varvec{M}}} \sqrt{{\varvec{M}}}\). Therefore:

This completes the proof of the sufficiency. The necessity is obvious. \(\square \)

Lemma 9

For any \({\varvec{M}} = \sum _{i=1}^{n} {\varvec{M}}_i \succeq {\varvec{O}}\), with each \({\varvec{M}}_i \succeq {\varvec{O}}\), we have \({\varvec{M}} {\varvec{x}} = {\varvec{0}} \iff {\varvec{M}}_i {\varvec{x}} = {\varvec{0}} \) for each \(i \in \left[ 1:n \right] \).

Proof

Because \({\varvec{M}} \succeq {\varvec{O}}\) and \({\varvec{M}}_i \succeq {\varvec{O}}\), by Lemma 8, it suffices to show \({\varvec{x}}^{\intercal } {\varvec{M}} {\varvec{x}} = 0 \iff {\varvec{x}}^{\intercal } {\varvec{M}}_i\, {\varvec{x}} = 0\), which is true because \( {\varvec{x}}^{\intercal } {\varvec{M}} {\varvec{x}} = \sum _{i=1}^{n} {\varvec{x}}^{\intercal } {\varvec{M}}_i\, {\varvec{x}} \) with each summand \({\varvec{x}}^{\intercal } {\varvec{M}}_i\, {\varvec{x}} \ge 0\). Therefore \({\varvec{x}}^{\intercal } {\varvec{M}} {\varvec{x}} = 0 \iff {\varvec{x}}^{\intercal } {\varvec{M}}_i\, {\varvec{x}} = 0\) for each \(i \in \left[ 1:n \right] \). \(\square \)

1.2 Proof of Theorem 2

We assume \(\mu _i\) being a tuning parameter \(\mu _i \ge 0\). Then the following chain of generalized equality holds:

Therefore \({\varvec{I}} - {\varvec{Q}}_i \succeq {\varvec{O}}\), and thus \( \varvec{{\mathcal {P}}}_{\mathrm {II}} = \sum _{i=1}^{n} \left( {\varvec{I}} - {\varvec{Q}}_i\right) \succeq {\varvec{O}} \).

Proof

We shall prove the result by showing: \((a)\iff (b)\), \((a)\iff (c)\), \((a)\iff (d)\), and \((c)\iff (e)\).

\((a)\iff (b)\). Because \(\varvec{{\mathcal {P}}}_{\mathrm {II}} = n {\varvec{I}} - \varvec{{\mathcal {Q}}}_{\mathrm {II}} \), thus \(\varvec{{\mathcal {P}}}_{\mathrm {II}} {\varvec{1}} = n {\varvec{1}} - \varvec{{\mathcal {Q}}}_{\mathrm {II}} {\varvec{1}} \), which means \(\varvec{{\mathcal {P}}}_{\mathrm {II}} {\varvec{1}} = {\varvec{0}} \iff \varvec{{\mathcal {Q}}}_{\mathrm {II}} {\varvec{1}} = n {\varvec{1}}\).

\((a)\iff (c)\). Note that \(\varvec{{\mathcal {P}}}_{\mathrm {II}} = \sum _{i=1}^{n} \left( {\varvec{I}} - {\varvec{Q}}_i\right) \succeq {\varvec{O}}\) with \({\varvec{I}} - {\varvec{Q}}_i \succeq {\varvec{O}}\). Therefore by Lemma 9, we have \( \varvec{{\mathcal {P}}}_{\mathrm {II}} {\varvec{1}} = {\varvec{0}} \iff \left( {\varvec{I}} - {\varvec{Q}}_i\right) {\varvec{1}} = {\varvec{0}} \iff {\varvec{Q}}_i {\varvec{1}} = {\varvec{1}} \).

\((a)\iff (d)\). Applying a translation \({\varvec{t}}\) to \({\varvec{S}}\), we have:

Sufficiency: if \(\varvec{{\mathcal {P}}}_{\mathrm {II}} {\varvec{1}} = {\varvec{0}}\), then \(- 2 \mathbf {tr}\left( {\varvec{S}} \varvec{{\mathcal {P}}}_{\mathrm {II}} {\varvec{1}} {\varvec{t}}^{\intercal } \right) = \mathbf {tr}\left( {\varvec{t}} {\varvec{1}}^{\intercal } \varvec{{\mathcal {P}}}_{\mathrm {II}} {\varvec{1}} {\varvec{t}}^{\intercal }\right) = 0 \).

Necessity: note that \(\varvec{{\mathcal {P}}}_{\mathrm {II}}\) is positive semidefinite, so \(\mathbf {tr}\left( {\varvec{t}} {\varvec{1}}^{\intercal } \varvec{{\mathcal {P}}}_{\mathrm {II}} {\varvec{1}} {\varvec{t}}^{\intercal }\right) \ge 0\), which means \(\mathbf {tr}\left( {\varvec{S}} \varvec{{\mathcal {P}}}_{\mathrm {II}} {\varvec{1}} {\varvec{t}}^{\intercal }\right) \le 0\). Since Eq. (34) holds for any \({\varvec{t}}\), flip the sign of \({\varvec{t}}^{\intercal }\), we obtain \(\mathbf {tr}\left( {\varvec{S}} \varvec{{\mathcal {P}}}_{\mathrm {II}} {\varvec{1}} {\varvec{t}}^{\intercal }\right) \ge 0\). Therefore:

Note that \({\varvec{1}}^{\intercal } \varvec{{\mathcal {P}}}_{\mathrm {II}} {\varvec{1}}\) is a scalar, hence:

for any \({\varvec{t}}\), if and only if \( {\varvec{1}}^{\intercal } \varvec{{\mathcal {P}}}_{\mathrm {II}} {\varvec{1}} = 0 \iff \varvec{{\mathcal {P}}}_{\mathrm {II}} {\varvec{1}} ={\varvec{0}} \) by Lemma 8.

\((c)\iff (e)\). Let \( \varvec{\bar{Q}}_i = \varvec{{\mathcal {B}}}_i^{\intercal } \left( \varvec{{\mathcal {B}}}_i \varvec{{\mathcal {B}}}_i^{\intercal } \right) ^{-1} \varvec{{\mathcal {B}}}_i \), and \(\left( \varvec{{\mathcal {B}}}_i^{\intercal }\right) ^{\dagger } = \left( \varvec{{\mathcal {B}}}_i \varvec{{\mathcal {B}}}_i^{\intercal } \right) ^{-1} \varvec{{\mathcal {B}}}_i \). Applying the Woodbury matrix identity, the matrix \({\varvec{Q}}_i\) can be expanded as:

Then we write:

with \({\varvec{I}} - \varvec{\bar{Q}}_i \succeq {\varvec{O}}\) and \({\varvec{{\Delta }}_i} \succeq {\varvec{O}}\). Noting that \(\mu _i \ge 0\), by Lemma 9, we have \(\left( {\varvec{I}} - {\varvec{Q}}_i\right) {\varvec{1}} = {\varvec{0}}\) if and only if \(\left( {\varvec{I}} - \varvec{\bar{Q}}_i\right) {\varvec{1}} = {\varvec{0}}\) and \( {\varvec{{\Delta }}_i} {\varvec{1}} = {\varvec{0}} \). In other words, \({\varvec{Q}}_i {\varvec{1}} = {\varvec{1}}\) if and only if \(\varvec{\bar{Q}}_i {\varvec{1}} = {\varvec{1}}\) and \({\varvec{{\Delta }}_i} {\varvec{1}} = {\varvec{0}}\).

Note that \(\varvec{\bar{Q}}_i = \varvec{{\mathcal {B}}}_i^{\intercal } \left( \varvec{{\mathcal {B}}}_i \varvec{{\mathcal {B}}}_i^{\intercal } \right) ^{-1} \varvec{{\mathcal {B}}}_i\) is the orthogonal projection into the range space of \(\varvec{{\mathcal {B}}}_i^{\intercal }\), which states \(\varvec{\bar{Q}}_i {\varvec{1}} = {\varvec{1}} \iff {\varvec{1}} \in \mathbf {Range}\left( \varvec{{\mathcal {B}}}_i^{\intercal } \right) \iff \exists {\varvec{x}}\) such that \(\varvec{{\mathcal {B}}}_i^{\intercal } {\varvec{x}} = {\varvec{1}}\). Therefore such an \({\varvec{x}}\) is \({\varvec{x}} = \left( \varvec{{\mathcal {B}}}_i^{\intercal }\right) ^{\dagger } {\varvec{1}}\).

Since \({\varvec{Z}}_i \left( \varvec{{\mathcal {B}}}_i \varvec{{\mathcal {B}}}_i^{\intercal }\right) ^{-1} {\varvec{Z}}_i^{\intercal } \succeq {\varvec{O}}\), we know \({\varvec{I}} + \mu _i {\varvec{Z}}_i \left( \varvec{{\mathcal {B}}}_i \varvec{{\mathcal {B}}}_i^{\intercal }\right) ^{-1} {\varvec{Z}}_i^{\intercal }\) is exactly positive definite. Therefore \({\varvec{{\Delta }}_i} {\varvec{1}} = {\varvec{0}} \iff {\varvec{1}}^{\intercal } {\varvec{{\Delta }}_i} {\varvec{1}} = 0\) happens if and only if \({\varvec{Z}}_i \left( \varvec{{\mathcal {B}}}_i^{\intercal }\right) ^{\dagger } {\varvec{1}} = {\varvec{0}}\) which is \({\varvec{Z}}_i {\varvec{x}} = {\varvec{0}}\). \(\square \)

Proof of Theorem 3

Proof

\((b) \iff (c)\). We define \(\varvec{\bar{Q}}_i = \varvec{\Gamma }_i \varvec{{\mathcal {B}}}_i^{\intercal } \left( \varvec{{\mathcal {B}}}_i \varvec{\Gamma }_i \varvec{\Gamma }_i \varvec{{\mathcal {B}}}_i^{\intercal } \right) ^{-1} \varvec{{\mathcal {B}}}_i \varvec{\Gamma }_i\). By the Woodbury matrix identity:

Then \({\varvec{P}}_i\) can be written as: \( {\varvec{P}}_i = \left( \varvec{\Gamma }_i - \varvec{\bar{Q}}_i \right) + \mu _i {\varvec{{\Delta }}_i} \) with \(\varvec{\Gamma }_i - \varvec{\bar{Q}}_i \succeq {\varvec{O}}\) and \(\varvec{{\Delta }}_i \succeq {\varvec{O}}\). Thus \({\varvec{P}}_i \succeq {\varvec{O}}\). By Lemma 9, \({\varvec{P}}_i {\varvec{1}} = {\varvec{0}}\) if and only if \(\left( \varvec{\Gamma }_i - \varvec{\bar{Q}}_i \right) {\varvec{1}} = {\varvec{0}}\) and \(\varvec{{\Delta }}_i {\varvec{1}} = {\varvec{0}}\)

Note that \(\varvec{\bar{Q}}_i\) is the orthogonal projection to \(\mathbf {Range}\left( \varvec{\Gamma }_i \varvec{{\mathcal {B}}}_i^{\intercal } \right) \). Thus we have:

if and only if there exists \({\varvec{x}}\) such that \(\varvec{\Gamma }_i \varvec{{\mathcal {B}}}_i^{\intercal } {\varvec{x}} = \varvec{\Gamma }_i {\varvec{1}}\). Such an \({\varvec{x}}\) is \({\varvec{x}} = \left( \varvec{\Gamma }_i \varvec{{\mathcal {B}}}_i^{\intercal }\right) ^{\dagger } {\varvec{1}}\).

Note that \(\varvec{{\Delta }}_i \succeq {\varvec{O}}\) and \({\varvec{I}} + \mu _i {\varvec{Z}}_i \left( \varvec{{\mathcal {B}}}_i \varvec{\Gamma }_i \varvec{\Gamma }_i \varvec{{\mathcal {B}}}_i^{\intercal }\right) ^{-1} {\varvec{Z}}_i^{\intercal } \succ {\varvec{O}}\), thus we have:

\((a) \iff (b)\). Since \({\varvec{P}}_i \succeq {\varvec{O}}\), the proof is immediate by Lemma 9.

This completes the proof. \(\square \)

Rigid Procrustes Analysis with Two Shapes

The pairwise Procrustes problem aims to solve for a similarity transformation (a scale factor \(s>0\), a rotation/reflection \({\varvec{R}} \in \mathrm {O}(3)\), and a translation \({\varvec{t}} \in {\mathbb {R}}^3\)), that minimizes the registration error between the two point-clouds \({\varvec{D}}_1\) and \({\varvec{D}}_2\) with known correspondences:

This problem admits a closed-form solution (Horn 1987; Horn et al. 1988; Kanatani 1994).

Let \( \varvec{\bar{D}}_1 = {\varvec{D}}_1 - \frac{1}{m} {\varvec{D}}_1 {\varvec{1}} {\varvec{1}}^{\intercal } \), and \( \varvec{\bar{D}}_2 = {\varvec{D}}_2 - \frac{1}{m} {\varvec{D}}_2 {\varvec{1}} {\varvec{1}}^{\intercal } \). Let the SVD of \(\varvec{\bar{D}}_1 \varvec{\bar{D}}_2^{\intercal }\) be \( \varvec{\bar{D}}_1 \varvec{\bar{D}}_2^{\intercal } = \varvec{\bar{U}} \varvec{\bar{\Sigma }} \varvec{\bar{V}}^{\intercal } \). Then the optimal rotation is \({\varvec{R}}^{\star } = \varvec{\bar{V}} \varvec{\bar{U}}^{\intercal }\) if \({\varvec{R}}^{\star } \in \mathrm {O}(3)\). If we require \({\varvec{R}}^{\star } \in \mathrm {SO}(3)\), the optimal rotation is:

The optimal scale factor and the optimal translation are then given by:

Reprints and permissions