Given a set of three points in coordinate frame A which are not
co-linear, and the same three physical points in coordinate frame B,
generate a transformation matrix to convert between frame A and frame
B.
Notes:
1) Each point is identified in each frame; that is P1a may be (1,2,3),
while P1b is (1,4,2) - the  point is the same physical point in space,
only the frame (and hence the coordinates) has changed.

2)in reality for this application, the frames may/may not be orthogonal.

3)the unit basis for each frame may / may not be identical.

Questions:
A) Is three points sufficient to solve this in 3d? Do we need three
points plus the origin for each system?

B) What are the steps required to generate the transformation matrix
given the coordinates of the points in both frames?  Please provide
Psuedocode if possible

C) Will this solution work for non-orthogonal frames?

Additional background:
The application relates to viewing a sample under a high resolution
microscope. The microscope has a "stage", or movable platform, that
can move the sample in x,y,z,rotation, and tilt axes.  However, once
the system is aligned, the only axes which will be moved are the x and
Y axes, which are (mostly) perpindicular to the field of view (when
the stage is at zero degrees tilt).

The stage is positioned by the stage electronics in units of mm - i.e.
1.23456mm. This could be considered frame A.  The sample (typically a
very flat semiconductor wafer) will have a different origin, and the
points we are aligning to could either be in mm, or could be in some
other unit, such as cell or grid units (I believe the units would be
irrelevent, as long as we can convert between the two).

Once aligned, the user could then enter the sample coordinate, the
transformation calculated, the correct stage cooridnates generated,
and we will drive to the correct sample feature.

The alignment process consists of the user driving to a known feature
on the sample (lets say P1b) and centers the point on the screen.  The
user clicks a button, and enters the sample coordinates (in terms of
the sample).  The stage electronics will already know the stage
position (we will call it P1a).  So now we have a single point, with
coordinates given with respect to two different frames.  This
procedure would continue for as many points as are required (I presume
three, perhaps four?)

The axes of motion of X and Y are not purely orthogonal - they are
close, but any answer which assumes perfect orthogonality most likely
will not produce acceptable results.
Thank you for your assistance.

Hi, eleceng-ga:

As I sketched in the Comment posted below, coordinates of 4 points
not lying in a common plane are required in both coordinate frames
to uniquely define a conversion between those coordinate systems.

One of the points might be an origin in one or the other coordinate
system, but this is not essential as we will see from the method of
solving for the conversion coefficients.

If you have only 3 points whose coordinates are known in both frames,
and these 3 points are not colinear, so that there is exactly one
plane passing through all 3 points, then you can still define a
conversion between the coordinate systems that is known to be valid
in that plane.  However it will not be possible to determine what
the conversion is outside of the plane without further information.

As I'll show below, the conversion is stated most simply in terms
of a formula that uses both a matrix M and a vector offset.  I'll
explain how the formula is derived, and how to compute the matrix
and the vector offset.  Then I'll illustrate the computation with
an example using Excel to do the "heavy lifting".

*  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *

Let four points in space be assigned three-dimensional coordinates
according to two different coordinate frames, respectively:

  point A:   (ax1,ay1,az1)   and   (ax2,ay2,az2)
  
  point B:   (bx1,by1,bz1)   and   (bx2,by2,bz2)
  
  point C:   (cx1,cy1,cz1)   and   (cx2,cy2,cz2)
  
  point D:   (dx1,dy1,dz1)   and   (dx2,dy2,dz2)
  
[All coordinate entries are assumed to be real numbers.]
  
The units (scales of measurement) may be different between the two
coordinate systems, and even from one axis to another within the same
coordinate system.  Moreover the coordinate axes within a coordinate
system need not be orthogonal (at right angles, perpendicular) to
one another.   We require only that in each coordinate system the
coordinates obey the laws of vector arithmetic (linear algebra),
e.g. addition of coordinates corresponds to vector addition (which
is therefore commutative), and scalar multiplication is identified
with coordinate-wise multiples (and is thus "compatible" with vector
addition in the sense of the distributive law).

*  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *

Geometrically the condition that the 4 points don't lie in one commmon
plane means that they form the vertices of a tetrahedron (though not
necessarily a regular tetrahedron).  Points in this solid tetrahedron
can be expressed by an application of coefficients to the coordinates
of the points in either system, where the coefficients are nonnegative
and sum to 1.  Furthermore points outside of the tetrahedron can be
expressed in a similar way by using coefficients that sum to 1 but are
not all nonnegative.

