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 |