How to find a transformation (non-uniform scaling and similarity) that maps one set of points to another?

Question

I have 2 set of points: $p_i$ and $q_i$ and I need to find and affine transformation that maps each point $p_i$ onto $q_i$. However there is a constraint that angles between unit vectors should be preserved, i.e. the transformation can be factored into scaling and rotation. Number of points is greater than $N$, so the problem is overdetermined, and I use least-squares formulation:

$$ \min_{R,S} \sum \|R S p_i - q_i\|^2 , \quad R \in \operatorname{SO}(N), S \in \operatorname{diag}(\mathbb{R}^N)$$

The problem is similar to How to solve an overdetermined system of point mappings via rotation and translation, except it doesn't take scaling into account. Following the answer to that question I came up with the following: $$ \begin{align} \|R S p_i-q_i\|^2&=(R S p_i-q_i)^\top(R S p_i-q_i) \\ &=p_i^\top S R^\top R S p_i+q_i^\top q_i-2q_i^\top R S p_i \\ &= p_i^\top S^2 p_i+q_i^\top q_i-2q_i^\top R S p_i\ \end{align} $$ $$ \def\argmin{\operatorname{argmin}} \argmin_{R,S}\sum \|R S p_i - q_i\|^2=\argmin_{R,S}\sum_i p_i^\top S^2 p_i-2q_i^\top R S p_i $$ $$ \def\tr{\operatorname{tr}} \begin{align} \sum_i p_i^\top S^2 p_i-2q_i^\top R S p_i &=\sum_i \tr p_i^\top S^2 p_i-2\tr q_i^\top R S p_i \\ &=\sum_i \tr p_i p_i^\top S^2 -2 \tr p_iq_i^\top R S \\ &=\tr\left(\left(\sum_ip_ip_i^\top\right)S^2 + \left(-2\sum_ip_iq_i^\top\right) R S\right) \\ &=:\tr C S^2+B R S \end{align} $$ This helped to simplify the problem to $N\times N$ matrices, so the number of points doesn't matter now. However I'm stuck at this point and can't figure out, how this expression can be minimized w.r.t $R$ and $S$.

I was able to solve this problem for 2-dimensional points by writing out matrices R and S in terms of their elements and then solving minimization problem by the method of Lagrange multipliers.

$$ R=\begin{pmatrix} r_1 & -r_2 \\ r_2 & r_1 \\ \end{pmatrix}, S=\begin{pmatrix} s_x & 0 \\ 0 & s_y \end{pmatrix} $$ Then the problem becomes to minimize $$ \tr{ \begin{pmatrix} c_{11} & c_{12} \\ c_{21} & c_{22} \\ \end{pmatrix} \begin{pmatrix} s_x^2 & 0 \\ 0 & s_y^2 \\ \end{pmatrix}+ \begin{pmatrix} b_{11} & b_{12} \\ b_{21} & b_{22} \\ \end{pmatrix} \begin{pmatrix} r_1 & -r_2 \\ r_2 & r_1 \\ \end{pmatrix} \begin{pmatrix} s_x & 0 \\ 0 & s_y \\ \end{pmatrix} } $$ with the constraint $r_1^2+r_2^2=1$. Solving it resulted in $$ k=\frac{b_{11}^2c_{22}+b_{22}^2c_{11}-b_{12}^2c_{22}-b_{21}^2c_{11}}{b_{12}b_{11}c_{22}-b_{21}b_{22}c_{11}}\\ r_1^2=\frac{2}{k^2+4\mp k\sqrt{k^2+4}}\\ r_2^2=\frac{2}{k^2+4\pm k\sqrt{k^2+4}}\\ s_x=-\frac{b_{11}}{2c_{11}}r_1-\frac{b_{12}}{2c_{11}}r_2\\ s_y=\frac{b_{21}}{2c_{22}}r_2-\frac{b_{22}}{2c_{22}}r_1 $$

However I'm not sure that I made no mistakes and I don't see any way to generalize this method to more than 2 dimensions.

I have also tried to find the transformation using pseudoinverses, but in that case I was able to get only affine transformations, and failed to add constraints, because they happen to be non-linear.

Answer 1

Thank you very much for your question and partial answer. I have encountered the same problem and only require the solution for 2D points, so your answer was very helpful to me. There is plenty of information available for finding the rigid transformation between corresponding points (SVD solutions in the links you posted / the Kabsch algorithm), which can optionally include a uniform scale component. Your post was the only information that I was able to find regarding solving for nonuniform scale, and crucially without allowing for shear.

I followed and agree with all your derivation steps up to and including the definition of the minimization problem, i.e. this point from your answer:

