Computing mesh normals on the GPU's other pipeline
While building MeshMaker I noticed I was constantly recomputing normals on the CPU and submitting them to the GPU. At the time I created the OpenCLPlayground to try OpenCL on the AMD Radeon in my iMac.
I never used that approach. The bandwidth for copying the results back was too slow and the whole thing was pretty complex. Now in 2026 I decided to revisit it using WebGPU, CUDA, and Slang.
Here it is, running in your browser.
That is a cloth, simulated and shaded entirely on the GPU. It moves every frame, so its normals are rebuilt from scratch each frame, and the slider runs that past a million triangles. Drag the resolution up and the GPU stays smooth. Flip the normals over to the CPU and watch the same work fall apart.
What a normal is, and why recompute it
A triangle’s normal is the direction it faces, the cross product of two of its edges. Smooth shading needs one per vertex instead, so a shared corner reads as rounded rather than faceted.
The cross product’s length is twice the triangle’s area, so summing the raw face normals already weights each one by its area, and a single normalize gives the vertex normal. That is the whole computation, the same handful of lines on a CPU core or ten thousand GPU lanes.
Here they are on a static shape.
The orange lines are the per-vertex normals, and smoothing them makes the faceted icosahedron read as a rounded ball. Its normals never change though, so you would compute them once and forget about them, which is exactly why a static mesh is the wrong demo for the GPU.
The algorithm, on the CPU first
Here is the reference version in plain TypeScript. It scatters each area-weighted face normal into its three vertices, then normalizes the accumulated vertex normals in a second pass.
export function computeVertexNormals(positions: Float32Array, indices: Uint32Array): Float32Array {
const normals = new Float32Array(positions.length); // zero-initialized
// Walk every triangle once. Its face normal is the cross product of two
// edges. We leave it un-normalized on purpose: the cross product's length is
// twice the triangle's area, so longer faces pull their shared vertices more.
for (let t = 0; t < indices.length; t += 3) {
const a = indices[t] * 3;
const b = indices[t + 1] * 3;
const c = indices[t + 2] * 3;
const e1x = positions[b] - positions[a];
const e1y = positions[b + 1] - positions[a + 1];
const e1z = positions[b + 2] - positions[a + 2];
const e2x = positions[c] - positions[a];
const e2y = positions[c + 1] - positions[a + 1];
const e2z = positions[c + 2] - positions[a + 2];
const nx = e1y * e2z - e1z * e2y;
const ny = e1z * e2x - e1x * e2z;
const nz = e1x * e2y - e1y * e2x;
// Scatter: add this face's normal to all three of its vertices.
normals[a] += nx; normals[a + 1] += ny; normals[a + 2] += nz;
normals[b] += nx; normals[b + 1] += ny; normals[b + 2] += nz;
normals[c] += nx; normals[c + 1] += ny; normals[c + 2] += nz;
}
// Normalize each accumulated vertex normal back to unit length.
for (let v = 0; v < normals.length; v += 3) {
const len = Math.hypot(normals[v], normals[v + 1], normals[v + 2]);
if (len > 0) {
normals[v] /= len;
normals[v + 1] /= len;
normals[v + 2] /= len;
}
}
return normals;
}A small harness imports that exact function and the same mesh generator the demos use, runs them on a bare icosahedron, and checks the result. The output below is captured at build time.
// The bare icosahedron: 12 vertices, 20 faces, every vertex on the unit sphere.
const mesh = icosphere(0);
const vertexCount = mesh.positions.length / 3;
const triangleCount = mesh.indices.length / 3;
console.log(`mesh: ${vertexCount} vertices, ${triangleCount} triangles`);
const normals = computeVertexNormals(mesh.positions, mesh.indices);
// Every vertex normal should come out unit length.
let maxLenError = 0;
for (let v = 0; v < normals.length; v += 3) {
const len = Math.hypot(normals[v], normals[v + 1], normals[v + 2]);
maxLenError = Math.max(maxLenError, Math.abs(len - 1));
}
console.log(`all normals unit length: ${maxLenError < 1e-6} (max error ${maxLenError.toExponential(1)})`);
// On a sphere sampled symmetrically, each vertex normal points the same way as
// the vertex itself. So the computed normal should match the normalized
// position. This is the invariant that tells us the averaging is right.
let maxDirError = 0;
for (let v = 0; v < normals.length; v += 3) {
const len = Math.hypot(mesh.positions[v], mesh.positions[v + 1], mesh.positions[v + 2]);
for (let k = 0; k < 3; k++) {
maxDirError = Math.max(maxDirError, Math.abs(normals[v + k] - mesh.positions[v + k] / len));
}
}
console.log(`normals match normalized position: ${maxDirError < 1e-6} (max error ${maxDirError.toExponential(1)})`);
// Show the first three, rounded so the output is stable across machines.
const round = (x: number) => (Math.abs(x) < 1e-7 ? 0 : Number(x.toFixed(4)));
for (let v = 0; v < 3; v++) {
const i = v * 3;
console.log(` normal[${v}] = (${round(normals[i])}, ${round(normals[i + 1])}, ${round(normals[i + 2])})`);
}mesh: 12 vertices, 20 triangles
all normals unit length: true (max error 3.1e-8)
normals match normalized position: true (max error 2.7e-8)
normal[0] = (-0.5257, 0.8507, 0)
normal[1] = (0.5257, 0.8507, 0)
normal[2] = (-0.5257, -0.8507, 0)
The same thing on the GPU
The CPU version scatters: it loops over triangles and adds into vertices. On the GPU that is a trap. Thousands of triangle threads would write the same vertex at once, and that is a data race. The obvious fix is an atomic add, except WGSL has no atomics on floats, so you would have to integer-encode each component and decode it in a second pass.
So I flipped the loop. Instead of one thread per triangle scattering outward, the GPU runs one thread per vertex that gathers inward. Each thread walks only the triangles touching its own vertex and writes its own slot exactly once. No contention, no atomics. It needs a small vertex-to-triangle table built once on the CPU, which is cheap.
The bindings declare the buffers the kernel reads and the one it writes.
struct Counts {
vertexCount : u32,
triangleCount : u32,
};
@group(0) @binding(0) var<storage, read> positions : array<f32>;
@group(0) @binding(1) var<storage, read> indices : array<u32>;
@group(0) @binding(2) var<storage, read> adjStart : array<u32>;
@group(0) @binding(3) var<storage, read> adjTris : array<u32>;
@group(0) @binding(4) var<storage, read_write> normals : array<f32>;
@group(0) @binding(5) var<uniform> counts : Counts;And the kernel itself reads almost line for line like the CPU loop, just inside out.
@compute @workgroup_size(64)
fn main(@builtin(global_invocation_id) gid : vec3<u32>) {
let v = gid.x;
if (v >= counts.vertexCount) { return; }
// Walk only the triangles incident to this vertex, summing their
// area-weighted normals. One thread, one output slot, no atomics.
var n = vec3<f32>(0.0, 0.0, 0.0);
for (var i = adjStart[v]; i < adjStart[v + 1u]; i = i + 1u) {
n = n + faceNormal(adjTris[i]);
}
if (dot(n, n) > 0.0) { n = normalize(n); }
normals[3u * v] = n.x;
normals[3u * v + 1u] = n.y;
normals[3u * v + 2u] = n.z;
}Why the cloth, and not the icosahedron
A cloth is the cleanest example I know. It is a grid of particles pulled on by gravity and wind, held together by distance constraints, and it never stops moving. Because it moves, the normals have to be rebuilt every frame instead of computed once, so now the cost is real and recurring. The whole thing runs on the GPU as a sequence of compute passes.
First a Verlet integration step moves every free particle.
@compute @workgroup_size(64)
fn integrate(@builtin(global_invocation_id) gid : vec3<u32>) {
let i = gid.x;
if (i >= sim.vertexCount) { return; }
if (pinned(i)) { return; }
let p = getPos(i);
let prev = vec3<f32>(prevPos[3u * i], prevPos[3u * i + 1u], prevPos[3u * i + 2u]);
let dt = sim.gravity.w;
// Gentle out-of-plane gust that varies across the sheet, so a firmly hanging
// cloth ripples into travelling curtain folds instead of staying flat. Kept
// positive and weaker than gravity, so the cloth billows but never blows away.
let t = sim.wind.w;
let gust = 0.5 + 0.5 * sin(p.x * 6.0 + t * 2.0) * cos(p.y * 2.0 + t * 0.7);
let accel = sim.gravity.xyz + sim.wind.xyz * gust;
// Verlet integration: velocity is implied by (p - prev).
var next = p + (p - prev) * sim.rest.w + accel * dt * dt;
// Push out of the interactive sphere.
if (sim.sphere.w > 0.0) {
let d = next - sim.sphere.xyz;
let dist = length(d);
if (dist < sim.sphere.w) {
next = sim.sphere.xyz + d / max(dist, 1e-5) * sim.sphere.w;
}
}
prevPos[3u * i] = p.x; prevPos[3u * i + 1u] = p.y; prevPos[3u * i + 2u] = p.z;
setPos(i, next);
}Then the distance constraints pull neighbours back toward their rest length. Each one is a small spring between two particles. The spring has a rest length. When the two particles move too far apart or too close, it slides both of them back along the line between them until the gap returns to that rest length. This is the distance constraint from Position Based Dynamics.
In one step the spring shares the correction between its two ends.
Here is the current distance, is the rest length, and points along the spring.
Solving one spring is easy. Solving all of them at once on the GPU brings back the write race from the normals. Two springs that share a particle would move it at the same time. The fix is to run the springs in batches. No two springs in a batch share a particle, so a whole batch can run at once. The batches run one after another. On a regular grid the split is just even and odd, eight batches in all. The trick has a name, graph coloring.
@compute @workgroup_size(64)
fn solve(@builtin(global_invocation_id) gid : vec3<u32>) {
let i = gid.x;
if (i >= sim.vertexCount) { return; }
let c = i % sim.cols;
let r = i / sim.cols;
let diag = sqrt(sim.rest.x * sim.rest.x + sim.rest.y * sim.rest.y);
if (AXIS == 0u) { // structural, horizontal
if (c % 2u == COLOR && c + 1u < sim.cols) { solvePair(i, i + 1u, sim.rest.x); }
} else if (AXIS == 1u) { // structural, vertical
if (r % 2u == COLOR && r + 1u < sim.rows) { solvePair(i, i + sim.cols, sim.rest.y); }
} else if (AXIS == 2u) { // shear, down-right diagonal
if (r % 2u == COLOR && c + 1u < sim.cols && r + 1u < sim.rows) {
solvePair(i, i + sim.cols + 1u, diag);
}
} else { // shear, down-left diagonal
if (r % 2u == COLOR && c > 0u && r + 1u < sim.rows) {
solvePair(i, i + sim.cols - 1u, diag);
}
}
}One catch shows up only at scale. Distance constraints pass tension between neighbours one step at a time, so on a small sheet a handful of solver passes reach all the way to the pins, but on a million-particle sheet they fall far short and the lower half sags away as if it were heavier. The fix is a long-range attachment: cap each particle’s distance to its nearest pin directly, using the taut-path length to that pin. That holds the whole sheet up in a single pass, no matter how many particles there are, so the cloth behaves the same at every resolution.
// Long-range attachment: cap each particle's distance to its nearest pin. The
// pin position is fixed, so this needs no propagation and holds the whole sheet
// up in a single pass, which is what keeps the drape the same at every
// resolution instead of letting a high-res sheet sag away.
@compute @workgroup_size(64)
fn applyLRA(@builtin(global_invocation_id) gid : vec3<u32>) {
let v = gid.x;
if (v >= sim.vertexCount || pinned(v)) { return; }
let pinPos = vec3<f32>(lra[4u * v], lra[4u * v + 1u], lra[4u * v + 2u]);
let maxDist = lra[4u * v + 3u];
let p = getPos(v);
let d = p - pinPos;
let len = length(d);
if (len > maxDist) { setPos(v, pinPos + d / len * maxDist); }
}And then the normals, recomputed from the freshly moved positions with the very same gather kernel from above.
@compute @workgroup_size(64)
fn computeNormals(@builtin(global_invocation_id) gid : vec3<u32>) {
let v = gid.x;
if (v >= sim.vertexCount) { return; }
var n = vec3<f32>(0.0, 0.0, 0.0);
for (var i = adjStart[v]; i < adjStart[v + 1u]; i = i + 1u) {
n = n + faceNormal(adjTris[i]);
}
if (dot(n, n) > 0.0) { n = normalize(n); }
normals[3u * v] = n.x;
normals[3u * v + 1u] = n.y;
normals[3u * v + 2u] = n.z;
}One buffer, two pipelines
This is the part I was chasing back in the MeshMaker days.
When the compute passes finish, the new positions and normals are already sitting in GPU buffers. The render pipeline takes those same buffers as its vertex inputs and draws them in place. There is no map, no copy, no trip through CPU memory between simulating the cloth and drawing it.
@fragment
fn fs(in : VertexOut, @builtin(front_facing) front : bool) -> @location(0) vec4<f32> {
var N = normalize(in.worldNormal);
if (!front) { N = -N; } // light the back face too
let L = normalize(camera.lightDir.xyz);
let diffuse = max(dot(N, L), 0.0);
let ambient = 0.24;
// Cool front, warm back, so a fold that turns the cloth over is obvious.
let base = select(vec3<f32>(0.85, 0.42, 0.28), vec3<f32>(0.34, 0.52, 0.86), front);
return vec4<f32>(base * (ambient + diffuse * 0.9), 1.0);
}The buffers carry two usages at once, storage for the compute passes and vertex data for the render pass, and the GPU just reads them where they lie. That hand-off is the whole point of the general-purpose pipeline: you compute the data where it will be used, instead of carrying it back to the CPU.
The bandwidth you do not spend
Switch the demo’s normals from the GPU to the CPU and the frame rate collapses, especially as you raise the resolution. The simulation still runs on the GPU, but now every frame has to copy the positions back to the CPU, recompute the normals in JavaScript, and upload them again. The reading on the demo shows the bytes moving and the milliseconds spent.
On a discrete card those buffers live across the PCIe bus, so pulling them back to the CPU is a real transfer, plus the latency of waiting for the queue to drain. At a million triangles that is megabytes a frame, in each direction. Keeping the result on the GPU avoids the whole thing.
The same kernel, four ways
WebGPU is the one option that runs live on this page, so it got the demo. Outside the browser the field is crowded.
The kernel itself barely changes across them. Here is that per-vertex gather in CUDA, then OpenCL (the dialect MeshMaker’s original playground used), then Slang.
__global__ void vertexNormals(const float* pos, const unsigned* indices,
const unsigned* adjStart, const unsigned* adjTris,
float* normals, unsigned vertexCount) {
unsigned v = blockIdx.x * blockDim.x + threadIdx.x;
if (v >= vertexCount) return;
float nx = 0.f, ny = 0.f, nz = 0.f;
for (unsigned i = adjStart[v]; i < adjStart[v + 1]; ++i) {
unsigned t = adjTris[i];
unsigned a = indices[3 * t], b = indices[3 * t + 1], c = indices[3 * t + 2];
float e1x = pos[3*b] - pos[3*a], e1y = pos[3*b+1] - pos[3*a+1], e1z = pos[3*b+2] - pos[3*a+2];
float e2x = pos[3*c] - pos[3*a], e2y = pos[3*c+1] - pos[3*a+1], e2z = pos[3*c+2] - pos[3*a+2];
nx += e1y * e2z - e1z * e2y;
ny += e1z * e2x - e1x * e2z;
nz += e1x * e2y - e1y * e2x;
}
float len = sqrtf(nx*nx + ny*ny + nz*nz);
if (len > 0.f) { nx /= len; ny /= len; nz /= len; }
normals[3*v] = nx; normals[3*v+1] = ny; normals[3*v+2] = nz;
}__kernel void vertexNormals(__global const float* pos, __global const uint* indices,
__global const uint* adjStart, __global const uint* adjTris,
__global float* normals, uint vertexCount) {
uint v = get_global_id(0);
if (v >= vertexCount) return;
float3 n = (float3)(0.0f, 0.0f, 0.0f);
for (uint i = adjStart[v]; i < adjStart[v + 1]; ++i) {
uint t = adjTris[i];
uint a = indices[3*t], b = indices[3*t+1], c = indices[3*t+2];
float3 pa = (float3)(pos[3*a], pos[3*a+1], pos[3*a+2]);
float3 pb = (float3)(pos[3*b], pos[3*b+1], pos[3*b+2]);
float3 pc = (float3)(pos[3*c], pos[3*c+1], pos[3*c+2]);
n += cross(pb - pa, pc - pa);
}
float len = length(n);
if (len > 0.0f) n /= len;
normals[3*v] = n.x; normals[3*v+1] = n.y; normals[3*v+2] = n.z;
}RWStructuredBuffer<float> positions;
RWStructuredBuffer<uint> indices;
RWStructuredBuffer<uint> adjStart;
RWStructuredBuffer<uint> adjTris;
RWStructuredBuffer<float> normals;
uniform uint vertexCount;
float3 getPos(uint i) {
return float3(positions[3 * i], positions[3 * i + 1], positions[3 * i + 2]);
}
[shader("compute")]
[numthreads(64, 1, 1)]
void computeMain(uint3 tid : SV_DispatchThreadID) {
uint v = tid.x;
if (v >= vertexCount) return;
float3 n = float3(0.0);
for (uint i = adjStart[v]; i < adjStart[v + 1]; ++i) {
uint t = adjTris[i];
float3 a = getPos(indices[3 * t]);
float3 b = getPos(indices[3 * t + 1]);
float3 c = getPos(indices[3 * t + 2]);
n += cross(b - a, c - a);
}
if (dot(n, n) > 0.0) n = normalize(n);
normals[3 * v] = n.x;
normals[3 * v + 1] = n.y;
normals[3 * v + 2] = n.z;
}Set them next to the WGSL above and they are the same idea four times. Loop over the incident triangles, cross two edges, accumulate, normalize. What differs is the host code around the kernel, how you allocate and bind buffers, and whether the result can be handed to a renderer without a copy.
Slang is the interesting one. It is a language rather than a runtime: it moved from NVIDIA to the Khronos Group in late 2024, and one Slang source compiles to SPIR-V for Vulkan, MSL for Metal, HLSL for Direct3D, CUDA, and PTX. The benchmark below compiles it to PTX and launches it through the CUDA driver.
So I ran all of them, on the same mesh, on my own machine, and let them race. Each one compiles and runs the kernel two hundred times and reports its best. They all produce the same normals, which the harness checks with a shared checksum, so the only thing that varies is speed.
| Kernel | Ran on | ms / frame | Mverts/s | vs CPU |
|---|---|---|---|---|
| JavaScript | CPU (single thread, JS) | 13.12 | 20 | — |
| WGSL (WebGPU) | nvidia blackwell* | 0.023 | 11,398 | 570× |
| CUDA | NVIDIA GeForce RTX 5070 Ti | 0.024 | 10,923 | 547× |
| OpenCL | NVIDIA GeForce RTX 5070 Ti | 0.021 | 12,483 | 625× |
| Slang | NVIDIA GeForce RTX 5070 Ti | 0.021 | 12,483 | 625× |
npm run bench (captured 2026-05-30).
* WebGPU does not expose the specific device model, only the vendor and architecture, as a fingerprinting safeguard. It is the same physical GPU the native rows name in full.
The numbers are what you would expect and still fun to see. A single CPU thread is the slowest by a wide margin. The four GPU paths land within a hair of each other, because they are all compiling nearly the same kernel down to the same hardware. WGSL, the one that also runs in the browser demo above, is right there with native CUDA.
The whole table is reproducible. npm run bench generates the shared mesh, builds each kernel with whatever toolchains are installed, runs them, and writes the numbers this table reads. Anything it cannot build, it skips out loud rather than quietly dropping.
Then and now
The OpenCLPlayground in MeshMaker was a side experiment that never shipped. Setting up an OpenCL context, sharing a buffer with OpenGL, and keeping the two pipelines in step was enough friction that it stayed a playground.
The friction is mostly gone now. WebGPU puts the compute pipeline, the render pipeline, and the shared buffer behind one API, in the browser, with no driver to install. The thing I wanted years ago, computing mesh data on the general-purpose side and handing it straight to the renderer, is now a cloth you can throw a million triangles at from your phone.