Objective

Combining a base geometry, a distortion function, and a shading material within a WebGL context enables the generation of highly diverse, interactive visual outputs. The goal is to apply spatial distortions to 3D shapes dynamically. A key constraint is compatibility with p5.js, which lacks an API for passing custom vertex attributes to shaders. Therefore, the implementation must rely solely on standard vertex data: position, normal, texture coordinates, and transformation matrices.

The Normal Vector Problem

Applying arbitrary distortions to a mesh invalidates its pre-calculated surface normals, resulting in incorrect lighting and reflections. In a standard rendering pipeline, the vertex shader passes normals to the fragment shader, where they are interpolated across triangles to determine light interaction. While rotating a mesh simply applies the same rotation to its normals, deforming a mesh (e.g., twisting a plane with a sine wave) alters the surface curvature, meaning the original normals no longer correspond to the modified geometry.

Ineffective Approaches

Fragment Shader Derivatives

Ignoring vertex normals and computing them entirely in the fragment shader using screen-space derivatives is a tempting workaround. GLSL provides dFdx and dFdy to determine the rate of change of a variable across adjacent pixels. The cross product of these two tangent vectors yields a surface normal.

varying vec3 vWorldPos;
uniform vec3 uTint;

void main() {
  vec3 computedNorm = normalize(cross(dFdx(vWorldPos), dFdy(vWorldPos)));
  vec3 lightDir = vec3(0.0, 1.0, 0.0);
  float diffuse = max(dot(computedNorm, lightDir), 0.0);
  gl_FragColor = vec4(diffuse * uTint, 1.0);
}

However, this approach produces faceted shading. Because the mesh consists of discrete triangles, the screen-space derivatives are constant across a single triangle, preventing the smooth interpolation typically provided by varying vertex normals.

Manual Vertex Shader Derivatives

Computing the analytical derivatives of the distortion function directly in the vertex shader allows for smooth normal interpolation. Unfortunately, manually writing the mathematical derivatives for complex distortion functions is highly error-prone and tedious.

Core Solution: Automated Shader Code Generation

Since calculating derivatives follows strict algorithmic rules (like the chain rule), a computer can perform it reliably. By constructing an abstract syntax tree of the mathematical operations, it is possible to automatically generate both the evaluation code and its derivative code.

Automatic Differentiation Architecture

Each mathematical operation is modeled as an object capable of recursively generating its own GLSL code and derivative code. For instance, a sine operation applies the chain rule by multiplying the cosine of its argument by the argument's derivative.

class SineNode {
  constructor(arg) {
    this.arg = arg;
  }

  generateExpr() {
    return `sin(${this.arg.generateExpr()})`;
  }

  generateDeriv() {
    return `cos(${this.arg.generateExpr()}) * ${this.arg.generateDeriv()}`;
  }
}

By utilizing a builder pattern to construct these operation trees, the generated GLSL can be seamlessly injected into the vertex shader. This technique allows for complex procedural displacements with automatically computed partial derivatives.

import { buildAutoDiff } from 'custom-glsl-derivative';

const vertexSrc = `
  void main() {
    vec4 pos = vec4(aVertexPos, 1.0);
    float px = pos.x;
    float py = pos.y;
    ${buildAutoDiff((ctx) => {
      const px = ctx.var('px');
      const py = ctx.var('py');
      const tick = ctx.var('tick');
      let wave = ctx.const(0);
      for (let i = 0; i < 3; i++) {
        wave = wave.add(ctx.sin(
          ctx.add(wave.mult(0.5), px.mult(1.5), py.mult(2.8), tick.mult(0.002))
        ));
      }
      wave = wave.mult(0.1);
      wave.assign('zOffset');
      wave.derivative('dzx', px);
      wave.derivative('dzy', py);
    })}
    pos.z = zOffset;
    vec3 tx = vec3(1.0, 0.0, dzx);
    vec3 ty = vec3(0.0, 1.0, dzy);
    vTransformedNormal = uNormMat * normalize(cross(tx, ty));

    gl_Position = uProjMat * uModelViewMat * pos;
  }
`;

Applying Derivatives to Normal Calculation

Deforming a Plane

For a flat plane displaced along the Z-axis, the surface tangents are easily derived from the partial derivatives of the offset with respect to X and Y. The new normal is the normalized cross product of these two tangent vectors.

vec3 deformedPos = aPosition;
deformedPos.z += heightOffset;

vec3 tanX = vec3(1.0, 0.0, dOffset_dx);
vec3 tanY = vec3(0.0, 1.0, dOffset_dy);
vec3 recalculatedNorm = normalize(cross(tanX, tanY));

Deforming Arbitrary Meshes

Applying the same Z-displacement to a non-planar shape requires preserving the original curvature. This is achieved by calculating the rotation induced by the displacement on a hypothetical plane, and then applying that exact rotation to the mesh's original normal vector. The rotation axis is the cross product of the original plane normal and the newly computed normal, and the rotation angle is the arc cosine of their dot product.

mat4 computeRotationMatrix(vec3 rotAxis, float rotAngle) { /* ... */ }

vec3 tanX = vec3(1.0, 0.0, dOffset_dx);
vec3 tanY = vec3(0.0, 1.0, dOffset_dy);
vec3 warpedFlatNorm = normalize(cross(tanX, tanY));
vec3 flatBaseNorm = vec3(0.0, 0.0, 1.0);

float rotAngle = acos(dot(flatBaseNorm, warpedFlatNorm));
vec3 rotAxis = normalize(cross(flatBaseNorm, warpedFlatNorm));
mat4 rotMat = computeRotationMatrix(rotAxis, rotAngle);

vec3 finalNorm = (rotMat * vec4(aOriginalNorm, 0.0)).xyz;

3D Vector Displacement

When the distortion outputs a full 3D vector instead of a scalar Z-offset, partial derivatives must account for all three axes. By defining a composite direction variable, such as u = y + z, the Jacobian of the displacement function can be used to find the directional derivatives. The tangent vectors are formed by adding the respective partial derivatives to the basis vectors, and the same rotation matrix logic applies to transform the original normals correctly.

vec3 tanX = vec3(1.0, 0.0, 0.0) + dVec_dx;
vec3 tanYZ = vec3(0.0, 1.0, 1.0) + dVec_dy + dVec_dz;
vec3 warpedFlatNorm = normalize(cross(tanX, tanYZ));
vec3 flatBaseNorm = normalize(cross(vec3(1.0, 0.0, 0.0), vec3(0.0, 1.0, 1.0)));

float rotAngle = acos(dot(flatBaseNorm, warpedFlatNorm));
vec3 rotAxis = normalize(cross(flatBaseNorm, warpedFlatNorm));
mat4 rotMat = computeRotationMatrix(rotAxis, rotAngle);

vec3 finalNorm = (rotMat * vec4(aOriginalNorm, 0.0)).xyz;