How can I work out if a point x,y,z is contained by a cone ?

January 14, 2007 5:33 PM

How can I work out if a point x,y,z is contained by a cone ?

I'm exploring some 3d programming, but sadly geometry is not my strong point.

I can define a cone object :

P1(x,y,z) - point of cone
P2(x,y,z) - centre of the base of cone
R1 - radius of the base.

How do I test if a point P3(x,y,z) is inside the surface of the cone ?
If it is inside the cone, how do I determine the closest point on the surface of the cone to P3 ?

I'd like to write a generic function so I can reuse it in other programming languages.

I'd prefer the cone base to be an elipse, I guess I could add a second raduis and rotation around the x-axis, but I'm out of my depth already.
I've looked through wikipedia and wolfram's but the crazy notation they use are all but meaningless to me.

Thanks

posted by matholio to computers & internet (15 comments total)

ok, here's a sketch with spots where math needs to be filled in:

the trick is to move all 3 points (p1,p2,and the test point which i'll call tp) around and simplify the problem.

1) move them so that p1 is at the origin(0,0,0), i.e. subtract p1 from p2 and tp,

2) rotate the p1-p2 line so that it falls along an axis (x,y,or z) it doesnt matter which one, just pick one, and apply the same rotational transformations to tp.

3) then you can see if p3 at least falls within p1 and p2, a precondition for being in the cone

4) if it does fall in that region, then measure the distance of tp to the axis that p1 and p2 lie along

5) so if it's closer to the axis than the radius of teh cone would be at that point of the axis
posted by garethspor at 5:49 PM on January 14, 2007

It should be possible to figure out the angle of the cone, and the angle of the point relative to the cone origin and its center line. If the latter angle is greater than the former, then the point is outside. (To save on math, you can compare the respective tangents instead of doing arctangents to get back to the actual angle.)

As garethspor says, this requires normalizing the coordinate system to the cone itself.
posted by Steven C. Den Beste at 5:59 PM on January 14, 2007

I see. yes that's making sense.

I'm not clear how to find the angles mentioned though.

Do I use the origin, cone length and diameter to form a right angle triangle ?
posted by matholio at 6:58 PM on January 14, 2007

I dug this up on google groups, see if the last answer has what you need.
posted by Osmanthus at 7:11 PM on January 14, 2007

