Find the point on a line that's closest to an infinite cone

Question

I have two objects:

• A line $L$ defined by one of its points $l_0$ and a direction vector $d$. (I.e. $l_0 + t \cdot d$)
• An infinite cone $C$ defined by its apex $a$, a direction $e$ and a half-angle $\alpha < \frac{\pi}{2} $.

(Geogebra)

I want to find the point on the line that is closest to the cone. In the case that the two objects intersect (in more than one point), I would expect the algorithm to return the point that's furthest "inside" the cone.

So far I failed to figure this out myself. Approaches I thought about/tried:

• The point closest to the cone's axis is not necessarily the one closest to the cone itself. I found out by finding a counter-example.
• I read about ways to calculate the intersection(s) between a line and a cone, but that doesn't help me in the case that there is no intersection: then I have no clue.
• I thought about finding a formula for cone-point-distance and then substitute my line equation for the point, finding a minimum in that function in regards to $t$. But I haven't found anything that looks like I could analytically find a minimum.

With this, I don't know how to proceed, hence this question. Any help is greatly appreciated :)

Answer 1

animated geogebra figure here.

Fig. 1: Initial cone (C) in violet. "Enlarged cone" (C') in pale blue, tangent to given line (L). We have chosen to take point $B$ at the same altitude as apex $A$. The shortest distance point $S$ on $(L)$ is projected as $S'$ on $z$ axis.

The equation of the cone is

$$x^2+y^2=t^2(z-a)^2 \ \ \text{where} \ t:=\tan \alpha \tag{1}$$

The issue is equivalent to find a smallest distance segment between the given line and a certain generatrix $AS''$ of cone $(C)$. It is known that this smallest distance segment is both orthogonal to the line and the generatrix ; if one considers now any cone $(C')$ with the same apex and axis as the initial cone, with half-aperture angle $\alpha'$, the shortest distance segment from line $(L)$ to cone $(C')$ will be contained in the same plane.

The generic equation of such cones is evidently :

$$x^2+y^2=t'^2(z-a)^2 \ \ \text{where} \ t':=\tan (\alpha') \tag{2}$$

Besides, the parametric equations of the line can be taken in this way :

$$\begin{cases}x&=& x_0 + \lambda u \\ y&=& y_0 + \lambda v\\ z&=& a + \lambda w \end{cases}\tag{3}$$

Indeed, one can assume that the origin point on the line has the same ordinate as the apex of the cone (the only case where this is not possible is for horizontal lines ; a case deserving a - simple - separate treatment).

Now, consider among all cones $(C')$ the one which allows to transform the initial problem where the line doesn't cross $(C)$ into an equivalent tangency issue for cone $(C')$ : we are going to look for a value of $t'$ such that the generic point of the line (equations (3)) is a double point of $(C')$.

Let us express that the point $M$ belongs to the cone by plugging equations (3) into (2) :

$$(x_0+\lambda u)^2+(y_0+\lambda v)^2=t'^2(\lambda w)^2 $$

giving rise to the quadratic equation in $\lambda$:

$$\lambda^2(u^2+v^2-t'^2w^2)+2 \lambda(x_0u+y_0v)+(x_0^2+y_0^2)=0\tag{4}$$

A double root means annihilation of the discriminant of (4) :

$$(x_0u+y_0v)^2-(x_0^2+y_0^2)(u^2+v^2-t'^2w^2)=0,$$

a condition giving the looked-for value of $t'$ :

$$t'^2=\frac{1}{w^2}\left(u^2+v^2-\frac{(x_0u+y_0v)^2}{x_0^2+y_0^2}\right).\tag{5}$$

The double root condition gives the unique solution for (4) :

$$\lambda_S=- \frac{x_0u+y_0v}{u^2+v^2-t'^2w^2}=-\frac{x_0^2+y_0^2}{x_0u+y_0v}\tag{6}$$

(by taking (5) into account). It remains to plug this value of $\lambda$ into (4) to get the coordinates of the shortest distance point $S$ on the line (see figure).

Important remark : As asked by Hosam H, the shortest distance from $S$ to the given cone is obtained plainly by computing

$$d=AS \sin(\alpha'-\alpha)$$

where the half-aperture angle of the enlarged cone is computed in this way :

$$\alpha'=\arccos\left(\frac{\vec{AS}.\vec{AO}}{\|\vec{AS}\|\|\vec{AO}\|}\right)$$

Remark : as pointed out by Lukas, there is a case where this method is not applicable : when (4) is not a quadratic equation, i.e. when its leading term zero, i.e., when

$$(u^2+v^2-t'^2w^2)=0 \ \iff \ \text{using (5) : } \ x_0u+y_0v=0$$

Answer 2

Disclaimer: This is not the optimal solution of the problem. The optimal solution is in my other answer to this question.

The equation of the cone with vertex $\mathbf{V}$ and axis unit vector $\mathbf{a}$ and semi-vertical angle $\theta$ is given by

$ ( \mathbf{r} - \mathbf{V}) ^T Q_c (\mathbf{r} - \mathbf{V}) = 0 \hspace{30pt} (1) $

where $ r = [x, y, z]^T, Q_c = \cos(\theta)^2 I - {\mathbf{a a}}^T \hspace{30pt} (1.1) $

while the line is given by

$ \mathbf{P}(t) = \mathbf{\ell_0} + t \ \mathbf{d} \hspace{30pt} (2)$

If $\mathbf{r}$ and $\mathbf{P}$ are the pair of points that are at the minimum distance possible, then the vector $\mathbf{P} - \mathbf{r}$ is along the gradient (the normal vector) of the cone, and also $\mathbf{P} - \mathbf{r}$ is orthogonal to the line direction vector $\mathbf{d}$. Thus, we now have

$ \mathbf{P}- \mathbf{r} = k \ Q_c (\mathbf{r} - \mathbf{V}) \hspace{30pt} (3)$ for some $ k \in \mathbb{R} $

$ (\mathbf{P} - \mathbf{r}) \cdot \mathbf{d} = 0 \hspace{30pt} (4)$

The unknowns here at $x, y, z$ of $\mathbf{r}$ and $t$ of the point $\mathbf{P}$. Substitute what $\mathbf{P}$ is:

$ \mathbf{\ell_0} + t \mathbf{d} - \mathbf{r} = k \ Q_c (\mathbf{r} - \mathbf{V} ) \hspace{30pt}(3) $

$ ( \mathbf{\ell_0} + t \mathbf{d} - \mathbf{r} ) \cdot \mathbf{d} = 0 \hspace{30pt} (4) $

We need to eliminate $k$ from the vector equation $(3)$, and this can done by taking the ratio of the components of the respective vectors on either side of the equation, and then cross multiplying. This will give

If $ \mathbf{\ell_0} = [x_0, y_0, z_0 ]^T $ and $ \mathbf{d} = [d_x, d_y, d_z]^T $, and $\mathbf{V} =[V_x, V_y, V_z ]^T $ , then the above process leads to two scalar equations

$ (x_0 + t d_x - x) (Q_{21} (x - V_x) + Q_{22} (y - V_y) + Q_{23} (z - V_z) ) = \\ = (y_0 + t d_y - y) (Q_{11} (x - V_x) + Q_{12} (y - V_y) + Q_{13} (z - V_z) ) \hspace{30pt}(5)$

and,

$ (x_0 + t d_x - x) (Q_{31} (x - V_x) + Q_{32} (y - V_y) + Q_{33} (z - V_z) ) = \\ = (z_0 + t d_z - z) (Q_{11} (x - V_x) + Q_{12} (y - V_y) + Q_{13} (z - V_z) ) \hspace{30pt} (6) $

The above two equations (equations $(5)$ and $(6)$ ) are quadratic in $x,y,z,t$, and in addition to them we have the quadratic equation $(1)$ and the linear equation $(4)$.

These equations can be solved using the famous Newton-Raphson multivariate root finder, and with only $4$ variables, and quadratic and linear functions, it converges to the solution very quickly (less than $10$ iterations).

However, it should be noted that the Newton-Raphson method returns a root based on the initial guess, and there can be several roots to the above equations, one corresponding to the minimum, while the others corresponding to saddle points, it cannot be guaranteed that we'll get the correct root, and this limits the viability of using the Newton-Raphson method. The alternative is to use a solver as the ones available on Mathematica, or elsewhere. I happen to have the source code of a solver of $3$ quadratic equation in $3$ unknowns. In this case, we have $3$ quadratic equations and $1$ linear equation in the four variables $x,y,z,t$. The linear equation makes it possible to express the $x,y,z,t$ in terms of three parameters $u_1, u_2, u_3$. And then the three quadratic equation of $x,y,z,t$ are transformed into three quadratic equations in terms of $u_1, u_2, u_3$. And these three I can solve with my routine (that I have the source code for). The result of solving is a list of all possible solutions, which can be up to $8$ different solutions. Now all the solutions are known at once, and by evaluating the distance function at each of them we can pick the minimizing one.

Final remark: This solution I wrote several days ago, and today I found an alternative method that is much better than the one here. Please check my other solution, which only requires solving a simple scalar quadratic equation.

Answer 3

With Lagrange multipliers.

Giving the cone $(p-p_0)'A(p-p_0)=0$ and the line $q = q_0+\mu v$ with $\|v\|=1$ we need to find the point $p^*$ such that $\|q-p^*\|$ is minimum. Forming the Lagrangian

$$ L(p,\mu,\lambda) = \|q_0+\mu v-p\|^2+\lambda(p-p_0)'A(p-p_0) $$

the stationary points are determined by solving

$$ \cases{ (q_0+\mu v-p)'-\lambda(p-p_0)'A = 0\\ (q_0+\mu v-p)'v=0\\ (p-p_0)'A(p-p_0)=0 } $$

from the first equation we have

$$ q_0+\mu v-p -\lambda A'(p-p_0) = 0 \Rightarrow (\lambda A'+I)p=\lambda A'p_0+q_0+\mu v $$

from the second equation we have

$$ \mu = (p-q_0)'v $$

and substituting into the previous equation

$$ (\lambda A'+I)p=\lambda A' p_0+q_0+((p-q_0)'v)v $$

but

$$ ((p-q_0)'v)v = (v v')(p-q_0) $$

so we follow with

$$ (\lambda A'+I)p=\lambda A' p_0+q_0+(v v')(p-q_0)\Rightarrow p = (\lambda A'-v v'+I)^{-1}(\lambda A'p_0-(v v'-I)q_0)=\phi(\lambda) $$

now substituting $\phi(\lambda)$ into the conic we can determine $\lambda^*$ such that

$$ (\phi(\lambda^*)-p_0)'A(\phi(\lambda^*)-p_0)=0 $$

once we have $\lambda^*$ we have also $p^*$. Follows a plot for the case

$$ \cases{ v = \left(\frac{1}{\sqrt{3}},-\frac{1}{\sqrt{3}},\frac{1}{\sqrt{3}}\right)'\\ p_0 = (0,0,0)'\\ q_0 = (-8,5,0)'\\ A = \left( \begin{array}{ccc} \frac{19}{28} & -\frac{1}{7} & -\frac{3}{14} \\ -\frac{1}{7} & \frac{13}{28} & -\frac{3}{7} \\ -\frac{3}{14} & -\frac{3}{7} & \frac{3}{28} \\ \end{array} \right)\\ \lambda^* = 1.08342 } $$

Answer 4

Without loss of generality, you can place the apex at the origin and orient the coordinate system such that z-axis is along the cone axis.

Use homogeneous coordinates to define a ray as

$$ \pmatrix{x \\ y \\ z \\ 1} = \pmatrix{o_x\\o_y\\o_z\\1} + t\,\pmatrix{e_x\\e_y\\e_z\\0} \tag{1}$$

where $o_i$ are the coordinates of any point on the line, and $e_i$ are the components of the direction of the line.

In the same coordinate system, the equation of the cone can be described by the quadratic form

$$ \pmatrix{x \\ y \\ z \\ 1}^\top \begin{bmatrix} m & 0 & 0 & 0 \\ 0 & m & 0 & 0 \\ 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix} \pmatrix{x \\ y \\ z \\ 1} = 0 \tag{2} $$

where $m = \cot^2 \left( \psi/2 \right)$ is the conicity of the cone of included angle $\psi$.

Plug (1) into (2) to find the intersection points from the resulting quadratic

$$ t^2 \left( m\,(e_x^2 + e_y^2) - e_z^2 \right) +2 t \left( m (e_x o_x + e_y o_y) - e_z o_z \right) + \left( m (o_x^2+o_y^2)-o_z^2 \right) = 0 \tag{3}$$

If the line does not intersect, then the real part of the solution for $t$ will describe the closest point.

$$ t_{\rm closest} = \text{-} \frac{ m ( e_x o_x + e_y o_y) - e_z o_z}{m (e_x^2 + e_y^2)-e_z^2} $$

Answer 5

We have a right circular cone with apex $a$, and axis along unit vector $e$, with semi-vertical angle $\alpha$, and we have the line

$ P(t) = \ell_0 + t \ d $

The (signed) distance from $P(t)$ to the surface of the cone is simply

$ d = \| u \| \sin( \theta - \alpha) \hspace{20pt} (1)$

where the vector $u = P(t) - a = \ell_0 - a + t d $

And, $\theta$ is the acute angle between $u$ and the axis $e$, i.e.

$ \theta = \cos^{-1} \dfrac{|u. e|}{\| u \|} \hspace{20pt}(2)$

Expanding $(1)$,

$ d = \| u \| \bigg( \sin \theta \cos \alpha - \cos \theta \sin \alpha \bigg) $

Using $(2)$ this becomes

$ d = \| u \| \bigg( \cos \alpha \sqrt{ 1 - \left( \dfrac{u.e}{\|u\|} \right)^2 } - \sin \alpha \dfrac{|u.e|}{\|u \|} \bigg) $

And this simplifies to,

$ d = \cos \alpha \sqrt{ \| u \|^2 - (u.e)^2 } - \sin \alpha \ |u.e| $

Pluggin in $u$, this becomes

$ d = \cos \alpha \sqrt{ (\ell_0 - a + t d ) \cdot (\ell_0 - a + t d) - ( e \cdot (\ell_0 - a + t d ) )^2 } - \sin \alpha \ |e \cdot (\ell_0 - a + t d) | $

Differentiating $d$ with respect to $t$, and equating to $0$, gives us,

$ 0 = \cos \alpha \dfrac{ d \cdot (\ell_0 - a + t d) - (e \cdot d) (e \cdot (\ell_0 - a + t d ) }{\sqrt{ (\ell_0 - a + t d ) \cdot (\ell_0 - a + t d) - ( e \cdot (\ell_0 - a + t d ) )^2 }} - S_1 \sin \alpha \ (e \cdot d ) \hspace{20pt}(3) $ where $S_1 = \operatorname{sgn}( e \cdot (\ell_0 - a + t d) ) $

($S_1$ will be $+1$ if its argument is positive and $-1$ if it is negative)

Define the following constants

$ A = d \cdot d - (e \cdot d) ^ 2$

$ B = d \cdot (\ell_0 - a ) - (e \cdot d ) ( e \cdot (\ell_0 - a) ) $

$ C = (\ell_0 - a) \cdot (\ell_0 - a) - ( e \cdot (\ell_0 - a) )^2$

$ D = e \cdot d $

Then the above equation becomes

$ \cos \alpha \dfrac{ A t + B }{\sqrt{A t^2 + 2 B t + C}} \pm D \sin \alpha = 0 $

Hence,

$ \cos^2 \alpha (A t + B)^2 = D^2 \sin^2 \alpha ( A t^2 + 2 B t + C ) \hspace{20pt} (4) $

and this is a quadratic equation in $t$. To put it in standard form, re-arrange $(4)$ by combining like powers of $t^2 , t , 1$, you will get

$ a_2 t^2 + a_1 t + a_0 = 0 \hspace{20pt}(5)$

where

$a_2 = A^2 \cos^2 \alpha - A D^2 \sin^2 \alpha $

$a_1 = 2 A B \cos^2 \alpha - 2 B D^2 \sin^2 \alpha $

$a_0 = B^2 \cos^2 \alpha - C D^2 \sin^2 \alpha$

Note that $(5)$ and $(4)$ are exactly the same equation.

Now, the two solutions of $(5)$ are:

$ t_{1,2} = \dfrac{ - a_1 \pm \sqrt{ a_1^2 - 4 a_2 a_0 } }{2 a_2}$

To validate the two roots $t_1, t_2$ obtained, we have to substitute the roots in equation $(3)$. One of the two roots will be valid (that is, it will satisfy equation (3)) while the other will not. To find the minimum distance, just evaluate $d$ at the valid root. To find the closest point to the cone, substitute the valid root into the line equation.

I've tested this method against my other method, and found that both methods give the same results.

Final Remark:

In addition to the above, it is necessary to check first whether the line intersects the cone. Recall that

$ d = \cos \alpha \sqrt{ (\ell_0 - a + t d ) \cdot (\ell_0 - a + t d) - ( e \cdot (\ell_0 - a + t d ) )^2 } - \sin \alpha \ | e \cdot (\ell_0 - a + t d) | $

So, setting $d = 0$, we have to find if the following quadratic equation has real roots

$ \cos^2 \alpha ((\ell_0 - a + t d ) \cdot (\ell_0 - a + t d) - ( e \cdot (\ell_0 - a + t d ) )^2 ) = \sin^2 \alpha \ (e \cdot (\ell_0 - a + t d) )^2 $

In terms of $A,B,C,D$ defined above, this last equation is

$ \cos^2 \alpha \ ( A t^2 + 2 B t + C ) = \sin^2 \alpha \ ( D t + E )^2 \hspace{20pt}(6)$

where $ E = e \cdot (\ell_0 - a) $

Equation $(6)$ can be re-arranged as follows:

$b_2 t^2 + b_1 t + b_0 = 0 \hspace{20pt}(7)$

where

$b_2 = A \cos^2 \alpha - D^2 \sin^2 \alpha $

$b_1 = 2 B \cos^2 \alpha - 2 D E \sin^2 \alpha $

$b_0 = C \cos^2 \alpha - E^2 \sin^2 \alpha $

So that the solution of $(6)$ are

$ t_{3,4} = \dfrac{ - b_1 \pm \sqrt{ b_1^2 - 4 b_2 b_0 } }{ 2 b_2 } \hspace{20pt} (8)$

Provided that the discriminant

$ \Delta = b_1^2 - 4 b_2 b_0 \hspace{20pt}(9)$

Is non-negative, i.e. solutions will exist provided that $\Delta \ge 0 $.

If this equation has real roots then the line intersects the cone at values of $t$ equal to the two roots, and in this case the minimum distance predicted by taking the derivative does not apply, because it is a signed distance, meaning that if the line goes inside the cone, then its distance from the surface of the cone at these internal points will be negative. The most negative distance in this case is the one reported by taking the derivative and setting it to zero. So in the case, the negative distance is ignored and the minimum distance is taken to be $0$ corresponding to the crossings of the line with the cone.

Otherwise, if equation $(6)$ has no real roots, then the line does not cross the cone, and hence does not get inside the cone, and the distance $d$ will always be positive. In this case, the minimum (positive) distance and the closest point to the cone are determined by the valid root of $(5)$.