Can I calculate a transformation matrix given a set of points?

Question

I'm trying to deduct the 2D-transformation parameters from the result.

Given is a large number of samples in an unknown X-Y-coordinate system as well as their respective counterparts in WGS84 (longitude, latitude). Since the area is small, we can assume the target system to be flat, too.

Sadly I don't know which order of scale, rotate, translate was used, and I'm not even sure if there were 1 or 2 translations.

I tried to create a lengthy equation system, but that ended up too complex for me to handle. Basic geometry also failed me, as the order of transformations is unknown and I would have to check every possible combination order.

Is there a systematic approach to this problem?

Answer 1

Figuring out the scaling factor is easy, just choose any two points and find the distance between them in your X-Y space and your WGS84 space and the ratio of them is your scaling factor.

The rotations and translations is a little trickier, but not nearly as difficult when you learn that the result of applying any number of rotations or translations (in 2 dimensions only!) can be reduced to a single rotation about some unknown point by some unknown angle.

Suddenly you have N points to determine 3 unknowns, the axis of rotation (x and y coordinate) and the angle of rotation.

Calculating the rotation looks like this:

    Pr = R*(Pxy - Paxis_xy) + Paxis_xy

Pr is your rotated point in X-Y space which then needs to be converted to WGS84 space (if the axes of your coordinate systems are different).
R is the familiar rotation matrix depending on your rotation angle.
Pxy is your unrotated point in X-Y space.
Paxis_xy is the axis of rotation in X-Y space.

To actually find the 3 unknowns, you need to un-scale your WGS84 points (or equivalently scale your X-Y points) by the scaling factor you found and shift your points so that the two coordinate systems have the same origin.

First, finding the angle of rotation: take two corresponding pairs of points P1, P1' and P2, P2' and write out

     P1' = R(P1-A) + A
     P2' = R(P2-A) + A

where I swapped A = Paxis_xy for brevity. Subtracting the two equations gives:

     P2'-P1' = R(P2-P1)
     B = R * C
     Bx = cos(a) * Cx - sin(a) * Cy
     By = cos(a) * Cx + sin(a) * Cy
     By + Bx = 2 * cos(a) * Cx
     (By + Bx) / (2 * Cx) = cos(a)
     ...
     (By - Bx) / (2 * Cy) = sin(a)
     a = atan2(sin(a), cos(a)) <-- to get the right quadrant

And you have your angle, you can also do a quick check that cos(a) * cos(a) + sin(a) * sin(a) == 1 to make sure either you got all the calculations correct or that your system really is an orientation-preserving isometry (consists only of translations and rotations).

Now that we know a we know R and so to find A we do:

    P1` = R(P1-A) + A
    P1' - R*P1 = (I-R)A
    A = (inverse(I-R)) * (P1' - R*P1)

where the inversion of a 2x2 matrix is easy.

EDIT: There is an error in the above, or more specifically one case that needs to be treated separately.
There is one combination of translations and rotations that does not reduce to a single rotation and that is a single translation. You can think of it in terms of fixed points (how many points are unchanged after the operation).
A translation has no fixed points (all points are changed) and a rotation has 1 fixed point (the axis doesn't change). It turns out that two rotations leave 1 fixed point and a translation and a rotation leaves 1 fixed point, which (with a little proof that says the number of fixed points tells you the operation performed) is the reason that arbitrary combinations of these result in a single rotation.

What this means for you is that if your angle comes out as 0 then using the method above will give you A = 0 as well, which is likely incorrect. In this case you have to do A = P1' - P1.