From: robinson@mthcsc.wfu.edu (Stephen Robinson)
Subject: A Problem in Projective Geometry
Date: 6 Sep 1999 09:30:03 -0500
Newsgroups: sci.math.research
Keywords: Determine Euclidean motion from shadows of 5 moved points
I am interested in finding references on the following problem and its
relatives. It is relevant to some research that I am doing in medical
imaging. I have essentially solved the problem using elementary
methods, but I suspect that an expert in Projective Geometry might
know of a more elegant treatment. References to specific articles and
books would be wonderful.
The Problem: Suppose that 5 noncoplanar points in R^3 are given, call
them x1,x2,x3,x4 and x5. (I think of them as the corners of a
polyhedron.) An unknown rotation and shift, i.e. a rigid motion, move
the points to some position above the x-y plane. A light shines from
an unknown point, p. Assume that p is high enough above the x-y plane
so that 5 point shadows appear on the x-y plane, call them u1, u2, u3,
u4 and u5. Assume that the coordinates of u1,...,u5 are known, and
that they correspond to the points x1,...,x5 in the obvious way. Here
is the question: Under what conditions can I use the known
information, (x1,...,x5,u1,...,u5), to solve for the unknown rotation,
shift, and projection point p?
Thank You,
Steve
==============================================================================
From:
Subject: Re: A Problem in Projective Geometry
Date: Mon, 6 Sep 1999 00:03:05 -0500 (CDT)
Newsgroups: [missing]
To: robinson@mthcsc.wfu.edu
A dimension-counting argument suggests that you should be able to recoup
the data up to finite ambiguity, in generel. Indeed, in order to specify
rotation+shift+point_p requires 3+3+3=9 variables; the positions of the
five shadows u_i provides 10 equations. Thus the problem is actually over-
specified; there will be an equation in 10 variables which is to be
satisfied when the variables are equal to the coordinates (u_i)_x and (u_i)_y.
In a real setting you would have to use a least-squares criterion to
decide what 9 equations will be used to determine the variables.
dave
==============================================================================
From: Dave Rusin
Subject: Re: A Problem in Projective Geometry
Date: Tue, 21 Sep 1999 09:47:09 -0500 (CDT)
Newsgroups: [missing]
To: robinson@mthcsc.wfu.edu
A couple of days ago I responded to a request you had made earlier on
sci.math.research:
> Suppose that 5 noncoplanar points in R^3 are given, call
> them x1,x2,x3,x4 and x5. (I think of them as the corners of a
> polyhedron.) An unknown rotation and shift, i.e. a rigid motion, move
> the points to some position above the x-y plane. A light shines from
> an unknown point, p. Assume that p is high enough above the x-y plane
> so that 5 point shadows appear on the x-y plane, call them u1, u2, u3,
> u4 and u5. Assume that the coordinates of u1,...,u5 are known, and
> that they correspond to the points x1,...,x5 in the obvious way. Here
> is the question: Under what conditions can I use the known
> information, (x1,...,x5,u1,...,u5), to solve for the unknown rotation,
> shift, and projection point p?
I don't have the response I mailed you, but I'm pretty sure I was
talking nonsense then.
Well, I ran the calculations recently to see how this goes. I think it's
an interesting problem! I also think it's possible to analyze, at least
for any particular configuration of your (x1, ..., x5, u1, ..., u5).
If done intelligently, it seems one can make quite a bit of progress
simplifying the equations. However, the task of obtaining a full
symbolic solution to the general problem is beyond the capacity of
my computer (although not in principle impossible with today's tools).
Still, I'm pleased that the computations below describe the situations
fairly well and allow for accurate and rapid calculation.
I would be curious to know how the problem arose and how you might
use the solution.
================================
Let me show you how I analyzed this. I have to change to a different
notation, sorry.
We suppose we have 5 points which we rotate, translate, and then
project from a point P onto the xy plane, to get 5 points in the
plane. The problem is to use the first and last sets of points to
(1) determine the rotation, translation, and the projection point P.
(2) describe constraints on the sets of 5 initial and 5 final points
which allow a solution to exist.
I will show that the right answer to (2) seems to be, there is a single
polynomial in the initial data which must vanish in order for the
problem to be solvable; if that polynomial vanishes, the solutions to (1)
(not unique!) can be computed.
I think it's easier to ignore the initial coordinates [x_i, y_i, z_i] of the
five initial points and instead focus on the coordinates [X_i, Y_i, Z_i]
of the points after the rotation and translation. Except in degenerate
cases, five points XYZ will be the image of the original points xyz
under a unique (computable) Euclidean transformation iff the C(5,2)=10
interpoint distances among the XYZ match those of the xyz.
Here's the easy way to write equations which state that the interpoint
distances are correct: number the points with subscripts 0, 1, 2, 3, 4
just to highlight one of them as a reference point, and then let
V_i = [X_i - X_0, Y_i - Y_0, Z_i - Z_0] be the displacement vector
to the i-th point (so V_0 = [0,0,0] ). Then the ten distances can
be computed from the ten inner products < V_i , V_j > (1 <= i <= j <= 4)
and vice-versa. So we compute these numbers
IP_{i,j} = < v_i, v_j >
analogously from the initial five points, and then have equations which
state that {X_0, Y_0, ..., Y_4, Z_4} be coordinates of five points which
are a Euclidean motion away from the initial points:
(*) < V_i, V_j > = IP_{i,j}, (1 <= i <= j <= 4).
Observe that this appears to be ten equations in the 15 variables, but
in fact they are not all independent equations. Indeed, the ten numbers
IP_{i,j} are not completely arbitrary: they correspond to the
innerproducts among five actual points in 3-space iff the matrix IP
has signature (3, 0, 1), and so in particular the equation det(IP)=0
gives a relation among the ten numbers IP_{i,j}. In principle you
could use it to eliminate any one of the ten values IP_{i,j} in favor
of the other 9. I won't do so explicitly.
OK, so now I have equations which state that the XYZ's are the image of
the xyz's. I also need equations which say they project to the original
five given points in the plane.
So let the point of projection be, say, [p,q,r]. (These are unknowns, too).
Well, the [X_i,Y_i,Z_i] are consistent with the projection data iff
for each i it's true that [X_i,Y_i,Z_i], [p,q,r], and [r_i,s_i,0]
(the assumed image) are collinear; we will simply ignore the physical
constraint of "betweenness". Well, any three points A B C are collinear
iff the differences A-B and B-C are dependent vectors so we find
we need [X_i-p, Y_i-q, Z_i-r] * r = [p-r_i, q-s_i,r]*(Z_i-r), i.e.
r X_i - p Z_i = r_i (r- Z_i)
r Y_i - q Z_i = s_i (r- Z_i)
r Z_i - r Z_i = 0
or, r V_i = Z_i [p,q,r] + (r-Z_i) [r_i,s_i,0] - r [X0, Y0, Z0] .
Note also that if we assume for simplicity that r0=s0=0 then
Z0 [p,q,r] = r [X0, Y0, Z0]. In this case the equations which state that
the points XYZ project to the images rs0 may be written
(**) r V_i = (Z_i-Z_0) [p,q,r] + (r-Z_i) [r_i,s_i,0] .
If we combine (*) and (**) we obtain equations which state that the
five points XYZ form a polyhedron which is congruent to the one formed by
the xyz, and which projects to the given points in the plane. We may write
these conditions as 10 equations of the form
< r V_i, r V_j > = r^2 IP_{i,j}, (1 <= i <= j <= 4)
where r V_i is an abbreviation for the vector
r V_i = (Z_i-Z_0) [p,q,r] + (r-Z_i) [r_i,s_i,0] .
which means that the only unknowns now are Z_0, Z_1, ..., Z_4, p, q, r
Now, as I say, this will really be just 9 independent equations if the
numbers IP_{i,j} correspond to a real polyhedron. But that's still a
set of 9 equations in 8 unknowns, and so we expect they cannot
always be solved. More precisely, if we carry along the expressions
IP_{i,j} and r_i, s_j as formal variables, we should be able to use
the 9 equations to eliminate the 8 unknowns, leaving one (polynomial)
equation in these parameters, which must be satisfied if the original
problem is to be solvable.
================================
I'm sorry to say it exhausts my computing resources to describe that
polynomial explicitly, and I am using what are probably close to optimal
tools short of something approaching a supercomputer. But I will illustrate
what I have been doing with a numerical example.
Suppose you begin with five points with coordinates
[0,0,0], [2,0,0], [1,1,0], [1,2,1], [1,-1,3]
Use the identity Euclidean map, that is, don't move them at all. Then
project them down to the xy-plane from the point [p,q,r]=[0,0,4].
You will get images respectively at
[0,0], [2,0], [1,1], [4/3,8/3], [4,-4].
Now hold all these data fixed except suppose we replace the coordinate
of the second image with variables [r1, s1]. Which combinations of
these variables (besides the current [2,0] ) will lead to a solvable
problem? How can we find the solutions in those cases?
Well, if the problem is solvable, there will be five points XYZ in
space forming four vectors relative to [X0, Y0, Z0] whose inner products
agree with those of the first five points: the values are
IP11:=4;IP22:=2;IP33:=6;IP44:=11;IP12:=2;IP13:=2;IP14:=2;IP23:=3;IP24:=0;IP34:=2;
(You can check that det(IP)=0 but that the upper 3x3 block of IP is
positive definite). I use these 10 values to describe 10 equations "(**)"
in the 10 unknowns Z0 Z1 Z2 Z3 Z4 p q r r1 s1. Here they are, in case you
want to feed them into a CAS:
r^2*(-4)+(p*(Z1-Z0)+r1*(r-Z1))^2+(q*(Z1-Z0)+s1*(r-Z1))^2+r^2*(Z1-Z0)^2,
(r-Z2+p*(Z2-Z0))*(p*(Z1-Z0)+r1*(r-Z1))+(r-Z2+q*(Z2-Z0))*(q*(Z1-Z0)+s1*(r-Z1)) -2*r^2+r^2*(Z1-Z0)*(Z2-Z0),
r^2*(-2)+r^2*(Z1-Z0)*(Z3-Z0)+(p*(Z1-Z0)+r1*(r-Z1))*(r*4/3+Z3*(-4/3)+p*(Z3-Z0))+(q*(Z1-Z0)+s1*(r-Z1))*(r*8/3+Z3*(-8/3)+q*(Z3-Z0)),
r^2*(-2)+r^2*(Z1-Z0)*(Z4-Z0)+(p*(Z1-Z0)+r1*(r-Z1))*(r*4+Z4*(-4)+p*(Z4-Z0))+(q*(Z1-Z0)+s1*(r-Z1))*(r*(-4)+Z4*4+q*(Z4-Z0)),
r^2*(-2)+(r-Z2+p*(Z2-Z0))^2+(r-Z2+q*(Z2-Z0))^2+r^2*(Z2-Z0)^2,
r^2*(-3)+r^2*(Z2-Z0)*(Z3-Z0)+(r-Z2+p*(Z2-Z0))*(r*4/3+Z3*(-4/3)+p*(Z3-Z0))+(r-Z2+q*(Z2-Z0))*(r*8/3+Z3*(-8/3)+q*(Z3-Z0)),
r^2*(Z2-Z0)*(Z4-Z0)+(r-Z2+p*(Z2-Z0))*(r*4+Z4*(-4)+p*(Z4-Z0))+(r-Z2+q*(Z2-Z0))*(r*(-4)+Z4*4+q*(Z4-Z0)),
r^2*(-6)+r^2*(Z3-Z0)^2+(r*4/3+Z3*(-4/3)+p*(Z3-Z0))^2+(r*8/3+Z3*(-8/3)+q*(Z3-Z0))^2,
r^2*(-2)+r^2*(Z3-Z0)*(Z4-Z0)+(r*4+Z4*(-4)+p*(Z4-Z0))*(r*4/3+Z3*(-4/3)+p*(Z3-Z0))+(r*(-4)+Z4*4+q*(Z4-Z0))*(r*8/3+Z3*(-8/3)+q*(Z3-Z0)),
r^2*(-11)+r^2*(Z4-Z0)^2+(r*4+Z4*(-4)+p*(Z4-Z0))^2+(r*(-4)+Z4*4+q*(Z4-Z0))^2
Now, as I said we "simply" eliminate the 8 variables Z0,Z1,Z2,Z3,Z4,p,q,r to
get a relation between r1 and s1. This is something of a challenge in
general but I had luck with the Magma CAS by eliminating the Z's first,
then p and q; one of the remaining relations was then of the form r^4*(...)
where (...) involved only r1, s1; this is the relation we want.
(This took about 3 hours on my 450 MHz Linux box, with 128Meg RAM).
So here's the answer: with the data as above, there is a solution iff the
last image [r1, s1] lies on the curve described by this equation:
r1*s1*(-2687257018368) + r1^2*1963764744192 + r1^3*(-1947469873152) + s1^2*723492274176 + r1^4*(-2459327528960) + s1^3*(-1322212392960) + r1^5*3791991717888 + s1^4*(-87147413504) + r1^6*(-257282429952) + s1^5*2891539955712 + r1^7*(-2013847068160) + s1^6*(-1195480671232) + r1^8*1392725716800 + s1^7*316476828672 + r1^9*(-387303480000) + s1^8*(-40732774336) + r1^10*40752250000 + s1^9*1380391440 + s1^10*(-13627193) + r1*s1^2*6529587609600 + r1^2*s1*(-6050518401024) + r1*s1^3*(-10784675725312) + r1^3*s1*8161537556480 + r1*s1^4*2670617149440 + r1^4*s1*18872985796608 + r1*s1^5*(-4367018932224) + r1^5*s1*(-40417541036032) + r1*s1^6*1235964359168 + r1^6*s1*31106916575232 + r1*s1^7*(-201147114048) + r1^7*s1*(-12120526209600) + r1*s1^8*11109277504 + r1^8*s1*2402261388800 + r1*s1^9*(-103010360) + r1^9*s1*(-192900552500) + r1^2*s1^2*26374186991616 + r1^2*s1^3*(-11409178656768) + r1^3*s1^2*(-83057697914880) + r1^2*s1^4*14562607744000 + r1^3*s1^3*60163334238208 + r1^4*s1^2*82270699262976 + r1^2*s1^5*(-3971959486464) + ...
Whew! This seems like a mess already! But it's what I predicted: one
equation in the variables representing the "given" data which must be
satisfied in order for the problem to be solvable.
Now, one may check that if r1=2, s1=0 then the equation above is indeed
satisfied, but of course we knew the problem could be solved with those
data. Let's try another, but nearby solution. If r1=1.9, say, we compute
the polynomial has one solution near s1=0, namely s1=0.02323370156.
So let me show how to find the projection point [p,q,r] and the
Euclidean motion so that the images of the original points are
[0,0], [1.9, 0.02323370156], [1,1], [4/3,8/3], [4,-4].
With the current values of [r1,s1], my set of 10 equations should be
consistent; I should be able to use them to solve for the remaining
variables. Now, how to do that is a matter of preference: I tried,
for example, to eliminate all variables but r and r1. The resulting
mess has hundreds of lines! But it's just a polynomial of degree 20 in
r, once the current value r1=1.9 is substituted. We expect a root
close to the previous value r=4, and indeed r=3.7939490 is a root.
Similarly we eliminate all variables but r1, s1, r, q and get a
quartic in q; the root closest to the previous value q=0 is about
q=-0.0437014. Finally we eliminate only the Z's and find a quadratic
in p with a root near 0.9070586.
So it looks like we should be able to project some set of 5 points from
[ 0.9070586, -0.0437014, 3.7939490] instead of [0,0,4] to
get the five shadows we seek. Now, where exactly are those five points?
We may return to the 10 equations, which are now simple quadratics in
Z0 through Z4. It's easy to solve these e.g. by Newton's method.
(In fact, I can use N.M. to improve the accuracy of the previous values:
p = 0.907058662040766975040745284008,
q = -0.0437014142167977395067955677103,
r = 3.79394906249537291231162065101,
Z0 = 0.152327220168302272678120874843,
Z1 = -0.311754736981796463400985375212,
Z2 = -0.0469752384346046997031360265532,
Z3 = 0.957918193869912383126156797325,
Z4 = 2.80401245861823457049326132692,
s1 = 0.0232337015668761989194584211035,
r1 = 1.9
Looking back near the beginning of this letter, we find equations to
locate the points. First, from Z0 [p,q,r] = r [X0, Y0, Z0] we find one
of the points to be [0.0364184448, -0.001754613685, 0.1523272201683].
Working backwards with (**), we may reconstruct the vectors V_i = XYZ_i-XYZ_0:
(**) r V_i = (Z_i-Z_0) [p,q,r] + (r-Z_i) [r_i,s_i,0] .
We add these to [X0, Y0, Z0] to get the other [Xi, Yi, Zi]:
[X1,Y1,Z1]=[1.9815915449982959, 0.02873386487224762, -0.31175473698179646]
[X2,Y2,Z2]=[1.0011507643985590, 1.01292271508651231, -0.04697523843460434]
[X3,Y3,Z3]=[1.2257050398347134, 1.98233726371603986, 0.95791819386991238]
[X4,Y4,Z4]=[1.7140847432177094,-1.07599908648922958, 2.80401245861823457]
Finally, we may compute the Euclidean motion carrying the initial five
points to these. If the motion is a rotation followed by a translation,
then since [x0, y0, z0] was at the origin, the translation will have to
be simply [X0,Y0,Z0] again, and the rotation is simply the one which
carries the other [xi,yi,zi] to V_i. It is sufficient to solve the
linear system A [ xi, yi, zi] = V_i , i=1, 2, 3.
I get roughly A =
[ 0.97258655 -0.0078542305 0.2324085059 ]
[ 0.01524424 0.9994330895 -0.0300185409 ]
[-0.23204098 0.0327385199 0.9721549123 ]
So if you rotate the original five points by this matrix, translate
by [X0,Y0,Z0], and then project from [p,q,r], you should get five
shadows where I indicated -- four of them where the first "trivial" set of
shadows was, and the other moved slightly away to have r1=1.9 instead of
r1=2. As I promised, I could compute everything as long as the shadow
point (r1,s1) stayed on that complicated curve.
I hope this is of some use for your project.
dave
==============================================================================
From: Stephen Robinson
Subject: Re: A Problem in Projective Geometry
Date: Wed, 22 Sep 1999 14:43:02 -0400
Newsgroups: [missing]
To: Dave Rusin
Dear Dave,
Thanks for your response. I am in the process of writing up my results
and I'll send them when I am finished. It is a theorem that 5 points in
general position lead to a solution that can be solved for precisely.
Also, the problem can still be handled if 4 of the points are in a
plane, although then there are singular cases to be avoided.
One of my reasons for posting the problem is that I suspected that a
Projective Geometer somewhere had already solved similar problems, and
had probably done so more elegantly.
I will look over your note below this week to compare results.
Best Regards,
Steve
Dave Rusin wrote:
[letter above was quoted --djr]