Analytical Ray-Cone Intersection in GLSL

Analytically computing the intersection between a ray and a cone in a GLSL shader requires framing the cone as an implicit geometric surface, substituting the parametric ray equation into this implicit form to derive a quadratic equation, solving for the distance parameter \(t\), and validating the result against height and directional bounds. This process allows real-time ray marchers and ray tracers to render sharp, mathematically exact conical primitives without relying on polygonal approximations or step-based signed distance fields.

Geometric Formulation

A parametric ray is defined by its origin \(\mathbf{O}\) and normalized direction \(\mathbf{D}\):

\[\mathbf{P}(t) = \mathbf{O} + t\mathbf{D}, \quad t \ge 0\]

An infinite double cone with apex \(\mathbf{V}\), normalized axis vector \(\mathbf{A}\), and half-angle \(\theta\) satisfies the condition that the angle between any vector on the surface \((\mathbf{P} - \mathbf{V})\) and the axis \(\mathbf{A}\) equals \(\theta\):

\[\frac{(\mathbf{P} - \mathbf{V}) \cdot \mathbf{A}}{\|\mathbf{P} - \mathbf{V}\|} = \cos\theta\]

Squaring both sides eliminates the square root from the vector magnitude and yields the implicit quadratic equation for an infinite double cone:

\[((\mathbf{P} - \mathbf{V}) \cdot \mathbf{A})^2 = \|\mathbf{P} - \mathbf{V}\|^2 \cos^2\theta\]

Deriving the Quadratic Equation

Let \(\mathbf{\Delta} = \mathbf{O} - \mathbf{V}\). Substituting \(\mathbf{P}(t) - \mathbf{V} = \mathbf{\Delta} + t\mathbf{D}\) into the implicit equation:

\[(t(\mathbf{D} \cdot \mathbf{A}) + (\mathbf{\Delta} \cdot \mathbf{A}))^2 = (t^2(\mathbf{D} \cdot \mathbf{D}) + 2t(\mathbf{D} \cdot \mathbf{\Delta}) + (\mathbf{\Delta} \cdot \mathbf{\Delta})) \cos^2\theta\]

Since \(\mathbf{D}\) is normalized, \(\mathbf{D} \cdot \mathbf{D} = 1\). Expanding and rearranging terms into the standard quadratic form \(a t^2 + b t + c = 0\) gives:

\[a = (\mathbf{D} \cdot \mathbf{A})^2 - \cos^2\theta\] \[b = 2 \left( (\mathbf{D} \cdot \mathbf{A})(\mathbf{\Delta} \cdot \mathbf{A}) - (\mathbf{D} \cdot \mathbf{\Delta})\cos^2\theta \right)\] \[c = (\mathbf{\Delta} \cdot \mathbf{A})^2 - (\mathbf{\Delta} \cdot \mathbf{\Delta})\cos^2\theta\]

Boundary and Height Constraints

Because the implicit equation describes a double cone extending infinitely in both directions along \(\mathbf{A}\), the solutions must be filtered:

1. Single Sheet Constraint: For a cone extending only in the positive direction of \(\mathbf{A}\), the projection of the hit point onto the axis must be positive: \((\mathbf{P}(t) - \mathbf{V}) \cdot \mathbf{A} > 0\).
2. Height Truncation: For a finite cone of height \(h\), the projection must satisfy \(0 \le (\mathbf{P}(t) - \mathbf{V}) \cdot \mathbf{A} \le h\).
3. Caps (Optional): If the cone is closed at its base, a standard ray-plane intersection can be evaluated at distance \(h\) along \(\mathbf{A}\), bounded by radius \(r = h \tan\theta\).

Surface Normal Calculation

The gradient of the implicit function provides the surface normal at intersection point \(\mathbf{P}\):

\[\nabla f(\mathbf{P}) = 2 ((\mathbf{P} - \mathbf{V}) \cdot \mathbf{A})\mathbf{A} - 2(\mathbf{P} - \mathbf{V})\cos^2\theta\]

Normalizing the outward-facing vector gives:

\[\mathbf{N} = \text{normalize}\left( (\mathbf{P} - \mathbf{V}) - \frac{(\mathbf{P} - \mathbf{V}) \cdot \mathbf{A}}{\cos^2\theta}\mathbf{A} \right)\]

GLSL Implementation

struct Cone {
    vec3 apex;
    vec3 axis;      // Normalized axis direction
    float cosTheta; // cos(half-angle)
    float height;   // Maximum height along axis
};

bool intersectCone(
    vec3 ro, 
    vec3 rd, 
    Cone cone, 
    out float hitT, 
    out vec3 hitNormal
) {
    vec3 delta = ro - cone.apex;

    float dDotA = dot(rd, cone.axis);
    float deltaDotA = dot(delta, cone.axis);
    float dDotDelta = dot(rd, delta);
    float deltaDotDelta = dot(delta, delta);
    float cos2 = cone.cosTheta * cone.cosTheta;

    float a = dDotA * dDotA - cos2;
    float b = 2.0 * (dDotA * deltaDotA - dDotDelta * cos2);
    float c = deltaDotA * deltaDotA - deltaDotDelta * cos2;

    float discriminant = b * b - 4.0 * a * c;
    if (discriminant < 0.0) {
        return false;
    }

    float sqrtDisc = sqrt(discriminant);
    float t0 = (-b - sqrtDisc) / (2.0 * a);
    float t1 = (-b + sqrtDisc) / (2.0 * a);

    if (t0 > t1) {
        float temp = t0;
        t0 = t1;
        t1 = temp;
    }

    // Test intersections against height bounds
    float tCandidates[2] = float[2](t0, t1);
    for (int i = 0; i < 2; ++i) {
        float t = tCandidates[i];
        if (t > 0.0) {
            vec3 p = ro + t * rd;
            float projection = dot(p - cone.apex, cone.axis);

            if (projection >= 0.0 && projection <= cone.height) {
                hitT = t;
                
                // Analytical normal computation
                vec3 tipToPoint = p - cone.apex;
                vec3 normal = normalize(tipToPoint - (projection / cos2) * cone.axis);
                
                // Orient normal against incoming ray
                hitNormal = (dot(rd, normal) < 0.0) ? normal : -normal;
                return true;
            }
        }
    }

    return false;
}