That is, for any real numbers such that ra + rb + rc + rd = 1, these
expressions give the same point in both coordinate systems:

1st coordinate system:
  
  ra*(ax1,ay1,az1) + rb*(bx1,by1,bz1) + rc(cx1,cy1,cz1) + rd*(dx1,dy1,dz1)

2nd coordinate system:
  
  ra*(ax2,ay2,az2) + rb*(bx2,by2,bz2) + rc(cx2,cy2,cz2) + rd*(dx2,dy2,dz2)

Expressions of these kinds are called "barycentric coordinates" and are
discussed in greater detail here for the two-dimensional case (a plane):

[The uses of homogeneous barycentric coordinates in plane euclidean geometry]
(by Paul Yiu)
http://www.math.fau.edu/yiu/barycentric.pdf

and in brief detail here for the general n-dimensional case:

[Barycentric Coordinates]
http://www.cap-lore.com/MathPhys/bcc.html

If the problem we needed to solve was to convert coordinates for only one
point from one coordinate frame to another, we could approach it by first
solving for ra, rb, rc, rd satisfying the condition of summing to 1 and of
providing the right coordinates of the point in the given coordinate frame,
then combining the alternative coordinates of A,B,C,D using the _same_
coefficients ra, rb, rc, rd to obtain this point's coordinates in the other
frame.

However if we need to convert a large number of points from one coordinate
frame to the other, it will pay us to invest a bit of effort up front in
deriving a formula that makes the conversion more convenient.

*  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *

As a convenience I will label the coordinates that we want to convert from
as being those of the "known" coordinate frame, while the target coordinates
that we want to convert _to_ are those of the "unknown" coordinate frame.
Of course this choice of terms is relative to whichever direction we want the
conversion to go, and doesn't mean anything intrinsic about either frame.

The mathematical name for the conversion from one system of coordinates to
another as outlined here is an "affine" mapping.  If the two coordinate
frames shared a common origin, then this would be the special case of a
"linear" mapping.

The general formula for such affine mappings in three dimensions is:

   (qx,qy,qz)  =  (px,py,pz) * M  + (Ox,Oy,Oz)

where (px,py,pz) are the coordinates in the "known" frame and (qx,qy,qz)
are the coordinates in the "unknown" frame.  M is a 3x3 matrix (to be
determined) and (ox,oy,oz) is a vector "offset" (also to be determined).

The special case of a linear mapping could be identified by a zero offset:

   (Ox,Oy,Oz)  =  (0,0,0)

In the general case of an affine mapping, the offset (ox,oy,oz) expresses
the origin of the "known" coordinate frame in terms of the "unknown" frame,
and thus the two origins are _not_ the same point if (ox,oy,oz) is nonzero.

With the 4 points A,B,C,D (not lying in a common plane) we have exactly the
right amount of data to uniquely determine the nine entries of the matrix M
together with the three offset coordinates.

One way to do this is to set up a system of 12 simultaneous linear equations
in the 12 unknowns (nine entries of M together with three offset coordinates).

However a simple trick allows us to "eliminate" the three offset coordinates
in the first step, and solve a system of 9 linear equations in 9 unknowns.

But let's start from the beginning.  If we plug in the known coordinates of
any single point, say A, to the formula above, and assuming the 1st coordinate
system is the "known" from while the 2nd coordinate system is the "unknown"
frame, then:

   (ax2,ay2,az2)  =  (ax1,ay1,az1) * M  +  (Ox,Oy,Oz)
   
This single "vector" equation is the equivalent of three "scalar" equations:

    ax1*M(1,1) + ay1*M(2,1) + az1*M(3,1) + Ox  =  ax2
    
    ax1*M(1,2) + ay1*M(2,2) + az1*M(3,2) + Oy  =  ay2
    
    ax1*M(1,3) + ay1*M(2,3) + az1*M(3,3) + Oz  =  az2

