3
$\begingroup$

Looking for alternatives or improvements for accomplishing the following:

Given the position of points A, B, and C in both the x-y and u-v planes, determine transformation functions to map values between the x-y and u-v planes.
Points A and B are on the v axis. Point C is on the u axis.

enter image description here

Purpose: The transform functions will be used to compute the uv coordinates of a point P given its xy coordinates, or its xy coordinates given its uv coordinates.

Below is what I have developed so far.

  1. Recommendations for improvements are welcome.
  2. Can this the result is expressed as an AffineTransform? This implementation uses (r + v).m but the AffineTransform uses r.m + v, perhaps there is a way to revise the implementation?
(*example data*)
{Axy, Bxy, Cxy, 
   Pxy} = {{0.2, 0.8}, {0.1, 0.15}, {0.8, 0.25}, {0.6, 0.7}};
{Auv, Cuv} = {{0, 75},  {100, 0}};

(*the unit vectors for the u-v plane, w.r.t. the x-y coordinate*)
Vxy = Normalize[ Axy - Bxy];
Uxy = Normalize[{Last[Vxy], - 
    First[Vxy] }]; (*axis orthogonal to Vxy*)
(*the origin of the u-v plane in x-y coordinates is the intersection \
of the v and u axis*)
UVoxy = First@(  
   Axy + k1 Vxy /. Solve[ Axy + k1 Vxy == Cxy + k2 Uxy , {k1, k2} ] );

(*scale factors for u-v axis*)
Uscale =   First@Cuv / EuclideanDistance[UVoxy, Cxy];
Vscale =   Last@Auv / EuclideanDistance[UVoxy, Axy];
mUVscale = {{Uscale, 0}, {0, Vscale}};
(*rotation from x-y to u-v*)
r1 = RotationMatrix[{{1, 0}, Uxy}] ;
(*combined rotation and scaling*)
mXYtoUV = r1 . mUVscale;
mUVtoXY = DiagonalMatrix[(1 / Diagonal@mUVscale)] . Transpose[r1]  ;

