Verify if point is inside a cone in 3D space

Question

Consider:

• X(x1,y1,z1) the point I need to verify if it is inside a cone.
• M(x2,y2,z2) the vertex of the cone. (the top point of the cone)
• N(x3,y3,z3) the point in the middle of the cone's base.

I found out that if a point X is on the cone, it needs to verify this equation:

cos(alfa) * ||X-M|| * ||N|| = dot(X-M,N)

Where dot is the scalar product of 2 vectors, and alfa is the angle between these 2 vectors.

Based on the formula, I calculated that:

X-M = (x1-x2,y1-y2,z1-z2)

So,

cos(alfa)
  * Math.sqrt((x1-x2)^2+(y1-y2)^2+(z1-z2)^2)
  * Math.sqrt(x3^2 + y3^2+z3^2)
= x3(x1-x2) + y3(y1-y2) + z3(z1-z2)

Unfortunatelly the above calculations seem to give me wrong results. What am I doing wrong?

Also I suspect that to check if X is inside the cone, I have to put <= instead of = in the formula. Is this correct?

The usage of this is: I develop a game where a machine gun has to start firing when an object is in its 'view'. This view will be a cone. The cone's vertex would be in the machine gun, the base of the cone will be at some known distance ahead. Any object entering this cone, the machine gun will shoot it.

Answer 1

I totally agree with Tim: we need "angle" (aperture) of cone to get the answer.

Let's do some coding then! I'll use some terminology from here.

Result-giving function:

/**
 * @param x coordinates of point to be tested 
 * @param t coordinates of apex point of cone
 * @param b coordinates of center of basement circle
 * @param aperture in radians
 */
static public boolean isLyingInCone(float[] x, float[] t, float[] b, 
                                    float aperture){

    // This is for our convenience
    float halfAperture = aperture/2.f;

    // Vector pointing to X point from apex
    float[] apexToXVect = dif(t,x);

    // Vector pointing from apex to circle-center point.
    float[] axisVect = dif(t,b);

    // X is lying in cone only if it's lying in 
    // infinite version of its cone -- that is, 
    // not limited by "round basement".
    // We'll use dotProd() to 
    // determine angle between apexToXVect and axis.
    boolean isInInfiniteCone = dotProd(apexToXVect,axisVect)
                               /magn(apexToXVect)/magn(axisVect)
                                 >
                               // We can safely compare cos() of angles 
                               // between vectors instead of bare angles.
                               Math.cos(halfAperture);


    if(!isInInfiniteCone) return false;

    // X is contained in cone only if projection of apexToXVect to axis
    // is shorter than axis. 
    // We'll use dotProd() to figure projection length.
    boolean isUnderRoundCap = dotProd(apexToXVect,axisVect)
                              /magn(axisVect)
                                <
                              magn(axisVect);
    return isUnderRoundCap;
}

Below are my fast implementations of basic functions, required by the upper code to manipulate vectors.

static public float dotProd(float[] a, float[] b){
    return a[0]*b[0]+a[1]*b[1]+a[2]*b[2];
}

static public float[] dif(float[] a, float[] b){
    return (new float[]{
            a[0]-b[0],
            a[1]-b[1],
            a[2]-b[2]
    });
}

static public float magn(float[] a){
    return (float) (Math.sqrt(a[0]*a[0]+a[1]*a[1]+a[2]*a[2]));
}

Have fun!

Answer 2

You need to check whether the angle between your difference vector (X-M) and your center vector (N) is less than or equal to the angle of your cone (which you haven't specified in the question). This will tell you if the position vector (X) is inside the infinite cone, and you can then also check for distance (if you want). So,

float theta = PI/6; //half angle of cone
if (acos(dot(X-M, N)/(norm(X-M)*norm(N)) <= theta) doSomething();

For performance, you could also normalize N (convert it to a vector with length 1) and store the length separately. You could then compare norm(X-M) to the length, giving you a round-bottom cone (for which I'm sure a name exists, but I don't know it).

Edit: Forgot the inverse cosine, the dot product is equal to norm(U)*norm(V)*cos(Angle) so we have to invert that operation to compare the angles. In this case, the acos should be fine because we want positive and negative angles to compare equally, but watch out for that.

Edit: Radians.