finding out 2D transformation given a list of sample points

Question

I have a bunch of 2D points that are rotated by $\theta$ and translated by $(\Delta_x, \Delta_y)$. I.e.

$$ x'=x \cos(\theta)-y \sin(\theta)+\Delta_x \\ y'=x \sin(\theta)+y \cos(\theta)+\Delta_y $$

My input is a list of $(x, y, x', y')$, which is noisy. I want to find $\theta$, $\Delta_x$ and $\Delta_y$ which give minimal error. This sounds like a regression problem, but I don't know how to handle the trigonometric functions.

If the trigonometric functions are too hard to handle, I can simplify the problem by assuming $\theta$ to be small, so that we can use the following approximation: $$ x'=x - y \theta+\Delta_x \\ y'=x \theta+y+\Delta_y $$

Even with this simplification, I still cannot apply standard least-square linear regression. There are two equations and $\theta$ is shared between them.

Answer 1

This is almost least squares regression except for the constraint on the parameter. In fact, by viewing $(x,y)\sim x + iy = z$ and $(x^\prime, y^\prime) = x^\prime + i y^\prime = z^\prime$ as complex numbers, we may write the model as

$$z_j^\prime = \alpha + \beta z_j + \varepsilon_j$$

where $\alpha = \Delta_x + i\Delta_y$ and $\beta = \exp(i\theta)$.

Assuming the (complex) random errors $\varepsilon_j$ have zero mean, a common variance, and are independent, it is easy to derive the conclusion that the fit must pass through the "point of means" $(\bar z, \bar z^\prime)$. This determines $\hat \alpha$, reducing the problem to the form

$$z_j^\prime - \bar z^\prime = \beta(z_j - \bar z) + \varepsilon_j$$

subject to the constraint $|\beta|^2 = 1$. This is readily solved using Lagrange multipliers. However, even a univariate minimizer will have no problems with this: ask it to minimize the sum of squares

$$\sum_j |(z_j^\prime - \bar z^\prime) - \exp(i\theta)(z_j - \bar z)|^2$$

Initialize the solution by fitting the unconstrained model (using ordinary regression) to find $\hat\beta$. Rescale $\hat\beta$ to unit length, then apply the univariate minimizer.

All this can be done entirely without using complex numbers, simply by translating the foregoing back into 2D coordinates. This is the method used in the R code below. Only two lines of code are needed:

lambda <- function(theta) sum((xy.prime.shift - rotate(theta) %*% xy.shift)^2)
fit <- optimize(lambda, interval=c(0, 2*pi))

Everything else is preliminary processing or post-processing to apply the fit.

From left to right, the figure shows the original points $z$, the target points $z^\prime$, the graph of the sum of squares against $\theta$ (the optimal value of $\theta$ is marked in red), and the comparison of the fitted values (gray circles) to the target values (red disks). The black arrows show the residual displacements.

By minimizing the sum of squares, this example assumed the components of the errors were uncorrelated and of equal variance. It is readily generalized to the case where the errors may be correlated and/or have unequal variances, without making any fundamental change to the overall approach.

Those familiar with the Gauss-Markov theorem and Maximum Likelihood estimation will have no difficulties obtaining standard errors for the estimates and testing hypotheses. This is not necessary in many applications: often the purpose of this exercise is to find the optimal transformation between two sets of points. However, such auxiliary information would be very useful in estimating the accuracy of the fit when it is applied to new starting points (that is, to compute prediction intervals). Despite the potential value of such an approach, I have never seen anyone do that.

---

#
# Preliminaries.
#
set.seed(17)
rotate <- function(theta) matrix(cos(theta)*c(1,0,0,1)+sin(theta)*c(0,-1,1,0), 2, 2)
#
# Generate data.
#
n <- 10
xy <- matrix(rnorm(n*2), nrow=2)                # Original points
delta <- c(1,-2)                                # True mean displacement
theta <- 1                                      # Rotation (radians)
beta <- rotate(theta)                           # Rotation matrix
epsilon <- matrix(rnorm(n*2, sd=0.5), nrow=2)   # Errors
xy.prime <- beta %*% xy + delta + epsilon       # New (target) points
#
# Find an initial rotation by conducting an unconstrained fit.
#
xy.shift <- xy - rowMeans(xy)                   # Recentered original points
xy.prime.shift <- xy.prime - rowMeans(xy.prime) # Recentered target points
fit.0 <- lm(t(xy.prime.shift) ~ t(xy.shift))    # Unconstrained fit
b.0 <- t(coef(fit.0)[-1, ])                     # Estimate of beta
theta.0 <- atan2(b.0[1,2], b.0[1,1])            # The corresponding angle
delta.hat <- rowMeans(xy.prime - xy)            # The estimated displacment
#
# Fit the data.
#
lambda <- function(theta) sum((xy.prime.shift - rotate(theta) %*% xy.shift)^2)
fit <- optimize(lambda, interval=c(0, 2*pi))
#
# Extract the estimated rotation angle and compute the fitted values.
#
theta.hat <- fit$minimum
xy.pred <- rotate(theta.hat) %*% xy.shift + rowMeans(xy.prime)
#
# Plot things.
#
#pairs(t(rbind(xy, xy.prime))) # Not too revealing!
par(mfrow=c(1,4))
plot(t(xy), pch=".", asp=1, type="b", col="Gray", xlab="X", ylab="Y", 
     main="Original points")
text(xy[1, ], xy[2, ], 1:n)

plot(t(xy.prime), pch=".", asp=1, type="b", col="Gray", xlab="X", ylab="Y", 
     main="Target points")
text(xy.prime[1, ], xy.prime[2, ], 1:n)

x <- seq(0, 2*pi, length.out=101)
plot(x, sapply(x, lambda), type="l", xlab="theta", ylab="Objective",
     main="Objective function")

points(theta.hat, fit$objective, pch=16, col="Red")
plot(t(cbind(xy.prime, xy.pred)), type="n", xlab="X", ylab="Y", main="Fit")
lines(t(xy.prime), col="Gray")
arrows(xy.pred[1, ], xy.pred[2, ], xy.prime[1, ], xy.prime[2, ], length=0.1)
points(t(xy.pred), col="Gray")
points(t(xy.prime), pch=16, col="Red")