Applying the same principle to each of the three remaining points B,C,D
will give 3 more vector equations, equiv. to 9 more scalar equations:

    bx1*M(1,1) + by1*M(2,1) + bz1*M(3,1) + Ox  =  bx2
    
    bx1*M(1,2) + by1*M(2,2) + bz1*M(3,2) + Oy  =  by2
    
    bx1*M(1,3) + by1*M(2,3) + bz1*M(3,3) + Oz  =  bz2

    cx1*M(1,1) + cy1*M(2,1) + cz1*M(3,1) + Ox  =  cx2
    
    cx1*M(1,2) + cy1*M(2,2) + cz1*M(3,2) + Oy  =  cy2
    
    cx1*M(1,3) + cy1*M(2,3) + cz1*M(3,3) + Oz  =  cz2

    dx1*M(1,1) + dy1*M(2,1) + dz1*M(3,1) + Ox  =  dx2
    
    dx1*M(1,2) + dy1*M(2,2) + dz1*M(3,2) + Oy  =  dy2
    
    dx1*M(1,3) + dy1*M(2,3) + dz1*M(3,3) + Oz  =  dz2

My notation was chosen so that in this simultaneous linear system, the
small letters refer to "given" values and the capital letters are "to be
determined".

If the top three equations are subtracted from the respective every third
equation below it, the offset coordinates are eliminated as variables in
this system:

    (bx1-ax1)*M(1,1) + (by1-ay1)*M(2,1) + (bz1-az1)*M(3,1)  =  bx2-ax2
    
    (cx1-ax1)*M(1,1) + (cy1-ay1)*M(2,1) + (cz1-az1)*M(3,1)  =  cx2-ax2
    
    (dx1-ax1)*M(1,1) + (dy1-ay1)*M(2,1) + (dz1-az1)*M(3,1)  =  dx2-ax2


    (bx1-ax1)*M(1,2) + (by1-ay1)*M(2,2) + (bz1-az1)*M(3,2)  =  by2-ay2
    
    (cx1-ax1)*M(1,2) + (cy1-ay1)*M(2,2) + (cz1-az1)*M(3,2)  =  cy2-ay2
    
    (dx1-ax1)*M(1,2) + (dy1-ay1)*M(2,2) + (dz1-az1)*M(3,2)  =  dy2-ay2


    (bx1-ax1)*M(1,3) + (by1-ay1)*M(2,3) + (bz1-az1)*M(3,3)  =  bz2-az2
    
    (cx1-ax1)*M(1,3) + (cy1-ay1)*M(2,3) + (cz1-az1)*M(3,3)  =  cz2-az2
    
    (dx1-ax1)*M(1,3) + (dy1-ay1)*M(2,3) + (dz1-az1)*M(3,3)  =  dz2-az2
    
It should be noticed that this is not simply a system of nine equations,
but actually three "decoupled" systems of 3 equations in 3 unknowns, and
moreover sharing the same "matrix of coefficients" on the left-hand side.

A compacted "matrix" formulation of the system can then be stated:

      P * M  = Q

where:


        / (bx1-ax1) (by1-ay1) (bz1-az1) \
       |                                 |
   P = |  (cx1-ax1) (cy1-ay1) (cz1-az1)  |
       |                                 |
        \ (dx1-ax1) (dy1-ay1) (dz1-az1) /


        / (bx2-ax2) (by2-ay2) (bz2-az2) \
       |                                 |
   Q = |  (cx2-ax2) (cy2-ay2) (cz2-az2)  |
       |                                 |
        \ (dx2-ax2) (dy2-ay2) (dz2-az2) /


The geometric condition on A,B,C,D not lying in a common plane amounts
to the matrix algebra condition that both P and Q are invertible, so:

         M  = P¯¹ * Q

gives an explicit expression for M, and returning to the original three
equations one easily solves for the offset (Ox,Oy,Oz):


   (ax2,ay2,az2)  =  (ax1,ay1,az1) * M  +  (Ox,Oy,Oz)

  
   (Ox,Oy,Oz)  =  (ax2,ay2,az2)  -  (ax1,ay1,az1) * M


*  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *

Let's do an example, using Excel to do the matrix arithmetic for us.

We are supposed to have two coordinate frames, and four points whose
coordinates are given for both coordinate frames (no plane containing
all four of the points), say:

            1st "known" frame        2nd "unknown" frame

   A:        (0.1, 1.2, 3.5)           (1.2, 0.6, -0.8)
   
   B:        (2.1, 1.1, 2.0)           (1.0, 1.0, -0.5)
   
   C:        (1.8, 0.2, 0.9)           (-0.3, 1.5, 1.3)
   
   D:        (3.1, 1.1, 0.0)           (-1.0, 0.2, 2.3)

