Calculating SDF Normals with Central Differences in GLSL

In raymarching, Signed Distance Fields (SDFs) mathematically describe 3D volumes rather than explicit polygonal geometry. Because there are no predefined polygonal mesh faces or vertex attributes, lighting and shading require dynamically deriving the surface normal at the point of ray intersection. This article explains how the surface normal corresponds to the gradient of the distance field, how the central differences method numerically approximates this gradient in GLSL, and how to implement the calculation efficiently in a shader.

Understanding the SDF Gradient

A Signed Distance Function, \(f(p)\), takes a point \(p = (x, y, z)\) in 3D space and returns the shortest Euclidean distance to the surface of a geometry. On the surface, the value of \(f(p)\) is zero.

The surface normal \(\vec{n}\) at any point on the boundary of an SDF is defined by the spatial gradient vector:

\[\nabla f(p) = \left( \frac{\partial f}{\partial x}, \frac{\partial f}{\partial y}, \frac{\partial f}{\partial z} \right)\]

Because an SDF represents distance, moving directly away from the surface along its perpendicular direction yields the steepest increase in distance. Therefore, computing the partial derivatives along each axis and normalizing the resulting vector provides the unit surface normal:

\[\vec{n} = \frac{\nabla f(p)}{\|\nabla f(p)\|}\]

The Central Differences Method

While analytical derivation is possible for simple primitives, complex scenes combining multiple primitives via boolean operations, smooth blending, or domain repetition make symbolic differentiation impractical. Instead, finite difference approximations are used.

The central differences technique evaluates the function at a small offset \(\epsilon\) both ahead of and behind the evaluation point along each axis:

\[\frac{\partial f}{\partial x} \approx \frac{f(x + \epsilon, y, z) - f(x - \epsilon, y, z)}{2\epsilon}\]

\[\frac{\partial f}{\partial y} \approx \frac{f(x, y + \epsilon, z) - f(x, y - \epsilon, z)}{2\epsilon}\]

\[\frac{\partial f}{\partial z} \approx \frac{f(x, y, z + \epsilon) - f(x, y, z - \epsilon)}{2\epsilon}\]

Central differences are preferred over forward differences (\([f(p + \epsilon) - f(p)] / \epsilon\)) because the truncation error is second-order, \(O(\epsilon^2)\), rather than first-order, \(O(\epsilon)\). This symmetry yields significantly higher accuracy and smoother normals across sharp transitions. Because the resulting gradient vector is normalized, the common scaling factor of \(2\epsilon\) in the denominator cancels out and can be omitted from the shader calculation.

GLSL Implementation

The standard central differences implementation requires six evaluations of the scene's distance function (two taps along each of the three cardinal axes).

// Function representing the scene's signed distance field
float map(vec3 p);

vec3 calcNormal(vec3 p) {
    const float eps = 0.001; // Small step size
    const vec2 h = vec2(eps, 0.0);

    return normalize(vec3(
        map(p + h.xyy) - map(p - h.xyy),
        map(p + h.yxy) - map(p - h.yxy),
        map(p + h.yyx) - map(p - h.yyx)
    ));
}

Choosing the Epsilon Value

The step size eps directly affects the visual quality of the normals:

• Too Large: High-frequency geometric details, sharp edges, and crevice features will appear overly smoothed or beveled.
• Too Small: 32-bit floating-point precision limits will introduce severe noise and banding artifacts into the normal estimation.

An eps value between 0.0005 and 0.002 works well for scenes modeled within unit-scale coordinates. For scenes viewed across vast distances, scaling eps proportionally with the camera distance to the intersection point helps mitigate aliasing and floating-point errors.