Point inside right circular cone

Question

I have already posted this on SO here with code however I'm really thinking that my line of thought (mathematical that is) is incorrect and that's why my code isn't working as intended so I decided to ask here since this is a mathematical problem.

I have a right circular cone with its tip ($\vec{tip} = (x_{tip},y_{tip},z_{tip})$) being the first point of the set of points and the center of the base being a point lying at the same $x$, $y$ coordinates but with $z$ equal to $0$ ($\vec{base} = (x_{tip},y_{tip},0)$). I have also a given angle between the cone's axis (tip to base) and the slant.

For every point in the set I have to check if it's inside the cone or not. For simplicity reasons I mark a point as valid if it is either really inside or at least lying on the cone's surface. A point is invalid if it lies outside.

In order to that I check two things:

• Is the point above the cone? - due to the range of values I have and the way I define my cone's base (at $z = 0$) a point cannot go below the cone's base (however it can lie on it in which case it is a valid one). This check is pretty simple. The height of the cone is already given by its tip (the first point from the set) so I just have to do a one-dimensional test comparing $z$ values. If the point is above the given height of the cone it is clearly outside.

$$ checkHeightOkay(z) : \begin{cases} yes & z_{point}\leq z_{tip} \\ no & else \end{cases} $$

• Is the point on the inside of/lying on the circumference of the slice of the cone at the given height or not? - this test is a two-dimensional one. After making sure that the point is not above the cone (that is it's $z$ coordinate is okay) I calculate the slice of the cone at the height of the point (not that in my case the cone goes along the $z$ and not the $x$ axis as in the image below):

with the following formula for the radius:

$$\frac{r_{slice}}{z_{point}} = \frac{R_{cone}}{h_{cone}} \Rightarrow r_{slice} = \frac{R_{cone} \cdot z_{point}}{h_{cone}}$$

For measuring if the point's x and y coordinates are part of this slice I use Euclidean distance:

$$dist(\vec{axis_{cone}}, \vec{point}) = \sqrt[]{(x_{tip} - x_{point})^2 + (y_{tip} - y_{point})^2}$$

And as a last step I check if $dist$ is $\le$ (inside) or $> r_{slice}$ (outside).

It seems however that this doesn't work since some points that are supposed to be valid are marked as invalid and vice versa (for more details and code you can check the link at the beginning of this post).

The algorithm is pretty straight forwards and I have really no idea what I'm doing wrong here.

Answer 1

One mistake, that I can spot, is that the relation:

$$\frac{r_{slice}}{z_{point}} = \frac{R_{cone}}{h_{cone}}$$

should be changed to

$$\frac{r_{slice}}{z_{tip}-z_{point}} = \frac{R_{cone}}{h_{cone}}.$$

I think the rest is fine.

Answer 2

This is not a direct answer to the stated question, but shows how to solve the general problem, if others stumble onto this particular question looking for a general, efficient solution.

---

If the apex of your cone is at $$\vec{a} = (x_a, y_a, z_a)$$ and the center of the base is at $$\vec{b} = (x_b, y_b, z_b)$$ and has radius $r$, then the height of the cone is $$h = \left\lVert \vec{a} - \vec{b} \right\rVert = \sqrt{(x_a-x_b)^2 + (y_a-y_b)^2 + (z_a-z_b)^2}$$ and the unit normal (unit axis, from apex towards center of base) vector is $$\hat{n} = \frac{\vec{b} - \vec{a}}{\left\lVert \vec{a} - \vec{b} \right\rVert} = (x_n, y_n, z_n) = \left ( \frac{x_b - x_a}{h}, \frac{y_b - y_a}{h}, \frac{z_b - z_a}{h} \right )$$ The constant $$C = \frac{r^2}{h} \; \; \; \; \left ( = r \tan\left(\frac{\theta}{2}\right)\right)$$ is useful later for checking point inclusion. If you happen to need it, the apex angle $\theta$ fulfills $$\cos\left(\frac{\theta}{2}\right) = \frac{r}{\sqrt{h^2 + r^2}}$$ i.e $$\theta = 2 \arccos\left(\frac{r}{\sqrt{h^2 + r^2}}\right) = 2 \arctan\left(\frac{r}{h}\right) = 2 \arctan\left(r, h\right)$$

To determine if point $\vec{p} = (x, y, z)$ is within the cone defined by apex $\vec{a} = (x_a, y_a, z_a)$, axis unit normal $\hat{n} = (x_n, y_n, z_n)$, height of the cone $h$ (and radius $r$ at the base), and $C$, the "cone shape constant" ($C = r\tan\left(\frac{\theta}{2}\right)$), first calculate the height $Y$ of the point $\vec{p}$ (from apex towards base): $$Y = \hat{n} \cdot \left ( \vec{p} - \vec{a} \right ) = x_n (x - x_a) + y_n (y - y_a) + z_n (z - z_a)$$

If $Y \lt 0$, point $\vec{p}$ is above the apex. If $Y \gt h$, the point is below the base. (Remember, $Y = 0$ at the apex, increasing to $Y = h$ at the base.)

Next, calculate the squared distance $X^2$ between the point $\vec{p}$ and the axis of the cone: $$X^2 = \left\lVert \vec{p} - \vec{a} - y \hat{n} \right\rVert^2 = \left( x - x_a - Y x_n \right)^2 + \left( y - y_a - Y y_n \right)^2 + \left( z - z_a - Y z_n \right)^2$$ The point $\vec{p}$ is within the cone if $$X^2 \le C Y$$

In pseudocode, one can calculate the features of the right circular cone first:

# (xA, yA, zA) = apex
# (xB, yB, zB) = center of base
#            r = radius at base

h = sqrt((xA-xB)*(xA-xB) + (yA-yB)*(yA-yB) + (zA-zB)*(zA-zB))

xN = (xB-xA)/h
yN = (yB-yA)/h
zN = (zB-zA)/h

C = r*r/h

Y0 = xN*xA + yN*yA + zN*zA

# Remember variables
#    xA,yA,zA      apex
#    xN,yN,zN      axis unit normal
#    Y0            height at apex
#    h             height of cone
#    C             related to cone apex angle

For each right circular cone, the above precalculations use 1 square root operation, 7 multiplications, 4 divisions, 11 additions or multiplications, and requires 9 variables stored per cone.

For each point x,y,z, the determination is then

coneY = x*xN + y*yN + z*zN - Y0
if (coneY < 0) then:
    # point is above the apex
else if (coneY > h) then:
    # point is below the base
end if

tempX = X - xA - coneY*xN
tempY = Y - yA - coneY*yN
tempZ = Z - zA - coneY*zN

coneR2 = tempX*tempX + tempY*tempY + tempZ*tempZ

if (coneR2 > C * coneY) then:
    # point is outside the cone
else:
    # point is inside the cone, or on the surface of the cone
end if

which requires a maximum of 10 multiplications and 11 additions or subtractions, and up to three comparisons per point (using precalculated cone details).

---

If the axis of a right circular cone is parallel to the $z$ axis, and passes through origin, with apex at $(0,0,h)$ and base at origin, radius $r$ at base, then point $(x,y,z)$ is within the cone if and only if $$\begin{cases} z \ge 0 \\ z \le h \\ x^2 + y^2 \le (z - h)^2 \frac{r^2}{h^2} \end{cases}$$