3
$\begingroup$

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}}]

enter image description here

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?

ADDENDUM:

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.

$\endgroup$
  • 1
    $\begingroup$ Try FindGeometricTransform . It 's not necessary to require a scaling point and/or a rotationpoint , that is the task ogf the fitting procedure. $\endgroup$ – Ulrich Neumann Jan 9 at 11:43
  • $\begingroup$ 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. $\endgroup$ – JimB Jan 9 at 14:28
  • $\begingroup$ 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$. $\endgroup$ – JimB Jan 9 at 14:32
  • 1
    $\begingroup$ I answered the question about what the error is in this question of yours. $\endgroup$ – Carl Lange Jan 9 at 17:57
  • $\begingroup$ Please see this follow up question: mathematica.stackexchange.com/questions/189592/… $\endgroup$ – mrz Jan 17 at 10:26
4
$\begingroup$

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*)
$\endgroup$
  • $\begingroup$ Please see the addendum. $\endgroup$ – mrz Jan 9 at 12:33
  • $\begingroup$ 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? $\endgroup$ – mrz Jan 9 at 13:11
  • 1
    $\begingroup$ Angle is around -.38 Degree. $\endgroup$ – Ulrich Neumann Jan 9 at 13:19
  • 1
    $\begingroup$ 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! $\endgroup$ – Ulrich Neumann Jan 9 at 15:19
  • 1
    $\begingroup$ @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... $\endgroup$ – Ulrich Neumann Jan 15 at 11:54
3
$\begingroup$

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}]

enter image description here

$\endgroup$
  • $\begingroup$ Please see the addendum. $\endgroup$ – mrz Jan 9 at 12:33
  • $\begingroup$ Please read also my last comment to Ulrich Neumann. $\endgroup$ – mrz Jan 9 at 13:44
2
$\begingroup$

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

enter image description here

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]
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – mrz Jan 9 at 13:32
  • $\begingroup$ Please read also my last comment to Ulrich Neumann. $\endgroup$ – mrz Jan 9 at 13:44
  • 1
    $\begingroup$ @mrz Included the plotting script $\endgroup$ – Cesareo Jan 9 at 14:43
  • 1
    $\begingroup$ 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. $\endgroup$ – Carl Lange Jan 10 at 18:10
  • 1
    $\begingroup$ 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. $\endgroup$ – Carl Lange Jan 10 at 23:01
2
$\begingroup$

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]]

enter image description here

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – mrz Jan 17 at 10:33
  • 1
    $\begingroup$ @mrz, the estimated angle is 359.607 Degrees (6.27632/Degree). $\endgroup$ – kglr Jan 17 at 10:36
  • $\begingroup$ Great. How do you receive this value from theta? $\endgroup$ – mrz Jan 17 at 10:38
  • 1
    $\begingroup$ @mrz, theta / Degree transforms theta (in radians) to degrees. $\endgroup$ – kglr Jan 17 at 10:40
  • $\begingroup$ 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? $\endgroup$ – mrz Jan 17 at 10:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.