Cylinders Through Five Points: Computational Algebra and Geometry

Daniel Lichtblau
Wolfram Research, Inc.
100 Trade Centre Dr.
Champaign, IL 61820
danl@wolfram.com

ICMS, Castro Urdiales, Spain
September 1-3, 2006

Abstract

We address the following question: Given five points in (ℝ^3), determine a right circular cylinder containing those points. We obtain algebraic equations for the axial line and radius parameters and show that these give six solutions in the generic case. An even number (0, 2, 4, or 6) will be real valued and hence correspond to actual cylinders in (ℝ^3). We will investigate computational and theoretical matters related to this problem. In particular we will show how exact and numeric Gröbner bases, equation solving, and related symbolic-numeric methods may be used to advantage. We will also discuss some applications.

Outline of the Problem

Given five points in (ℝ^3), we are to determine all right circular cylinders containing those points.

Questions of importance:

How do we know there are finitely many in generic case? How many are there? (Depends...Are we working in real or complex space?)

How do we find the axial line and radius parameters?

Given the cylinder parameters, how do we obtain its implicit equation?

Reversing this, how can one obtain parameters from the implicit form?

How might we display them graphically?

Given six or more points, how do we find the coordinates of a cylinder in (ℝ^3) that "best" fits those points?

Given five points chosen with random uniform distribution in a cube, what is the expected probability that one lies inside the convex hull of the other four. Related to "no real cylinder" case . Also to an old recently solved problem in integral geometry.

To what extent can computational methods be used to prove enumerative geometry or other types of results related to this problem?

Basics

Terminology

We say "cylinder" for any solution, and "real cylinder" for the real valued solutions.

We call configurations "generic" if they do not have multiple solutions and if all sufficiently small perturbations give same number of solutions. Usually we assume this of our configurations.

Easy to show

There are finitely many solutions (plausible, because we get one equation for each data point and require five parameters to describe a cylinder).

If points are real valued then complex solutions pair off.

One expects six solutions. Reason: Take five "random" points. Solve for cylinder parameters. You "always" get six solutions. This is the Shape Lemma at work for you.

Related Work

Already known

Number of solutions in generic case is indeed six. First shown 1977 (Bottema and Veldkamp). Various other proofs appearing in recent years. We will show two simple computational proofs at end of this report.

Related to this presentation

This is a companion to "Cylinders Through Five Points: Complex and Real Enumerative Geometry", which is really part 2 but was presented already at ADG 2006. The focus here is more on computational methods and related problems.

Computing cylinders through five points

Our set up

Parametrize the axis line as {y=a x+b,z=c x+d} . So we need to solve for a, b, c, d, and a radius r. One might argue that this only captures "generic" cases. But it avoids issues with double counting if we allow most general form of direction vector (because its negative gives same cylinder).

Project points onto axis. For (Subscript[j, th]) point ((Subscript[x, j]) ,(Subscript[y, j]),(Subscript[z, j])) we need length of orthogonal projection (Subscript[perp, j]). 
Work with equation  ((‖(Subscript[perp, j])‖)^2)=(r^2).
After some minor algebra, we have five polynomials of form

(b^2)+(b^2) (c^2)-2 a b c d+(d^2)+(a^2) (d^2)-(r^2)-(a^2) (r^2)-(c^2) (r^2)+2 a b (Subscript[x, j])+2 c d (Subscript[x, j])+(a^2) (Subscript[x, j]^2)+(c^2) (Subscript[x, j]^2)-2 b (Subscript[y, j])-2 b (c^2) (Subscript[y, j])+2 a c d (Subscript[y, j])-2 a (Subscript[x, j]) (Subscript[y, j])+(Subscript[y, j]^2)+(c^2) (Subscript[y, j]^2)+2 a b c (Subscript[z, j])-2 d (Subscript[z, j])-2 (a^2) d (Subscript[z, j])-2 c (Subscript[x, j]) (Subscript[z, j])-2 a c (Subscript[y, j]) (Subscript[z, j])+(Subscript[z, j]^2)+(a^2) (Subscript[z, j]^2)=0

Computing cylinders through five points

An example

Points: (7,9,8), (8,-4,-10), (-4,1,4), (-9,-9,-10), and (-7,-10,-10).

This has two real valued solutions:

How to find these solutions? Find all six using a (perhaps numeric) polynomial solver. Can be done e.g. by homotopy continuation or reduction to an eigensystem. In Mathematica NSolve uses the latter approach. The actual code for finding cylinder parameters is quite simple. As I need it later I show it below.

perp[vec1_,vec_,offset_]:=vec1-offset-Projection[vec1-offset,vec,Dot]

