Coordinate transformation can this be done with AffineTransform?

Question

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.

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

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.

Answer 1

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.

Answer 2

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.

Answer 3

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;

Answer 4

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

(* 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

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