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]