GLSL Tangent Space from Screen Derivatives

This article details the mathematical formulation and GLSL implementation for dynamically reconstructing tangent-space vectors (the tangent and bitangent) within fragment shaders. By evaluating pixel-quad screen derivatives via standard GLSL functions, developers can perform tangent-space normal mapping on arbitrary meshes without storing or interpolating explicit tangent attributes per vertex.

The Mathematical Foundation

Surface mapping assumes a parameterized two-dimensional manifold where surface positions \(P \in \mathbb{R}^3\) correspond to texture coordinates \((u, v) \in \mathbb{R}^2\). The tangent vector \(T\) and bitangent vector \(B\) represent the directional derivatives of position along the principal texture axes:

\[T = \frac{\partial P}{\partial u}, \quad B = \frac{\partial P}{\partial v}\]

A differential change in surface position \(dP\) relative to texture coordinates is represented by the total derivative:

\[dP = \frac{\partial P}{\partial u} du + \frac{\partial P}{\partial v} dv = T \, du + B \, dv\]

On rasterized geometry, fragment shaders do not have direct access to continuous analytical derivatives with respect to \(u\) and \(v\). Instead, modern hardware evaluates screen-space derivatives along pixel axes (\(x\) and \(y\)) across \(2 \times 2\) fragment quads:

\[\frac{\partial P}{\partial x} = T \frac{\partial u}{\partial x} + B \frac{\partial v}{\partial x}\]

\[\frac{\partial P}{\partial y} = T \frac{\partial u}{\partial y} + B \frac{\partial v}{\partial y}\]

Solving the Linear System

Let:

• \(dp_1 = \frac{\partial P}{\partial x}\) (dFdx(P))
• \(dp_2 = \frac{\partial P}{\partial y}\) (dFdy(P))
• \(du_1 = \frac{\partial u}{\partial x}\), \(dv_1 = \frac{\partial v}{\partial x}\) (dFdx(uv))
• \(du_2 = \frac{\partial u}{\partial y}\), \(dv_2 = \frac{\partial v}{\partial y}\) (dFdy(uv))

In matrix notation:

\[\begin{bmatrix} dp_1 & dp_2 \end{bmatrix} = \begin{bmatrix} T & B \end{bmatrix} \begin{bmatrix} du_1 & du_2 \\ dv_1 & dv_2 \end{bmatrix}\]

Multiplying both sides by the inverse of the \(2 \times 2\) Jacobian matrix of UV coordinates yields the tangent and bitangent vectors:

\[\begin{bmatrix} T & B \end{bmatrix} = \begin{bmatrix} dp_1 & dp_2 \end{bmatrix} \begin{bmatrix} du_1 & du_2 \\ dv_1 & dv_2 \end{bmatrix}^{-1}\]

The determinant \(D\) of the UV coordinate matrix is:

\[D = du_1 dv_2 - du_2 dv_1\]

When \(D \neq 0\), the matrix inverse is:

\[\begin{bmatrix} du_1 & du_2 \\ dv_1 & dv_2 \end{bmatrix}^{-1} = \frac{1}{D} \begin{bmatrix} dv_2 & -du_2 \\ -dv_1 & du_1 \end{bmatrix}\]

Multiplying out the terms produces the explicit equations for the unnormalized tangent and bitangent vectors:

\[T = \frac{dp_1 dv_2 - dp_2 dv_1}{D}\]

\[B = \frac{dp_2 du_1 - dp_1 du_2}{D}\]

Because the vectors will be normalized during shading, the scalar division by \(D\) can be ignored, preserving vector direction while requiring only a sign correction to account for UV parameterization mirroring:

\[T = (dp_1 dv_2 - dp_2 dv_1) \cdot \text{sgn}(D)\]

\[B = (dp_2 du_1 - dp_1 du_2) \cdot \text{sgn}(D)\]

Orthonormalization with the Geometric Normal

Due to perspective projection distortion, non-uniform scaling, or non-orthogonal UV mappings, the derived \(T\) and \(B\) vectors may not be strictly orthogonal to the interpolated surface normal \(N\). The Gram-Schmidt process resolves this by projecting \(T\) onto the plane orthogonal to \(N\):

\[T_{\perp} = \text{normalize}(T - N (N \cdot T))\]

The final bitangent can then either be evaluated using the same process from \(B\) to accommodate non-orthogonal UV coordinates, or computed via a cross-product if an orthogonal tangent frame is assumed:

\[B_{\perp} = \text{cross}(N, T_{\perp}) \cdot \text{sgn}(D)\]

GLSL Implementation

The following function implements this derivation to compute a \(3 \times 3\) tangent-to-world (TBN) matrix:

mat3 cotangent_frame(vec3 N, vec3 p, vec2 uv) {
    // Evaluate screen-space partial derivatives
    vec3 dp1 = dFdx(p);
    vec3 dp2 = dFdy(p);
    vec2 duv1 = dFdx(uv);
    vec2 duv2 = dFdy(uv);

    // Solve the linear system
    vec3 dp2perp = cross(dp2, N);
    vec3 dp1perp = cross(N, dp1);
    
    vec3 T = dp2perp * duv1.x + dp1perp * duv2.x;
    vec3 B = dp2perp * duv1.y + dp1perp * duv2.y;

    // Construct an orthonormal basis
    float invmax = inversesqrt(max(dot(T, T), dot(B, B)));
    return mat3(T * invmax, B * invmax, N);
}

This formulation projects the derivatives directly using the geometric normal via vector cross products, avoiding explicit matrix determinant inversion while producing an identical orthonormal basis. Normal maps can then be transformed directly to world space by multiplying the sampled normal vector with the resulting matrix.