Compute the matrices P and Q as above by subtracting the coordinates of
A from the respective rows of coordinates given by B,C,D in each frame:


        /  2.0  -0.1  -1.5  \
       |                     |
   P = |   1.7  -1.0  -2.6   |
       |                     |
        \  3.0  -0.1  -3.5  /


        / -0.2   0.4   0.3  \
       |                     |
   Q = |  -1.5   0.9   2.1   |
       |                     |
        \ -2.2  -0.4   3.1  /


Excel has a built-in functions for matrix inversion and multiplication
that can be used within what is called an "array formula" that computes
a whole rectangle of cells at once.  Using this I found approximately:


        / 0.983471074   0.666115702  -1.360330579  \
       |                                            |
   M = | -0.70661157   -1.673553719   0.995867769   |
       |                                            |
        \ 1.491735537   0.733057851  -2.080165289  /


using M  = P¯¹ * Q, and that approximately:


   (Ox,Oy,Oz) = (-3.271487603, -0.024049587, 5.421570248)


using (Ox,Oy,Oz)  =  (ax2,ay2,az2)  -  (ax1,ay1,az1) * M.

I then verified, using additional array formulas that the conversions
from the 1st "known" coordinate frame to the 2nd "unknown" coordinate
frame give numerical results accurate to the display of Excel.

*  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *

TIPS on Excel Array Formulas:

I put the matrix P into cells A1:C3 and matrix Q into cells E1:G3.

To compute M  = P¯¹ * Q, I highlighted the cells I1:K3, typed this
formula into the entry box along the top of the spreadsheet:

  = MMULT(MINVERSE(A1:C3),E1:G3)

and then hit Ctrl-Shift-Enter.  This last is what assigns the cells
as a whole to the array formula, which is then displayed with curly
braces like this:

  {= MMULT(MINVERSE(A1:C3),E1:G3)}

to distinguish it from a regular formula.  But the curly braces are
not entered manually; Excel just displays the formula that way after
you hit Ctrl-Shift-Enter.

To compute the offset, I put (ax1,ay1,az1) in cells A5:C5, and I put
(ax2,ay2,az2) in cells E5:G5.  I then created the following array
formula on cells I5:K5:

  {= E5:G5 - MMULT(A5:C5,I1:K3)}
  
which again displays with curly braces that I did not actually type.
This last is the Excel equivalent of:

  (Ox,Oy,Oz)  =  (ax2,ay2,az2)  -  (ax1,ay1,az1) * M
  
and its matrix multiplication references the value of M computed by my
first array formula above.


Please let me know if some part of my Answer would benefit from further
Clarification.


regards, mathtalk-ga

---

Request for Answer Clarification by eleceng-ga on 18 Aug 2005 11:07 PDT

Thank you for your detailed answer - I am working through the info right now.

It occured to me that there is another scenario that will exist:

Typically the user coordinates will be 2d; that is, all points of
reference will be constrained to a flat plane.  However, the actual
sample may be slightly tilted, resulting in an actual z axis shift as
you traverse the sample(z axis vector is parallel to the view).  I
think that the user would enter x and Y user coordinates, and some
arbitrary z coordinate (0?) would be assumed.  The actual Z
coordinates of the stage may or may not change, depending on how flat
the sample is, how far apart the reference points are, and the
referencing methodology:

A)We can drive X and Y only, correcting for the change in Z (focus)
manually as we move about (depending on how far away from 0 degrees
the tilt is). X and Y errors are corrected automaticaly, as they are a
function of the reference points we pick - if tilt were at 45 degrees,
the Y axis would be compressed by 0.707, or 0.707 inches stage travel
= 1 inch sample travel (NOT in user coordinates but physically).  This
effectively reduces the problem to strictly a 2d affine
transformation?

B) We can drive X and Y, but correct the Z (focus) by driving the Z
axis until the image is once again in focus.  This gives us a 3d
change in stage coordinates, but only a 2d change in user coordinates.
 However, a change in Z is not guarenteed - the sample may be
perfectly flat....

As the sample is traversed to pick the alignment points, will this
method work if all the user reference points are in the same plane?
How could I modify it so that it would work?

Thanks again for your help......

---

Clarification of Answer by mathtalk-ga on 18 Aug 2005 18:02 PDT