Then the problem becomes to minimize $$ \mathrm{tr}{ \begin{pmatrix} c_{11} & c_{12} \\ c_{21} & c_{22} \\ \end{pmatrix} \begin{pmatrix} s_x^2 & 0 \\ 0 & s_y^2 \\ \end{pmatrix}+ \begin{pmatrix} b_{11} & b_{12} \\ b_{21} & b_{22} \\ \end{pmatrix} \begin{pmatrix} r_1 & -r_2 \\ r_2 & r_1 \\ \end{pmatrix} \begin{pmatrix} s_x & 0 \\ 0 & s_y \\ \end{pmatrix} } $$ with the constraint $r_1^2+r_2^2=1$.

In your comment regarding how you solved this, I believe that I followed the same steps, using the Lagrange multiplier method. I will enumerate these steps in a bit more detail than you did, but I think our derivation proceeded the same way, basically until the final solution, and I obtained a different answer. First, the partial derivatives plus constraint equations that I used (following the Lagrange method):

$$ \frac{\partial}{\partial s_x} = b_{11} r_{1} + b_{12} r_{2} + 2 c_{11} s_{x} = 0 \\ \frac{\partial}{\partial s_y} = - b_{21} r_{2} + b_{22} r_{1} + 2 c_{22} s_{y} = 0 \\ \frac{\partial}{\partial r_1} = b_{11} s_{x} + b_{22} s_{y} - 2 \lambda r_{1} = 0 \\ \frac{\partial}{\partial r_2} = b_{12} s_{x} - b_{21} s_{y} - 2 \lambda r_{2} = 0 \\ r_{1}^{2} + r_{2}^{2} = 1 $$

I combined the $\frac{\partial}{\partial r_1}$ and $\frac{\partial}{\partial r_2}$ equations in order to eliminate $\lambda$, and then in the resulting equation substituted for $s_x$ and $s_y$. One can see in this equation that the $r_1^{2}$ and $r_2^{2}$ terms share common coefficients, which can be divided out. The resulting coefficient from the ${r_1}{r_2}$ cross term is how I arrived at what you defined as $k$, but with the sign flipped:

$$ r_{1}^{2} + \frac{- b_{11}^{2} c_{22} + b_{12}^{2} c_{22} + b_{21}^{2} c_{11} - b_{22}^{2} c_{11}}{b_{11} b_{12} c_{22} - b_{21} b_{22} c_{11}}r_{1} r_{2} - r_{2}^{2} = 0 $$

So, if one substitutes this coefficient with $k$ and then use the constraint equation to substitute $r_2$, the result is a quadratic equation, but with an square root in the cross term, meaning this is actually a quartic equation:

$$ 2 r_{1}^{2} + k r_{1} \sqrt{1 - r_{1}^{2}} - 1 = 0 $$

where

$$ k=\frac{- b_{11}^{2} c_{22} + b_{12}^{2} c_{22} + b_{21}^{2} c_{11} - b_{22}^{2} c_{11}}{b_{11} b_{12} c_{22} - b_{21} b_{22} c_{11}} $$

Instead of solving this myself, I just plugged it into a solver (sympy), and obtained this result:

$$ r_{1} = \{ - \frac{\sqrt{- \frac{2 k}{\sqrt{k^{2} + 4}} + 2}}{2}, \frac{\sqrt{- \frac{2 k}{\sqrt{k^{2} + 4}} + 2}}{2}, - \frac{\sqrt{\frac{2 k}{\sqrt{k^{2} + 4}} + 2}}{2}, \frac{\sqrt{\frac{2 k}{\sqrt{k^{2} + 4}} + 2}}{2} \} \\ r_{2} = \{ \frac{\sqrt{\frac{2 k}{\sqrt{k^{2} + 4}} + 2}}{2}, \frac{\sqrt{\frac{2 k}{\sqrt{k^{2} + 4}} + 2}}{2}, \frac{\sqrt{- \frac{2 k}{\sqrt{k^{2} + 4}} + 2}}{2}, \frac{\sqrt{- \frac{2 k}{\sqrt{k^{2} + 4}} + 2}}{2} \} $$

As far as I can tell, these do not simplify to your solution, so basically just this final part is where my solution diverges from yours (essentially the solution to the quartic equation).

EDIT: This was not correct, the solutions can be simplified to each other, so the answer is indeed the same as user502144's answer. I will leave my answer here, just in case any later readers want more details on solving the minimization problem.

The other scale unknowns are then solved as before (same as in your solution):

$$ s_x = - \frac{b_{11} r_{1} + b_{12} r_{2}}{2 c_{11}} \\ s_y = \frac{b_{21} r_{2} - b_{22} r_{1}}{2 c_{22}} $$

I coded the result, which requires a method of handling the multiple solutions to the quartic. For this I eliminated any solutions where $s_x$ and $s_y$ have different signs, and then I take as my final solution the one that produces the greatest variance accounted for in the fit. If random source points are scaled and rotated (in that order) and then translated to create target points, the algorithm is able to fit the transformation perfectly. The algorithm also works if the points are not scaled at all (i.e., also produces a correct rigid-only fit).