Identify all the right circular cones passing through six arbitrary points

geometrylinear algebraquadricssystem identification

I have this interesting question. Given $6$ arbitrary points, I want to identify all the possible circular cones passing through them.

The equation of a right circular cone whose vertex is at $\mathbf{r_0}$ and whose axis is along the unit vector $\mathbf{a}$, and whose semi-vertical angle is $\theta$ is given by

$$ (\mathbf{r} – \mathbf{r_0})^T ( (\cos \theta)^2 \ \mathbf{I} – \mathbf{aa}^T ) (\mathbf{r} – \mathbf{r_0}) = 0 $$

where $\mathbf{r}$ is the position vector of a point on the surface of the cone. Counting the number of parameters, we have $3$ for $\mathbf{r_0}$, $2$ for $\mathbf{a}$ and $1$ for $\theta$. We therefore need at least $6$ points on the surface of the cone to specify it.

Question: What is the procedure for extracting the parameters of a right circular cone passing through $6$ arbitrary points?

What I have tried: I have parameterized the axis unit vector as

$ \mathbf{a} = ( \sin t \cos s, \sin t \sin s , \cos t ) $

and written

$\mathbf{r_0} = (x, y, z) $

Now define the functions

$ f_1 = (\mathbf{r_1} – \mathbf{r_0})^T ( (\cos \theta)^2 \ \mathbf{I} – \mathbf{aa}^T ) (\mathbf{r_1} – \mathbf{r_0}) $

$ f_2 = (\mathbf{r_2} – \mathbf{r_0})^T ( (\cos \theta)^2 \ \mathbf{I} – \mathbf{aa}^T ) (\mathbf{r_2} – \mathbf{r_0}) $

$ f_3 = (\mathbf{r_3} – \mathbf{r_0})^T ( (\cos \theta)^2 \ \mathbf{I} – \mathbf{aa}^T ) (\mathbf{r_3} – \mathbf{r_0}) $

$ f_4 = (\mathbf{r_4} – \mathbf{r_0})^T ( (\cos \theta)^2 \ \mathbf{I} – \mathbf{aa}^T ) (\mathbf{r_4} – \mathbf{r_0}) $

$ f_5 = (\mathbf{r_5} – \mathbf{r_0})^T ( (\cos \theta)^2 \ \mathbf{I} – \mathbf{aa}^T ) (\mathbf{r_5} – \mathbf{r_0}) $

$ f_6 = (\mathbf{r_6} – \mathbf{r_0})^T ( (\cos \theta)^2 \ \mathbf{I} – \mathbf{aa}^T ) (\mathbf{r_6} – \mathbf{r_0}) $

Now using the multivariate Newton-Raphson method I could iterate to find a solution for the parameter vector $(t, s, x, y, z , \theta )$ that will make $f_1, f_2, f_3, f_4, f_5, f_6$ all zero.

This works, but it's iterative and at best converges to one of the possible right circular cones.

Is there a way to generate all the possible right circular cones that are solutions, i.e. passing the $6$ given points?

Best Answer

Not really sure whether this answers the question, but here's how I would approach it. Basically, we solve a system of six polynomials in the (free) online version of Mathematica. The lines given below as preformatted text can be copied and pasted into Mathematica.

For a right circular cone with vertex $(a,b,c),$ semi-vertical angle $\theta$,and whose axis points in the direction $(s,t,1)$ the equation is

$$(s (x-a)+t (y-b)+(z-c))^2-\cos^2(\theta)(s^2+t^2+1)((x-a)^2+(y-b)^2+(z-c)^2)=0.$$

(a derivation can be found at http://www.mndcollegerajur.org/uploads/department/cone.pdf, Section 9.4.1 on pg 8.)

Note that the direction vector $(s,t,1)$ isn't fully general .. but it makes the Mathematica code simpler for this exposition. Also, for Mathematica we set $k=\cos^2(\theta)$ so the master expression is

f=(s(x-a)+t(y-b)+1(z-c))^2-(s^2+t^2+1)((x-a)^2+(y-b)^2+(z-c)^2)k

We pick six points $(0,0,0),(-1,0,1),\dots,(0,2,2)$ on the cone $x^2+y^2=z^2$ to create the expressions $f_1,f_2...f_6:$

f1=f /. {x->0,y->0,z->0}
f2=f /. {x->-1,y->0,z->1}
f3=f /. {x->1,y->0,z->1}
f4=f /. {x->-2,y->0,z->2}
f5=f /. {x->2,y->0,z->2}
f6=f /. {x->0,y->2,z->2}

resulting in, e.g. the following expression for $f_1:$

enter image description here

The corresponding system of polynomials can then be solved using:

sol=Solve[f1==0&&f2==0&&f3==0&&f4==0&&f5==0&&f6==0,{s,t,a,b,c,k}]

yielding two solutions:

enter image description here

We recover the two functions $g_1,g_2$ using

g1=f /. sol[[1]] 
g2=f /. sol[[2]]

giving

$$ \begin{align} g_1(x,y,z) &= (-2 y + z)^2 + (-x^2 - y^2 - z^2)/2 \\ g_2(x,y,z) &= z^2 + (-x^2 - y^2 - z^2)/2 \end{align} $$

We recognize $g_2$ as the cone that we started with to create sample points. But $g_1$ is a bit of a surprise. The vertex is at the origin, but the axis is not vertical. Here's a picture (using ContourPlot3D[g1==0,{x,-4,4},{y,-4,4},{z,-4,4}]):

enter image description here

The example used here is a simple one. I've played with this a bit, modifying the points. You can get solutions with complex parameters, leading to cones that are in $\mathbb C^3$ but don't show in $\mathbb R^3$. As well, there can be more than two solutions.

I believe there are plenty of computer algebra systems available, so it shouldn't be difficult to make this work elsewhere.

If OP would like to try out a different set of six points I'd be happy to plug them in.

A note on how the system of equations was solved:

Although OP was ostensibly a question about fitting right circular cones, the comments suggest that it's really about solving a system of polynomial equations. I've used the Mathematica Solve function as a black box to do this, but I'm unaware of a way to coax Mathematica into revealing the steps it took.

My understanding, however, is that computer algebra systems use so-called Gröbner bases to solve these systems. The systems are transformed into systems that are larger and more complicated but that disentangle variables so that the systems are easier to solve.

A good way to learn more is to use a search term such as "using grobner bases for solving".

This leads to places such as:

Cook, Solving systems of polynomial equations and mathSE, Using Gröbner bases for solving polynomial equations. The latter link asks how this method can be used to solve the intersection of two conics and has many useful links.

My feeling is that even if there was a way to get Mathematica to "show its work", the result would be pages and pages of very hard to follow manipulations. This just seems to be the nature of these methods .. beyond very simple examples they are not amenable to pencil and paper and are best left to computers.

Related Question