# Fitting of 2d data points with a function considering scaling, rotation and translation

I have the following set of 2d data points:

data1=
{
{21.557, 801.607}, {5.84689, 800.425}, {50.9284, 770.49},
{46.4516, 750.192}, {32.9808, 671.931}, {48.8067, 673.198},
{3.59394, 671.167}, {18.1513, 671.949}, {64.1628, 670.801},
{13.1805, 652.588}, {55.6619, 651.298}, {26.9262, 650.35},
{41.4876, 650.752}, {5.45129, 635.602}, {20.3858, 633.391},
{64.1931, 632.506}, {33.9168, 631.006}, {58.7559, 613.401},
{36.0045, 612.007}, {23.5348, 608.289}, {54.6781, 598.251},
{26.4914, 548.723}, {65.0549, 531.442}, {82.9996, 514.631},
{74.4132, 479.425}, {58.3295, 458.015}, {27.1816, 413.334}
}


I want to apply ScalingTransform, TranslationTransform and RotationTransform to find the best fit to transform data1 into data2, whereby:

data2=
{
{1530.03, 790.2}, {1514.13, 789.}, {1559.17, 758.9},
{1554.5, 738.5}, {1540.5, 660.237}, {1556.15, 661.154},
{1511.34, 659.395}, {1525.63, 660.167}, {1572.13, 658.656},
{1520.66, 640.844}, {1562.55, 639.132}, {1533.79, 638.607},
{1548.37, 638.933}, {1512.62, 623.985}, {1526.88, 621.69},
{1571.44, 620.556}, {1540.44, 618.794}, {1565.69, 601.532},
{1543.06, 600.093}, {1530.22, 596.423}, {1560.9, 586.053},
{1532.93, 536.587}, {1571.9, 519.25}, {1590.15, 501.882},
{1580.39, 467.111}, {1564.73, 445.615}, {1532.8, 400.935}
}


The corresponding points of data1 that should be transformed into data2 are already sorted and at the same position of the lists.

Here are plots of the two data sets:

plot1 = ListPlot[data1, PlotRange -> {{1, 91}, {300, 900}},
PlotStyle -> Red, Frame -> True,
FrameLabel -> {{"y", ""}, {"x", "data1"}},
BaseStyle -> {FontWeight -> "Bold", FontSize -> 15,
FontFamily -> "Calibri"}, ImageSize -> Large];

plot2 = ListPlot[data2, PlotRange -> {{1510, 1600}, {300, 900}},
PlotStyle -> Blue, Frame -> True,
FrameLabel -> {{"y", ""}, {"x", "data2"}},
BaseStyle -> {FontWeight -> "Bold", FontSize -> 15,
FontFamily -> "Calibri"}, ImageSize -> Large];

GraphicsColumn[{plot1, plot2}, ImageSize -> Large,
Spacings -> {{0, 0}, {0, 50}}]


I use the following naming:

s = ScalingTransform[{sx, sy}, {psx, psy}];
t = TranslationTransform[{vecx, vecy}];
r = RotationTransform[theta, {prx, pry}];


The combined transformation for each point {x, y} of data1 is:

combinedTransformation = s.t.r;


and finally :

combinedTransformation[{x, y}] =

{sx (prx (-Cos[theta]) + prx + pry Sin[theta]) + psx (-sx) + psx +
sx x Cos[theta] - sx y Sin[theta] + sx vecx,
sy (-(prx Sin[theta]) + pry (-Cos[theta]) + pry) + psy (-sy) + psy +
sy x Sin[theta] + sy y Cos[theta] + sy vecy}


The fitting parameters are: sx, sy, vecx, vecy, theta.

The scaling is centered at the point {psx, psy} and the 2d rotation is around the point {prx, pry}.

I would set {psx, psy} = {1, 1} and {prx, pry} = {1, 1}.

How can I transform data1 best into data2 and how can I obtain the fitting paramaters and their errors?

I already tried the same as what is proposed below by Ulrich Neumann and Carl Lange.

• The problem with FindGeometricTransform is, it is not described how the error is obtained - I need this for a paper. See this question.

• Second FindGeometricTransform does not give me the rotation angle and scaling factor in x and y separately, which are not exactly the same. FindGeometricTransform shows only the transformation function (or matrix) which is not enough for me.

