Ray Infinite Cylinder Intersection in GLSL

Calculating the analytical intersection between a 3D ray and an infinite cylinder involves solving a quadratic equation derived from the geometric definition of a cylinder. In GLSL, this process projects the ray's origin and direction vectors onto the plane perpendicular to the cylinder's axis, reducing the 3D problem into a 2D circle-ray intersection to analytically solve for the ray distance parameter \(t\).

Mathematical Formulation

A 3D ray is defined as:

\[\mathbf{r}(t) = \mathbf{o} + t\mathbf{d}\]

Where:

• \(\mathbf{o}\) is the ray origin.
• \(\mathbf{d}\) is the normalized ray direction.
• \(t \ge 0\) is the distance along the ray.

An infinite cylinder can be defined by a point on its axis \(\mathbf{p}\), a normalized axis vector \(\mathbf{a}\), and a radius \(r\). Any point \(\mathbf{x}\) lies on the cylinder's surface if its perpendicular distance to the cylinder's axis equals \(r\):

\[\|(\mathbf{x} - \mathbf{p}) - ((\mathbf{x} - \mathbf{p}) \cdot \mathbf{a})\mathbf{a}\|^2 = r^2\]

Substituting the ray equation \(\mathbf{x} = \mathbf{o} + t\mathbf{d}\) and defining \(\Delta\mathbf{p} = \mathbf{o} - \mathbf{p}\) yields a quadratic equation in terms of \(t\):

\[A t^2 + 2B t + C = 0\]

Where:

• \(\mathbf{d}_\perp = \mathbf{d} - (\mathbf{d} \cdot \mathbf{a})\mathbf{a}\)
• \(\Delta\mathbf{p}_\perp = \Delta\mathbf{p} - (\Delta\mathbf{p} \cdot \mathbf{a})\mathbf{a}\)
• \(A = \mathbf{d}_\perp \cdot \mathbf{d}_\perp\)
• \(B = \mathbf{d}_\perp \cdot \Delta\mathbf{p}_\perp\)
• \(C = \Delta\mathbf{p}_\perp \cdot \Delta\mathbf{p}_\perp - r^2\)

The discriminant is:

\[\Delta = B^2 - A C\]

• If \(\Delta < 0\), the ray misses the cylinder completely.
• If \(\Delta \ge 0\), the intersections are given by \(t = \frac{-B \pm \sqrt{\Delta}}{A}\).

GLSL Implementation

The following GLSL function calculates the entry and exit intersection distances. It returns vec2(t_near, t_far). If the ray does not intersect the cylinder, it returns vec2(-1.0).

// Intersects a ray with an infinite cylinder
// ro: Ray origin
// rd: Ray direction (must be normalized)
// pa: Point on cylinder axis
// axis: Normalized axis direction of the cylinder
// radius: Cylinder radius
vec2 intersectInfiniteCylinder(vec3 ro, vec3 rd, vec3 pa, vec3 axis, float radius) {
    vec3 dp = ro - pa;

    // Project direction and origin onto plane perpendicular to cylinder axis
    vec3 d_proj = rd - axis * dot(rd, axis);
    vec3 dp_proj = dp - axis * dot(dp, axis);

    float a = dot(d_proj, d_proj);
    float b = dot(d_proj, dp_proj);
    float c = dot(dp_proj, dp_proj) - radius * radius;

    float discriminant = b * b - a * c;

    if (discriminant < 0.0) {
        return vec2(-1.0); // No intersection
    }

    float sqrt_d = sqrt(discriminant);
    float t0 = (-b - sqrt_d) / a;
    float t1 = (-b + sqrt_d) / a;

    return vec2(t0, t1);
}

Computing the Surface Normal

To shade the cylinder at the intersection point \(\mathbf{hit} = \mathbf{o} + t\mathbf{d}\), project the vector from the axis to the hit point to get the surface normal:

vec3 getCylinderNormal(vec3 hitPoint, vec3 pa, vec3 axis) {
    vec3 d = hitPoint - pa;
    return normalize(d - axis * dot(d, axis));
}