solveCylinders[pts_List,vec_,offset_,prec_:Automatic]:=Module[{exprs,k,perps},perps=Table[perp[pts[[k]],vec,offset],{k,5}];exprs=(Numerator[Together[#1.#1-rsqr]]&)/@perps;NSolve[exprs,{a,b,c,d,rsqr},WorkingPrecision->prec]]

The implicit equation

The orthodox way

We start with a parametrization: Find two unit vectors pairwise orthogonal and orthogonal to the direction of the axis. Say P is on the axis, v is a direction vector, (Subscript[w, 1]) and (Subscript[w, 2]) are the perps. A point on the cylinder is parametrized by (t,θ) as P+t v + (Subscript[w, 1]) cos(θ)+(Subscript[w, 2]) sin(θ) . To make it algebraic we consider the sine and cosine terms as algebraic variables with the usual trig identity linking them. Use a Gröbner basis computation to eliminate these parameters, obtaining an implicit relation satisfied by the cylinder parameters (a,b,c,d,r).

Result:

(b^2)+(b^2) (c^2)-2 a b c d+(d^2)+(a^2) (d^2)-(r^2)-(a^2) (r^2)-(c^2) (r^2)+(2 a b+2 c d) x+((a^2)+(c^2)) (x^2)+(-2 b-2 b (c^2)+2 a c d) y-2 a x y+(1+(c^2)) (y^2)+(2 a b c-2 d-2 (a^2) d) z-2 c x z-2 a c y z+(1+(a^2)) (z^2)

The smart way

Step 1. Use the formulation we described for finding the distance from a point to the axial line:  ((‖(Subscript[perp, j])‖)^2)-(r^2) gives exactly the polynomial we seek (that is, the same as the one above).

Step 2. Feel foolish.
Applies only to those of us, like myself, who did it the hard way first.

An example

We take cylinder with parameters:

a=3,b=2,c=4,d=-1,r=(Sqrt[21])

Implicit polynomial defining the cylinder:

Parameters from implicit form

The easiest way

Step 1. Take general implicit equation and specific one for cylinder at hand. Equate coefficients.

Step 2. This gives equations in the parameters. Solve them.

Our example

In Mathematica one might use SolveAlways to automatically equate coefficients and solve. But pretty much any symbolic math program can do this in some way.

SolveAlways[(b^2)+(b^2) (c^2)-2 a b c d+(d^2)+(a^2) (d^2)-(r^2)-(a^2) (r^2)-(c^2) (r^2)+(2 a b+2 c d) x+((a^2)+(c^2)) (x^2)+(-2 b-2 b (c^2)+2 a c d) y-2 a x y+(1+(c^2)) (y^2)+(2 a b c-2 d-2 (a^2) d) z-2 c x z-2 a c y z+(1+(a^2)) (z^2)==-420+4 x+25 (x^2)-92 y-6 x y+17 (y^2)+68 z-8 x z-24 y z+10 (z^2)/.(r^2)->rsqr,{x,y,z}]

Best fit to overdetermined cylinders

General idea

Step 1. Pick five points,

Step 2. Solve for cylinder parameters, obtain candidate solutions. Discard complex ones.

Step 3. Form sum of squares of distances of all points to the remaining candidates. Take the one with the smallest sum of squares.

Step 4. Do a (nonlinear) least squares minimization, using the candidate's values as start points.

Refinements

Might try several sets of five, use the most promising candidate.

Might make effort to choose five points not too "close" to one another, in attempt to reduce ill conditioning.

Applications

Geometric tolerancing (metrology)

Ftting object to point cloud in scene reconstruction

First step to fitting peptides and other biomacromolecules to a helix

Graphing cylinders through five points

Goals

Good view of the cylinder

See how it hits the points

An example with considerable symmetry

We work with two regular tetrahedra glued along a face in the horizontal plane.

dpoints={{1,0,0},{-1/2,(Sqrt[3])/2,0},{-1/2,-(Sqrt[3])/2,0},{0,0,(Sqrt[2])},{0,0,-(Sqrt[2])}};

vec={a,c,1};offset={b,d,0};

First we solve for the parameter values. This configuration gives six real cylinders, all with radius 9/10.

solns=solveCylinders[dpoints,vec,offset,Infinity];
FullSimplify[{a,b,c,d,rsqr}/.solns]

Graphing cylinders ...

Plotting

We try to give an idea of how one solution cylinder "fits" the set of points. The radius is shrunk slightly to give a better view.

{p1,p2,p3,p4,p5}=dpoints;
top1=(vec*2+offset)/.solns[[1]];
bot1=(-vec*2+offset)/.solns[[1]];
Graphics3D[{{Red,Sphere[#,0.04]&/@{p1,p2,p3},Green,Sphere[p4,0.04],Blue,Sphere[p5,0.04]},{White,Opacity[0.1],Polygon[{{p1,p2,p4},{p2,p3,p4},{p3,p1,p4}}],Polygon[{{p1,p2,p5},{p2,p3,p5},{p3,p1,p5}}]},{Yellow,Opacity[0.6],Cylinder[{top1,bot1},0.70]}},Axes->True,AxesLabel->{"x","y",None}]

Graphing cylinders ...

Here we shrink the radius considerably and change the orientation in order to better see how the axis cuts through the double tetrahedra.

Graphics3D[{{Red,Sphere[#,0.04]&/@{p1,p2,p3},Green,Sphere[p4,0.04],Blue,Sphere[p5,0.04]},{White,Opacity[0.1],Polygon[{{p1,p2,p4},{p2,p3,p4},{p3,p1,p4}}],Polygon[{{p1,p2,p5},{p2,p3,p5},{p3,p1,p5}}]},{Yellow,Opacity[0.6],Cylinder[{top1,bot1},0.05]}},Axes->True,AxesLabel->{"x","y",None}]

Graphing cylinders ...

Another sort of plot as a parametrized surface:

Significance of the double tetrahedra configuration example

Related to problem of finding cylinders of given fixed radius through four given points. There are 12 solutions to the equations that result from this problem. All solutions can be real. In the example above all cylinders had the same radius of 9/10. Consider the top tetrahedron to give the four points. Gluing another onto each face gives six cylinders that go through those four points and all have radius of 9/10. Seems like 24 cylinders, but they pair off for a total of 12.

Numbers of real cylinders that arise

Basics

Already saw number is even.

Already saw it can be as large as 6.

Using pseudorandom sets of five points in unit cube will show the other possibilities all arise.

Obvious on reflection: We get no real cylinders whenever one point is in the (tetrahedral) hull of the other four. Reason: All projections onto planes keep it inside the hull of the projected quadrilateral, hence planar quadratic through the five projected points cannot be elliptical, hence not a circle.

Numbers of real cylinders...

Interesting to study

How often do such random examples have one point enclosed by the other four? Related to expected volume of a tetrahedron with points picked randomly in unit cube (generalization of a problem of Sylvester). Only recently solved, using symbolic calculus among other things. It is 3977/216000-(π^2)/2160.

How do other cases of no real cylinders arise. This is discussed in the companion to this talk.

More generally, can we (either algebraically or geometrically) classify the cases of 0, 2, 4, or 6 real solutions? Little seems to be known about this.
Suggested by referee: look at discriminant varieties.
Me: Good idea. What's a discriminat variety?
Okay, that wasn't quite my response. But I will say the computations appear to be daunting. One problem is that it takes extra polynomials and variables just to enforce the condition that no pair of points coincides.

Proving the generic number of cylinders is six

To simplify the computations we reformulate so as to have two equations in two variables. Idea:

Without loss of generality we have one point at the origin, another at (1,0,0), and a third in the x y coordinate plane.

Project these onto the set of planes through the origin, parametrized generically by a normal vector (a,b,1).

In each such plane these three points determine a circle, and we get one equation for each of the remaining two points in order that they project onto the same circle (which is the condition that the five be cocylindrical).

So our points are (0,0,0), (1,0,0), ((Subscript[x, 2]),(Subscript[y, 2]),0), ((Subscript[x, 3]),(Subscript[y, 3]),(Subscript[z, 3])), and ((Subscript[x, 4]),(Subscript[y, 4]),(Subscript[z, 4])) and our variables are (a,b) where direction vector is (a,b,1). After some linear algebra we obtain the polynomials below.

(-(Subscript[x, 3]) (Subscript[y, 2])-(b^2) (Subscript[x, 3]) (Subscript[y, 2])+(Subscript[x, 3]^2) (Subscript[y, 2])+(b^2) (Subscript[x, 3]^2) (Subscript[y, 2])+(Subscript[x, 2]) (Subscript[y, 3])+(b^2) (Subscript[x, 2]) (Subscript[y, 3])-(Subscript[x, 2]^2) (Subscript[y, 3])-(b^2) (Subscript[x, 2]^2) (Subscript[y, 3])+2 a b (Subscript[x, 2]) (Subscript[y, 2]) (Subscript[y, 3])-2 a b (Subscript[x, 3]) (Subscript[y, 2]) (Subscript[y, 3])-(Subscript[y, 2]^2) (Subscript[y, 3])-(a^2) (Subscript[y, 2]^2) (Subscript[y, 3])+(Subscript[y, 2]) (Subscript[y, 3]^2)+(a^2) (Subscript[y, 2]) (Subscript[y, 3]^2)-b (Subscript[x, 2]) (Subscript[z, 3])-(b^3) (Subscript[x, 2]) (Subscript[z, 3])+b (Subscript[x, 2]^2) (Subscript[z, 3])+(b^3) (Subscript[x, 2]^2) (Subscript[z, 3])+a (Subscript[y, 2]) (Subscript[z, 3])+a (b^2) (Subscript[y, 2]) (Subscript[z, 3])-2 a (b^2) (Subscript[x, 2]) (Subscript[y, 2]) (Subscript[z, 3])-2 a (Subscript[x, 3]) (Subscript[y, 2]) (Subscript[z, 3])+b (Subscript[y, 2]^2) (Subscript[z, 3])+(a^2) b (Subscript[y, 2]^2) (Subscript[z, 3])-2 b (Subscript[y, 2]) (Subscript[y, 3]) (Subscript[z, 3])+(a^2) (Subscript[y, 2]) (Subscript[z, 3]^2)+(b^2) (Subscript[y, 2]) (Subscript[z, 3]^2),-(Subscript[x, 4]) (Subscript[y, 2])-(b^2) (Subscript[x, 4]) (Subscript[y, 2])+(Subscript[x, 4]^2) (Subscript[y, 2])+(b^2) (Subscript[x, 4]^2) (Subscript[y, 2])+(Subscript[x, 2]) (Subscript[y, 4])+(b^2) (Subscript[x, 2]) (Subscript[y, 4])-(Subscript[x, 2]^2) (Subscript[y, 4])-(b^2) (Subscript[x, 2]^2) (Subscript[y, 4])+2 a b (Subscript[x, 2]) (Subscript[y, 2]) (Subscript[y, 4])-2 a b (Subscript[x, 4]) (Subscript[y, 2]) (Subscript[y, 4])-(Subscript[y, 2]^2) (Subscript[y, 4])-(a^2) (Subscript[y, 2]^2) (Subscript[y, 4])+(Subscript[y, 2]) (Subscript[y, 4]^2)+(a^2) (Subscript[y, 2]) (Subscript[y, 4]^2)-b (Subscript[x, 2]) (Subscript[z, 4])-(b^3) (Subscript[x, 2]) (Subscript[z, 4])+b (Subscript[x, 2]^2) (Subscript[z, 4])+(b^3) (Subscript[x, 2]^2) (Subscript[z, 4])+a (Subscript[y, 2]) (Subscript[z, 4])+a (b^2) (Subscript[y, 2]) (Subscript[z, 4])-2 a (b^2) (Subscript[x, 2]) (Subscript[y, 2]) (Subscript[z, 4])-2 a (Subscript[x, 4]) (Subscript[y, 2]) (Subscript[z, 4])+b (Subscript[y, 2]^2) (Subscript[z, 4])+(a^2) b (Subscript[y, 2]^2) (Subscript[z, 4])-2 b (Subscript[y, 2]) (Subscript[y, 4]) (Subscript[z, 4])+(a^2) (Subscript[y, 2]) (Subscript[z, 4]^2)+(b^2) (Subscript[y, 2]) (Subscript[z, 4]^2))

The generic number...

Now we must count solutions to this system.

PROOF 1. We form a Gröbner basis with respect to a degree based term ordering for the polynomials. Looking at the head terms we find that there are 6 monomials in (a,b) that are not reducible with respect to this basis and hence 6 solutions to the system.

PROOF 2. We compute the resultant of the pair of polynomials with respect to one of the two variables. We obtain a polynomial of degree 6 in the other (with large symbolic coefficients). This means there are at most 6 solutions. As we already know there are at least that many, this suffices to show that there are generically six solutions.

The generic number...

Remarks about these proofs.

Both rely on having a fairly simple formulation of the problem.

With more than two equations, a "standard" resultant would not likely suffice unless used in iterated fashion.

For even slightly more complicated problems, a Gröbner basis computation involving symbolic coefficients is likely to bog down.

A related Gröbner basis approach is to use the Bezout bound and subtract off the count of solutions at infinity. Works for this problem but often fails because of presence of dimensional component(s) at infinity.

Another Gröbner basis approach that seems promising is to work in a mixed algebra with symbolic perturbation variables. Idea is to show we have six solutions not just at a given configuration, but at all configurations in a neighborhood thereof. This gives the generic count.

Summary and open questions

What we did thus far

Discussed computational methods for finding cylinders through a given set of five points in (ℝ^3).

Covered several related problems and computational approaches thereto.

Combined geometric reasoning with Gröbner bases and other tools to study problems from enumerative and computational geometry.

Summary and open questions

Open questions

Describe the configuration space variety for which we have multiplicity.

Related: Describe the configuration space variety on which the number of real solutions changes.

Easy to show: if five points are coplanar then generic number of (complex) solutions is 4. Of these either 0 or 2 are real (depends on whether the quadratic containing the five is a hyperbola or an ellipse). (Picture below for case of 2 real solutions).
Is coplanarity a necessary condition that we have fewer than six solutions?

Easy to show: If points are collinear or lie on two parallel lines, then there are infinitely many solutions.
Are these necessary conditions or are there other configurations wherein we have an infinite solution set?