(*test: original xy points to uv*)
testUV =  (# - UVoxy) . mXYtoUV & /@ {Axy, Bxy, Cxy, Pxy} ;

(*test: computed uv points to xy*)
testXY01 = (# . mUVtoXY + UVoxy) & /@ testUV;

(*test: original uv points to xy*)
testXY02 = (# . mUVtoXY + UVoxy) & /@ {Auv,  Cuv} ;

(*report: xy to uv*)
Framed@Labeled[
  Grid[{ { "pt", Style["UV data", Bold],  Style["UV computed", Bold]}
    , { Column[{"A", "B", "C", "P" }]
     , {Auv, {"-", "-"}, Cuv, {"-", "-"}} // MatrixForm 
      , testUV // MatrixForm
     } 
    }
   , Frame -> All], Style["xy to uv", 16, Bold], Top]

(*report: uv to xy*)
Framed@Labeled[
  Grid[{ { "pt", Style["xy data", Bold],  Style["xy computed", Bold]}
    , { Column[{"A", "B", "C", "P" }]
     , {Axy, Bxy, Cxy, Pxy} // MatrixForm 
      , testXY01 // MatrixForm
     } 
    }
   , Frame -> All], Style["uv to xy", 16, Bold], Top]

results

results

Finally, here is the code to display a sketch of the system

(*sketch of the system*)
Graphics[{AbsolutePointSize[10]
  , Point[{Axy, Bxy, Cxy, Pxy}]
  , AbsoluteThickness[1]
  , AbsoluteDashing[10]
  , Darker@Blue
  , Arrowheads[{-0.05, +0.05}]
  , Arrow[{UVoxy - 0.35 Vxy, UVoxy + 0.65 Vxy}]
  , Arrowheads[+0.05]
  , Arrow[{UVoxy, UVoxy + 0.95 Uxy}]
  , Darker@Green
  , Text[Style["A", Bold, 18], Axy + {0.03, 0.03}, {-1, 0}]
  , Text[Style["B", Bold, 18], Bxy + {0.03, 0.03}, {-1, 0}]
  , Text[Style["C", Bold, 18], Cxy + {0.03, 0.03}, {-1, 0}]
  , Darker@Red
  , Text[Style["P", Bold, 18], Pxy + {0.03, 0.03}, {-1, 0}]
  , Darker@Blue
  , Text[Style["u", Bold, Italic, 14], UVoxy + 0.95 Uxy , {-1, 0}]
  , Text[Style["v", Bold, Italic, 14], UVoxy + 0.65 Vxy , {-1, -1}]
  , Text[Style["-v", Bold, Italic, 14], UVoxy - 0.35 Vxy , {0.5, 1}]
  }
 , Axes -> True
 , PlotRange -> {{-0.1, 1.2}, {-0.1, 1.1} }
 , Frame -> True
 ]

EDIT #1 and #2 Explaining the motivation for the question may help address any mistakes or ambiguities in the problem statement.

The transformations described above will be used in a program where the user imports a scanned image of a plot, like the one shown below, and the plot is displayed in a Mathematica Graphics object. It is likely that the axis of the plot and the axis of the Graphics object will not align. Rather than align the axes, the intent is to have transformation functions that will allow mapping values between the Graphics coordinate system (the xy axis) and the scanned image coordinate system (the uv axis).

As a first step, the system will have the user identify points A, B, and C with locator points; the Mathematica graphics object will "know" the xy values of those locator points. The system will also have the user manually enter the uv coordinates of the A,B,C points; the user can read those values from the plot. (The intent in this step is for the user to enter the minimum amount of information required to define the uv coordinate system. In my first implementation, only the "uv coordinates" of points A and C were required.)

With that information, the transformation functions can be computed. Those functions will be used to overlay the results of other calculations on the plot or to select features of the plot (with different locator points).

For example, in a second step, the user will also place locator points on points D and E in the sketch below. The Mathematica Graphics object will know the "xy coordinates" of these two points; a transform function will be used to compute their "uv coordinates" of points D and E. Then the "uv coordinates" of A,B,C, D and E will be known.

In a third step, the system will compute the "uv coordinates" of F and G, using the "uv coordinates" of A, B, C, D and E.

Then in a fourth step, the system will display the F and G points on the the Mathematica Graphics object. This will employ the other transform function which will convert the "uv coordinates" of the F and G points, computed in step 3, into "xy coordinates" required by the Mathematica Graphics object.

Below is an example of an scanned image that has coordinates that do not align with the horizontal and vertical and points A through G.

use of transform function

$\endgroup$
4
  • $\begingroup$ As stated, this is ambiguous. The xy origin is mapped to the uv origin, Projection[C-B,A-B]+B. But where should xy $(1,0)$ be mapped? To $C$? Your transform is unique only up to two independent scalars; one for the x axis and one for the y axis. I'm not sure what choice your code makes. $\endgroup$
    – Adam
    Apr 17, 2021 at 21:21
  • $\begingroup$ the xy origin is offset from the uv origin. the code computes the uv origin in xy coordinates, by computing the intersection of the u and v axis, with the statement: UVoxy = First@( Axy + k1 Vxy /. Solve[ Axy + k1 Vxy == Cxy + k2 Uxy , {k1, k2} ] ); in the sketch the x-y coordinates are the frame coordinates. $\endgroup$
    – user6546
    Apr 17, 2021 at 22:59
  • $\begingroup$ Just to be clear: the data of the problem consists of the coordinates of A, B and C. The "uv axes" are the inferred from the points. Then what you are asking for is an affine transformation that outputs the coordinates of A, B and C in the "uv axes". Is this correct? $\endgroup$
    – A.G.
    Apr 18, 2021 at 12:09
  • $\begingroup$ Thanks for your comment. Am asking for two affine transformations. One transform will be used to determine the "uv coordinates" of a point P given only its "xy coordinates." The other transform will be used to determine the "xy coordinates" of a point F given only its "uv coordinates." Will update the question. $\endgroup$
    – user6546
    Apr 18, 2021 at 13:11

4 Answers 4

3
$\begingroup$

Try FindGeometricTransform to describe the transformation {x,y}<->{u,v}.

Therefore it's necessary to know the three points A,B,C (I added Buv ).

{Axy, Bxy, Cxy} = {{0.2, 0.8}, {0.1, 0.15}, {0.8, 0.25}}
{Auv, Buv, Cuv} = {{0, 75}, {0, -45.378 }, {100, 0}, {0, -45.378 }}

trafo = FindGeometricTransform[  {Auv, Buv, Cuv}, {Axy, 
Bxy, Cxy}] [[2]]

Last line {0,0,1} of the transformationmatrix

TransformationMatrix[trafo]
(*{{146.067, -22.4719, -11.236}, 
{39.2312, 179.161,-76.1753}, 
{0., 0.,1.}}*)  

indicates an affine transformation!

trafo[{Axy, Cxy, Bxy}] // Chop
(*{{0, 75.}, {100., 0}, {0, -45.378}}*)

trafo[ {0.6, 0.7}]
(*{60.6742, 72.7764}*)

The inverse transfomation follows to InverseFunction[trafo].

This approach works in the same way if more than 3 points are known.

$\endgroup$
2
  • $\begingroup$ Thanks! FindGeometricTransform works well for this question. $\endgroup$
    – user6546
    Apr 18, 2021 at 18:18
  • $\begingroup$ You're welcome. It even works for distorted images of coordinate systems! $\endgroup$ Apr 19, 2021 at 6:17
1
$\begingroup$

To clear up ambiguities, the following code finds the transformation from the points in xy space $(0,1),(1,0)$ and $(0,0)$ to the points in uv space $a,b$ and $c$ projected onto $\overline{ab}$ respectively.

With[{a = {ax, ay}, b = {bx, by}, c = {cx, cy}, 
m = {{mx, mxy}, {myx, myy}}}, 
With[{v = Projection[c - b, a - b, Dot] + b}, 
m /. Simplify@
First@Solve[{a == m.{0, 1} + v, c == m.{1, 0} + v}, Flatten@m]]]

The result is a matrix m that can be used to make an affine transform, i.e. AffineTransform[{%,v}], if it's still in the inner With.

The math isn't too bad to do without Solve. Consider

With[{v=Projection[c - b, a - b, Dot] + b},m=AffineTransform[{
RotationMatrix[ArcTan @@ (a - b)].
ScalingMatrix[Norm[# - v] & /@ {a, c}], v}]]

This latter code doesn't do negative scaling, so it is broken when the points are oriented the wrong way (c has to be to the right etc.). You could entirely remove the ScalingMatrix...

With regards to the end goal, I did something similar here Detect and correct sheared squares in image, though my code was nasty.


Also I really enjoy doing interactive things with MMA, I'm interested in developing this further with Manipulate etc.

$\endgroup$
1
$\begingroup$

Assume a,b,c are given as:

{a, b, c} = {{0.2, 0.8}, {0.1, 0.15}, {0.8, 0.25}}; 

The origin of the uv system is on a-b: orig= b+ x(a-b), where we need to determine x. This is done by using the fact that (c-orig) is perpendicular to (a-b):

orig = b + x (a - b) /. Solve[(b + x (a - b) - c).(a - b) == 0, x][[1]];

With the origin, it is easy to determine unit vectors along u and v:

{eu, ev} = {(a - orig)/Norm[a - orig], (c - orig)/Norm[c - orig]};

To map coordinates x/y to u/v we need first to subtract orig from x/y und then project onto eu and ev. To map from u/v to x/y we need to multiply eu by u and ev by v, add both together and finally add orig:

xy2uv[p : {x_, y_}] := {eu, ev}.(p - orig);
uv2xy[p : {u_, v_}] := p.{eu, ev} + orig;

As a test, the u/v coordinates of c should be {0,v}:

t=xy2uv[c]
(*{-1.38778*10^-17, 0.676654}*)  

And if we back map t we should again get the coordinates of c:

uv2xy[t]
(*{0.8, 0.25}*)

Addendum

If the u/v coordinates have a different scale, the transformation is:

xy2uv[p : {x_, y_}] := {eu, ev}.(p - orig) scaleuv;
uv2xy[p : {u_, v_}] := (p/scaleuv).{eu, ev} + orig;
$\endgroup$
2
  • $\begingroup$ thanks for the clear explanations. how to extend this so that the scales of the two coordinate systems are different? In the example data, the value of C in the uv system is {100,0}. Also, this statement places c on the v axis, when it should be on the u axis. {eu, ev} = {(a - orig)/Norm[a - orig], (c - orig)/Norm[c - orig]}; $\endgroup$
    – user6546
    Apr 18, 2021 at 13:50
  • $\begingroup$ You are right, nobody is perfect. eu and ev should be reversed: . {eu, ev} = {(c - orig)/Norm[c - orig],(a - orig)/Norm[a - orig], }. To change the scale of the u/v coordinates you have to multiply/ divide the u/v coordinates as shown in the addendum to my answer. $\endgroup$ Apr 18, 2021 at 16:21
1
$\begingroup$

Here is a solution that uses AffineTransform and Solves for coefficients. You can specify a scale.

{a, b, c, p} = {{0.2, 0.8}, {0.1, 0.15}, {0.8, 0.25}, {0.6, 0.7}};
(* Specify scale here; for example, u and v's lengths are 1/2 *)
scale = 1/2; 
(* Define the parameters *)
ClearAll[u0, ux, uy, v0, vx, vy, sol, UV];
m := {{ux, uy}, {vx, vy}};
V := {u0, v0};
UV = AffineTransform[{m, V}];
(* Find values of the parameters *)
sol = FindInstance[
   (* Call w the origin of the uv coordinates *)
   (* a and b are on the wv axis, c is on the wu axis *)
   UV[a][[1]] == 0 && UV[b][[1]] == 0 && UV[c][[2]] == 0 &&
    (* normalize / scale *)
    ux^2 + uy^2 == 1/scale^2 && 
    vx^2 + vy^2 == 1/scale^2 &&
    (* u and v are orthogonal: scalar product is 0 *)
    First@m . Last@m == 0,
       {u0, ux, uy, v0, vx, vy}];
(* assign the results to parameters *)
{u0, ux, uy, v0, vx, vy} = {u0, ux, uy, v0, vx, vy} /. First@sol;
TableForm[{UV[{x, y}], UV[a], UV[b], UV[c], UV[p]},
 TableHeadings -> {{"Formula", "a", "b", "c", "p"}, {"u", "v"}}]

enter image description here

(* Inverse transformation *)
XY = InverseFunction[UV];
w = XY[{0, 0}]; (* origin of the uv axes *)
TableForm[{XY[{u, v}], w},
  TableHeadings -> {{"Formula", "w"}, {"x", "y"}}]

"Just checking -- the following should = {x,y}:"
XY[UV[{x, y}]] // Simplify

enter image description here

(* sketch of the system *)
Graphics[{
  AbsolutePointSize[5],
  {Point[{a, b, c}]},
  {Blue, Point[p]},
  {Green, AbsolutePointSize[10], Point[w]},
  (* Extend both u and v axes a bit *)
  {Black, Dashed, 
   Line[{a - (b - a), b + (b - a)}]},
  {Black, Dashed, Line[{w - (c - w), c + (c - w)}]},
  {Black, Arrow[{w, XY[{1, 0}]}]},
  {Black, Arrow[{w, XY[{0, 1}]}]},
  {Blue, Point[{0, 0}]},
  {Blue, Arrow[{{0, 0}, {1, 0}}]},
  {Blue, Arrow[{{0, 0}, {0, 1}}]}
  }, Axes -> True, PlotRange -> {{-0.1, 1.5}, {-0.1, 1.5}}, 
 AspectRatio -> 1, Frame -> True]

enter image description here

$\endgroup$
2
  • $\begingroup$ thanks! very helpful. especially seeing how the constraints can be used to create the parameters of the Affine transformation. The scaling is unclear to me, especially if the u and v axes are not the same. Will look for other references for that, though recommendations are welcome. $\endgroup$
    – user6546
    Apr 19, 2021 at 22:01
  • $\begingroup$ @user6546 Glad you find it helpful! A good starting point is Wolfram's page on affine transformations. (mathworld.wolfram.com/AffineTransformation.html). A good way to see all this is to interpret the affine transform not so much as a transformation that changes objects as a way of obtaining coordinate in two different affine frames (en.wikipedia.org/wiki/Affine_space#Coordinates). (btw, the Wikipedia page is in need of explanations :) ) $\endgroup$
    – A.G.
    Apr 19, 2021 at 23:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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