Hi, eleceng-ga:

The same sort of calculations can be applied to the case of
determining a mapping from two "user" coordinates to three "physical"
coordinates.  More details on this below, but first a couple of
general remarks I should add to my original Answer.

1.  I used "row" vectors in the write-up mainly because they are easy
to type out in this text-based medium.  You could convert everything
to column vectors just by taking the transposes in the formulas
(which, if eleceng stands for "electrical engineer", you will probably
remember also reverses the order of matrix multiplications).

2.  This calculation assumes the measured or given coordinates are
exact.  In many application one has the same sort of need to develop a
conversion based on coordinates that have errors of measurement.  The
specific case of input coordinates known exactly and the output
coordinates containing measurement errors is widely addressed by a
least-squares fit, the precise generalization of the fitting of a
straight line to data points whose y-coordinates are subject to
experimental error.

Now back to the two-dimensional "user" coordinates case:  As I briefly
noted in at the start of my Answer, the availability of three points
(not all in a common line) is enough to allow you to define a mapping
valid on the plane containing those three points.  The mapping can in
fact output three "physical" coordinates even though this output
depends functionally on just two "user" coordinates, and the affine
mapping will work very much the same as before.

Let's call (u,v) the user coordinates, and (x,y,z) the physical coordinates.

Suppose there are three points A,B,C (not all in one line) for which
we know both kinds of coordinates:

   A:   (ua,va)     (xa,ya,za)

   B:   (ub,vb)     (xb,yb,zb)

   C:   (uc,vc)     (xc,yc,zc)

Then the mapping from user coordinates to physical coordinates will be:

   (x,y,z)  =  (u,v)*M  +  (Ox,Oy,Oz)

where we need to determine the 2x3 matrix M and the vector (Ox,Oy,Oz)
from the information contained in the three pairs of corresonding
coordinates.

If we proceed as before, the equivalent calculation is now:

        /  (ub-ua)  (vb-va)  \
   P =  |                    |
        \  (uc-ua)  (vc-va)  /


        /  (xb-xa)  (yb-ya)  (zb-za)  \
   Q =  |                             |
        \  (xc-xa)  (yc-ya)  (zc-za)  /


and:   M  =  P¯¹ * Q 

    (Ox,Oy,Oz) = (xa,ya,za) - (ua,ub) * M


If all the Z-coordinates in the "physical" coordinate frame are zero,
then naturally the third column of Q (and thus the third column of M)
will be zeroes, and so also will Oz = 0.


regards, mathtalk-ga

---

Request for Answer Clarification by eleceng-ga on 19 Aug 2005 09:39 PDT

Hello Again,

Concerning the system converting between 3d and 2d - 

I was able to generate the 2x3 matrix for conversion from "user" to
"stage" coordinates, but then to go back, wouldn't I need to invert
the non-square matrix?  Any ideas on the best way to handle this?

Thanks again for you help!

---

Clarification of Answer by mathtalk-ga on 19 Aug 2005 17:14 PDT

In the context you are talking about, namely the z-coordinate of the
physical frame represents a slight tilt (possibly none at all), just
discard the z-coordinate for the purpose of going "backwards".  The x-
& y-coordinates should suffice to give an invertible 2x2 matrix for
determining the "user" u- & v-coordinates.

Thanks for the generous rating (and tip)!

regards, mathtalk-ga

Answer speaks for itself.....

A coordinate frame usually consists of a choice of origin together
with a basis (or set of coordinates) relative to that origin.

If two coordinate frames share a common origin, then the conversion
from one to the other is a matrix multiplication.  Thus, in three
dimensions you have to determine the nine entries of a three by three
matrix, and three "independent" points with the respective coordinates
in both coordinate frames is exactly sufficient to determine them.

But if the coordinate frames do not share a common origin, then the
conversion is an "affine" mapping, rather than a linear one, and
besides a matrix multiplication there is a "shift" or vector addition
to pin down.  In essence one needs the coordinates in both systems of
a fourth point.

The four points should satisfy:

 - No more than two in any line.

 - No more than three in any plane.

The equations for the matrix & vector entries are then a simple linear
system defined by the coordinates of the points.  This works even if
the mapping is not "orthogonal", which means (or should mean) angle
preserving in the case of a linear map and (more generally) distance
preserving in the case of an affine map.

regards, mathtalk-ga