From: rusin@vesuvius.math.niu.edu (Dave Rusin) Subject: Re: "center" of 3 skewed lined Date: 5 Aug 2000 22:02:09 GMT Newsgroups: sci.math,geometry.forum,geometry.research,alt.cad,comp.graphics.algorithms Summary: An answer exists. It may not be the answer you want. In article <398A3DF6.15E486A@lems.brown.edu>, Frederic Leymarie wrote: >I have a problem in geometry (or is it algebra?!) It's both, but I'll treat it as algebra since someone else has commented on the geometry. > Given 3 (mutually) skew lines* in Euclidean 3D space, > find their "center", i.e., the unique point along the > 1D loci of points at equal minimal distance from the > 3 lines . > >I can derive an iterative algorithm to find that point, but I wonder >if there is a geometric construct for this locus, or maybe >an algebraic formula. No geometric construct in the strict sense, nor an algebraic formula if you mean using radicals and so on, but we can describe the coordinates as roots of equations. I will give a procedure and an example. We can do most of this as an inner-product-space calculation, that is, we need to know what '0', '+', and '.' mean, but we don't need coordinates (yet). First, how far is a point from a line? Let u1 be a unit vector in that line (i.e. the vector of length 1 joining two points in the line). Let pi be the plane u1^\perp, and let u0 be the point where pi meets the line (that is, u0 is really the vector joining the origin to that point.) Now given any other point P in space, we compute its distance to the line by projecting down to pi. We want P - x u1 to lie in pi so we need (P - x u1) . u1 = 0 or simply x = (P . u1). Thus P - (P.u1) u1 lies in pi and its distance to the line is its distance to u0, which is the square root of || (P - (P.u1) u1) - u0 ||^2 = (P - (P.u1)u1 - u0).(P - (P.u1)u1 - u0) = (P.P) - (P.u1)^2 + (u0.u0) - 2 (P.u0) (Exercise: find the Pythagorean Theorem hidden in this formula.) Now given another line we may likewise compute vectors v0, v1 (unique up to the sign of v1) with v1.v1=1, v1.v0=0, and the square of the distance from P to that line given by (P.P) - (P.v1)^2 + (v0.v0) - 2 (P.v0). In particular, P is equidistant between the lines iff (P.u1)^2 - (P.v1)^2 + 2 (P.u0) - 2 (P.v0) = (u0.u0) - (v0.v0) which is (assuming P, u1, and v1 are in 'general position') a quadratic equation in P, meaning that the set of such points is one of the quadric surfaces. Finally, given a third line, described by a third pair of vectors w0, w1, we have two equations to describe the locus of points equidistant between the lines: (P.u1)^2 - (P.v1)^2 + 2 (P.u0) - 2 (P.v0) = (u0.u0) - (v0.v0) (P.u1)^2 - (P.w1)^2 + 2 (P.u0) - 2 (P.w0) = (u0.u0) - (w0.w0) As an intersection of two quadrics, this locus should be an elliptic curve, and in particular is not rationally parameterizable except in degenerate cases. I ran some computations with the specific values of the parameters below; in that case, anyway, the curve turns out indeed to have genus 1, with two singular points (where y=1). Now we can do some computations. We can choose coordinates to do so in a reasonable way. Use the Gram-Schmidt process to form an orthonormal basis [u1, u2, u3] from the ordered basis [u1, v1, w1], and then express every vector as linear combinations of these three: u1 = u1 v1 = a u1 + b u2 for some a and b with a^2 + b^2 = 1, and w1 = c u1 + d u2 + e u3 where c^2 + d^2 + e^2 = 1. Also (thanks to orthogonality) we will be able to write u0 = f u2 + g u3 v0 = bh u1 - ah u2 + k u3 w0 = l u1 + m u2 + n u3 where c l + d m + e n = 0. Finally, write P = x u1 + y u2 + z u3. We would like to know the relations on x,y, and z which characterize P being in that curve. We can now treat this as a coordinate problem in R^3 if we wish. We have two quadratic equations in x,y,z with coefficients which will be polynomials in the quantities a,b,c,d,e,f,g,h,k,l,m,n (which are essentially arbitrary except for the constraints a^2 + b^2 = 1, c^2 + d^2 + e^2 = 1, and c l + d m + e n = 0. ) The quadratic equations are, say, eq1 and eq2 : x^2 - (a x + b y)^2 + 2 (f y + g z) - 2(b h x - a h y + k z) = f^2+g^2-h^2-k^2 x^2 - (c x + d y + e z)^2 + 2 (f y + g z) - 2(l x + m y + n z) = f^2+g^2-l^2-m^2-n^2 Note that the first of these two is linear in z, so we may easily eliminate z from our set of equations. (The solution set of these equations is then a certain algebraic curve in the x-y plane, but it's not particularly interesting.) Next, recall that the original goal was to choose the point on this geometric locus to minimize distance to the lines. Of course it's sufficient to minimize the square of the distance, which is (P.P) - (P.u1)^2 - 2 (P.u0) + (u0.u0) = y^2+z^2 - 2 (f y + g z) + f^2+g^2 We simply treat this as a Lagrange Multipliers problem. We set the gradient of the 5-variable function F = dist2 - p*eq1 - q*eq2 equal to zero. This gives five equations in five unknowns. It's easy to see two are linear in p,q and another as we noted is linear in z, so we easily reduce to two equations in two unknowns x,y. These last can be solved using elimination or resultants. Well, these equations are pretty ugly. I looked at them with Maple, and you may wish to follow along at home: I set eq1:=x^2-(a*x+b*y)^2+2*(f*y+g*z)-2*(b*h*x-a*h*y+k*z)-(f^2+g^2-h^2-k^2); eq2:=x^2-(c*x+d*y+e*z)^2+2*(f*y+g*z)-2*(l*x+m*y+n*z)-(f^2+g^2-l^2-m^2-n^2); dd:=y^2+z^2-2*(f*y+g*z)+f^2+g^2; F:= dd - p*eq1 - q*eq2: vars:=[x,y,z,p,q]: grads:=[seq(diff(F,v),v=vars)]: to get the five equations. Now reduce a bit: readlib(eliminate): solve({grads[1],grads[2]},{p,q}):numer(factor(subs(%,grads))): xyz:=[%[3]/2, %[4],%[5]]; myz:=solve(xyz[2],z);numer(factor(subs(z=%,xyz))):xy:=[%[1],%[3]]: At this point, length(xy); shows the answer to be 32K long! One can proceed further but the answer is so long as to be almost useless [*]. Instead, I found it helpful to try a test configuration: test:={a=7/25,b=24/25,c=4/5,d=3/13,e=36/65,f=1/4,g=5/7,h=2/3,k=-2, l=-3,m=-4,n=6}; xy0:=numer(factor(subs(test,xy))); myx:=eliminate({op(xy0)},x); A component of this last calculation shows the condition y must satisfy: 1008862006740136425979835953961899574014479125446656*y^12 +37363815126157608565211009320552646255906028836093952*y^11 +773536246701798013730116055554696879439611556807049216*y^10 +20049285816202817845316456235484146912035641817655410688*y^9 +276071515918302450197467384192708489550322364194333130752*y^8 +2870384151743365199225264094665258034042792724641651621888*y^7 +36375406621070048498826133497000515055637340125370567884800*y^6 +179013277572518950236754011134614971709679247839555198386176*y^5 -1048948766504179978078497617619824784229769544216796296069120*y^4 -15755781346593974832290961864310372806700962712160733978580992*y^3 -63925751753070871588766372541615867533094983303946818894289216*y^2 -103282269729291020590568338936076562094018880020498083867757200*y -58447101458117825021105998736019932504464176567473747313825625 must vanish. Assuming these computations to representative of the general case, we see that y is given as a root of a polynomial of degree 12, generically irreducible over Q. For these particular data, we find there are four real roots to the equation for y, leading to four critical points on the curve of points equidistant from the three lines. (I'd bet you could find other data which would lead to polynomials with anywhere from 2 to 12 real roots, but I didn't bother to go experimenting.) It is necessary to do the computations with rather great precision, but this is what I believe are there coordinates, and the squares of the distances of these critical points to the lines: x y z dist^2 - 2.80046939 - 3.74799272 1.32226182 16.35358079 3.12534520 6.91858531 7.52003607 90.78826797 -30.84546726 -25.06149229 17.70321989 929.29552650 26.02925268 -12.98343296 -112.35873822 12960.63249016 So we have found the "center" of the set of three skew lines to be near (-2.800, -3.748, 1.322), but we observe that this point is defined algebraically to be in an extension field of degree 12 over the field of definition of the original lines, and in particular (12 not being a power of 2) we see the center cannot be described with the usual Euclidean geometric constructions. Of course I'm being cavalier here about genericity: I have been assuming various vectors to be in 'general position', I have been ignoring singular points on the curve of equidistant points (where the method of Lagrange Multipliers can fail), and I have simply assumed the existence of a "center" where a minimum exists. I imagine that these situations can be dealt with more carefully should there be need to do so. dave [*] Actually Maple failed for some reason. I did manage to get it to treat l,m,n as variables, with the other 9 parameters fixed as in test. The computations took several minutes, and the answer has something like 1386 terms. I asked it to treat h and k as variables too, with the other seven parameters still fixed. After about an hour it quit, having found nothing. Of course treating all nine parameters as variables seems hopeless with Maple. (This is Maple V release 5.1) I also tried to get Magma to do the elimination. It would almost surely eventually succeed. Here is typical input to Magma V2.5-1: Q:=RationalField(); R:=PolynomialRing(Q,8); a:=7/25;b:=24/25;c:=4/5;d:=3/13;e:=36/65;f:=1/4;g:=5/7;h:=2/3;k:=-2; I:=ideal; J:=EliminationIdeal(I,{y,l,m,n}); This runs for a long time and uses something like 40M of RAM. I use a Linux box with 128M RAM and another 128M swap space, and it probably cannot even handle a problem with five or so of the parameters treated as variables. Computing the degree-12(?) polynomial with explicit coefficients as polynomials in a,b,c,d,e,f,g,h,k,l,m,n would require far more resources, I imagine. Even if one managed to do so: would you really want to see the answer?