Marschner Hair Shading: GLSL Cross-Fiber Scattering

The Marschner hair shading model evaluates cross-fiber scattering by decomposing light interaction into distinct azimuthal paths that occur across the hair fiber's circular or elliptical cross-section. In real-time rendering using GLSL, evaluating these cross-fiber paths requires projecting incident light and viewing vectors onto the normal plane perpendicular to the hair strand, calculating relative azimuthal angles, applying Bravais' law to determine effective refraction indices, and summing the contributions of surface reflection (\(R\)), transmission (\(TT\)), and internal reflection (\(TRT\)). This article details the mathematical formulation of cross-fiber scattering and its shader implementation in GLSL.

The Cross-Section Coordinate System

Hair geometry is defined per-fragment by a longitudinal tangent vector \(\mathbf{u}\) running along the fiber. Cross-fiber scattering occurs purely in the plane orthogonal to \(\mathbf{u}\).

To evaluate azimuthal scattering:

  1. The incident light vector \(\mathbf{\omega}_i\) and exitant view vector \(\mathbf{\omega}_r\) are projected onto the plane normal to \(\mathbf{u}\).
  2. The normalized projections \(\mathbf{\omega}_{i,\perp}\) and \(\mathbf{\omega}_{r,\perp}\) define the 2D azimuthal plane.
  3. The relative azimuthal angle \(\phi\) is calculated as the difference between the exitant azimuth \(\phi_r\) and the incident azimuth \(\phi_i\):

\[\phi = \phi_r - \phi_i\]

Because the light ray strikes the cylinder at an oblique angle \(\theta_d = (\theta_r - \theta_i) / 2\), the relative refractive index of the hair medium \(\eta\) is adjusted using Bravais' Law to determine the equivalent 2D index of refraction \(\eta'\):

