Calculating new positions for points on a 2D plane that has moved

Question

I have a 2D plane, with few points on it. Coordinates for all points are known. The black plane is basically a piece of paper and the points are holes in it. I then take the black plane and rotate it by 180 degrees on the A axis, then also move it around a bit (how much is not known) and then rotate it on the vertical axis (how much is not known) to end up in the final position that can be seen at the bottom of the drawing (the plane is still 2D and flat). I now manually measure the new coordinates for the red and blue points, so now I know their coordinates and also how much they moved on the A,B axes relative to their previous position.

With all this information, how do I calculate the new coordinates of the remaining two points? (In case this can't be done with just knowing new coordinates of two points, I can add a third one)

Excuse the crude drawing.

Could this be solved by using something like a 2d matrix?

I take the positions of the red, blue and one of the black points

|x1,x2,x3|
|y1,y2,y3|

|6,10,7|
|7,4,5|

And have a matrix representing a triangle.

Then after I moved and rotated the the triangle I have these known values:

|x1,x2,x3|
|y1,y2,y3|

|7,12,?|
|4,6,?|

How to calculate those two missing numbers in this new matrix, so the shape of the triangle remains unchanged ?

Answer 1

As mentioned by the OP, the transformation done to the points consists on a 2D reflection and a translation. Thus, we seek to find a reflection matrix $R$, and a translation vector $t$.

If we denote, for an observation $i$, the coordinates of a point before the transfomation as $p_i$ and the points after the transformation as $q_i$, then $R$ and $t$ should satisfy the following: $$ q_i = R\,p_i + t $$ However, since the coordinates of the points are measurements, they are subjected to noise. Hence, we can't get a perfect transformation and we will have some error. Thereby, the previous expression will present some error $\varepsilon_i$: $$ \varepsilon_i = q_i - (R\,p_i + t ) $$

If we choose to minimize the sum of squared errors, i.e.: $$ R, t = \arg\min_{R,\,t} \sum_{i=1}^n \Vert\varepsilon_i\Vert^2 $$ corresponding to the $n$ observations (2 in the OP's problem), then a closed-form solution is derived in "Least-squares estimation of transformation parameters between two point patterns, S. Umeyama", at equations 34-43.

However, since we are interested in a reflection, instead of a rotation, then equations 39 and 43 should be reversed, to ensure a reflection instead of a rotation.

Following this procedure we would reach the following transformation: $$ R = \begin{bmatrix} 0.96561576 & -0.25997347\\ -0.25997347 & -0.96561576 \end{bmatrix}, \qquad % t = \begin{bmatrix} 3.20492804 \\ 12.39067446 \end{bmatrix} $$ Graphically:

In case it is useful, here is the python code that I used to get the previous result (I'm sorry that I'm not used to c#). I've tried to comment each step. For clarity, the matrix/vector dimensions are indicated in a comment to the right of them.

import matplotlib.pyplot as plt
import numpy as np

# measured prior to the transformation.
X = np.array([
    [6.0, 10.0], 
    [7.0, 4.0]
    ]) # (2, 2)
# measured after the transformation.
Y = np.array([[7.0, 12.0], 
              [4.0, 6.0]]) # (2, 2)

# get centers of mass in the prior and post configuration.
centroid_x = X.mean(axis=1, keepdims=True) # (2, 1)
centroid_y = Y.mean(axis=1, keepdims=True) # (2, 1)

# sample cross-covariance.
n = 2 # we have only two samples (points)
cov_xy = (Y - centroid_y) @ (X - centroid_x).T / n # (2, 2)

# svd.
U, d, Vt = np.linalg.svd(cov_xy) # (2, 2), (2,), (2, 2) 

# some useful values for later.
UV_det = np.linalg.det(U @ Vt) # float
cov_det = np.linalg.det(cov_xy) # float
cov_rank = np.linalg.matrix_rank(cov_xy) # int

# ensure non-degenerate case.
if cov_rank < 1:
    raise ValueError

# ensure reflection instead of rotation.
S = np.eye(2) # (2, 2)
if (cov_rank==1 and UV_det>0) or (cov_rank>1 and cov_det>0):
    S[1, 1] = -1 # set bottom-right element of S to -1.
    
# estimate transformation.
R = U @ S @ Vt # (2, 2)
t = centroid_y - R.dot(centroid_x) # (2, 1)

# visualize and test with the other points.
X_other = np.array([
    [7.0, 9.0], 
    [5.0, 6.0]
    ]) # (2, 2)
fig, axes = plt.subplots(ncols=2)

# plot points prior to the transformation and their coords.
axes[0].scatter(X[0], X[1], c=['r', 'b'], s=50)
axes[0].scatter(X_other[0], X_other[1], c='k', s=50)

# plot points after the transformation and their coords.
Y_ = R.dot(X) + t
Y_other = R.dot(X_other) + t

axes[1].scatter(Y_[0], Y_[1], c=['r', 'b'], s=50)
axes[1].scatter(Y[0], Y[1], c=['r', 'b'], s=50, alpha=0.5)
axes[1].scatter(Y_other[0], Y_other[1], c='k', s=50)

# annotations.
for x in np.concatenate((X.T, X_other.T)):
    axes[0].annotate(f'{tuple(x.astype(int))}', x)
    
for y in np.concatenate((Y_.T, Y_other.T)):
    axes[1].annotate(f'{tuple(np.round(y, 1))}', y)

# style.
axes[0].set(
    xlim=(1, 13),
    ylim=(1, 9),
    xticks=np.arange(1, 14, 2),
    yticks=np.arange(1, 10, 2),
    aspect='equal',
    title='before'
    )

axes[1].set(
    xlim=(1, 13),
    ylim=(1, 9),
    xticks=np.arange(1, 14, 2),
    yticks=np.arange(1, 10, 2),
    aspect='equal',
    title='after'
    )
# axes[0].grid(True)
# axes[1].grid(True)

---

Edit

When doing a reflection or not is known before hand, this alternative code can be used

import matplotlib.pyplot as plt
import numpy as np

# specify if a reflection is done.
reflection_done = False

# measured prior to the transformation.
X = np.array([
    [2.0, 4.0], 
    [2.0, 1.0]
    ]) # (2, 2)
# measured after the transformation.
Y = np.array([[3.0, 5.0], 
              [2.0, 1.0]]) # (2, 2)

# get centers of mass in the prior and post configuration.
centroid_x = X.mean(axis=1, keepdims=True) # (2, 1)
centroid_y = Y.mean(axis=1, keepdims=True) # (2, 1)

# sample cross-covariance.
n = 2 # we have only two samples (points)
cov_xy = (Y - centroid_y) @ (X - centroid_x).T / n # (2, 2)

# svd.
U, d, Vt = np.linalg.svd(cov_xy) # (2, 2), (2,), (2, 2) 

# some useful values for later.
UV_det = np.linalg.det(U @ Vt) # float
cov_det = np.linalg.det(cov_xy) # float
cov_rank = np.linalg.matrix_rank(cov_xy) # int

# ensure non-degenerate case.
if cov_rank < 1:
    raise ValueError

# ensure the correct transformation.
S = np.eye(2) # (2, 2)
if reflection_done: 
    if (cov_rank==1 and UV_det>0) or (cov_rank>1 and cov_det>0):
        S[1, 1] = -1 # set bottom-right element of S to -1.
else:
    if (cov_rank==1 and UV_det<0) or (cov_rank>1 and cov_det<0):
        S[1, 1] = -1
    
# estimate transformation.
R = U @ S @ Vt # (2, 2)
t = centroid_y - R.dot(centroid_x) # (2, 1)

# visualize and test with the other points.
X_other = np.array([
    [3.0, 7.0, 9.0], 
    [4.0, 5.0, 6.0]
    ]) # (2, 2)
fig, axes = plt.subplots(ncols=2)

# plot points prior to the transformation and their coords.
axes[0].scatter(X[0], X[1], c=['r', 'b'], s=50)
axes[0].scatter(X_other[0], X_other[1], c='k', s=50)

# plot points after the transformation and their coords.
Y_ = R.dot(X) + t
Y_other = R.dot(X_other) + t

axes[1].scatter(Y_[0], Y_[1], c=['r', 'b'], s=50)
axes[1].scatter(Y[0], Y[1], c=['r', 'b'], s=50, alpha=0.5)
axes[1].scatter(Y_other[0], Y_other[1], c='k', s=50)

# annotations.
for x in np.concatenate((X.T, X_other.T)):
    axes[0].annotate(f'{tuple(x.astype(int))}', x)
    
for y in np.concatenate((Y_.T, Y_other.T)):
    axes[1].annotate(f'{tuple(np.round(y, 1))}', y)

# style.
get_xlim = lambda axes: (
    min(axes[0].dataLim.xmin, axes[1].dataLim.xmin) - 1, 
    max(axes[0].dataLim.xmax, axes[1].dataLim.xmax) + 1)
get_ylim = lambda axes: (
    min(axes[0].dataLim.ymin, axes[1].dataLim.ymin) - 1, 
    max(axes[0].dataLim.ymax, axes[1].dataLim.ymax) + 1)

axes[0].set(
    xlim=get_xlim(axes),
    ylim=get_ylim(axes),
    # xticks=np.arange(1, 14, 2),
    # yticks=np.arange(1, 10, 2),
    aspect='equal',
    title='before'
    )

axes[1].set(
    xlim=get_xlim(axes),
    ylim=get_ylim(axes),
    # xticks=np.arange(1, 14, 2),
    # yticks=np.arange(1, 10, 2),
    aspect='equal',
    title='after'
    )
# axes[0].grid(True)
# axes[1].grid(True)