Take the vector P2-P1 and the vector P3-P1. Normalize them both to unit length. Their dot product is then cos(angle P3-P1-P2). If this number is greater than or equal to the cosine of the half-angle at the apex of the cone, then the point is inside the cone. (If it's exactly equal, then P3 is on the cone.)

To find the nearest point on the cone: I could be wrong, but the nearest point (call it P4) should be one for which the line P3-P4 intersects the cone's surface at a right angle. (Visualize a sphere centered on P3 and tangent to the inside of the cone.) Of course, P4 is on the surface of the cone, so that means that you know the value of (P4-P1)*(P2-P1)/|P2-P1|; it has to be the cosine of the half-angle, from the first paragraph. And you know that P4-P3 is perpendicular to P4-P1, which means their dot product is zero. Since you're working in three dimensions, there's still one degree of freedom to get rid of --- these two constraints will still leave you with a closed path of possible points on the surface of the cone --- you could choose the P4 that's closest to P3, but the math on that might be bothersome; or you could maybe do something with cross-and-dot-products to constrain the answer to one in which P1,P2,P3,P4 are all coplanar. (Though you'll still end up with two solutions in that last case.)
posted by hattifattener at 7:41 PM on January 14, 2007

Personally, I'd recommend you want to represent your cone's origin and orientation as a subspace represented by a 4x4 matrix. You can use that matrix to map the worldspace point into the space occupied by the cone. To construct such a matrix, each row in the upper-left 3x3 is one of the basis vectors of the cone's coordinate system, the right column is the offset to the cone's origin.

Having transformed the point into cone-space, you just need to check whether the ratio of the distance of the point from the centreline to the height above the base is small enough for the cone's slope.

So, imagine a cone with the long axis along x, 2 units radius in y and 1 units radius in z and origin at the world origin. The subspace matrix is this:

Mcone=
| 0 0.5 0 0 |
| 0 0 1 0 |
| 1 0 0 0 |
| 0 0 0 1 |

Given a point in world space, P_w, P_c=M_coneP_w.

r² = P_c.x² + P_c.y²
h² = P_c.z²
InCone == (P_c.z >= 0 && P_c.z <= coneheight && h²slope² <= 1 - r²)

Slope is the slope of the cone side. 0 = parallel, i.e. a cylinder. 1 = 45 degrees, etc.
posted by polyglot at 7:46 PM on January 14, 2007

no offence to hattifattener (wtf is a hatti!?), cosines and normalisation suck because they are expensive and not so precise. Much better to keep to linear operations, no divides and no square roots.

The formulae for mine (I coulda blundered on them, they're off the top of my head) are based on the geometry of a right triangle that is a radial slice of one side of your cone. The baseline is the cone's baseline, length one (remember, cone size is defined by the matrix used upfront to map to conespace. the radius of the cone in cone space is one), the vertical is the cone's centreline and the hypotenuse is a line along the curved part of the cone.

h² is the height above the baseline squared, r² is the distance along the baseline squared. Reason for the squares is so that we don't have to square root both halves of the equation, which would suck. It's all high-school geometry after that.
posted by polyglot at 7:52 PM on January 14, 2007

I mostly agree with polyglot; however, I would put a normalization in as a preprocess, then try to avoid keeping around a full basis for the cone (when you only really need the axis of rotation.) I would derive x and y by taking the projected z (i.e., the point on the center of rotation closest to the candidate point,) back into world space, and subtracting it from the original candidate point in world space. You do need the length (here encoded in the normalization) to do the inverse transform, though.

Something like (warning! untested code!):
// Preprocess:
delta = p2 - p1;
deltaNorm = delta / length(delta);
deltaBegin = p1 dot deltaNorm;
slopeSquared = Radius * Radius / (coneheight * coneheight);

// Inner loop:
cp = candidate point
cpProj = cp dot deltaNorm;
cpProjNorm = cpProj - deltaBegin;
if ( cpProjNorm < 0 || cpprojnorm> coneheight ) return false;
cpDeproj = cp - (deltaNorm * cpProj);
cpRadiusSquared = slopeSquared * cpProjNorm * cpProjNorm;
return ( cpDeproj dot cpDeproj < cpRadiusSquared );
posted by blenderfish at 11:00 PM on January 14, 2007

polyglot: You don't actually need to do any cosine operations in my method. It happens that some numbers you end up with are cosines of interesting angles, but the only operations are multiplication, addition, and square root (for the normalization). I suspect you'll need a sqrt to get the transformation matrix from the cone parameters, and I also suspect that if you crunch through the algebra our techniques are the same.

I hadn't noticed the possible elliptical-cone requirement, though. That makes your "unit-cone transform" technique come out ahead on clarity. The "nearest point on cone" problem is a lot harder for an elliptical cone, though!
posted by hattifattener at 11:35 PM on January 14, 2007

Upon further reading of the question, I see the postet wanted to potentially deal with elliptical cones, in which case I think polyglot's approach is better.
A couple thoughts, though.
1> It's slightly more general than is required (i.e., it handles cylinders and cones with the top chopped off parallel to the base.) Putting the apex of the cone at the
origin (not the center of the base,) and assuming a slope of 1 would change the last check to (r^2 < h^2)br> 2> There is an implicit normalization (or matrix inversion) to get from the data the question provided (P1, P2, and friends), which are essentially a matrix
going from cone to world space, to your M[cone] matrix, which goes the other way. So, you still have that little bit of potential numerical nastiness. (Though, it
isn't happening in the inner loop, thankfully.)
3> In actual code, it would be good to wait to calculate r^2 until after the coneheight check.

On preview: agreeing with hattifattner.

So, moving on to the 'closest point to the surface' problem...
That's kind of esoteric. I'm kind of wondering what the inquisitor is planning on using that for. (Physics penetration?)
posted by blenderfish at 11:54 PM on January 14, 2007

hattifattener: I would say problem with your approach is that it has a normalization in the inner loop; when you find P3-P1.

Typically, I assume (since I'm coming at this with a games perspective,) when doing collision detection that the geometry will go through a setup phase at the beginning of a frame, then get checked against many times before the next frame. So, the more work one can offload to the former part of the process, the better. Of course, we don't know whether the poster wants to move the cone each time he checks it (in which case, it's a wash,) and we don't even know if performance is important.

Didn't see a problem with your use of cosine though; you
were using it to name an algebraic quantity, not actually calling cos() or (even more scary) acos(), any more that polyglot is calling tan() when he says "1 = 45 degrees."
posted by blenderfish at 12:09 AM on January 15, 2007

I had written down a bunch of maths for finding the nearest point on cone based on using the same right-triangle representing a slice of the cone and forming a smaller right-triangle with its baseline being the cone-wall and the ascender leading to the point in question. Analysis of that leads you to the nearest point in that radial slice, but this is the nearest point in general only if the cone does not have an elliptical base. If you want a solution for circular cones only, I can write it down for you, or hopefully you could figure it out from that description.

Consider a highly elliptical cone with a point on the long radial axis but near the cone's primary axis. The nearest point in that radial slice is further away than the two equally nearest points that are slightly off the short axis. d'oh. So you can't solve for the nearest point by analysing a 2D slice of the cone. It also means that you cannot do the analysis in the cone-space I defined earlier because that loses the elliptical nature of the cone.

My next approach would be to go for a numerical (Newton/Raphson) approach solving for a min distance. Represent the cone as a parametric surface over theta and altitude. You can compute the cartesian position of any point on the cone from those two parameters and therefore express the distance from your test point to the parametric point. Solve that using a 2-dimensional (theta, h) solver; a simple slope-following algorithm should be sufficient since I don't think there are any local minima to avoid on a simple shape like the cone.

That function would be damn interesting to look at :)
posted by polyglot at 12:57 AM on January 15, 2007

Thanks for all your comments.
Well, thanks for making them, I'm completely unable to use most of them, I'm just so out of my depth.

My use for this was in an alife project I'm toying with.
I wanted to define a sensor for an agent, I wanted the sensor to have varying capabilities, rather like eyes, ears and noses do. (wide, narrow, depth, angle)

*Closest cone surface* was, for an idea to have a fall-off function to simulate reduced accuracy on the periphery, but tbh, I hadn't really thought it through. I guess distance from the center line is easier and probably more useful.

I'm encouraged by the presumably non-trivial effort of the posters, and somewhat embarrassed that I don't actually get what's been written.

Probably the most poignant phrase was "It's all high-school geometry after that.". So it's back to the classroom for me (google), and I know I'm going to find a lot of info, much of which I expect I youthfully thought 'what would I ever need this for?'.

All that said, I'll probably opt for a rougher, simpler solution and just test vertical and horizontal angles and distance away from origin.
That will give me a square cone with a curved base right ?

Thank you all.
posted by matholio at 5:19 PM on January 15, 2007

Well if we've aimed the explanations at the wrong level, that's at least partly our fault for assuming too much. When I say "high school maths", I mean what's taught at year 12 (final year of school, whatever you have over there) Pure Mathematics here in Australia, which people only do if they're planning on a degree in the science and engineering field. It ain't for everyone.

If you're looking for information on this sort of stuff, I have to recommend Wikipedia because it has the most awesome mathematical explanations. Start with Linear Algebra, see the List of Topics (scroll down for link) and then follow up with Affine transformations, subspaces and homogeneous coordinates. Once you've got the basic geometry down, have a look into Quadric surfaces, of which the cone is a special case. Browse around there lots and soak up as much as you can and try googling for those terms.

You'd also do well to read tutorials & documentation for 3D programming, they often contain good information on geometry for people who don't know geometry yet. Reading about raytracing will also tell you a lot since most raytracers can perform intersections with arbitrary quadric surfaces and the equations to solve that are well known & widely described.

Angles from centre, i.e. atan2(x, z) and atan2(y, z) - assuming line of sight is z axis - and applying a threshold to those will you the points in a rectangular pyramid with the apex of the pyramid at the origin.

Good luck... geometry and computer graphics can be awesome fun and produce spectacular results. Of course you should feel free to ask more geometry questions here :)
posted by polyglot at 9:52 PM on January 15, 2007

Okay, hopefully this can clear up some things, or at least introduce some terms for you to research on your own:
(This is hatti's approach, broken down.)

The "dot product" of P1 and P2 is (P1.x * P2.x + P1.y * P2.y+ P1.z * P1.z) In goes two 3-d vectors, out comes a 1-d "scalar."

The "length" or "magnitude" of a vector is the square root of the dot product of the vector with itself.

The "difference" between two vectors is just the componentwise difference. So, if we have D=P1-P2, D.x = P1.x - P2.x, D.y = P1.y - P2.y, and D.z = P1.z - P2.z

The distance between two points is the magnitude of the difference between those points.

When you take the dot product of any two vectors P1 and P2, the quantity you get (remember, it's a 1d scalar) will happen to be the magnitude of P1 times the magnitude of P2 times the cosine of the angle between those vectors. So, the dot product of two perpendicular vectors is 0.

To 'scale' a vector by a scalar, multiply each component by that scalar. So, if we say Q = P1 * s, then Q.x = P1.x * s, Q.y = P1.y * s, and Q.z = P1.z * s.

To 'normalize' a vector, scale it by the multiplicative inverse of (that is one over) it's magnitude. A normalized vector has a magnitude of One, or is said to be 'unit-lengthed'.

The dot product, then, of two normalized vectors is just the cosine of the angle between them (since the magnitudes are one, and multiplying by one doesn't change the value.)


Okay, that being out of the way, lets say you have your sensor with P1 being the location of the 'eye', and P2 being a point in the middle of the 'eye's vision. If you take P2-P1 (subtract each element in P2 by P1, and put the 3 resulting numbers in a vector,) then you normalize that quantity, you will have a unit-length vector that points in the direction the eye is looking.

Now, if you do the same thing with your candidate point (we'll call it CP,) that is, take CP-P1, and normalize it, you'll have the direction from the center of your eye to your candidate point.

So, the question is probably then, for a sensor, how far from 'straight ahead' is the CP object in the eye's vision. So, what is the angle between 'straight ahead' and 'toward the CP.' Well, since both of your normalized vectors have length one, taking the dot product of those two vectors will give you the cosine of the angle between the two. 1 means dead-on-looking-straight-at-it, and 0 means completely-to-the-side-of-the-eye. Negative numbers mean your CP is behind the eye's vision cone.

Now, you can either compare that value with the cosine of some kind of threshold angle, or you can take the arc-cosine (a.k.a. inverse cosine) to get back the exact angle, to do whatever you want with.

For distance, you just keep around the magnitude of CP-P1 that you had to calculate anyway to normalize it.
posted by blenderfish at 10:47 PM on January 15, 2007