Can I calculate the transformation matrix with many points taken into account?

Question

I am trying to subtract 2D conversion options from the result.

A large number of samples are indicated in the unknown coordinate system XY, as well as their corresponding analogues in WGS84 (longitude, latitude). Since the area is small, we can assume that the target system is also flat.

Unfortunately, I don’t know what scaling order, rotation, translation was used, and I’m not even sure if there were 1 or 2 translations.

I tried to create a long system of equations, but it turned out to be too complicated for me. The basic geometry also did not help me, since the order of the transforms is unknown, and I would have to check every possible order of combinations.

Is there a systematic approach to this problem?

Answer 1

It’s easy to figure out the scale factor, just select any two points and find the distance between them in the XY space and the WGS84 space, and the ratio between them is your scaling factor.

Rotations and translations are a bit more complicated, but not so difficult when you learn that the result of applying any number of rotations or translations (in only 2 dimensions!) Can be reduced to one rotation around some unknown point at some unknown angle.

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

The rotation calculation is as follows:

Pr = R*(Pxy - Paxis_xy) + Paxis_xy 

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

To actually find the 3 unknowns, you need to expand your WGS84 points (or, respectively, scale your XY points) using the scaling factor you found, and shift the points so that the two coordinate systems have the same origin.

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

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

where I abbreviated A = Paxis_xy for short. Subtraction of 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 own angle, you can also quickly check that cos(a) * cos(a) + sin(a) * sin(a) == 1 , to make sure that you correctly corrected all the calculations or that your system really is an orientation-preserving isometry (consists of translations and rotations only).

Now that we know a , we know R and therefore find a :

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

where the inversion of the 2x2 matrix is ​​simple .

EDIT . In the case above, there is an error, or more specifically, one case that needs to be handled separately.
There is one combination of translations and rotations that cannot be reduced to one rotation, and this is one translation. You can think about it in terms of fixed points (how many points have not changed after the operation).
The translation has no fixed points (all points are changed), and rotation has 1 fixed point (the axis does not change). It turns out that two turns leave 1 fixed point, and translation and rotation leave 1 fixed point, which (with a little proof, which indicates that the number of fixed points informs you of the operation performed), is the reason that arbitrary combinations of them lead to one rotation.

This means that if your angle is expressed as 0 , then using the above method will also give you A = 0 , which is probably incorrect. In this case, you must do A = P1' - P1 .

Answer 2

If I understand the question correctly, you have n points (X ₁ , Y ₁ ), ..., (X _n , Y _n ), corresponding points, say, (x ₁ , y ₁ ), ..., (x _n , y _n ) in a different coordinate system, and the former are presumably obtained from the latter by rotation, scaling, and translation.

Note that this data does not define a fixed rotation / scaling point or the order in which "should" operations should be performed. On the other hand, if you know this in advance or choose them arbitrarily, you will find a rotation, translation, and scaling factor that transforms the data as intended.

For example, you can select any point, for example, p ₀ = [X ₁ , Y ₁ ] ^T (column vector) as a fixed point of rotation and scaling and subtract its coordinates from the coordinates of two other points to get p ₂ = [X _2- X ₁ , Y ₂ -Y ₁ ] ^T and p ₃ = [X ₃ -X ₁ , Y ₃ -Y ₁ ] ^T. Also take the column vectors q ₂ = [x ₂ -x ₁ , y ₂ -y ₁ ] ^T q ₃ = [x ₃ -x ₁ , y <sub>3sub> -y _{1sub>] ^T. Now [p 2 p 3 ] = A * [q 2 q 3 ], where A is unknwon 2x2 representing the scaling of the company. You can solve it (unless you are lucky and choose degenerate points), because A = [p 2 p 3 ] * [q 2 q 3 ] ^-1 where ^-1 denotes the inverse matrix (from the matrix 2x2 [q 2 q 3sub> ]). Now, if the transformation between coordinate systems is really a roto-scaling translator, all points must satisfy P k = A * (Q k -q 0 ) + p 0 , where P k = [X k , Y k ] ^T Q k = [x k , y k ] ^T q 0 = [x 1 , y 1 ] ^T and k = 1, ..., n.}

If you want, you can easily determine the scaling and rotation parameter from components A or combine b = -A * q ₀ + p ₀ to get P _k = A * Q _k + b.

The above method responds poorly to noise or selects degenerate points. If necessary, this can be eliminated by using, for example, basic component analysis, which is also only a few lines of code if MATLAB or some other linear algebra tools are available.