\[\eta' = \frac{\sqrt{\eta^2 - \sin^2 \theta_d}}{\cos \theta_d}\]

Evaluating the Three Scattering Paths

Light entering the fiber is split into distinct paths denoted by \(p\):

1. The \(R\) Path (Direct Surface Reflection, \(p=0\))

Light reflects directly off the cuticle surface without penetrating the cortex. The deflection angle satisfies specular reflection off the outer cylinder:

  • It produces a uncolored, sharp highlight.
  • Evaluated azimuthally, the distribution \(N_R(\phi)\) peaks when the hair surface acts as a mirror normal to the half-angle. In GLSL, it is standard to approximate this with a wrapped Gaussian or cosine distribution centered at \(\phi = 0\).

2. The \(TT\) Path (Forward Transmission, \(p=1\))

Light refracts into the hair, traverses the cortex, and refracts out the back:

  • It generates a strong forward-scattering transmission colored by cortex pigments (eumelanin and pheomelanin).
  • The deflection angle centers around \(\phi = \pi\).
  • The ray attenuation follows the Beer-Lambert law based on the path length \(l\) through the cylinder: \(T = \exp(-\sigma_a l)\), where \(\sigma_a\) is the absorption coefficient.

3. The \(TRT\) Path (Internal Reflection and Caustics, \(p=2\))

Light refracts into the fiber, reflects off the back inner interface, and refracts out near the original entry side:

  • It produces a secondary, colored highlight shifted in position relative to the \(R\) path.
  • \(TRT\) scattering produces caustics due to stationary points in the deflection function:

\[\frac{d\Theta(p, \gamma)}{d\gamma} = 0\]

where \(\gamma\) is the impact parameter angle. These caustic peaks create bright, sharp glints along the fiber flanks.

GLSL Implementation

Evaluating the exact analytical roots of the \(TRT\) cubic deflection function per-pixel is computationally expensive. Real-time GLSL implementations compute \(N_p(\phi)\) using pre-integrated Gaussian approximations parameterized by roughness and azimuthal deflection angles.

// Computes cross-fiber azimuthal scattering (N components)
vec3 evaluateCrossFiberScattering(
    vec3 L,          // Vector toward light
    vec3 V,          // Vector toward camera
    vec3 T,          // Hair tangent
    vec3 absorption, // Hair pigment absorption coefficient
    float eta        // Hair IOR (typically 1.55)
) {
    // 1. Longitudinal angles to evaluate Bravais' index
    float sinThetaI = dot(L, T);
    float sinThetaR = dot(V, T);
    float cosThetaI = sqrt(max(0.0, 1.0 - sinThetaI * sinThetaI));
    float cosThetaR = sqrt(max(0.0, 1.0 - sinThetaR * sinThetaR));
    
    float thetaD = 0.5 * (asin(sinThetaR) - asin(sinThetaI));
    float cosThetaD = cos(thetaD);

    // Bravais effective refractive index
    float etaPrime = sqrt(eta * eta - (1.0 - cosThetaD * cosThetaD)) / cosThetaD;

    // 2. Azimuthal projections
    vec3 Lp = normalize(L - sinThetaI * T);
    vec3 Vp = normalize(V - sinThetaR * T);

    float cosPhi = clamp(dot(Lp, Vp), -1.0, 1.0);
    float phi = acos(cosPhi); // Relative azimuth in [0, PI]

    // 3. Fresnel reflectance at normal incidence
    float f0 = (eta - 1.0) / (eta + 1.0);
    f0 = f0 * f0;
    float F = f0 + (1.0 - f0) * pow(1.0 - max(0.0, dot(Lp, normalize(Lp + Vp))), 5.0);

    // 4. Azimuthal Path Distributions (Np approximations)
    // R Path: Focused backward reflection
    float sR = 0.15; // Azimuthal roughness parameter
    float NR = (1.0 / (sqrt(2.0 * 3.14159265) * sR)) * exp(-0.5 * (phi * phi) / (sR * sR));

    // TT Path: Centered around PI (forward scattering)
    float sTT = 0.25;
    float phiTT = phi - 3.14159265;
    float NTT = (1.0 / (sqrt(2.0 * 3.14159265) * sTT)) * exp(-0.5 * (phiTT * phiTT) / (sTT * sTT));

    // TRT Path: Caustic peaks shifted from 0
    float sTRT = 0.2;
    float phiTRT = phi - 0.75; // Caustic shift
    float NTRT = (1.0 / (sqrt(2.0 * 3.14159265) * sTRT)) * exp(-0.5 * (phiTRT * phiTRT) / (sTRT * sTRT));

    // 5. Internal absorption (Beer-Lambert)
    // Approximate mean cross-sectional path lengths
    float pathLengthTT = 2.0 * cos(asin(clamp(1.0 / etaPrime, -1.0, 1.0)));
    float pathLengthTRT = 4.0 * cos(asin(clamp(1.0 / etaPrime, -1.0, 1.0)));

    vec3 T_TT = exp(-absorption * pathLengthTT);
    vec3 T_TRT = exp(-absorption * pathLengthTRT);

    // 6. Combine paths with Fresnel attenuation
    vec3 colorR   = vec3(F * NR);
    vec3 colorTT  = (1.0 - F) * (1.0 - F) * T_TT * NTT;
    vec3 colorTRT = (1.0 - F) * (1.0 - F) * F * T_TRT * NTRT;

    return colorR + colorTT + colorTRT;
}

Energy Conservation and Normalization

When evaluating these paths in a complete shader, the azimuthal functions \(N_p(\phi)\) are multiplied by the longitudinal scattering functions \(M_p(\theta)\) such that the bidirectional scattering distribution function (BSDF) is:

\[S = S_R + S_{TT} + S_{TRT} = M_R N_R + M_{TT} N_{TT} + M_{TRT} N_{TRT}\]

To maintain energy balance, the integral of each \(N_p(\phi)\) over the interval \([-\pi, \pi]\) must evaluate to unity prior to multiplying by the internal transmission and Fresnel attenuation factors. When scale variations or rough cross-sections are simulated, the variance of the Gaussian distributions is widened dynamically without altering the integral sum.