• Try FindGeometricTransform . It 's not necessary to require a scaling point and/or a rotationpoint , that is the task ogf the fitting procedure. Commented Jan 9, 2019 at 11:43
• How do you define "the error"? Is that the mean distance between each point in data2 to the nearest point in the transformed data1? Or the square root of the mean of the square of those distances? Or something else? If the former, then (using @CarlLange 's code) data1Transformed = transform@data1; data2Nearest = Flatten[Nearest[data2, #] & /@ data1Transformed, 1]; Mean[Norm[#] & /@ (data1Transformed - data2)]might do it.
– JimB
Commented Jan 9, 2019 at 14:28
• If by "their errors" you mean the errors in the individual parameters, you'd need to specify a probabilistic model that generates the transformation parameters. Much like in a linear regression you need not just $y=a+bx$ but $y=a+bx+error$.
– JimB
Commented Jan 9, 2019 at 14:32
• I answered the question about what the error is in this question of yours. Commented Jan 9, 2019 at 17:57
– mrz
Commented Jan 17, 2019 at 10:26

Try FindGeometricTransform

trafo = FindGeometricTransform[data2, data1 ];
F = TransformationMatrix[trafo[[2]]]


F[[{1, 2}, 3]] is the offset. Matrix

T= F[[{1, 2}, {1,2}]]


describes rotation and scaling .

S = MatrixPower[ Transpose[T].T , 1/2]  (* scaling matrix*)
(*{{0.970832, -0.00629071}, {-0.00629071, 1.00107}}*)
R = Inverse[Transpose[T]].S             (* rotation matrix *)
(*{{0.999918, 0.0128058}, {-0.0128058, 0.999918}}*)
T - R.S // Chop                         (*T==R.S*)


The scaling factors are given by the eigenvalues of S. The rotation angle can be obtained by

J = #.# &[Flatten[RotationMatrix[\[CurlyPhi]] - R]];
NMinimize[{J, 0 <= \[CurlyPhi] <= 2 Pi  }, \[CurlyPhi]]
(*{4.36514*10^-15, {\[CurlyPhi] -> 6.27038}}*)
\[CurlyPhi]/Degree /. %[[2]]  (* angle in degree*)
(*359.266*)

– mrz
Commented Jan 9, 2019 at 12:33
• Thank you very much. What is the angle in degree (should be in the order of 0.5 degree)? The scaling factors should be about 1,004 in y and about 1.001 in x (data2 are larger by this factors). Do you understand what the error of FindGeometricTransform exactly means?
– mrz
Commented Jan 9, 2019 at 13:11
• Angle is around -.38 Degree. Commented Jan 9, 2019 at 13:19
• That's because the general mapping used by GeometericTransform maps vectors x->y in the form y=(A.x+b)/(c.x+d) and describes the central projection completly. The affine case c=0,d=1 is only a rough approximation! Commented Jan 9, 2019 at 15:19
• @mrz Usually the mean of nondiagonal elements of S describe shear: (S[[1,2]]+S[2,1])/2, (S[[1,3]]+S[3,1])/2,(S[[2,3]]+S[3,2])/2... Commented Jan 15, 2019 at 11:54

Per Ulrich Neumann's comment, FindGeometricTransform will do the job very nicely.

We get the transform by doing

transform = FindGeometricTransform[data2, data1][[2]]


This gives us a TransformationFunction, in this case:

$$\text{TransformationFunction}\left[\left( \begin{array}{ccc} 0.970671 & 0.00652924 & 1502.57 \\ -0.0187224 & 1.00107 & -12.4938 \\ -0.0000212516 & -\text{2.8535460791719293\grave{ }*{}^{\wedge}-7} & 1. \\ \end{array} \right)\right]$$

Now we can apply that TransformationFunction to our data and plot the result:

ListPlot[{transform@data1, data2}]


– mrz
Commented Jan 9, 2019 at 12:33
– mrz
Commented Jan 9, 2019 at 13:44

First change the function combinedTransformation to

combinedTransformation[{x_, y_}] =
{sx (prx (-Cos[theta]) + prx + pry Sin[theta]) + psx (-sx) + psx +
sx x Cos[theta] - sx y Sin[theta] + sx vecx,
sy (-(prx Sin[theta]) + pry (-Cos[theta]) + pry) + psy (-sy) + psy +
sy x Sin[theta] + sy y Cos[theta] + sy vecy}


and then try

v   = Map[combinedTransformation, data1] - data2;
err = Sum[v[[k]].v[[k]], {k, 1, Length[v]}];
sol = NMinimize[err, {prx, pry, psx, psy, sx, sy, vecx, vecy, x, y, theta}]


and the result

The plot was produced as

plot1 = ListPlot[Map[combinedTransformation, data1] /. sol[[2]],
PlotStyle -> Red, Frame -> True,
FrameLabel -> {{"y", ""}, {"x", "data1 and data2"}},
BaseStyle -> {FontWeight -> "Bold", FontSize -> 15,
FontFamily -> "Calibri"}, ImageSize -> Large]

plot2 = ListPlot[data2, PlotStyle -> Blue, Frame -> True,
FrameLabel -> {{"y", ""}, {"x", "data1 and data2"}},
BaseStyle -> {FontWeight -> "Bold", FontSize -> 15,
FontFamily -> "Calibri"}, ImageSize -> Large]

Show[plot1, plot2]

• Thanks a lot for this solution. You get sx -> 1.00395, sy -> 1.00219, which is what I expect. Theta (theta -> -0.00683116) is probably in radian which is 0.39 degree, also what I expect. Only the translation surprises me a little bit: vecx -> 1513.99, vecy -> -7.23374, compared to the results of Ulrich Neumann and Carl Lange (vecx ca. 1503 and vecy ca. -13). I trust their vecy value because of this estimation vecy = Mean[data2[[All, 2]] - data1[[All, 2]]]=-11.95. On the other hand: vecx=Mean[data2[[All, 1]] - data1[[All,1]]]=1507.11 is between your result and their.
– mrz
Commented Jan 9, 2019 at 13:32
– mrz
Commented Jan 9, 2019 at 13:44
• @mrz Included the plotting script Commented Jan 9, 2019 at 14:43
• Consider using FindMinimum instead of NMinimize: {er, ru} = FindMinimum[ Mean[MapThread[ EuclideanDistance, {combinedTransformation /@ data1, data2}]], {prx, pry, psx, psy, sx, sy, vecx, vecy, x, y, theta}, MaxIterations -> 200] gives roughly the same error as FindGeometricTransform. Commented Jan 10, 2019 at 18:10
• Well, the mean euclidean distance is minimized to the same amount as the FindGeometricTransform one. I don't know that the values are strange, I suppose you have better knowledge. You could consider giving the various variables default conditions if you have a decent idea of what the values might be. It's unclear to me that the message is specifically an error message - it looks like it's informational. Commented Jan 10, 2019 at 23:01

An alternative approach is to use NonlinearModelFit after reorganizing your data:

ClearAll[trans, model]
trans[sx_, sy_, tx_, ty_, θ_] := Composition[ScalingTransform[{sx, sy}, {1, 1}],
TranslationTransform[{tx, ty}], RotationTransform[θ, {1, 1}]]
model[sx_, sy_, tx_, ty_, θ_][x_] := Module[{h, v},
Flatten@Transpose@CoefficientArrays[trans[sx, sy, tx, ty, θ][{h, v}], {h, v}]. Array[x, 6]]

designmat = ArrayFlatten[{{#, 0}, {0, #}}] &@(Prepend[#, 1] & /@ data1);
response = Join @@ Transpose[data2];
nlm = NonlinearModelFit[Join[designmat, List /@ response, 2],
{model[sx, sy, tx, ty, θ][x], 0 <= θ <= 2 Pi}, {sx, sy, tx, ty, {θ, Pi}},
Array[x, 6]];

Row[{ListPlot[{data2, Transpose[Partition[#, Length[#]/2] &@nlm["PredictedResponse"]]},
PlotStyle -> {Directive[PointSize[Medium], Blue],
Directive[PointSize[.03], Opacity[.4], Red]}, ImageSize -> 400],
MapAt[Style[#, 16] &, nlm["ParameterTable"], {1}]}, Spacer[10]]


• This is a very interesting solution. Thanky you. Only th angle I do not understand: In the two solutions below both get Ulrich Neumann and Cesareo get an angle of about 0.38 to 0.39 degree.
– mrz
Commented Jan 17, 2019 at 10:33
• @mrz, the estimated angle is 359.607 Degrees (6.27632/Degree).
– kglr
Commented Jan 17, 2019 at 10:36
• Great. How do you receive this value from theta?
– mrz
Commented Jan 17, 2019 at 10:38
• @mrz, theta / Degree transforms theta (in radians) to degrees.
– kglr
Commented Jan 17, 2019 at 10:40
• All values which you receive are similar as of Ulrich Neumann and Cesareo except for the translation. Ulrich Neumann got {1502.57, -13.1302} and Cesareo got {1513.99, -7.23374}. You got {1496.06, -13.01}. My approximate estimation for the vertical translation is = Mean[data2[[All, 2]] - data1[[All, 2]]]=-11.95 and for the horizontal translation = Mean[data2[[All, 1]] - data1[[All,1]]]=1507.1.Which values of the there solutions are most reliable?
– mrz
Commented Jan 17, 2019 at 10:58