Closest Point on a 3D Triangle in GLSL

This article explains how to determine the closest point on an arbitrary 3D triangle to a query point within an OpenGL Shading Language (GLSL) shader. The process relies on identifying which Voronoi feature region of the triangle—a vertex, an edge, or the interior face—the query point projects onto, and then computing the corresponding clamped coordinates. Below is a detailed mathematical breakdown followed by an optimized, division-safe GLSL implementation.

The Algorithm

A 3D triangle defined by vertices \(A\), \(B\), and \(C\) splits the surrounding space into seven distinct geometric regions:

  1. Three vertex regions (\(A\), \(B\), \(C\))
  2. Three edge regions (\(AB\), \(BC\), \(CA\))
  3. One interior face region

Rather than projecting the point onto the triangle's plane and clipping against edge lines, the standard approach (derived from Christer Ericson's Real-Time Collision Detection) computes scalar products to determine which region contains the closest point. This method is computationally efficient, completely eliminates square roots, and defers divisions until the exact target feature is identified.

GLSL Implementation

vec3 closestPointOnTriangle(vec3 p, vec3 a, vec3 b, vec3 c) {
    vec3 ab = b - a;
    vec3 ac = c - a;
    vec3 ap = p - a;

    // Check if P is in vertex region outside A
    float d1 = dot(ab, ap);
    float d2 = dot(ac, ap);
    if (d1 <= 0.0 && d2 <= 0.0) return a;

    // Check if P is in vertex region outside B
    vec3 bp = p - b;
    float d3 = dot(ab, bp);
    float d4 = dot(ac, bp);
    if (d3 >= 0.0 && d4 <= d3) return b;

    // Check if P is in edge region of AB
    float vc = d1 * d4 - d3 * d2;
    if (vc <= 0.0 && d1 >= 0.0 && d3 <= 0.0) {
        float v = d1 / (d1 - d3);
        return a + v * ab;
    }

    // Check if P is in vertex region outside C
    vec3 cp = p - c;
    float d5 = dot(ab, cp);
    float d6 = dot(ac, cp);
    if (d6 >= 0.0 && d5 <= d6) return c;

    // Check if P is in edge region of AC
    float vb = d5 * d2 - d1 * d6;
    if (vb <= 0.0 && d2 >= 0.0 && d6 <= 0.0) {
        float w = d2 / (d2 - d6);
        return a + w * ac;
    }

    // Check if P is in edge region of BC
    float va = d3 * d6 - d5 * d4;
    if (va <= 0.0 && (d4 - d3) >= 0.0 && (d5 - d6) >= 0.0) {
        float w = (d4 - d3) / ((d4 - d3) + (d5 - d6));
        return b + w * (c - b);
    }

    // P is inside the face region; project onto the plane using barycentric coordinates
    float denom = 1.0 / (va + vb + vc);
    float v = vb * denom;
    float w = vc * denom;
    return a + ab * v + ac * w;
}

Breakdown of the Steps

  1. Vertex Tests: Dot products between the relative edge vectors and the relative point vectors check if the point lies outside the angle subtended by adjacent edges. If both dot products are non-positive at vertex \(A\), \(A\) is the closest point. The same logic applies to \(B\) and \(C\).
  2. Edge Tests: Quantities va, vb, and vc represent scaled normal components of the sub-triangles formed by the point and edge pairs. When a normal component is zero or negative and the point lies between the endpoints along that edge, the closest point lies directly on that segment. A linear interpolation parameter clamps the projection along that edge.
  3. Face Projection: If the point fails all vertex and edge region conditions, it must reside inside the triangle's prism. The unnormalized barycentric coordinates va, vb, and vc are normalized by their sum (denom), reconstructing the exact orthogonal projection on the triangle's surface.