fix(cloth-icm): keep drop cloth free of intersections

The drop scene used the tangled-start settings, so layers stayed
intersecting. It now has its own grid (35x35, gap 0.035) and
re-checks proximity after 5 of 10 substeps.
After landing, ICM on drops from ~260 to ~1 hit per frame.

ROADMAP.md

@@ -19,6 +19,7 @@ Open ideas for the playground frontend. Tick a box when done.
 
 ## Features
 
+- [ ] **Cloth ICM drop scene**: hits mostly solved, cloth settling and CCD open. With ICM on, drop scene 35x35 (spacing 0.145, same 4.9 m cloth), gap 0.035, one mid-frame proximity re-detection: ~0 intersections after landing (avg 1.5/frame and 82 % of frames at 0 after frame 150; 0 at the marks 100-450 except 11 at 600), 9-10 ms/frame headless; was 263 at frame 900 on 50x50. ICM off: avg 148 (max 356) after frame 150 (125 at frames 100-240), 68 avg over 1800 frames, bounded but reaches the 3 m/s speed cap at times. Finer grids or a gap >= 0.04 at this size make ICM off blow up (thousands of hits, 60-100 ms frames). Cloth still does not settle (fold flutter beside the sphere, max speed 0.6-1.0 m/s with ICM on, also with collisions off). Next: Bridson 2002 CCD + velocity repulsion, then damping of the fold flutter.
 - [ ] **Shareable demo state**: keep demo settings in query params (fractal zoom/center, sorting algorithm, seeds). Opening the link restores the view.
 - [ ] **New algorithm demos** (pattern: route + `RouterConstants` + `AlgorithmsService` entry + i18n + information panel):
   - [x] Boids / flocking (GPU path like cloth)
@@ -30,4 +31,5 @@ Open ideas for the playground frontend. Tick a box when done.
   - [x] Reaction-diffusion (Gray-Scott) shader (reuse the fractal shader setup)
   - [x] **Lenia** (continuous Game of Life, shader ping-pong like Gray-Scott)
   - [x] **Fluid (Stable Fluids)** (velocity, pressure and dye shader ping-pong like Lenia)
+  - [x] **Cloth self-collision (ICM)** (Volino and Magnenat-Thalmann 2006, CPU: BVH with normal-cone culling, contour tracing, intersection contour minimization)
   - [x] **Ant colony** (ACO on TSP, CPU + Canvas 2D)

src/app/app.routes.ts

@@ -16,6 +16,7 @@ export const routes: Routes = [
   { path: RouterConstants.FRACTAL3d.PATH, loadComponent: () => import('./pages/algorithms/fractal3d/fractal3d.component').then(m => m.Fractal3dComponent), data: { seo: 'ALGORITHMS' } },
   { path: RouterConstants.PENDULUM.PATH, loadComponent: () => import('./pages/algorithms/pendulum/pendulum.component').then(m => m.default), data: { seo: 'ALGORITHMS' } },
   { path: RouterConstants.CLOTH.PATH, loadComponent: () => import('./pages/algorithms/cloth/cloth.component').then(m => m.ClothComponent), data: { seo: 'ALGORITHMS' } },
+  { path: RouterConstants.CLOTH_ICM.PATH, loadComponent: () => import('./pages/algorithms/cloth-icm/cloth-icm.component').then(m => m.ClothIcmComponent), data: { seo: 'ALGORITHMS' } },
   { path: RouterConstants.BOIDS.PATH, loadComponent: () => import('./pages/algorithms/boids/boids.component').then(m => m.BoidsComponent), data: { seo: 'ALGORITHMS' } },
   { path: RouterConstants.PHYSARUM.PATH, loadComponent: () => import('./pages/algorithms/physarum/physarum.component').then(m => m.PhysarumComponent), data: { seo: 'ALGORITHMS' } },
   { path: RouterConstants.GRAY_SCOTT.PATH, loadComponent: () => import('./pages/algorithms/gray-scott/gray-scott.component').then(m => m.GrayScottComponent), data: { seo: 'ALGORITHMS' } },

src/app/constants/RouterConstants.ts

@@ -55,6 +55,11 @@
     LINK: '/algorithms/cloth',
   };
 
+  static readonly CLOTH_ICM = {
+    PATH: 'algorithms/cloth-icm',
+    LINK: '/algorithms/cloth-icm',
+  };
+
   static readonly BOIDS = {
     PATH: 'algorithms/boids',
     LINK: '/algorithms/boids',

src/app/constants/UrlConstants.ts

@@ -21,6 +21,8 @@
   static readonly DOUBLE_PENDULUM_WIKI = 'https://de.wikipedia.org/wiki/Doppelpendel'
   static readonly CLOTH_SIMULATION_WIKI = 'https://en.wikipedia.org/wiki/Cloth_modeling'
   static readonly XPBD_WIKI = 'https://www.emergentmind.com/topics/extended-position-based-dynamics-xpbd'
+  static readonly ICM_PAPER = 'https://doi.org/10.1145/1141911.1142007'
+  static readonly KDOP_WIKI = 'https://en.wikipedia.org/wiki/Bounding_volume'
   static readonly BOIDS_WIKI = 'https://de.wikipedia.org/wiki/Boids'
   static readonly PARTICLE_LIFE_WIKI = 'https://particle-life.com'
   static readonly PHYSARUM_WIKI = 'https://en.wikipedia.org/wiki/Physarum_polycephalum'

src/app/pages/algorithms/algorithms.service.ts

@@ -78,6 +78,13 @@ export class AlgorithmsService {
       routerLink: RouterConstants.CLOTH.LINK,
       icon: 'texture'
     },
+    {
+      id: 'clothIcm',
+      title: 'ALGORITHM.CLOTH_ICM.TITLE',
+      description: 'ALGORITHM.CLOTH_ICM.DESCRIPTION',
+      routerLink: RouterConstants.CLOTH_ICM.LINK,
+      icon: 'layers'
+    },
     {
       id: 'boids',
       title: 'ALGORITHM.BOIDS.TITLE',

src/app/pages/algorithms/cloth-icm/cloth-icm-renderer.ts

@@ -0,0 +1,275 @@
+import {
+  ArcRotateCamera, Color3, Color4, CreateLineSystemVertexData, LinesMesh, Mesh, MeshBuilder, Scene, ShaderMaterial, StandardMaterial, Vector3, VertexBuffer, VertexData
+} from '@babylonjs/core';
+import {CLOTH_FRAGMENT_SHADER_GLSL, CLOTH_VERTEX_SHADER_GLSL} from '../cloth/cloth-glsl.shader';
+import {ClothIcmScene} from './cloth-icm-scenes';
+import {ClothIcmSolver} from './cloth-icm-solver';
+import {ClothBvh} from './collision/cloth-bvh';
+import {ClothTopology} from './collision/cloth-topology';
+import {DOP_SIZE} from './collision/dop18';
+import {Contour} from './collision/intersection-contours';
+import {IntersectionSet} from './collision/intersection-set';
+
+const CONTOUR_COLOR = new Color3(0.1, 1, 0.2);
+const CULLED_COLOR = new Color4(0.2, 0.8, 1, 1);
+const TESTED_COLOR = new Color4(0.6, 0.6, 0.6, 1);
+const OBSTACLE_COLOR = new Color3(0.35, 0.55, 0.9);
+const GROUND_SIZE = 8;
+/** Lines of group 1 are drawn after the cloth with a cleared depth buffer, so they are never hidden by it. */
+const OVERLAY_RENDERING_GROUP = 1;
+
+interface BvhBox {
+  node: number;
+  corners: Vector3[];
+  color: Color4;
+}
+
+/** Babylon side of the page: cloth mesh, obstacles, contour lines and BVH boxes. The solver never touches it. */
+export class ClothIcmRenderer {
+  private scene: Scene | null = null;
+  private clothMesh: Mesh | null = null;
+  private clothMaterial: ShaderMaterial | null = null;
+  private obstacleMaterial: StandardMaterial | null = null;
+  private obstacles: Mesh[] = [];
+  private contourLines: LinesMesh | null = null;
+  private bvhLines: LinesMesh | null = null;
+  private vertexBuffer = new Float32Array(0);
+  private culledRegion = new Uint8Array(0);
+  private bvhBoxes: BvhBox[] = [];
+  private bvhPaths: Vector3[][] = [];
+  private bvhPathColors: Color4[][] = [];
+  private bvhShownDepth = -1;
+
+  init(scene: Scene, topology: ClothTopology, sceneDef: ClothIcmScene): void {
+    this.dispose();
+    this.scene = scene;
+    this.vertexBuffer = new Float32Array(topology.vertexCount * 3);
+    this.createClothMesh(scene, topology);
+    this.createObstacles(scene, sceneDef);
+    this.setupCamera(scene);
+  }
+
+  update(solver: ClothIcmSolver, showBvh: boolean, bvhDepth: number): void {
+    if (!this.scene || !this.clothMesh) {
+      return;
+    }
+    this.updateClothPositions(solver.positions);
+    this.updateContourLines(this.scene, solver.intersections, solver.contours);
+    this.updateBvhLines(this.scene, solver.bvh, showBvh, bvhDepth);
+  }
+
+  dispose(): void {
+    this.contourLines?.dispose();
+    this.bvhLines?.dispose();
+    this.clothMesh?.dispose();
+    this.clothMaterial?.dispose();
+    this.obstacles.forEach(mesh => mesh.dispose());
+    this.obstacleMaterial?.dispose();
+    this.contourLines = null;
+    this.bvhLines = null;
+    this.bvhBoxes = [];
+    this.bvhPaths = [];
+    this.bvhPathColors = [];
+    this.clothMesh = null;
+    this.clothMaterial = null;
+    this.obstacles = [];
+    this.obstacleMaterial = null;
+    this.scene = null;
+  }
+
+  private createClothMesh(scene: Scene, topology: ClothTopology): void {
+    const {gridWidth, gridHeight, vertexCount} = topology;
+    const uvs = new Float32Array(vertexCount * 2);
+    for (let vertex = 0; vertex < vertexCount; vertex++) {
+      uvs[vertex * 2] = (vertex % gridWidth) / (gridWidth - 1);
+      uvs[vertex * 2 + 1] = 1 - Math.floor(vertex / gridWidth) / (gridHeight - 1);
+    }
+
+    const vertexData = new VertexData();
+    vertexData.positions = new Float32Array(vertexCount * 3);
+    vertexData.indices = topology.triangles;
+    vertexData.uvs = uvs;
+    const mesh = new Mesh('clothIcm', scene);
+    vertexData.applyToMesh(mesh, true);
+    // The bounding box is not refit when the positions change.
+    mesh.alwaysSelectAsActiveMesh = true;
+
+    const material = new ShaderMaterial('clothIcmMat', scene, {
+      vertexSource: CLOTH_VERTEX_SHADER_GLSL,
+      fragmentSource: CLOTH_FRAGMENT_SHADER_GLSL
+    }, {
+      attributes: ['position', 'uv'],
+      uniforms: ['viewProjection']
+    });
+    material.backFaceCulling = false;
+    mesh.material = material;
+    this.clothMesh = mesh;
+    this.clothMaterial = material;
+  }
+
+  /** The scene has no light, so the obstacles are unlit wireframes. */
+  private createObstacles(scene: Scene, sceneDef: ClothIcmScene): void {
+    const material = new StandardMaterial('clothIcmObstacle', scene);
+    material.wireframe = true;
+    material.disableLighting = true;
+    material.emissiveColor = OBSTACLE_COLOR;
+    this.obstacleMaterial = material;
+
+    const {sphere, groundY} = sceneDef;
+    if (sphere) {
+      const mesh = MeshBuilder.CreateSphere('clothIcmSphere', {diameter: sphere.r * 2, segments: 24}, scene);
+      mesh.position.set(sphere.x, sphere.y, sphere.z);
+      mesh.material = material;
+      this.obstacles.push(mesh);
+    }
+    if (groundY !== null) {
+      const mesh = MeshBuilder.CreateGround('clothIcmGround', {width: GROUND_SIZE, height: GROUND_SIZE, subdivisions: 16}, scene);
+      mesh.position.y = groundY;
+      mesh.material = material;
+      this.obstacles.push(mesh);
+    }
+  }
+
+  private setupCamera(scene: Scene): void {
+    const camera = scene.activeCamera as ArcRotateCamera | null;
+    if (camera) {
+      camera.alpha = Math.PI / 4;
+      camera.beta = Math.PI / 2.5;
+      camera.radius = 15;
+    }
+  }
+
+  private updateClothPositions(positions: Float32Array): void {
+    const buffer = this.vertexBuffer;
+    for (let vertex = 0; vertex < buffer.length / 3; vertex++) {
+      buffer[vertex * 3] = positions[vertex * 4];
+      buffer[vertex * 3 + 1] = positions[vertex * 4 + 1];
+      buffer[vertex * 3 + 2] = positions[vertex * 4 + 2];
+    }
+    this.clothMesh?.updateVerticesData(VertexBuffer.PositionKind, buffer);
+  }
+
+  /**
+   * The hit points come from the detection, which runs on the predicted positions before the solver moves the cloth,
+   * so the lines can be a little off the drawn cloth (up to about 0.05 m).
+   * The mesh is created once and gets new geometry every frame. Disposing and creating it again would drop its shader
+   * (reference counted) and recompile it; the line shader compiles asynchronously, so a mesh whose point layout changes
+   * every frame (the drop scene) was never ready to draw.
+   */
+  private updateContourLines(scene: Scene, set: IntersectionSet, contours: Contour[]): void {
+    if (contours.length === 0) {
+      this.contourLines?.setEnabled(false);
+      return;
+    }
+    const vertexData = CreateLineSystemVertexData({lines: contours.map(contour => this.contourPath(set, contour))});
+    if (!this.contourLines) {
+      const mesh = new LinesMesh('clothIcmContours', scene);
+      mesh.color = CONTOUR_COLOR;
+      mesh.renderingGroupId = OVERLAY_RENDERING_GROUP;
+      mesh.alwaysSelectAsActiveMesh = true;
+      this.contourLines = mesh;
+    }
+    vertexData.applyToMesh(this.contourLines, true);
+    this.contourLines.setEnabled(true);
+  }
+
+  private contourPath(set: IntersectionSet, contour: Contour): Vector3[] {
+    const path = contour.points.map(hit => new Vector3(set.points[hit * 3], set.points[hit * 3 + 1], set.points[hit * 3 + 2]));
+    if (contour.closed) {
+      path.push(path[0].clone());
+    }
+    return path;
+  }
+
+  /** The boxes are built once per shown level; every frame only their corners and colors are rewritten in place. */
+  private updateBvhLines(scene: Scene, bvh: ClothBvh, showBvh: boolean, depth: number): void {
+    if (!showBvh) {
+      this.bvhLines?.dispose();
+      this.bvhLines = null;
+      this.bvhBoxes = [];
+      return;
+    }
+    this.markCulledRegions(bvh);
+    if (!this.bvhLines || this.bvhShownDepth !== depth) {
+      this.rebuildBvhLines(scene, bvh, depth);
+    }
+    for (const box of this.bvhBoxes) {
+      box.color.copyFrom(this.culledRegion[box.node] === 1 ? CULLED_COLOR : TESTED_COLOR);
+      this.moveBox(bvh, box);
+    }
+    if (this.bvhLines) {
+      MeshBuilder.CreateLineSystem('clothIcmBvh', {lines: this.bvhPaths, colors: this.bvhPathColors, instance: this.bvhLines});
+    }
+  }
+
+  private rebuildBvhLines(scene: Scene, bvh: ClothBvh, depth: number): void {
+    this.bvhLines?.dispose();
+    this.bvhLines = null;
+    this.bvhBoxes = [];
+    this.bvhPaths = [];
+    this.bvhPathColors = [];
+    this.bvhShownDepth = depth;
+    for (let node = 0; node < bvh.nodeCount; node++) {
+      const isShown = bvh.depth[node] === depth || (bvh.depth[node] < depth && bvh.isLeaf(node));
+      if (isShown) {
+        this.addBox(node);
+      }
+    }
+    if (this.bvhPaths.length === 0) {
+      return;
+    }
+    const mesh = MeshBuilder.CreateLineSystem('clothIcmBvh', {lines: this.bvhPaths, colors: this.bvhPathColors, updatable: true}, scene);
+    mesh.renderingGroupId = OVERLAY_RENDERING_GROUP;
+    mesh.alwaysSelectAsActiveMesh = true;
+    this.bvhLines = mesh;
+  }
+
+  /** A node is marked when its own self-test or the one of an ancestor was culled. Parents have smaller ids than children. */
+  private markCulledRegions(bvh: ClothBvh): void {
+    if (this.culledRegion.length !== bvh.nodeCount) {
+      this.culledRegion = new Uint8Array(bvh.nodeCount);
+    }
+    for (let node = 0; node < bvh.nodeCount; node++) {
+      const parent = bvh.parent[node];
+      const parentCulled = parent >= 0 ? this.culledRegion[parent] : 0;
+      this.culledRegion[node] = parentCulled | bvh.culled[node];
+    }
+  }
+
+  /** A box is 6 polylines over 8 shared corners (bottom 0-3, top 4-7); all share one color object. */
+  private addBox(node: number): void {
+    const corners = Array.from({length: 8}, () => new Vector3());
+    const color = new Color4();
+    const bottom = corners.slice(0, 4);
+    const top = corners.slice(4);
+    const paths = [
+      [...bottom, bottom[0]],
+      [...top, top[0]],
+      ...bottom.map((corner, index) => [corner, top[index]])
+    ];
+    for (const path of paths) {
+      this.bvhPaths.push(path);
+      this.bvhPathColors.push(path.map(() => color));
+    }
+    this.bvhBoxes.push({node, corners, color});
+  }
+
+  /** Axis-aligned part of the 18-DOP (axes 0 to 2). */
+  private moveBox(bvh: ClothBvh, box: BvhBox): void {
+    const base = box.node * DOP_SIZE;
+    const {bounds} = bvh;
+    const {corners} = box;
+    const minX = bounds[base];
+    const minY = bounds[base + 1];
+    const minZ = bounds[base + 2];
+    const maxX = bounds[base + 9];
+    const maxY = bounds[base + 10];
+    const maxZ = bounds[base + 11];
+    for (let index = 0; index < 4; index++) {
+      const x = index === 1 || index === 2 ? maxX : minX;
+      const z = index >= 2 ? maxZ : minZ;
+      corners[index].set(x, minY, z);
+      corners[index + 4].set(x, maxY, z);
+    }
+  }
+}

src/app/pages/algorithms/cloth-icm/cloth-icm-scenes.spec.ts

@@ -0,0 +1,71 @@
+import {ClothBvh} from './collision/cloth-bvh';
+import {detectCollisions} from './collision/collision-detection';
+import {buildClothTopology} from './collision/cloth-topology';
+import {IntersectionSet} from './collision/intersection-set';
+import {ClothIcmSceneId} from './cloth-icm.model';
+import {applyContactFriction, initialPositions, projectObstacles, sceneOf} from './cloth-icm-scenes';
+
+function intersectionsAtStart(id: ClothIcmSceneId): number {
+  const topology = buildClothTopology(30, 30);
+  const positions = initialPositions(id, topology, 0.1, 1);
+  const bvh = new ClothBvh(topology);
+  bvh.refit(positions);
+  const set = new IntersectionSet();
+  detectCollisions(bvh, positions, topology, 0, set, []);
+  return set.count;
+}
+
+describe('cloth icm scenes', () => {
+  it('starts the tangled scene with intersections', () => {
+    expect(intersectionsAtStart('tangled')).toBeGreaterThan(0);
+  });
+
+  it('starts the drop scene without intersections', () => {
+    expect(intersectionsAtStart('drop')).toBe(0);
+  });
+
+  it('pins only the top row of the tangled scene', () => {
+    const topology = buildClothTopology(6, 6);
+    const positions = initialPositions('tangled', topology, 0.1, 2);
+    const pinned = Array.from({length: topology.vertexCount}, (_, vertex) => positions[vertex * 4 + 3] === 0);
+    expect(pinned.filter(Boolean).length).toBe(6);
+    expect(pinned.slice(0, 6).every(Boolean)).toBe(true);
+  });
+
+  it('pushes vertices out of the sphere and above the ground', () => {
+    const scene = sceneOf('drop');
+    const sphere = scene.sphere!;
+    const groundY = scene.groundY!;
+    const positions = new Float32Array([sphere.x, sphere.y + 0.1, sphere.z, 1, 0, groundY - 10, 5, 1, 0, groundY + 3, 5, 1]);
+    projectObstacles(positions, scene, 0.1);
+    const distance = Math.hypot(positions[0] - sphere.x, positions[1] - sphere.y, positions[2] - sphere.z);
+    expect(distance).toBeCloseTo(sphere.r + 0.1);
+    expect(positions[5]).toBeCloseTo(groundY + 0.1);
+    expect(positions[9]).toBe(groundY + 3);
+  });
+
+  it('removes tangential motion of a vertex that rests on the ground', () => {
+    const scene = sceneOf('drop');
+    const groundY = scene.groundY!;
+    const previous = new Float32Array([0, groundY + 0.1, 0, 1]);
+    const positions = new Float32Array([1, groundY + 0.1, 2, 1]);
+    applyContactFriction(positions, previous, scene, 0.1, 0.5);
+    expect(positions[0]).toBeCloseTo(0.5);
+    expect(positions[2]).toBeCloseTo(1);
+    expect(positions[1]).toBeCloseTo(groundY + 0.1);
+  });
+
+  it('leaves a vertex in the air alone when friction is applied', () => {
+    const scene = sceneOf('drop');
+    const previous = new Float32Array([0, 3, 0, 1]);
+    const positions = new Float32Array([1, 3, 2, 1]);
+    applyContactFriction(positions, previous, scene, 0.1, 0.5);
+    expect(Array.from(positions)).toEqual([1, 3, 2, 1]);
+  });
+
+  it('leaves pinned vertices alone', () => {
+    const positions = new Float32Array([0, 1, 0, 0]);
+    projectObstacles(positions, sceneOf('drop'), 0.1);
+    expect(positions[1]).toBe(1);
+  });
+});

src/app/pages/algorithms/cloth-icm/cloth-icm-scenes.ts

@@ -0,0 +1,218 @@
+import {ClothIcmSceneId} from './cloth-icm.model';
+import {ClothTopology} from './collision/cloth-topology';
+
+interface Sphere {
+  x: number;
+  y: number;
+  z: number;
+  r: number;
+}
+
+export interface ClothIcmScene {
+  gravityY: number;
+  sphere: Sphere | null;
+  groundY: number | null;
+}
+
+const SCENES: Record<ClothIcmSceneId, ClothIcmScene> = {
+  tangled: {gravityY: 0, sphere: null, groundY: null},
+  drop: {gravityY: -9.81, sphere: {x: 0, y: -0.3, z: 0, r: 0.7}, groundY: -1}
+};
+
+const DROP_HEIGHT = 0.9;
+// The back half hangs one spacing behind the front half, so the fold edges have their rest length.
+const FOLD_GAP_IN_SPACINGS = 1;
+const BUMP_HEIGHT_IN_SPACINGS = 1.5;
+const BUMP_RADIUS_IN_SPACINGS = 8;
+const BUMP_MARGIN_IN_SPACINGS = 2;
+// Keeps the back half off the grid lines of the front half; hits through a vertex or an edge are dropped by the detection.
+const BACK_OFFSET_X_IN_SPACINGS = 0.13;
+const BACK_OFFSET_Y_IN_SPACINGS = 0.07;
+
+export function sceneOf(id: ClothIcmSceneId): ClothIcmScene {
+  return SCENES[id];
+}
+
+/** Vec4 layout [x, y, z, invMass]. */
+export function initialPositions(id: ClothIcmSceneId, topology: ClothTopology, spacing: number, invMass: number): Float32Array {
+  return id === 'tangled' ? tangledPositions(topology, spacing, invMass) : dropPositions(topology, spacing, invMass);
+}
+
+/**
+ * Hangs from the pinned top row. The rows below the middle are folded back up behind the front half
+ * and a ridge across the back half pokes through the front half: open contours from border to border at frame 0.
+ */
+function tangledPositions(topology: ClothTopology, spacing: number, invMass: number): Float32Array {
+  const {gridWidth, gridHeight} = topology;
+  const positions = new Float32Array(topology.vertexCount * 4);
+  const frontRows = Math.floor(gridHeight / 2);
+  const top = ((gridHeight - 1) * spacing) / 2;
+  const frontEndY = top - (frontRows - 1) * spacing;
+  const halfBackHeight = ((gridHeight - frontRows - 1) * spacing) / 2;
+  const backCenterY = frontEndY + halfBackHeight;
+  // Small grids need a narrower ridge so that it stays inside the back half.
+  const ridgeRadius = Math.min(BUMP_RADIUS_IN_SPACINGS * spacing, halfBackHeight - BUMP_MARGIN_IN_SPACINGS * spacing);
+
+  const profile = bumpProfile(BUMP_HEIGHT_IN_SPACINGS * spacing, ridgeRadius);
+  for (let row = 0; row < gridHeight; row++) {
+    for (let column = 0; column < gridWidth; column++) {
+      const base = (row * gridWidth + column) * 4;
+      const isFront = row < frontRows;
+      const x = (column - (gridWidth - 1) / 2 + (isFront ? 0 : BACK_OFFSET_X_IN_SPACINGS)) * spacing;
+      let y = top - row * spacing;
+      let z = 0;
+      if (!isFront) {
+        y = frontEndY + (row - frontRows + BACK_OFFSET_Y_IN_SPACINGS) * spacing;
+        const [ridgeY, ridgeZ] = ridgePoint(profile, y - backCenterY);
+        y = backCenterY + ridgeY;
+        z = -FOLD_GAP_IN_SPACINGS * spacing + ridgeZ;
+      }
+      positions[base] = x;
+      positions[base + 1] = y;
+      positions[base + 2] = z;
+      positions[base + 3] = row === 0 ? 0 : invMass;
+    }
+  }
+  return positions;
+}
+
+interface BumpProfile {
+  amplitude: number;
+  radius: number;
+  /** Arc length of the profile from the center to radius `i * radius / BUMP_TABLE_STEPS`. */
+  arcLength: Float64Array;
+}
+
+const BUMP_TABLE_STEPS = 512;
+
+/** Smooth bump with compact support: exactly flat beyond `radius`. */
+function bumpHeight(distance: number, amplitude: number, radius: number): number {
+  const falloff = Math.max(0, 1 - (distance * distance) / (radius * radius));
+  return amplitude * falloff ** 3;
+}
+
+function bumpProfile(amplitude: number, radius: number): BumpProfile {
+  const arcLength = new Float64Array(BUMP_TABLE_STEPS + 1);
+  const step = radius / BUMP_TABLE_STEPS;
+  for (let index = 1; index <= BUMP_TABLE_STEPS; index++) {
+    const rise = bumpHeight(index * step, amplitude, radius) - bumpHeight((index - 1) * step, amplitude, radius);
+    arcLength[index] = arcLength[index - 1] + Math.hypot(step, rise);
+  }
+  return {amplitude, radius, arcLength};
+}
+
+/**
+ * Bends the flat offset `u` (across the ridge, from its center line) into the ridge profile. The profile radius with
+ * the same arc length as `u` is used, so edges across the ridge keep their rest length, and edges along it are
+ * untouched. A ridge is a developable shape: no edge is stretched, the cloth has no stress that could flatten it.
+ * The flat rest of the sheet is shifted inward by the extra arc length.
+ */
+function ridgePoint(profile: BumpProfile, u: number): [number, number] {
+  const {arcLength, radius, amplitude} = profile;
+  const flatDistance = Math.abs(u);
+  const total = arcLength[BUMP_TABLE_STEPS];
+  let projected: number;
+  if (flatDistance >= total) {
+    projected = radius + (flatDistance - total);
+  } else {
+    let index = 1;
+    while (arcLength[index] < flatDistance) {
+      index++;
+    }
+    const share = (flatDistance - arcLength[index - 1]) / (arcLength[index] - arcLength[index - 1]);
+    projected = ((index - 1 + share) * radius) / BUMP_TABLE_STEPS;
+  }
+  return [Math.sign(u) * projected, bumpHeight(projected, amplitude, radius)];
+}
+
+function dropPositions(topology: ClothTopology, spacing: number, invMass: number): Float32Array {
+  const {gridWidth, gridHeight} = topology;
+  const positions = new Float32Array(topology.vertexCount * 4);
+  for (let row = 0; row < gridHeight; row++) {
+    for (let column = 0; column < gridWidth; column++) {
+      const base = (row * gridWidth + column) * 4;
+      positions[base] = (column - (gridWidth - 1) / 2) * spacing;
+      positions[base + 1] = DROP_HEIGHT;
+      positions[base + 2] = (row - (gridHeight - 1) / 2) * spacing;
+      positions[base + 3] = invMass;
+    }
+  }
+  return positions;
+}
+
+/** Pushes free vertices out of the sphere and above the ground, both fattened by `thickness`. */
+export function projectObstacles(positions: Float32Array, scene: ClothIcmScene, thickness: number): void {
+  const vertexCount = positions.length / 4;
+  const {sphere, groundY} = scene;
+  for (let vertex = 0; vertex < vertexCount; vertex++) {
+    const base = vertex * 4;
+    if (positions[base + 3] <= 0) {
+      continue;
+    }
+    if (sphere) {
+      pushOutOfSphere(positions, base, sphere, thickness);
+    }
+    if (groundY !== null && positions[base + 1] < groundY + thickness) {
+      positions[base + 1] = groundY + thickness;
+    }
+  }
+}
+
+function pushOutOfSphere(positions: Float32Array, base: number, sphere: Sphere, thickness: number): void {
+  const dx = positions[base] - sphere.x;
+  const dy = positions[base + 1] - sphere.y;
+  const dz = positions[base + 2] - sphere.z;
+  const distance = Math.hypot(dx, dy, dz);
+  const radius = sphere.r + thickness;
+  if (distance >= radius || distance < 1e-9) {
+    return;
+  }
+  const scale = radius / distance;
+  positions[base] = sphere.x + dx * scale;
+  positions[base + 1] = sphere.y + dy * scale;
+  positions[base + 2] = sphere.z + dz * scale;
+}
+
+const CONTACT_TOLERANCE = 1e-4;
+
+/**
+ * Removes a share of the tangential motion of the vertices that rest on an obstacle. Without it the cloth slides
+ * off the sphere without any resistance and piles up on the ground, where the folds cannot be resolved.
+ */
+export function applyContactFriction(positions: Float32Array, previous: Float32Array, scene: ClothIcmScene, thickness: number, friction: number): void {
+  const vertexCount = positions.length / 4;
+  const {sphere, groundY} = scene;
+  for (let vertex = 0; vertex < vertexCount; vertex++) {
+    const base = vertex * 4;
+    if (positions[base + 3] <= 0) {
+      continue;
+    }
+    if (groundY !== null && positions[base + 1] <= groundY + thickness + CONTACT_TOLERANCE) {
+      positions[base] -= friction * (positions[base] - previous[base]);
+      positions[base + 2] -= friction * (positions[base + 2] - previous[base + 2]);
+    }
+    if (sphere) {
+      slowOnSphere(positions, previous, base, sphere, thickness, friction);
+    }
+  }
+}
+
+function slowOnSphere(
+  positions: Float32Array, previous: Float32Array, base: number, sphere: Sphere,
+  thickness: number, friction: number
+): void {
+  const nx = positions[base] - sphere.x;
+  const ny = positions[base + 1] - sphere.y;
+  const nz = positions[base + 2] - sphere.z;
+  const distance = Math.hypot(nx, ny, nz);
+  if (distance > sphere.r + thickness + CONTACT_TOLERANCE || distance < 1e-9) {
+    return;
+  }
+  const dx = positions[base] - previous[base];
+  const dy = positions[base + 1] - previous[base + 1];
+  const dz = positions[base + 2] - previous[base + 2];
+  const normalPart = (dx * nx + dy * ny + dz * nz) / (distance * distance);
+  positions[base] -= friction * (dx - normalPart * nx);
+  positions[base + 1] -= friction * (dy - normalPart * ny);
+  positions[base + 2] -= friction * (dz - normalPart * nz);
+}

src/app/pages/algorithms/cloth-icm/cloth-icm-solver.spec.ts

@@ -0,0 +1,145 @@
+import {ClothIcmConfig, ClothIcmSceneId, SCENE_SETTINGS} from './cloth-icm.model';
+import {ClothIcmSolver} from './cloth-icm-solver';
+
+const DT = 1 / 60;
+const TANGLED_SPACING = 0.1;
+
+function config(gridSize: number, thickness = 0.01, proximityRedetectEvery = 0, spacing = TANGLED_SPACING): ClothIcmConfig {
+  return {gridWidth: gridSize, gridHeight: gridSize, spacing, particleInvMass: 1 / (spacing * spacing), thickness, proximityRedetectEvery};
+}
+
+/** Mean hits per frame over frames 100..240 and the share of those frames with 0 hits, drop scene with the page settings. */
+function dropHitsAfterLanding(icmEnabled: boolean, proximityRedetectEvery: number): {average: number; zeroShare: number} {
+  const {gridSize, spacing, thickness} = SCENE_SETTINGS.drop;
+  const solver = new ClothIcmSolver(config(gridSize, thickness, proximityRedetectEvery, spacing), 'drop');
+  const first = 100;
+  const last = 240;
+  let sum = 0;
+  let zeroFrames = 0;
+  for (let frame = 1; frame <= last; frame++) {
+    solver.step({dt: DT, icmEnabled, icmStep: 0.25});
+    if (frame > first) {
+      sum += solver.stats.intersections;
+      zeroFrames += solver.stats.intersections === 0 ? 1 : 0;
+    }
+  }
+  return {average: sum / (last - first), zeroShare: zeroFrames / (last - first)};
+}
+
+function run(
+  sceneId: ClothIcmSceneId, gridSize: number, frames: number, icmEnabled: boolean, thickness?: number, proximityRedetectEvery?: number
+): {solver: ClothIcmSolver; firstContours: number} {
+  const solver = new ClothIcmSolver(config(gridSize, thickness, proximityRedetectEvery), sceneId);
+  let firstContours = 0;
+  for (let frame = 0; frame < frames; frame++) {
+    solver.step({dt: DT, icmEnabled, icmStep: 0.25});
+    if (frame === 0) {
+      firstContours = solver.contours.length;
+    }
+  }
+  return {solver, firstContours};
+}
+
+function runFrames(solver: ClothIcmSolver, frames: number, icmStep: number, dt: number, icmEnabled = true): number {
+  for (let frame = 0; frame < frames; frame++) {
+    solver.step({dt, icmEnabled, icmStep});
+  }
+  return solver.stats.intersections;
+}
+
+describe('ClothIcmSolver', () => {
+  it('untangles the tangled scene with ICM on', () => {
+    const {solver, firstContours} = run('tangled', 30, 150, true);
+    expect(firstContours).toBeGreaterThan(0);
+    expect(solver.contours.length).toBe(0);
+  });
+
+  it('keeps the tangled scene tangled with ICM off', () => {
+    const {solver, firstContours} = run('tangled', 30, 150, false);
+    expect(firstContours).toBeGreaterThan(0);
+    expect(solver.contours.length).toBeGreaterThan(0);
+  });
+
+  it('keeps the dropped cloth finite and above the ground', () => {
+    const {gridSize, spacing, thickness, proximityRedetectEvery} = SCENE_SETTINGS.drop;
+    const solver = new ClothIcmSolver(config(gridSize, thickness, proximityRedetectEvery, spacing), 'drop');
+    runFrames(solver, 300, 0.25, DT);
+    const groundY = solver.scene.groundY ?? -Infinity;
+    for (let vertex = 0; vertex < solver.topology.vertexCount; vertex++) {
+      const y = solver.positions[vertex * 4 + 1];
+      expect(Number.isFinite(y)).toBe(true);
+      expect(y).toBeGreaterThanOrEqual(groundY);
+    }
+    expect(solver.positions.every(Number.isFinite)).toBe(true);
+  });
+
+  describe('drop scene after landing', () => {
+    const {proximityRedetectEvery} = SCENE_SETTINGS.drop;
+    const icmOn = dropHitsAfterLanding(true, proximityRedetectEvery);
+
+    it('keeps the dropped page grid almost free of intersections', () => {
+      expect(icmOn.average).toBeLessThanOrEqual(5);
+      expect(icmOn.zeroShare).toBeGreaterThanOrEqual(0.5);
+    });
+
+    it('shows the ICM effect', () => {
+      expect(dropHitsAfterLanding(false, proximityRedetectEvery).average).toBeGreaterThanOrEqual(3 * icmOn.average);
+    });
+
+    it('re-detects proximity inside the frame', () => {
+      expect(proximityRedetectEvery).toBeGreaterThan(0);
+      expect(dropHitsAfterLanding(true, 0).average).toBeGreaterThanOrEqual(3 * icmOn.average);
+    });
+  });
+
+  it('keeps the BVH culled flags of the frame-start detection', () => {
+    const {gridSize, spacing, thickness, proximityRedetectEvery} = SCENE_SETTINGS.drop;
+    const solver = new ClothIcmSolver(config(gridSize, thickness, proximityRedetectEvery, spacing), 'drop');
+    let framesWithMismatch = 0;
+    for (let frame = 0; frame < 60; frame++) {
+      solver.step({dt: DT, icmEnabled: true, icmStep: 0.25});
+      const flaggedNodes = solver.bvh.culled.reduce((sum, flag) => sum + flag, 0);
+      if (flaggedNodes !== solver.stats.culledNodes) {
+        framesWithMismatch++;
+      }
+    }
+    expect(framesWithMismatch).toBe(0);
+  });
+
+  it('reports collision stats', () => {
+    const {solver} = run('tangled', 30, 1, true);
+    expect(solver.stats.intersections).toBeGreaterThan(0);
+    expect(solver.stats.testedPairs).toBeGreaterThan(0);
+    expect(solver.stats.collisionMs).toBeGreaterThanOrEqual(0);
+  });
+
+  it.each([0.05, 0.3, 0.5])('resolves the tangled page grid at icmStep %s', icmStep => {
+    const solver = new ClothIcmSolver(config(50), 'tangled');
+    expect(runFrames(solver, 1, icmStep, DT)).toBeGreaterThan(0);
+    runFrames(solver, 299, icmStep, DT);
+    expect(solver.stats.intersections).toBe(0);
+  });
+
+  it('keeps the tangled page grid unchanged at rest with ICM off', () => {
+    const solver = new ClothIcmSolver(config(50), 'tangled');
+    const start = runFrames(solver, 1, 0.25, DT, false);
+    runFrames(solver, 299, 0.25, DT, false);
+    expect(solver.stats.intersections).toBe(start);
+  });
+
+  it('gives the same result at 144 fps as at 60 fps', () => {
+    const solver = new ClothIcmSolver(config(50), 'tangled');
+    runFrames(solver, 360, 0.25, 1 / 144);
+    expect(solver.stats.intersections).toBe(0);
+  });
+
+  it('never moves the pinned top row', () => {
+    const {solver} = run('tangled', 30, 30, true);
+    const initial = new ClothIcmSolver(config(30), 'tangled');
+    for (let vertex = 0; vertex < 30; vertex++) {
+      for (let axis = 0; axis < 3; axis++) {
+        expect(solver.positions[vertex * 4 + axis]).toBe(initial.positions[vertex * 4 + axis]);
+      }
+    }
+  });
+});

src/app/pages/algorithms/cloth-icm/cloth-icm-solver.ts

@@ -0,0 +1,224 @@
+import {ClothPhysicsParams, solveConstraints, updateVelocities} from '../cloth/strategies/cloth-cpu-physics';
+import {
+  CONTACT_FRICTION, ClothIcmConfig, ClothIcmSceneId, ClothIcmStats, DISTANCE_COMPLIANCE, emptyClothIcmStats, ICM_G0, ICM_SUBSTEPS, MAX_DT, MAX_SPEED, SIM_DT, VELOCITY_DAMPING
+} from './cloth-icm.model';
+import {applyContactFriction, ClothIcmScene, initialPositions, projectObstacles, sceneOf} from './cloth-icm-scenes';
+import {ClothBvh} from './collision/cloth-bvh';
+import {buildClothTopology, buildDistanceConstraints, ClothTopology} from './collision/cloth-topology';
+import {detectCollisions} from './collision/collision-detection';
+import {Contour, traceContours} from './collision/intersection-contours';
+import {IntersectionSet} from './collision/intersection-set';
+import {buildProximityConstraints, solveProximityConstraints} from './collision/proximity-constraints';
+import {accumulateIcmDisplacement, applyAveragedDisplacement} from './response/icm-response';
+
+const STEP_TOLERANCE = 1e-9;
+
+export interface ClothIcmFrameParams {
+  dt: number;
+  icmEnabled: boolean;
+  /** Maximum ICM displacement per frame as a fraction of the grid spacing. */
+  icmStep: number;
+}
+
+/** CPU cloth with self-collision handling. Pure TypeScript, no Babylon: one `step()` per frame. */
+export class ClothIcmSolver {
+  readonly topology: ClothTopology;
+  /** Frame-start detection; the renderer reads its culled flags and bounds. */
+  readonly bvh: ClothBvh;
+  readonly scene: ClothIcmScene;
+  /** Vec4 layout [x, y, z, invMass]. */
+  readonly positions: Float32Array;
+  readonly intersections = new IntersectionSet();
+  contours: Contour[] = [];
+  stats: ClothIcmStats = emptyClothIcmStats();
+
+  private readonly prevPositions: Float32Array;
+  private readonly velocities: Float32Array;
+  private readonly distanceConstraints: Float32Array;
+  private readonly icmDelta: Float32Array;
+  private readonly icmCount: Float32Array;
+  /** Vertices with an ICM displacement this step, and their positions before the proximity solve. */
+  private readonly icmVertices: Int32Array;
+  private icmVertexCount = 0;
+  private readonly proximitySnapshot: Float32Array;
+  private accumulator = 0;
+  private readonly proximityCandidates: number[] = [];
+  private proximityConstraints: Float32Array = new Float32Array(0);
+  private readonly midFrameIntersections = new IntersectionSet();
+  private readonly midFrameBvh: ClothBvh;
+  private readonly physics: ClothPhysicsParams;
+
+  constructor(private readonly config: ClothIcmConfig, sceneId: ClothIcmSceneId) {
+    this.topology = buildClothTopology(config.gridWidth, config.gridHeight);
+    this.bvh = new ClothBvh(this.topology);
+    this.midFrameBvh = new ClothBvh(this.topology);
+    this.scene = sceneOf(sceneId);
+    const vertexCount = this.topology.vertexCount;
+    this.positions = initialPositions(sceneId, this.topology, config.spacing, config.particleInvMass);
+    this.prevPositions = this.positions.slice();
+    this.velocities = new Float32Array(vertexCount * 4);
+    this.distanceConstraints = buildDistanceConstraints(this.topology, config.spacing);
+    this.icmDelta = new Float32Array(vertexCount * 3);
+    this.icmCount = new Float32Array(vertexCount);
+    this.icmVertices = new Int32Array(vertexCount);
+    this.proximitySnapshot = new Float32Array(vertexCount * 3);
+    this.physics = {
+      dt: 0, gravityY: this.scene.gravityY, compliance: DISTANCE_COMPLIANCE, numVertices: vertexCount,
+      windX: 0, windY: 0, windZ: 0, time: 0, elongation: 1
+    };
+  }
+
+  /** Advances by the frame time in steps of `SIM_DT`; a frame can run zero, one or two steps. */
+  step(params: ClothIcmFrameParams): void {
+    if (!(params.dt > 0)) {
+      return;
+    }
+    this.accumulator += Math.min(params.dt, MAX_DT);
+    while (this.accumulator >= SIM_DT - STEP_TOLERANCE) {
+      this.accumulator = Math.max(0, this.accumulator - SIM_DT);
+      this.simulateStep(params);
+    }
+  }
+
+  private simulateStep(params: ClothIcmFrameParams): void {
+    this.physics.dt = SIM_DT;
+    this.predictPositions();
+
+    const collisionStart = performance.now();
+    const counters = this.detect(this.bvh, this.intersections);
+    this.contours = traceContours(this.intersections);
+    this.prepareIcm(params);
+    this.buildProximity();
+    const detectMs = performance.now() - collisionStart;
+
+    const redetectMs = this.solve();
+    const collisionMs = detectMs + redetectMs;
+    applyContactFriction(this.positions, this.prevPositions, this.scene, this.config.thickness, CONTACT_FRICTION);
+    updateVelocities(this.positions, this.prevPositions, this.velocities, this.physics);
+    this.limitSpeed();
+
+    this.stats = {
+      intersections: this.intersections.count,
+      contours: this.contours.length,
+      testedPairs: counters.testedPairs,
+      culledNodes: counters.culledNodes,
+      collisionMs
+    };
+  }
+
+  /**
+   * Gravity and damping, then the predicted position. Written here and not taken from `integratePositions`:
+   * that one always adds a sine flutter force in y, even with the wind off.
+   */
+  private predictPositions(): void {
+    const {dt, gravityY, numVertices} = this.physics;
+    const {velocities} = this;
+    for (let vertex = 0; vertex < numVertices; vertex++) {
+      const base = vertex * 4;
+      if (this.positions[base + 3] <= 0) {
+        continue;
+      }
+      velocities[base] *= VELOCITY_DAMPING;
+      velocities[base + 1] = velocities[base + 1] * VELOCITY_DAMPING + gravityY * dt;
+      velocities[base + 2] *= VELOCITY_DAMPING;
+      for (let axis = 0; axis < 3; axis++) {
+        this.prevPositions[base + axis] = this.positions[base + axis];
+        this.positions[base + axis] += velocities[base + axis] * dt;
+      }
+    }
+  }
+
+  private limitSpeed(): void {
+    const {velocities} = this;
+    for (let vertex = 0; vertex < this.physics.numVertices; vertex++) {
+      const base = vertex * 4;
+      const speed = Math.hypot(velocities[base], velocities[base + 1], velocities[base + 2]);
+      if (speed > MAX_SPEED) {
+        const scale = MAX_SPEED / speed;
+        velocities[base] *= scale;
+        velocities[base + 1] *= scale;
+        velocities[base + 2] *= scale;
+      }
+    }
+  }
+
+  private detect(bvh: ClothBvh, intersections: IntersectionSet): {testedPairs: number; culledNodes: number} {
+    bvh.refit(this.positions);
+    return detectCollisions(bvh, this.positions, this.topology, this.config.thickness, intersections, this.proximityCandidates);
+  }
+
+  private buildProximity(): void {
+    this.proximityConstraints = buildProximityConstraints(this.positions, this.topology, this.proximityCandidates, this.config.thickness);
+  }
+
+  /** Only the proximity constraints follow the current positions; the frame-start hits and the ICM displacement stay. */
+  private redetectProximity(): number {
+    const start = performance.now();
+    this.midFrameBvh.coneCullingEnabled = this.bvh.coneCullingEnabled;
+    this.detect(this.midFrameBvh, this.midFrameIntersections);
+    this.buildProximity();
+    return performance.now() - start;
+  }
+
+  private prepareIcm(params: ClothIcmFrameParams): void {
+    this.icmDelta.fill(0);
+    this.icmCount.fill(0);
+    this.icmVertexCount = 0;
+    if (!params.icmEnabled || this.contours.length === 0) {
+      return;
+    }
+    const icmParams = {h0: params.icmStep * this.config.spacing, g0: ICM_G0};
+    accumulateIcmDisplacement(
+      this.positions, this.topology, this.bvh.triangleNormals, this.intersections, this.contours, icmParams, this.icmDelta, this.icmCount
+    );
+    for (let vertex = 0; vertex < this.icmCount.length; vertex++) {
+      if (this.icmCount[vertex] > 0) {
+        this.icmVertices[this.icmVertexCount++] = vertex;
+      }
+    }
+  }
+
+  /** The ICM displacement is spread over the iterations so that the distance constraints can react to it. Returns the re-detection time. */
+  private solve(): number {
+    const icmScale = 1 / ICM_SUBSTEPS;
+    const every = this.config.proximityRedetectEvery;
+    let redetectMs = 0;
+    for (let iteration = 0; iteration < ICM_SUBSTEPS; iteration++) {
+      if (every > 0 && iteration > 0 && iteration % every === 0) {
+        redetectMs += this.redetectProximity();
+      }
+      solveConstraints(this.positions, this.distanceConstraints, this.physics);
+      this.solveProximityBesideIcm();
+      applyAveragedDisplacement(this.positions, this.icmDelta, this.icmCount, icmScale);
+      projectObstacles(this.positions, this.scene, this.config.thickness);
+    }
+    return redetectMs;
+  }
+
+  /**
+   * Diploma 5.3 (p. 44-45): a proximity response that works against the ICM displacement of a vertex is dropped.
+   * Next to a contour, proximity pushes the vertices that already crossed to the wrong side and widens the intersection.
+   */
+  private solveProximityBesideIcm(): void {
+    const {positions, icmVertices, proximitySnapshot, icmDelta} = this;
+    for (let index = 0; index < this.icmVertexCount; index++) {
+      const vertex = icmVertices[index];
+      proximitySnapshot[index * 3] = positions[vertex * 4];
+      proximitySnapshot[index * 3 + 1] = positions[vertex * 4 + 1];
+      proximitySnapshot[index * 3 + 2] = positions[vertex * 4 + 2];
+    }
+    solveProximityConstraints(positions, this.proximityConstraints, this.config.thickness);
+    for (let index = 0; index < this.icmVertexCount; index++) {
+      const vertex = icmVertices[index];
+      const moveX = positions[vertex * 4] - proximitySnapshot[index * 3];
+      const moveY = positions[vertex * 4 + 1] - proximitySnapshot[index * 3 + 1];
+      const moveZ = positions[vertex * 4 + 2] - proximitySnapshot[index * 3 + 2];
+      const alongIcm = moveX * icmDelta[vertex * 3] + moveY * icmDelta[vertex * 3 + 1] + moveZ * icmDelta[vertex * 3 + 2];
+      if (alongIcm < 0) {
+        positions[vertex * 4] = proximitySnapshot[index * 3];
+        positions[vertex * 4 + 1] = proximitySnapshot[index * 3 + 1];
+        positions[vertex * 4 + 2] = proximitySnapshot[index * 3 + 2];
+      }
+    }
+  }
+}

src/app/pages/algorithms/cloth-icm/cloth-icm.component.html

@@ -0,0 +1,36 @@
+<mat-card class="w-full max-w-[1920px] p-5">
+  <mat-card-header>
+    <mat-card-title>{{ 'CLOTH_ICM.TITLE' | translate }}</mat-card-title>
+  </mat-card-header>
+  <mat-card-content>
+    <app-information [algorithmInformation]="algoInformation"/>
+    <div class="flex flex-col mb-4">
+      <div class="flex gap-4 mb-4 items-center flex-wrap mt-2.5 text-[0.9em]">
+        <mat-button-toggle-group class="rounded-sm overflow-hidden" [value]="sceneId()" (change)="setScene($event.value)" [aria-label]="'CLOTH_ICM.SCENE' | translate">
+          <mat-button-toggle value="tangled">{{ 'CLOTH_ICM.SCENE_TANGLED' | translate }}</mat-button-toggle>
+          <mat-button-toggle value="drop">{{ 'CLOTH_ICM.SCENE_DROP' | translate }}</mat-button-toggle>
+        </mat-button-toggle-group>
+        <button mat-raised-button color="primary" (click)="toggleIcm()">
+          {{ icmEnabled() ? ('CLOTH_ICM.ICM_OFF' | translate) : ('CLOTH_ICM.ICM_ON' | translate) }}
+        </button>
+        <button mat-raised-button color="primary" (click)="toggleCulling()">
+          {{ coneCulling() ? ('CLOTH_ICM.CULLING_OFF' | translate) : ('CLOTH_ICM.CULLING_ON' | translate) }}
+        </button>
+        <button mat-raised-button color="primary" (click)="toggleBvh()">
+          {{ showBvh() ? ('CLOTH_ICM.BVH_OFF' | translate) : ('CLOTH_ICM.BVH_ON' | translate) }}
+        </button>
+        <button mat-raised-button color="accent" (click)="restartSimulation()">
+          {{ 'CLOTH_ICM.RESTART_SIMULATION' | translate }}
+        </button>
+      </div>
+      <app-slider-field label="CLOTH_ICM.ICM_STEP" [min]="0.05" [max]="0.5" [step]="0.05" [(value)]="icmStep"/>
+      <app-slider-field label="CLOTH_ICM.BVH_DEPTH" [min]="0" [max]="maxBvhDepth" [step]="1" [(value)]="bvhDepth"/>
+      <p class="text-[0.9em] opacity-80">{{ 'CLOTH_ICM.HINT' | translate }}</p>
+    </div>
+    <app-metrics [metrics]="metrics()"/>
+    <app-babylon-canvas
+      [config]="renderConfig"
+      (sceneReady)="onSceneReady($event)"
+    />
+  </mat-card-content>
+</mat-card>

src/app/pages/algorithms/cloth-icm/cloth-icm.component.ts

@@ -0,0 +1,204 @@
+import {Component, computed, OnDestroy, signal} from '@angular/core';
+import {MatButton} from '@angular/material/button';
+import {MatButtonToggleModule} from '@angular/material/button-toggle';
+import {MatCard, MatCardContent, MatCardHeader, MatCardTitle} from '@angular/material/card';
+import {Observer, Scene} from '@babylonjs/core';
+import {TranslatePipe} from '@ngx-translate/core';
+import {UrlConstants} from '../../../constants/UrlConstants';
+import {BabylonCanvas, RenderConfig, SceneEventData} from '../../../shared/components/render-canvas/babylon-canvas.component';
+import {Metric, MetricsComponent} from '../../../shared/components/metrics/metrics';
+import {SliderFieldComponent} from '../../../shared/components/slider-field/slider-field';
+import {FpsCounter} from '../../../shared/fps-counter';
+import {Information} from '../information/information';
+import {AlgorithmInformation} from '../information/information.models';
+import {SOURCE_PATHS} from '../information/source-paths';
+import {ClothIcmConfig, ClothIcmSceneId, ClothIcmStats, emptyClothIcmStats, SCENE_SETTINGS} from './cloth-icm.model';
+import {ClothIcmRenderer} from './cloth-icm-renderer';
+import {ClothIcmSolver} from './cloth-icm-solver';
+
+const DEFAULT_ICM_STEP = 0.25;
+const MAX_BVH_DEPTH = 10;
+const STATS_EVERY_FRAMES = 10;
+
+@Component({
+  selector: 'app-cloth-icm',
+  imports: [
+    MatButton,
+    MatButtonToggleModule,
+    MatCard,
+    MatCardContent,
+    MatCardHeader,
+    MatCardTitle,
+    TranslatePipe,
+    BabylonCanvas,
+    SliderFieldComponent,
+    Information,
+    MetricsComponent
+  ],
+  templateUrl: './cloth-icm.component.html'
+})
+export class ClothIcmComponent implements OnDestroy {
+  private currentSceneData: SceneEventData | null = null;
+  private solver: ClothIcmSolver | null = null;
+  private renderer: ClothIcmRenderer | null = null;
+  private beforeRenderObserver: Observer<Scene> | null = null;
+  private afterRenderObserver: Observer<Scene> | null = null;
+  private frameCount = 0;
+
+  readonly maxBvhDepth = MAX_BVH_DEPTH;
+  readonly sceneId = signal<ClothIcmSceneId>('tangled');
+  readonly icmEnabled = signal(true);
+  readonly coneCulling = signal(true);
+  readonly showBvh = signal(false);
+  readonly bvhDepth = signal(4);
+  readonly icmStep = signal(DEFAULT_ICM_STEP);
+
+  private readonly fpsCounter = new FpsCounter();
+  private readonly stats = signal<ClothIcmStats>(emptyClothIcmStats());
+  readonly metrics = computed<Metric[]>(() => [
+    {label: 'METRICS.FPS', value: this.fpsCounter.fps()},
+    {label: 'METRICS.PARTICLES', value: SCENE_SETTINGS[this.sceneId()].gridSize ** 2},
+    {label: 'METRICS.INTERSECTIONS', value: this.stats().intersections},
+    {label: 'METRICS.CONTOURS', value: this.stats().contours},
+    {label: 'METRICS.TESTED_PAIRS', value: this.stats().testedPairs},
+    {label: 'METRICS.CULLED_NODES', value: this.stats().culledNodes},
+    {label: 'METRICS.COLLISION_MS', value: this.stats().collisionMs, unit: 'ms', fractionDigits: 1}
+  ]);
+
+  readonly renderConfig: RenderConfig = {
+    mode: '3D',
+    initialViewSize: 20
+  };
+
+  readonly algoInformation: AlgorithmInformation = {
+    title: 'CLOTH_ICM.EXPLANATION.TITLE',
+    entries: [
+      {
+        name: 'CLOTH_ICM.EXPLANATION.ICM_TITLE',
+        description: 'CLOTH_ICM.EXPLANATION.ICM',
+        link: UrlConstants.ICM_PAPER,
+        linkLabel: 'ACM',
+        translateName: true
+      },
+      {
+        name: 'CLOTH_ICM.EXPLANATION.BVH_TITLE',
+        description: 'CLOTH_ICM.EXPLANATION.BVH',
+        link: UrlConstants.KDOP_WIKI,
+        translateName: true
+      },
+      {
+        name: 'CLOTH_ICM.EXPLANATION.PROXIMITY_TITLE',
+        description: 'CLOTH_ICM.EXPLANATION.PROXIMITY',
+        link: UrlConstants.XPBD_WIKI,
+        linkLabel: 'XPBD',
+        translateName: true
+      }
+    ],
+    disclaimer: 'CLOTH_ICM.EXPLANATION.DISCLAIMER',
+    disclaimerBottom: '',
+    disclaimerListEntry: ['CLOTH_ICM.EXPLANATION.DISCLAIMER_1', 'CLOTH_ICM.EXPLANATION.DISCLAIMER_2', 'CLOTH_ICM.EXPLANATION.DISCLAIMER_3'],
+    howItWorks: {
+      intro: 'CLOTH_ICM.HOW_IT_WORKS.INTRO',
+      items: [
+        'CLOTH_ICM.HOW_IT_WORKS.ITEM_1',
+        'CLOTH_ICM.HOW_IT_WORKS.ITEM_2',
+        'CLOTH_ICM.HOW_IT_WORKS.ITEM_3',
+        'CLOTH_ICM.HOW_IT_WORKS.ITEM_4',
+        'CLOTH_ICM.HOW_IT_WORKS.ITEM_5'
+      ]
+    },
+    sources: SOURCE_PATHS.CLOTH_ICM
+  };
+
+  ngOnDestroy(): void {
+    this.disposeSimulation();
+  }
+
+  onSceneReady(event: SceneEventData): void {
+    this.currentSceneData = event;
+    this.createSimulation();
+  }
+
+  setScene(id: ClothIcmSceneId): void {
+    this.sceneId.set(id);
+    this.createSimulation();
+  }
+
+  toggleIcm(): void {
+    this.icmEnabled.update(value => !value);
+  }
+
+  toggleCulling(): void {
+    this.coneCulling.update(value => !value);
+  }
+
+  toggleBvh(): void {
+    this.showBvh.update(value => !value);
+  }
+
+  restartSimulation(): void {
+    this.createSimulation();
+  }
+
+  private createSimulation(): void {
+    if (!this.currentSceneData) {
+      return;
+    }
+    const {engine, scene} = this.currentSceneData;
+    this.disposeSimulation();
+
+    this.solver = new ClothIcmSolver(this.clothConfig(), this.sceneId());
+    this.renderer = new ClothIcmRenderer();
+    this.renderer.init(scene, this.solver.topology, this.solver.scene);
+    this.frameCount = 0;
+
+    this.beforeRenderObserver = scene.onBeforeRenderObservable.add(() => {
+      this.stepFrame(engine.getDeltaTime() / 1000);
+    });
+    this.afterRenderObserver = scene.onAfterRenderObservable.add(() => {
+      this.fpsCounter.tick();
+    });
+  }
+
+  private stepFrame(deltaSeconds: number): void {
+    if (!this.solver || !this.renderer) {
+      return;
+    }
+    this.solver.bvh.coneCullingEnabled = this.coneCulling();
+    this.solver.step({dt: deltaSeconds, icmEnabled: this.icmEnabled(), icmStep: this.icmStep()});
+    this.renderer.update(this.solver, this.showBvh(), this.bvhDepth());
+    this.frameCount++;
+    if (this.frameCount % STATS_EVERY_FRAMES === 0) {
+      this.stats.set(this.solver.stats);
+    }
+  }
+
+  private disposeSimulation(): void {
+    const scene = this.currentSceneData?.scene;
+    if (scene && this.beforeRenderObserver) {
+      scene.onBeforeRenderObservable.remove(this.beforeRenderObserver);
+    }
+    if (scene && this.afterRenderObserver) {
+      scene.onAfterRenderObservable.remove(this.afterRenderObserver);
+    }
+    this.beforeRenderObserver = null;
+    this.afterRenderObserver = null;
+    this.renderer?.dispose();
+    this.renderer = null;
+    this.solver = null;
+    this.fpsCounter.reset();
+    this.stats.set(emptyClothIcmStats());
+  }
+
+  private clothConfig(): ClothIcmConfig {
+    const {gridSize, spacing, thickness, proximityRedetectEvery} = SCENE_SETTINGS[this.sceneId()];
+    return {
+      gridWidth: gridSize,
+      gridHeight: gridSize,
+      spacing,
+      particleInvMass: 1 / (spacing * spacing),
+      thickness,
+      proximityRedetectEvery
+    };
+  }
+}

src/app/pages/algorithms/cloth-icm/cloth-icm.model.ts

@@ -0,0 +1,54 @@
+export interface ClothIcmConfig {
+  gridWidth: number;
+  gridHeight: number;
+  spacing: number;
+  particleInvMass: number;
+  /** Minimum distance kept by the proximity constraints. */
+  thickness: number;
+  /** Re-detect proximity pairs every this many substeps, 0 = only at frame start. */
+  proximityRedetectEvery: number;
+}
+
+export interface ClothIcmStats {
+  intersections: number;
+  contours: number;
+  testedPairs: number;
+  culledNodes: number;
+  collisionMs: number;
+}
+
+export type ClothIcmSceneId = 'tangled' | 'drop';
+
+export const ICM_SUBSTEPS = 10;
+export const VELOCITY_DAMPING = 0.98;
+/** The simulation always steps with this dt: damping and ICM step sizes are per step, so the result must not depend on the refresh rate. */
+export const SIM_DT = 1 / 60;
+/** Longest frame time that is simulated; slower frames lose time instead of piling up steps. */
+export const MAX_DT = 1 / 30;
+/** Speed limit in m/s: fast cloth tunnels through the thin proximity shell and creates new intersections. */
+export const MAX_SPEED = 3;
+/** Stabilizes the ICM displacement of short contours (g0). */
+export const ICM_G0 = 4;
+export const DISTANCE_COMPLIANCE = 1e-5;
+/** Share of the tangential motion removed per frame for vertices that rest on the sphere or the ground. */
+export const CONTACT_FRICTION = 0.8;
+
+export function emptyClothIcmStats(): ClothIcmStats {
+  return {intersections: 0, contours: 0, testedPairs: 0, culledNodes: 0, collisionMs: 0};
+}
+
+export interface SceneSettings {
+  gridSize: number;
+  spacing: number;
+  thickness: number;
+  proximityRedetectEvery: number;
+}
+
+/**
+ * Drop: a 35x35 cloth of the same size (wider spacing), gap 3.5x the tangled gap and one mid-frame re-detection:
+ * ICM on ends with ~0 hits after landing, ICM off stays bounded. A larger gap or a finer grid makes ICM off blow up.
+ */
+export const SCENE_SETTINGS: Record<ClothIcmSceneId, SceneSettings> = {
+  tangled: {gridSize: 50, spacing: 0.1, thickness: 0.01, proximityRedetectEvery: 0},
+  drop: {gridSize: 35, spacing: 0.145, thickness: 0.035, proximityRedetectEvery: 5}
+};

src/app/pages/algorithms/cloth-icm/collision/cloth-bvh.spec.ts

@@ -0,0 +1,115 @@
+import {initialPositions} from '../cloth-icm-scenes';
+import {ClothBvh} from './cloth-bvh';
+import {buildClothTopology} from './cloth-topology';
+import {DOP_AXES, DOP_SIZE} from './dop18';
+
+const SPACING = 0.1;
+
+function flatPositions(gridWidth: number, gridHeight: number): {positions: Float32Array; bvh: ClothBvh} {
+  const topology = buildClothTopology(gridWidth, gridHeight);
+  const positions = new Float32Array(topology.vertexCount * 4);
+  for (let vertex = 0; vertex < topology.vertexCount; vertex++) {
+    positions[vertex * 4] = (vertex % gridWidth) * SPACING;
+    positions[vertex * 4 + 1] = Math.floor(vertex / gridWidth) * SPACING;
+    positions[vertex * 4 + 3] = 1;
+  }
+  const bvh = new ClothBvh(topology);
+  bvh.refit(positions);
+  return {positions, bvh};
+}
+
+function tangled(gridSize: number): {positions: Float32Array; bvh: ClothBvh} {
+  const topology = buildClothTopology(gridSize, gridSize);
+  const positions = initialPositions('tangled', topology, SPACING, 1);
+  const bvh = new ClothBvh(topology);
+  bvh.refit(positions);
+  return {positions, bvh};
+}
+
+function countLeafPairs(bvh: ClothBvh, positions: Float32Array, margin = 0): number {
+  let pairs = 0;
+  bvh.selfPairs(positions, margin, () => {
+    pairs++;
+  });
+  return pairs;
+}
+
+describe('ClothBvh', () => {
+  it('builds a tree with one leaf per quad and children after their parent', () => {
+    const {bvh} = flatPositions(5, 4);
+    expect(bvh.nodeCount).toBe(2 * 4 * 3 - 1);
+    for (let node = 0; node < bvh.nodeCount; node++) {
+      if (!bvh.isLeaf(node)) {
+        expect(bvh.left[node]).toBeGreaterThan(node);
+        expect(bvh.right[node]).toBeGreaterThan(node);
+      }
+    }
+  });
+
+  it('refits bounds that contain every vertex of every leaf', () => {
+    const {positions, bvh} = tangled(8);
+    for (let node = 0; node < bvh.nodeCount; node++) {
+      if (!bvh.isLeaf(node)) {
+        continue;
+      }
+      const x0 = bvh.quadRange[node * 4];
+      const y0 = bvh.quadRange[node * 4 + 1];
+      const corners = [y0 * 8 + x0, y0 * 8 + x0 + 1, (y0 + 1) * 8 + x0, (y0 + 1) * 8 + x0 + 1];
+      for (const corner of corners) {
+        DOP_AXES.forEach((axis, index) => {
+          const projection = positions[corner * 4] * axis[0] + positions[corner * 4 + 1] * axis[1] + positions[corner * 4 + 2] * axis[2];
+          expect(projection).toBeGreaterThanOrEqual(bvh.bounds[node * DOP_SIZE + index] - 1e-5);
+          expect(projection).toBeLessThanOrEqual(bvh.bounds[node * DOP_SIZE + 9 + index] + 1e-5);
+        });
+      }
+    }
+  });
+
+  it('refits every inner node so it contains both child bounds', () => {
+    const {bvh} = tangled(8);
+    for (let node = 0; node < bvh.nodeCount; node++) {
+      if (bvh.isLeaf(node)) {
+        continue;
+      }
+      for (const child of [bvh.left[node], bvh.right[node]]) {
+        for (let axis = 0; axis < 9; axis++) {
+          expect(bvh.bounds[node * DOP_SIZE + axis]).toBeLessThanOrEqual(bvh.bounds[child * DOP_SIZE + axis]);
+          expect(bvh.bounds[node * DOP_SIZE + 9 + axis]).toBeGreaterThanOrEqual(bvh.bounds[child * DOP_SIZE + 9 + axis]);
+        }
+      }
+    }
+  });
+
+  it('culls the root of a flat cloth and reports no leaf pairs', () => {
+    const {positions, bvh} = flatPositions(10, 10);
+    let pairs = 0;
+    const culled = bvh.selfPairs(positions, 0, () => {
+      pairs++;
+    });
+    expect(bvh.culled[0]).toBe(1);
+    expect(culled).toBe(1);
+    expect(pairs).toBe(0);
+  });
+
+  it('does not cull the root of a cloth folded through itself and finds leaf pairs', () => {
+    const {positions, bvh} = tangled(30);
+    const pairs = countLeafPairs(bvh, positions);
+    expect(bvh.culled[0]).toBe(0);
+    expect(pairs).toBeGreaterThan(0);
+  });
+
+  it('tests fewer leaf pairs with the cone culling than without', () => {
+    const {positions, bvh} = flatPositions(10, 10);
+    const withCulling = countLeafPairs(bvh, positions);
+    bvh.coneCullingEnabled = false;
+    const withoutCulling = countLeafPairs(bvh, positions);
+    expect(bvh.culled.some(flag => flag === 1)).toBe(false);
+    expect(withCulling).toBeLessThan(withoutCulling);
+  });
+
+  it('finds more pairs with a bigger margin', () => {
+    const {positions, bvh} = flatPositions(10, 10);
+    bvh.coneCullingEnabled = false;
+    expect(countLeafPairs(bvh, positions, SPACING)).toBeGreaterThan(countLeafPairs(bvh, positions, 0));
+  });
+});

src/app/pages/algorithms/cloth-icm/collision/cloth-bvh.ts

@@ -0,0 +1,217 @@
+import {ClothTopology} from './cloth-topology';
+import {DOP_SIZE, dopFromTriangle, dopMerge, dopOverlap} from './dop18';
+import {borderProjectionIsSimple, computeTriangleNormals, coneMergeInto} from './normal-cone';
+
+/** A region whose normal cone is at least this wide (minus a tolerance) is never a height field. */
+const MAX_FLAT_CONE_HALF = Math.PI / 2 - 1e-3;
+
+/**
+ * Bounding volume hierarchy over the quads of the cloth grid. The tree shape comes from the grid and never changes,
+ * only the volumes are refit each frame. Children always have larger ids than their parent.
+ */
+export class ClothBvh {
+  readonly nodeCount: number;
+  /** -1 on leaves. */
+  readonly left: Int32Array;
+  readonly right: Int32Array;
+  readonly parent: Int32Array;
+  /** x0, y0, x1, y1 per node, in quads, end exclusive. */
+  readonly quadRange: Int32Array;
+  readonly depth: Uint8Array;
+  /** 18 per node. */
+  readonly bounds: Float32Array;
+  readonly coneAxis: Float32Array;
+  readonly coneHalf: Float32Array;
+  /** 3 per triangle, filled in refit. */
+  readonly triangleNormals: Float32Array;
+  /** Per frame: 1 = the self-test of this node was culled. Descendants are never visited. */
+  readonly culled: Uint8Array;
+  /** Switch for the Volino 1994 culling; off = every node is tested. */
+  coneCullingEnabled = true;
+
+  private readonly quadsX: number;
+  private readonly triangleScratch = new Float32Array(DOP_SIZE);
+  private readonly borderScratch: Uint32Array;
+  private readonly axisScratch = new Float32Array(3);
+  private nextNode = 0;
+
+  private positions: Float32Array = new Float32Array(0);
+  private margin = 0;
+  private onLeafPair: (quadA: number, quadB: number) => void = () => undefined;
+  private culledCount = 0;
+
+  constructor(private readonly topology: ClothTopology) {
+    const quadsX = topology.gridWidth - 1;
+    const quadsY = topology.gridHeight - 1;
+    this.quadsX = quadsX;
+    this.nodeCount = 2 * quadsX * quadsY - 1;
+    this.left = new Int32Array(this.nodeCount).fill(-1);
+    this.right = new Int32Array(this.nodeCount).fill(-1);
+    this.parent = new Int32Array(this.nodeCount).fill(-1);
+    this.quadRange = new Int32Array(this.nodeCount * 4);
+    this.depth = new Uint8Array(this.nodeCount);
+    this.bounds = new Float32Array(this.nodeCount * DOP_SIZE);
+    this.coneAxis = new Float32Array(this.nodeCount * 3);
+    this.coneHalf = new Float32Array(this.nodeCount);
+    this.triangleNormals = new Float32Array(topology.triangles.length);
+    this.culled = new Uint8Array(this.nodeCount);
+    this.borderScratch = new Uint32Array(2 * (quadsX + quadsY));
+    this.build(0, 0, quadsX, quadsY, 0, -1);
+  }
+
+  isLeaf(node: number): boolean {
+    return this.left[node] < 0;
+  }
+
+  /** Bottom-up: children have larger ids than parents, so one reverse pass is enough. */
+  refit(positions: Float32Array): void {
+    computeTriangleNormals(positions, this.topology.triangles, this.triangleNormals);
+    for (let node = this.nodeCount - 1; node >= 0; node--) {
+      if (this.isLeaf(node)) {
+        this.refitLeaf(positions, node);
+      } else {
+        this.refitInner(node);
+      }
+    }
+  }
+
+  /**
+   * Self-traversal. Calls `onLeafPair` once for every pair of different leaves whose volumes overlap within `margin`.
+   * Returns the number of nodes whose self-test was culled.
+   */
+  selfPairs(positions: Float32Array, margin: number, onLeafPair: (quadA: number, quadB: number) => void): number {
+    this.positions = positions;
+    this.margin = margin;
+    this.onLeafPair = onLeafPair;
+    this.culledCount = 0;
+    this.culled.fill(0);
+    this.selfNode(0);
+    return this.culledCount;
+  }
+
+  private build(x0: number, y0: number, x1: number, y1: number, depth: number, parent: number): void {
+    const node = this.nextNode++;
+    this.parent[node] = parent;
+    this.depth[node] = depth;
+    this.quadRange.set([x0, y0, x1, y1], node * 4);
+    const width = x1 - x0;
+    const height = y1 - y0;
+    if (width === 1 && height === 1) {
+      return;
+    }
+
+    this.left[node] = this.nextNode;
+    if (width >= height) {
+      const middle = x0 + Math.floor(width / 2);
+      this.build(x0, y0, middle, y1, depth + 1, node);
+      this.right[node] = this.nextNode;
+      this.build(middle, y0, x1, y1, depth + 1, node);
+    } else {
+      const middle = y0 + Math.floor(height / 2);
+      this.build(x0, y0, x1, middle, depth + 1, node);
+      this.right[node] = this.nextNode;
+      this.build(x0, middle, x1, y1, depth + 1, node);
+    }
+  }
+
+  private leafQuad(node: number): number {
+    return this.quadRange[node * 4 + 1] * this.quadsX + this.quadRange[node * 4];
+  }
+
+  private refitLeaf(positions: Float32Array, node: number): void {
+    const quad = this.leafQuad(node);
+    const offset = node * DOP_SIZE;
+    dopFromTriangle(positions, this.topology, quad * 2, this.bounds, offset);
+    dopFromTriangle(positions, this.topology, quad * 2 + 1, this.triangleScratch, 0);
+    dopMerge(this.bounds, offset, this.triangleScratch, 0, this.bounds, offset);
+    this.coneHalf[node] = coneMergeInto(this.triangleNormals, quad * 6, 0, this.triangleNormals, quad * 6 + 3, 0, this.coneAxis, node * 3);
+  }
+
+  private refitInner(node: number): void {
+    const left = this.left[node];
+    const right = this.right[node];
+    dopMerge(this.bounds, left * DOP_SIZE, this.bounds, right * DOP_SIZE, this.bounds, node * DOP_SIZE);
+    this.coneHalf[node] = coneMergeInto(
+      this.coneAxis, left * 3, this.coneHalf[left],
+      this.coneAxis, right * 3, this.coneHalf[right],
+      this.coneAxis, node * 3
+    );
+  }
+
+  /** Culling is per node, never per triangle: the region must be a height field over its own border. */
+  private selfTestIsCulled(node: number): boolean {
+    if (!this.coneCullingEnabled || this.coneHalf[node] >= MAX_FLAT_CONE_HALF) {
+      return false;
+    }
+    const borderLength = this.fillBorderVertices(node);
+    this.axisScratch[0] = this.coneAxis[node * 3];
+    this.axisScratch[1] = this.coneAxis[node * 3 + 1];
+    this.axisScratch[2] = this.coneAxis[node * 3 + 2];
+    return borderProjectionIsSimple(this.positions, this.borderScratch, this.axisScratch, borderLength);
+  }
+
+  /** Writes the closed border loop of the node into `borderScratch` and returns its length. */
+  private fillBorderVertices(node: number): number {
+    const x0 = this.quadRange[node * 4];
+    const y0 = this.quadRange[node * 4 + 1];
+    const x1 = this.quadRange[node * 4 + 2];
+    const y1 = this.quadRange[node * 4 + 3];
+    const gridWidth = this.topology.gridWidth;
+    const border = this.borderScratch;
+    let length = 0;
+    for (let x = x0; x < x1; x++) {
+      border[length++] = y0 * gridWidth + x;
+    }
+    for (let y = y0; y < y1; y++) {
+      border[length++] = y * gridWidth + x1;
+    }
+    for (let x = x1; x > x0; x--) {
+      border[length++] = y1 * gridWidth + x;
+    }
+    for (let y = y1; y > y0; y--) {
+      border[length++] = y * gridWidth + x0;
+    }
+    return length;
+  }
+
+  private selfNode(node: number): void {
+    if (this.isLeaf(node)) {
+      return;
+    }
+    if (this.selfTestIsCulled(node)) {
+      this.culled[node] = 1;
+      this.culledCount++;
+      return;
+    }
+    const left = this.left[node];
+    const right = this.right[node];
+    this.selfNode(left);
+    this.selfNode(right);
+    this.pairNodes(left, right);
+  }
+
+  private pairNodes(a: number, b: number): void {
+    if (!dopOverlap(this.bounds, a * DOP_SIZE, this.bounds, b * DOP_SIZE, this.margin)) {
+      return;
+    }
+    const aIsLeaf = this.isLeaf(a);
+    const bIsLeaf = this.isLeaf(b);
+    if (aIsLeaf && bIsLeaf) {
+      this.onLeafPair(this.leafQuad(a), this.leafQuad(b));
+      return;
+    }
+    const splitA = !aIsLeaf && (bIsLeaf || this.quadArea(a) >= this.quadArea(b));
+    if (splitA) {
+      this.pairNodes(this.left[a], b);
+      this.pairNodes(this.right[a], b);
+    } else {
+      this.pairNodes(a, this.left[b]);
+      this.pairNodes(a, this.right[b]);
+    }
+  }
+
+  private quadArea(node: number): number {
+    const base = node * 4;
+    return (this.quadRange[base + 2] - this.quadRange[base]) * (this.quadRange[base + 3] - this.quadRange[base + 1]);
+  }
+}

src/app/pages/algorithms/cloth-icm/collision/cloth-topology.spec.ts

@@ -0,0 +1,51 @@
+import {buildClothTopology, buildDistanceConstraints, trianglesShareVertex} from './cloth-topology';
+
+describe('cloth topology', () => {
+  const topology = buildClothTopology(4, 3);
+
+  it('has the expected counts', () => {
+    const quads = 3 * 2;
+    expect(topology.vertexCount).toBe(12);
+    expect(topology.triangles.length / 3).toBe(quads * 2);
+    // horizontal 3*3 + vertical 4*2 + diagonal 6
+    expect(topology.edges.length / 2).toBe(9 + 8 + 6);
+  });
+
+  it('gives border edges one triangle and inner edges two', () => {
+    const edgeCount = topology.edges.length / 2;
+    let borderEdges = 0;
+    for (let edge = 0; edge < edgeCount; edge++) {
+      expect(topology.edgeTriangles[edge * 2]).toBeGreaterThanOrEqual(0);
+      if (topology.edgeTriangles[edge * 2 + 1] === -1) {
+        borderEdges++;
+      }
+    }
+    expect(borderEdges).toBe(2 * (3 + 2));
+  });
+
+  it('lets both triangles of a quad share the diagonal edge', () => {
+    const first = Array.from(topology.triangleEdges.slice(0, 3));
+    const second = Array.from(topology.triangleEdges.slice(3, 6));
+    const shared = first.filter(edge => second.includes(edge));
+    expect(shared.length).toBe(1);
+    const edge = shared[0];
+    expect(topology.edgeTriangles[edge * 2]).toBe(0);
+    expect(topology.edgeTriangles[edge * 2 + 1]).toBe(1);
+  });
+
+  it('detects triangles that share a vertex', () => {
+    expect(trianglesShareVertex(topology, 0, 1)).toBe(true);
+    expect(trianglesShareVertex(topology, 0, 2)).toBe(true);
+    expect(trianglesShareVertex(topology, 0, 11)).toBe(false);
+  });
+
+  it('writes distance constraints as [a, b, rest, 1] with the second diagonal added', () => {
+    const constraints = buildDistanceConstraints(topology, 0.5);
+    expect(constraints.length / 4).toBe(23 + 6);
+    expect(Array.from(constraints.slice(0, 4))).toEqual([topology.edges[0], topology.edges[1], 0.5, 1]);
+    const diagonalRest = 0.5 * Math.SQRT2;
+    const last = constraints.length - 4;
+    expect(constraints[last + 2]).toBeCloseTo(diagonalRest);
+    expect(constraints[last + 3]).toBe(1);
+  });
+});

src/app/pages/algorithms/cloth-icm/collision/cloth-topology.ts

@@ -0,0 +1,103 @@
+export interface ClothTopology {
+  gridWidth: number;
+  gridHeight: number;
+  vertexCount: number;
+  /** 3 ids per triangle, same winding for all. */
+  triangles: Uint32Array;
+  /** 2 ids per edge (horizontal, vertical, one diagonal). */
+  edges: Uint32Array;
+  /** 2 per edge, -1 = border. */
+  edgeTriangles: Int32Array;
+  /** 3 per triangle. */
+  triangleEdges: Uint32Array;
+}
+
+/** Vertex id = y * gridWidth + x. Quad q = y * (gridWidth - 1) + x has triangles 2q and 2q + 1. */
+export function buildClothTopology(gridWidth: number, gridHeight: number): ClothTopology {
+  const vertexCount = gridWidth * gridHeight;
+  const triangles = buildTriangles(gridWidth, gridHeight);
+  const triangleCount = triangles.length / 3;
+  const edgeIndexByKey = new Map<number, number>();
+  const edgeList: number[] = [];
+  const edgeTriangleList: number[] = [];
+  const triangleEdges = new Uint32Array(triangleCount * 3);
+
+  for (let triangle = 0; triangle < triangleCount; triangle++) {
+    for (let corner = 0; corner < 3; corner++) {
+      const a = triangles[triangle * 3 + corner];
+      const b = triangles[triangle * 3 + (corner + 1) % 3];
+      const key = Math.min(a, b) * vertexCount + Math.max(a, b);
+      let edge = edgeIndexByKey.get(key);
+      if (edge === undefined) {
+        edge = edgeList.length / 2;
+        edgeIndexByKey.set(key, edge);
+        edgeList.push(a, b);
+        edgeTriangleList.push(triangle, -1);
+      } else {
+        edgeTriangleList[edge * 2 + 1] = triangle;
+      }
+      triangleEdges[triangle * 3 + corner] = edge;
+    }
+  }
+
+  return {
+    gridWidth,
+    gridHeight,
+    vertexCount,
+    triangles,
+    edges: Uint32Array.from(edgeList),
+    edgeTriangles: Int32Array.from(edgeTriangleList),
+    triangleEdges
+  };
+}
+
+function buildTriangles(gridWidth: number, gridHeight: number): Uint32Array {
+  const quadsX = gridWidth - 1;
+  const quadsY = gridHeight - 1;
+  const triangles = new Uint32Array(quadsX * quadsY * 6);
+  for (let y = 0; y < quadsY; y++) {
+    for (let x = 0; x < quadsX; x++) {
+      const quad = y * quadsX + x;
+      const v00 = y * gridWidth + x;
+      const v10 = v00 + 1;
+      const v01 = v00 + gridWidth;
+      const v11 = v01 + 1;
+      triangles.set([v00, v01, v11, v00, v11, v10], quad * 6);
+    }
+  }
+  return triangles;
+}
+
+export function trianglesShareVertex(topology: ClothTopology, a: number, b: number): boolean {
+  const {triangles} = topology;
+  for (let i = 0; i < 3; i++) {
+    const vertex = triangles[a * 3 + i];
+    for (let j = 0; j < 3; j++) {
+      if (vertex === triangles[b * 3 + j]) {
+        return true;
+      }
+    }
+  }
+  return false;
+}
+
+/** Layout of `solveConstraints`: [a, b, restLength, type] per edge, plus the second diagonal of every quad. */
+export function buildDistanceConstraints(topology: ClothTopology, spacing: number): Float32Array {
+  const {gridWidth, gridHeight, edges} = topology;
+  const constraints: number[] = [];
+  const addConstraint = (a: number, b: number): void => {
+    const dx = (a % gridWidth) - (b % gridWidth);
+    const dy = Math.floor(a / gridWidth) - Math.floor(b / gridWidth);
+    constraints.push(a, b, spacing * Math.hypot(dx, dy), 1);
+  };
+
+  for (let edge = 0; edge < edges.length / 2; edge++) {
+    addConstraint(edges[edge * 2], edges[edge * 2 + 1]);
+  }
+  for (let y = 0; y < gridHeight - 1; y++) {
+    for (let x = 0; x < gridWidth - 1; x++) {
+      addConstraint(y * gridWidth + x + 1, (y + 1) * gridWidth + x);
+    }
+  }
+  return Float32Array.from(constraints);
+}

src/app/pages/algorithms/cloth-icm/collision/collision-detection.ts

@@ -0,0 +1,50 @@
+import {ClothBvh} from './cloth-bvh';
+import {ClothTopology, trianglesShareVertex} from './cloth-topology';
+import {DOP_SIZE, dopFromTriangle, dopOverlap} from './dop18';
+import {IntersectionSet} from './intersection-set';
+
+export interface DetectionCounters {
+  /** Triangle pairs that went to the edge-triangle test. */
+  testedPairs: number;
+  culledNodes: number;
+}
+
+const triangleDops = new Float32Array(DOP_SIZE * 2);
+
+/**
+ * Broad phase (BVH self-traversal with normal-cone culling) and narrow phase (edge-triangle).
+ * Fills `set` with the hits and `proximityCandidates` with triangle pairs (2 ids each) that do not
+ * intersect but are closer than `margin`. Expects a refit BVH.
+ */
+export function detectCollisions(
+  bvh: ClothBvh, positions: Float32Array, topology: ClothTopology, margin: number,
+  set: IntersectionSet, proximityCandidates: number[]
+): DetectionCounters {
+  set.clear();
+  proximityCandidates.length = 0;
+  let testedPairs = 0;
+
+  const testTrianglePair = (triA: number, triB: number): void => {
+    if (trianglesShareVertex(topology, triA, triB)) {
+      return;
+    }
+    testedPairs++;
+    const hits = set.addTrianglePair(positions, topology, triA, triB);
+    if (hits > 0) {
+      return;
+    }
+    dopFromTriangle(positions, topology, triA, triangleDops, 0);
+    dopFromTriangle(positions, topology, triB, triangleDops, DOP_SIZE);
+    if (dopOverlap(triangleDops, 0, triangleDops, DOP_SIZE, margin)) {
+      proximityCandidates.push(triA, triB);
+    }
+  };
+
+  const culledNodes = bvh.selfPairs(positions, margin, (quadA, quadB) => {
+    testTrianglePair(quadA * 2, quadB * 2);
+    testTrianglePair(quadA * 2, quadB * 2 + 1);
+    testTrianglePair(quadA * 2 + 1, quadB * 2);
+    testTrianglePair(quadA * 2 + 1, quadB * 2 + 1);
+  });
+  return {testedPairs, culledNodes};
+}

src/app/pages/algorithms/cloth-icm/collision/dop18.spec.ts

@@ -0,0 +1,69 @@
+import {buildClothTopology} from './cloth-topology';
+import {DOP_SIZE, dopFromTriangle, dopMerge, dopOverlap} from './dop18';
+
+/** Vertex layout: [x, y, z, invMass]. */
+function vertices(...rows: number[][]): Float32Array {
+  return new Float32Array(rows.flat());
+}
+
+describe('18-DOP', () => {
+  const topology = buildClothTopology(2, 2);
+  const positions = vertices([0, 0, 0, 1], [1, 0, 0, 1], [0, 1, 0, 1], [1, 1, 0, 1]);
+
+  function triangleDop(offset: number, triangle: number, shift: number): Float32Array {
+    const moved = positions.slice();
+    for (let vertex = 0; vertex < 4; vertex++) {
+      moved[vertex * 4 + 2] += shift;
+    }
+    const out = new Float32Array(DOP_SIZE + offset);
+    dopFromTriangle(moved, topology, triangle, out, offset);
+    return out;
+  }
+
+  it('bounds the three vertices on every axis', () => {
+    const dop = triangleDop(0, 0, 0);
+    // triangle 0 = vertices 0, 2, 3: (0,0,0), (0,1,0), (1,1,0); axis 3 = x + y
+    expect(dop[0]).toBe(0);
+    expect(dop[9]).toBe(1);
+    expect(dop[3]).toBe(0);
+    expect(dop[12]).toBe(2);
+  });
+
+  it('overlaps with a margin but not without', () => {
+    const a = triangleDop(0, 0, 0);
+    expect(dopOverlap(a, 0, triangleDop(0, 0, 0), 0, 0)).toBe(true);
+    const apart = triangleDop(0, 0, 0.5);
+    expect(dopOverlap(a, 0, apart, 0, 0)).toBe(false);
+    expect(dopOverlap(a, 0, apart, 0, 0.5)).toBe(true);
+    expect(dopOverlap(a, 0, triangleDop(0, 0, 2), 0, 0.5)).toBe(false);
+  });
+
+  it('merges to a volume that contains both parts, in place and with offsets', () => {
+    const a = triangleDop(0, 0, 0);
+    const b = triangleDop(DOP_SIZE, 1, 3);
+    const merged = new Float32Array(DOP_SIZE);
+    dopMerge(a, 0, b, DOP_SIZE, merged, 0);
+    expect(merged[2]).toBe(0);
+    expect(merged[11]).toBe(3);
+    expect(dopOverlap(merged, 0, a, 0, 0)).toBe(true);
+    expect(dopOverlap(merged, 0, b, DOP_SIZE, 0)).toBe(true);
+  });
+
+  it('overlapsWithinMarginAlongADiagonalAxis', () => {
+    const flat = vertices([0, 0, 0, 1], [1, 0, 0, 1], [0, 1, 0, 1], [1, 1, 0, 1]);
+    const moved = flat.slice();
+    for (let vertex = 0; vertex < 4; vertex++) {
+      moved[vertex * 4] += 1.45;
+      moved[vertex * 4 + 1] += 1.45;
+    }
+    const a = new Float32Array(DOP_SIZE);
+    const b = new Float32Array(DOP_SIZE);
+    // Triangle 0 = (0,0,0), (0,1,0), (1,1,0). The moved copy starts at (1.45,1.45,0): closest gap is 0.45 * sqrt 2 = 0.64.
+    dopFromTriangle(flat, topology, 0, a, 0);
+    dopFromTriangle(moved, topology, 0, b, 0);
+    const distance = Math.hypot(0.45, 0.45);
+    const margin = 0.7;
+    expect(distance).toBeLessThan(margin);
+    expect(dopOverlap(a, 0, b, 0, margin)).toBe(true);
+  });
+});

src/app/pages/algorithms/cloth-icm/collision/dop18.ts

@@ -0,0 +1,59 @@
+import {ClothTopology} from './cloth-topology';
+
+export const DOP_AXIS_COUNT = 9;
+export const DOP_SIZE = DOP_AXIS_COUNT * 2;
+
+/** Not normalized: the slab bounds are plain dot products with these vectors. */
+export const DOP_AXES: readonly (readonly [number, number, number])[] = [
+  [1, 0, 0], [0, 1, 0], [0, 0, 1],
+  [1, 1, 0], [1, -1, 0],
+  [1, 0, 1], [1, 0, -1],
+  [0, 1, 1], [0, 1, -1]
+];
+
+/** Layout per node: [min0..min8, max0..max8]. Positions use the vec4 layout [x, y, z, invMass]. */
+export function dopFromTriangle(positions: Float32Array, topology: ClothTopology, triangle: number, out: Float32Array, offset: number): void {
+  for (let axis = 0; axis < DOP_AXIS_COUNT; axis++) {
+    out[offset + axis] = Infinity;
+    out[offset + DOP_AXIS_COUNT + axis] = -Infinity;
+  }
+  for (let corner = 0; corner < 3; corner++) {
+    const base = topology.triangles[triangle * 3 + corner] * 4;
+    const x = positions[base];
+    const y = positions[base + 1];
+    const z = positions[base + 2];
+    for (let axis = 0; axis < DOP_AXIS_COUNT; axis++) {
+      const direction = DOP_AXES[axis];
+      const projection = x * direction[0] + y * direction[1] + z * direction[2];
+      if (projection < out[offset + axis]) {
+        out[offset + axis] = projection;
+      }
+      if (projection > out[offset + DOP_AXIS_COUNT + axis]) {
+        out[offset + DOP_AXIS_COUNT + axis] = projection;
+      }
+    }
+  }
+}
+
+/** `out` may be the same array and offset as `a` or `b`. */
+export function dopMerge(a: Float32Array, offsetA: number, b: Float32Array, offsetB: number, out: Float32Array, offsetOut: number): void {
+  for (let axis = 0; axis < DOP_AXIS_COUNT; axis++) {
+    out[offsetOut + axis] = Math.min(a[offsetA + axis], b[offsetB + axis]);
+    const max = DOP_AXIS_COUNT + axis;
+    out[offsetOut + max] = Math.max(a[offsetA + max], b[offsetB + max]);
+  }
+}
+
+const AXIS_LENGTHS: readonly number[] = DOP_AXES.map(([x, y, z]) => Math.hypot(x, y, z));
+
+/** True when the slabs overlap on every axis; `margin` is a distance, scaled by the length of each (unnormalized) axis. */
+export function dopOverlap(a: Float32Array, offsetA: number, b: Float32Array, offsetB: number, margin: number): boolean {
+  for (let axis = 0; axis < DOP_AXIS_COUNT; axis++) {
+    const max = DOP_AXIS_COUNT + axis;
+    const axisMargin = margin * AXIS_LENGTHS[axis];
+    if (a[offsetA + axis] > b[offsetB + max] + axisMargin || b[offsetB + axis] > a[offsetA + max] + axisMargin) {
+      return false;
+    }
+  }
+  return true;
+}

src/app/pages/algorithms/cloth-icm/collision/edge-triangle.spec.ts

@@ -0,0 +1,48 @@
+import {EdgeTriangleHit, intersectEdgeTriangle, newEdgeTriangleHit} from './edge-triangle';
+
+/** Vertex layout: [x, y, z, invMass]. Triangle = vertices 0..2, edge = vertices 3 and 4. */
+function scene(edgeStart: number[], edgeEnd: number[]): Float32Array {
+  return new Float32Array([
+    0, 0, 0, 1,
+    4, 0, 0, 1,
+    0, 0, 4, 1,
+    ...edgeStart, 1,
+    ...edgeEnd, 1
+  ]);
+}
+
+function hit(edgeStart: number[], edgeEnd: number[]): {hit: boolean; out: EdgeTriangleHit} {
+  const out = newEdgeTriangleHit();
+  return {hit: intersectEdgeTriangle(scene(edgeStart, edgeEnd), 3, 4, 0, 1, 2, out), out};
+}
+
+describe('intersectEdgeTriangle', () => {
+  it('hits and reports the point, the edge parameter and the barycentrics', () => {
+    const result = hit([1, -1, 1], [1, 3, 1]);
+    expect(result.hit).toBe(true);
+    expect(result.out.t).toBeCloseTo(0.25);
+    expect([result.out.x, result.out.y, result.out.z]).toEqual([1, 0, 1]);
+    expect(result.out.u).toBeCloseTo(0.25);
+    expect(result.out.v).toBeCloseTo(0.25);
+  });
+
+  it('misses beside the triangle', () => {
+    expect(hit([3, -1, 3], [3, 1, 3]).hit).toBe(false);
+  });
+
+  it('misses when the segment ends before the plane', () => {
+    expect(hit([1, 1, 1], [1, 3, 1]).hit).toBe(false);
+    expect(hit([1, -3, 1], [1, -1, 1]).hit).toBe(false);
+  });
+
+  it('gives false for an edge parallel to the plane or inside it', () => {
+    expect(hit([1, 1, 1], [2, 1, 1]).hit).toBe(false);
+    expect(hit([1, 0, 1], [2, 0, 1]).hit).toBe(false);
+  });
+
+  it('counts a hit exactly on a vertex', () => {
+    const result = hit([0, -1, 0], [0, 1, 0]);
+    expect(result.hit).toBe(true);
+    expect(result.out.t).toBeCloseTo(0.5);
+  });
+});

src/app/pages/algorithms/cloth-icm/collision/edge-triangle.ts

@@ -0,0 +1,79 @@
+export interface EdgeTriangleHit {
+  /** Position of the hit along the edge, 0 at edgeA and 1 at edgeB. */
+  t: number;
+  u: number;
+  v: number;
+  x: number;
+  y: number;
+  z: number;
+}
+
+export function newEdgeTriangleHit(): EdgeTriangleHit {
+  return {t: 0, u: 0, v: 0, x: 0, y: 0, z: 0};
+}
+
+const EPS = 1e-6;
+const MIN_DETERMINANT = 1e-12;
+
+/**
+ * Moller-Trumbore with a segment range and inclusive borders (EPS).
+ * An edge parallel to the triangle plane or inside it gives false.
+ * Positions use the vec4 layout [x, y, z, invMass]; the other arguments are vertex ids.
+ */
+export function intersectEdgeTriangle(
+  positions: Float32Array, edgeA: number, edgeB: number, triA: number, triB: number, triC: number, out: EdgeTriangleHit
+): boolean {
+  const pa = edgeA * 4;
+  const pb = edgeB * 4;
+  const va = triA * 4;
+  const vb = triB * 4;
+  const vc = triC * 4;
+
+  const dx = positions[pb] - positions[pa];
+  const dy = positions[pb + 1] - positions[pa + 1];
+  const dz = positions[pb + 2] - positions[pa + 2];
+  const e1x = positions[vb] - positions[va];
+  const e1y = positions[vb + 1] - positions[va + 1];
+  const e1z = positions[vb + 2] - positions[va + 2];
+  const e2x = positions[vc] - positions[va];
+  const e2y = positions[vc + 1] - positions[va + 1];
+  const e2z = positions[vc + 2] - positions[va + 2];
+
+  const px = dy * e2z - dz * e2y;
+  const py = dz * e2x - dx * e2z;
+  const pz = dx * e2y - dy * e2x;
+  const determinant = e1x * px + e1y * py + e1z * pz;
+  if (Math.abs(determinant) < MIN_DETERMINANT) {
+    return false;
+  }
+  const inverse = 1 / determinant;
+
+  const tx = positions[pa] - positions[va];
+  const ty = positions[pa + 1] - positions[va + 1];
+  const tz = positions[pa + 2] - positions[va + 2];
+  const u = (tx * px + ty * py + tz * pz) * inverse;
+  if (u < -EPS || u > 1 + EPS) {
+    return false;
+  }
+
+  const qx = ty * e1z - tz * e1y;
+  const qy = tz * e1x - tx * e1z;
+  const qz = tx * e1y - ty * e1x;
+  const v = (dx * qx + dy * qy + dz * qz) * inverse;
+  if (v < -EPS || u + v > 1 + EPS) {
+    return false;
+  }
+
+  const t = (e2x * qx + e2y * qy + e2z * qz) * inverse;
+  if (t < -EPS || t > 1 + EPS) {
+    return false;
+  }
+
+  out.t = t;
+  out.u = u;
+  out.v = v;
+  out.x = positions[pa] + t * dx;
+  out.y = positions[pa + 1] + t * dy;
+  out.z = positions[pa + 2] + t * dz;
+  return true;
+}

src/app/pages/algorithms/cloth-icm/collision/intersection-contours.spec.ts

@@ -0,0 +1,50 @@
+import {traceContours} from './intersection-contours';
+import {IntersectionSet} from './intersection-set';
+
+function setOf(hitCount: number, segments: number[][]): IntersectionSet {
+  const set = new IntersectionSet();
+  set.count = hitCount;
+  set.segments = segments.flat();
+  return set;
+}
+
+describe('traceContours', () => {
+  it('walks one closed loop', () => {
+    const contours = traceContours(setOf(4, [[0, 1], [1, 2], [2, 3], [3, 0]]));
+    expect(contours.length).toBe(1);
+    expect(contours[0].closed).toBe(true);
+    expect([...contours[0].points].sort()).toEqual([0, 1, 2, 3]);
+  });
+
+  it('walks an open path from one border end to the other', () => {
+    const contours = traceContours(setOf(3, [[1, 0], [1, 2]]));
+    expect(contours.length).toBe(1);
+    expect(contours[0].closed).toBe(false);
+    const {points} = contours[0];
+    expect(points.length).toBe(3);
+    expect([points[0], points[2]].sort()).toEqual([0, 2]);
+    expect(points[1]).toBe(1);
+  });
+
+  it('finds two separate loops', () => {
+    const contours = traceContours(setOf(6, [[0, 1], [1, 2], [2, 0], [3, 4], [4, 5], [5, 3]]));
+    expect(contours.length).toBe(2);
+    expect(contours.every(contour => contour.closed)).toBe(true);
+  });
+
+  it('finds a loop next to an open path', () => {
+    const contours = traceContours(setOf(6, [[0, 1], [1, 2], [3, 4], [4, 5], [5, 3]]));
+    expect(contours.map(contour => contour.closed).sort()).toEqual([false, true]);
+  });
+
+  it('terminates on a hit with three segments', () => {
+    const contours = traceContours(setOf(4, [[0, 1], [0, 2], [0, 3]]));
+    const walked = contours.flatMap(contour => contour.points);
+    expect(new Set(walked).size).toBe(walked.length);
+    expect(walked.length).toBe(3);
+  });
+
+  it('returns nothing for an empty set', () => {
+    expect(traceContours(setOf(0, []))).toEqual([]);
+  });
+});

src/app/pages/algorithms/cloth-icm/collision/intersection-contours.ts

@@ -0,0 +1,65 @@
+import {IntersectionSet} from './intersection-set';
+
+export interface Contour {
+  /** Hit indices in walking order. Consecutive points are joined by a segment. */
+  points: number[];
+  /** Closed loops walk back to their first point; the first point is not repeated. */
+  closed: boolean;
+}
+
+/**
+ * Joins the segments of a set into contours. Hits with one segment (cloth border) start open paths,
+ * what is left are closed loops. A hit with three or more segments is walked through once and never loops forever.
+ */
+export function traceContours(set: IntersectionSet): Contour[] {
+  const neighbours = buildNeighbours(set);
+  const visited = new Uint8Array(set.count);
+  const contours: Contour[] = [];
+
+  for (let hit = 0; hit < set.count; hit++) {
+    if (neighbours[hit].length === 1 && !visited[hit]) {
+      contours.push(walk(hit, neighbours, visited, false));
+    }
+  }
+  for (let hit = 0; hit < set.count; hit++) {
+    if (neighbours[hit].length > 0 && !visited[hit]) {
+      contours.push(walk(hit, neighbours, visited, true));
+    }
+  }
+  return contours.filter(contour => contour.points.length > 1);
+}
+
+function buildNeighbours(set: IntersectionSet): number[][] {
+  const neighbours: number[][] = Array.from({length: set.count}, () => []);
+  for (let segment = 0; segment < set.segmentCount; segment++) {
+    const a = set.segments[segment * 2];
+    const b = set.segments[segment * 2 + 1];
+    neighbours[a].push(b);
+    neighbours[b].push(a);
+  }
+  return neighbours;
+}
+
+function walk(start: number, neighbours: number[][], visited: Uint8Array, closedWalk: boolean): Contour {
+  const points = [start];
+  visited[start] = 1;
+  let current = start;
+  let next = firstUnvisited(neighbours[current], visited);
+  while (next !== -1) {
+    points.push(next);
+    visited[next] = 1;
+    current = next;
+    next = firstUnvisited(neighbours[current], visited);
+  }
+  const closesBack = closedWalk && points.length > 2 && neighbours[current].includes(start);
+  return {points, closed: closesBack};
+}
+
+function firstUnvisited(candidates: number[], visited: Uint8Array): number {
+  for (const candidate of candidates) {
+    if (!visited[candidate]) {
+      return candidate;
+    }
+  }
+  return -1;
+}

src/app/pages/algorithms/cloth-icm/collision/intersection-set.spec.ts

@@ -0,0 +1,71 @@
+import {buildClothTopology} from './cloth-topology';
+import {traceContours} from './intersection-contours';
+import {IntersectionSet} from './intersection-set';
+
+/**
+ * Grid 4x2 (vertex id = y * 4 + x). Quad 0 lies flat in z = 0 (triangles 0 and 1).
+ * Quad 2 (triangles 4 and 5) stands on it and crosses z = 0 inside triangle 0.
+ */
+const topology = buildClothTopology(4, 2);
+const positions = new Float32Array([
+  0, 0, 0, 1, /* v0 */ 2, 0, 0, 1, /* v1 */ 0.6, 0.7, -1, 1, /* v2 */ 1.0, 1.3, -1, 1, /* v3 */
+  0, 2, 0, 1, /* v4 */ 2, 2, 0, 1, /* v5 */ 0.6, 0.7, 1, 1, /* v6 */ 1.0, 1.3, 1, 1 /* v7 */
+]);
+
+describe('IntersectionSet', () => {
+  it('skips adjacent triangles', () => {
+    const set = new IntersectionSet();
+    expect(set.addTrianglePair(positions, topology, 0, 1)).toBe(0);
+    expect(set.count).toBe(0);
+  });
+
+  it('turns two hits into one segment', () => {
+    const set = new IntersectionSet();
+    expect(set.addTrianglePair(positions, topology, 0, 4)).toBe(2);
+    expect(set.count).toBe(2);
+    expect(set.segmentCount).toBe(1);
+  });
+
+  it('stores a hit once when both triangles of its edge find it', () => {
+    const set = new IntersectionSet();
+    set.addTrianglePair(positions, topology, 0, 4);
+    set.addTrianglePair(positions, topology, 0, 5);
+    expect(set.count).toBe(3);
+    expect(set.segmentCount).toBe(2);
+  });
+
+  it('records the hit point on the pierced triangle', () => {
+    const set = new IntersectionSet();
+    set.addTrianglePair(positions, topology, 0, 4);
+    const triangles = Array.from(set.triangle);
+    expect(triangles).toEqual([0, 0]);
+    expect(set.points[2]).toBeCloseTo(0);
+    expect(set.points[5]).toBeCloseTo(0);
+  });
+
+  it('gives no segment for triangles that do not touch', () => {
+    const set = new IntersectionSet();
+    expect(set.addTrianglePair(positions, topology, 1, 4)).toBe(0);
+    expect(set.segmentCount).toBe(0);
+  });
+
+  it('feeds traceContours with one open path across two pairs', () => {
+    const set = new IntersectionSet();
+    set.addTrianglePair(positions, topology, 0, 4);
+    set.addTrianglePair(positions, topology, 0, 5);
+    const contours = traceContours(set);
+    expect(contours.length).toBe(1);
+    expect(contours[0].closed).toBe(false);
+    expect(contours[0].points.length).toBe(3);
+  });
+
+  it('clear() resets hits and segments', () => {
+    const set = new IntersectionSet();
+    set.addTrianglePair(positions, topology, 0, 4);
+    set.clear();
+    expect(set.count).toBe(0);
+    expect(set.segmentCount).toBe(0);
+    set.addTrianglePair(positions, topology, 0, 4);
+    expect(set.count).toBe(2);
+  });
+});

src/app/pages/algorithms/cloth-icm/collision/intersection-set.ts

@@ -0,0 +1,105 @@
+import {ClothTopology, trianglesShareVertex} from './cloth-topology';
+import {EdgeTriangleHit, intersectEdgeTriangle, newEdgeTriangleHit} from './edge-triangle';
+
+/**
+ * All edge-triangle hits of one frame and the contour segments between them.
+ * A hit is one edge piercing one triangle. The same hit is found from both triangles that share the edge,
+ * so it is stored once and can be an end of two segments.
+ */
+export class IntersectionSet {
+  count = 0;
+  edge: number[] = [];
+  triangle: number[] = [];
+  t: number[] = [];
+  u: number[] = [];
+  v: number[] = [];
+  /** 3 per hit. */
+  points: number[] = [];
+  /** 2 hit indices per segment. */
+  segments: number[] = [];
+
+  private readonly hitByKey = new Map<number, number>();
+  private readonly scratch: EdgeTriangleHit = newEdgeTriangleHit();
+  private readonly pairHits: {edge: number; triangle: number; hit: EdgeTriangleHit}[] = [];
+
+  clear(): void {
+    this.count = 0;
+    this.edge.length = 0;
+    this.triangle.length = 0;
+    this.t.length = 0;
+    this.u.length = 0;
+    this.v.length = 0;
+    this.points.length = 0;
+    this.segments.length = 0;
+    this.hitByKey.clear();
+  }
+
+  get segmentCount(): number {
+    return this.segments.length / 2;
+  }
+
+  /**
+   * Tests the 3 edges of each triangle against the other one. Adjacent triangles (shared vertex) are skipped.
+   * Exactly 2 hits give one segment. 1 or 3+ hits (a vertex or an edge lies on the other triangle) drop the pair.
+   * Returns the number of raw hits found, 0 = the triangles do not intersect.
+   */
+  addTrianglePair(positions: Float32Array, topology: ClothTopology, triA: number, triB: number): number {
+    if (trianglesShareVertex(topology, triA, triB)) {
+      return 0;
+    }
+    let found = 0;
+    found = this.collectHits(positions, topology, triA, triB, found);
+    found = this.collectHits(positions, topology, triB, triA, found);
+    if (found !== 2) {
+      return found;
+    }
+
+    const triangleCount = topology.triangles.length / 3;
+    const first = this.registerHit(this.pairHits[0], triangleCount);
+    const second = this.registerHit(this.pairHits[1], triangleCount);
+    this.segments.push(first, second);
+    return found;
+  }
+
+  /** Tests the edges of `owner` against `target`; keeps the hits in `pairHits` from index `found`. */
+  private collectHits(positions: Float32Array, topology: ClothTopology, owner: number, target: number, found: number): number {
+    const {triangles, triangleEdges, edges} = topology;
+    let total = found;
+    for (let corner = 0; corner < 3; corner++) {
+      const edge = triangleEdges[owner * 3 + corner];
+      const isHit = intersectEdgeTriangle(
+        positions, edges[edge * 2], edges[edge * 2 + 1],
+        triangles[target * 3], triangles[target * 3 + 1], triangles[target * 3 + 2], this.scratch
+      );
+      if (!isHit) {
+        continue;
+      }
+      if (total < this.pairHits.length) {
+        this.pairHits[total].edge = edge;
+        this.pairHits[total].triangle = target;
+        Object.assign(this.pairHits[total].hit, this.scratch);
+      } else {
+        this.pairHits.push({edge, triangle: target, hit: {...this.scratch}});
+      }
+      total++;
+    }
+    return total;
+  }
+
+  private registerHit(pairHit: {edge: number; triangle: number; hit: EdgeTriangleHit}, triangleCount: number): number {
+    const key = pairHit.edge * triangleCount + pairHit.triangle;
+    const existing = this.hitByKey.get(key);
+    if (existing !== undefined) {
+      return existing;
+    }
+    const index = this.count++;
+    this.hitByKey.set(key, index);
+    this.edge.push(pairHit.edge);
+    this.triangle.push(pairHit.triangle);
+    this.t.push(pairHit.hit.t);
+    this.u.push(pairHit.hit.u);
+    this.v.push(pairHit.hit.v);
+    this.points.push(pairHit.hit.x, pairHit.hit.y, pairHit.hit.z);
+    return index;
+  }
+}

src/app/pages/algorithms/cloth-icm/collision/normal-cone.spec.ts

@@ -0,0 +1,49 @@
+import {borderProjectionIsSimple, computeTriangleNormals, coneMerge} from './normal-cone';
+
+/** Vertex layout: [x, y, z, invMass]. */
+function loop(...points: number[][]): {positions: Float32Array; ids: Uint32Array} {
+  const positions = new Float32Array(points.flatMap(point => [point[0], point[1], point[2], 1]));
+  return {positions, ids: Uint32Array.from(points.map((_, index) => index))};
+}
+
+describe('normal cone', () => {
+  it('merges two normals 90 degrees apart into a cone of half 45 degrees', () => {
+    const merged = coneMerge([1, 0, 0], 0, [0, 1, 0], 0);
+    expect(merged.half).toBeCloseTo(Math.PI / 4);
+    expect(merged.axis[0]).toBeCloseTo(Math.SQRT1_2);
+    expect(merged.axis[1]).toBeCloseTo(Math.SQRT1_2);
+  });
+
+  it('keeps the half angle of a part that already covers the other', () => {
+    const merged = coneMerge([0, 0, 1], 1, [0, 0, 1], 0.2);
+    expect(merged.half).toBeCloseTo(1);
+  });
+
+  it('gives half PI for opposite axes', () => {
+    expect(coneMerge([0, 0, 1], 0, [0, 0, -1], 0).half).toBe(Math.PI);
+  });
+
+  it('computes unit triangle normals from the winding', () => {
+    const {positions} = loop([0, 0, 0], [1, 0, 0], [0, 1, 0]);
+    const normals = new Float32Array(3);
+    computeTriangleNormals(positions, Uint32Array.from([0, 1, 2]), normals);
+    expect(Array.from(normals)).toEqual([0, 0, 1]);
+  });
+});
+
+describe('borderProjectionIsSimple', () => {
+  it('accepts a flat loop', () => {
+    const {positions, ids} = loop([0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0]);
+    expect(borderProjectionIsSimple(positions, ids, [0, 0, 1])).toBe(true);
+  });
+
+  it('accepts a raised loop seen from above', () => {
+    const {positions, ids} = loop([0, 0, 0], [1, 0, 0.5], [1, 1, 0], [0, 1, 0.5]);
+    expect(borderProjectionIsSimple(positions, ids, [0, 0, 1])).toBe(true);
+  });
+
+  it('rejects a figure-eight loop', () => {
+    const {positions, ids} = loop([0, 0, 0], [1, 1, 0], [1, 0, 0], [0, 1, 0]);
+    expect(borderProjectionIsSimple(positions, ids, [0, 0, 1])).toBe(false);
+  });
+});

src/app/pages/algorithms/cloth-icm/collision/normal-cone.ts

@@ -0,0 +1,151 @@
+export interface NormalCone {
+  axis: [number, number, number];
+  half: number;
+}
+
+const MIN_AXIS_LENGTH = 1e-6;
+
+let projectedScratch = new Float64Array(0);
+
+/** Unit normal per triangle (vec4 positions, 3 floats per triangle in `out`). Degenerate triangles get (0, 1, 0). */
+export function computeTriangleNormals(positions: Float32Array, triangles: Uint32Array, out: Float32Array): void {
+  const triangleCount = triangles.length / 3;
+  for (let triangle = 0; triangle < triangleCount; triangle++) {
+    const a = triangles[triangle * 3] * 4;
+    const b = triangles[triangle * 3 + 1] * 4;
+    const c = triangles[triangle * 3 + 2] * 4;
+    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;
+    const length = Math.hypot(nx, ny, nz);
+    const base = triangle * 3;
+    if (length < MIN_AXIS_LENGTH * MIN_AXIS_LENGTH) {
+      out[base] = 0;
+      out[base + 1] = 1;
+      out[base + 2] = 0;
+    } else {
+      out[base] = nx / length;
+      out[base + 1] = ny / length;
+      out[base + 2] = nz / length;
+    }
+  }
+}
+
+function angleBetween(ax: number, ay: number, az: number, bx: number, by: number, bz: number): number {
+  const dot = ax * bx + ay * by + az * bz;
+  return Math.acos(Math.min(1, Math.max(-1, dot)));
+}
+
+/**
+ * Merges two cones with unit axes into the cone that contains both.
+ * Writes the new axis to `outAxis[outOffset..+2]` and returns the half angle (PI = every direction).
+ */
+export function coneMergeInto(
+  axisA: ArrayLike<number>, offsetA: number, halfA: number,
+  axisB: ArrayLike<number>, offsetB: number, halfB: number,
+  outAxis: Float32Array, outOffset: number
+): number {
+  const ax = axisA[offsetA];
+  const ay = axisA[offsetA + 1];
+  const az = axisA[offsetA + 2];
+  const bx = axisB[offsetB];
+  const by = axisB[offsetB + 1];
+  const bz = axisB[offsetB + 2];
+  const sx = ax + bx;
+  const sy = ay + by;
+  const sz = az + bz;
+  const length = Math.hypot(sx, sy, sz);
+  if (length < MIN_AXIS_LENGTH) {
+    outAxis[outOffset] = ax;
+    outAxis[outOffset + 1] = ay;
+    outAxis[outOffset + 2] = az;
+    return Math.PI;
+  }
+  const mx = sx / length;
+  const my = sy / length;
+  const mz = sz / length;
+  outAxis[outOffset] = mx;
+  outAxis[outOffset + 1] = my;
+  outAxis[outOffset + 2] = mz;
+  const reachA = angleBetween(mx, my, mz, ax, ay, az) + halfA;
+  const reachB = angleBetween(mx, my, mz, bx, by, bz) + halfB;
+  return Math.min(Math.PI, Math.max(reachA, reachB));
+}
+
+export function coneMerge(axisA: ArrayLike<number>, halfA: number, axisB: ArrayLike<number>, halfB: number): NormalCone {
+  const axis = new Float32Array(3);
+  const half = coneMergeInto(axisA, 0, halfA, axisB, 0, halfB, axis, 0);
+  return {axis: [axis[0], axis[1], axis[2]], half};
+}
+
+/**
+ * Volino 1994: projects the closed border loop onto the plane orthogonal to `axis` and checks that
+ * the projection has no crossing segments. Adjacent segments are skipped. O(B^2).
+ * Only the first `count` ids of `borderVertexIds` are used, so that callers can reuse one buffer.
+ */
+export function borderProjectionIsSimple(
+  positions: Float32Array, borderVertexIds: Uint32Array, axis: ArrayLike<number>, count = borderVertexIds.length
+): boolean {
+  const [ux, uy, uz, vx, vy, vz] = planeBasis(axis[0], axis[1], axis[2]);
+  if (projectedScratch.length < count * 2) {
+    projectedScratch = new Float64Array(count * 2);
+  }
+  const projected = projectedScratch;
+  for (let i = 0; i < count; i++) {
+    const base = borderVertexIds[i] * 4;
+    projected[i * 2] = positions[base] * ux + positions[base + 1] * uy + positions[base + 2] * uz;
+    projected[i * 2 + 1] = positions[base] * vx + positions[base + 1] * vy + positions[base + 2] * vz;
+  }
+
+  for (let i = 0; i < count; i++) {
+    const nextI = (i + 1) % count;
+    for (let j = i + 2; j < count; j++) {
+      const nextJ = (j + 1) % count;
+      if (nextJ === i) {
+        continue;
+      }
+      if (segmentsCross(projected, i, nextI, j, nextJ)) {
+        return false;
+      }
+    }
+  }
+  return true;
+}
+
+function planeBasis(x: number, y: number, z: number): [number, number, number, number, number, number] {
+  // Cross with the world axis that is least aligned with `axis`.
+  const useX = Math.abs(x) <= Math.abs(y) && Math.abs(x) <= Math.abs(z);
+  const useY = !useX && Math.abs(y) <= Math.abs(z);
+  const hx = useX ? 1 : 0;
+  const hy = useY ? 1 : 0;
+  const hz = useX || useY ? 0 : 1;
+  let ux = y * hz - z * hy;
+  let uy = z * hx - x * hz;
+  let uz = x * hy - y * hx;
+  const uLength = Math.hypot(ux, uy, uz);
+  ux /= uLength;
+  uy /= uLength;
+  uz /= uLength;
+  const vx = y * uz - z * uy;
+  const vy = z * ux - x * uz;
+  const vz = x * uy - y * ux;
+  return [ux, uy, uz, vx, vy, vz];
+}
+
+function orientation(p: Float64Array, a: number, b: number, c: number): number {
+  return (p[b * 2] - p[a * 2]) * (p[c * 2 + 1] - p[a * 2 + 1]) - (p[b * 2 + 1] - p[a * 2 + 1]) * (p[c * 2] - p[a * 2]);
+}
+
+function segmentsCross(p: Float64Array, a: number, b: number, c: number, d: number): boolean {
+  const o1 = orientation(p, a, b, c);
+  const o2 = orientation(p, a, b, d);
+  const o3 = orientation(p, c, d, a);
+  const o4 = orientation(p, c, d, b);
+  return o1 * o2 < 0 && o3 * o4 < 0;
+}

src/app/pages/algorithms/cloth-icm/collision/proximity-constraints.spec.ts

@@ -0,0 +1,83 @@
+import {buildClothTopology} from './cloth-topology';
+import {buildProximityConstraints, PROXIMITY_STRIDE, solveProximityConstraints} from './proximity-constraints';
+
+const H = 0.2;
+
+/**
+ * Grid 4x2: triangle 0 = vertices 0, 4, 5 (quad 0), triangle 5 = vertices 2, 7, 3 (quad 2). No shared vertex.
+ * Vertex layout [x, y, z, invMass]; the triangle 0 corners are placed explicitly, the rest of the grid is moved far away.
+ */
+const topology = buildClothTopology(4, 2);
+
+function positionsWith(triangle0: number[][], triangle5: number[][], invMass0 = [1, 1, 1], invMass5 = [1, 1, 1]): Float32Array {
+  const positions = new Float32Array(topology.vertexCount * 4);
+  for (let vertex = 0; vertex < topology.vertexCount; vertex++) {
+    positions.set([100 + vertex, 100, 100, 1], vertex * 4);
+  }
+  [0, 4, 5].forEach((vertex, index) => positions.set([...triangle0[index], invMass0[index]], vertex * 4));
+  [2, 7, 3].forEach((vertex, index) => positions.set([...triangle5[index], invMass5[index]], vertex * 4));
+  return positions;
+}
+
+const BASE_TRIANGLE = [[0, 0, 0], [4, 0, 0], [0, 0, 4]];
+// Far away from the base triangle unless a test places one corner above it.
+const FAR_TRIANGLE = [[50, 50, 50], [60, 50, 50], [50, 60, 50]];
+
+function distanceAbove(positions: Float32Array, vertex: number): number {
+  return positions[vertex * 4 + 1];
+}
+
+describe('proximity constraints', () => {
+  it('pushes a point 0.5 h above a triangle out to h', () => {
+    const positions = positionsWith(BASE_TRIANGLE, [[1, 0.5 * H, 1], ...FAR_TRIANGLE.slice(1)]);
+    const constraints = buildProximityConstraints(positions, topology, [0, 5], H);
+    expect(constraints.length).toBeGreaterThan(0);
+    expect(constraints.length % PROXIMITY_STRIDE).toBe(0);
+
+    solveProximityConstraints(positions, constraints, H);
+
+    const point = positions.slice(2 * 4, 2 * 4 + 3);
+    const triangleHeight = positions[1];
+    expect(point[1] - triangleHeight).toBeGreaterThanOrEqual(H - 1e-4);
+  });
+
+  it('leaves a point at 2 h untouched and builds no constraint for it', () => {
+    const positions = positionsWith(BASE_TRIANGLE, [[1, 2 * H, 1], ...FAR_TRIANGLE.slice(1)]);
+    const before = positions.slice();
+    const constraints = buildProximityConstraints(positions, topology, [0, 5], H);
+    expect(constraints.length).toBe(0);
+    solveProximityConstraints(positions, constraints, H);
+    expect(Array.from(positions)).toEqual(Array.from(before));
+  });
+
+  it('does not push beyond the thickness', () => {
+    const positions = positionsWith(BASE_TRIANGLE, [[1, 0.5 * H, 1], ...FAR_TRIANGLE.slice(1)]);
+    const constraints = buildProximityConstraints(positions, topology, [0, 5], H);
+    solveProximityConstraints(positions, constraints, H);
+    solveProximityConstraints(positions, constraints, H);
+    const gap = distanceAbove(positions, 2) - distanceAbove(positions, 0);
+    // Several constraints with slightly different normals share the point, so a small overshoot is fine.
+    expect(gap).toBeLessThan(H * 1.1);
+  });
+
+  it('keeps pinned vertices in place and moves the free point the whole way', () => {
+    const positions = positionsWith(BASE_TRIANGLE, [[1, 0.5 * H, 1], ...FAR_TRIANGLE.slice(1)], [0, 0, 0]);
+    const constraints = buildProximityConstraints(positions, topology, [0, 5], H);
+    solveProximityConstraints(positions, constraints, H);
+    expect(Array.from(positions.slice(0, 3))).toEqual([0, 0, 0]);
+    expect(distanceAbove(positions, 2)).toBeCloseTo(H, 4);
+  });
+
+  it('pushes two crossing edges apart to h', () => {
+    // Triangle 0 edge 0-4 runs along z at x=0,y=0; the other triangle's edge 2-7 runs along x at height 0.5 h, crossing over it.
+    const positions = positionsWith(
+      [[0, 0, -2], [0, 0, 2], [-2, 0, 0]],
+      [[-1, 0.5 * H, 0], [1, 0.5 * H, 0], [50, 50, 50]]
+    );
+    const constraints = buildProximityConstraints(positions, topology, [0, 5], H);
+    solveProximityConstraints(positions, constraints, H);
+    const edgeHeightA = (distanceAbove(positions, 0) + distanceAbove(positions, 4)) / 2;
+    const edgeHeightB = (distanceAbove(positions, 2) + distanceAbove(positions, 7)) / 2;
+    expect(edgeHeightB - edgeHeightA).toBeGreaterThanOrEqual(H - 1e-3);
+  });
+});

src/app/pages/algorithms/cloth-icm/collision/proximity-constraints.ts

@@ -0,0 +1,275 @@
+import {ClothTopology} from './cloth-topology';
+
+/** [type, i0, i1, i2, i3, w0, w1, w2, w3, nx, ny, nz] */
+export const PROXIMITY_STRIDE = 12;
+const POINT_TRIANGLE = 0;
+const EDGE_EDGE = 1;
+const MIN_DISTANCE = 1e-9;
+
+const closestWeights = new Float64Array(3);
+const closestParameters = new Float64Array(2);
+const solveIds = new Int32Array(4);
+const solveCoefficients = new Float64Array(4);
+
+let constraintBuffer = new Float32Array(PROXIMITY_STRIDE * 1024);
+let constraintLength = 0;
+
+/**
+ * Point-triangle and edge-edge pairs closer than `thickness`, for every candidate triangle pair (2 ids each).
+ * `n` points from the second primitive (triangle / edge 2) to the first (point / edge 1) at build time.
+ * The result is a view on one shared buffer: valid until the next call.
+ */
+export function buildProximityConstraints(
+  positions: Float32Array, topology: ClothTopology, candidatePairs: number[], thickness: number
+): Float32Array {
+  constraintLength = 0;
+  for (let pair = 0; pair < candidatePairs.length; pair += 2) {
+    const triA = candidatePairs[pair];
+    const triB = candidatePairs[pair + 1];
+    addPointTriangleConstraints(positions, topology, triA, triB, thickness);
+    addPointTriangleConstraints(positions, topology, triB, triA, thickness);
+    addEdgeEdgeConstraints(positions, topology, triA, triB, thickness);
+  }
+  return constraintBuffer.subarray(0, constraintLength);
+}
+
+function addPointTriangleConstraints(
+  positions: Float32Array, topology: ClothTopology, pointTriangle: number, triangle: number, thickness: number
+): void {
+  const {triangles} = topology;
+  const t0 = triangles[triangle * 3];
+  const t1 = triangles[triangle * 3 + 1];
+  const t2 = triangles[triangle * 3 + 2];
+  for (let corner = 0; corner < 3; corner++) {
+    const point = triangles[pointTriangle * 3 + corner];
+    closestPointOnTriangle(positions, point, t0, t1, t2, closestWeights);
+    const w0 = closestWeights[0];
+    const w1 = closestWeights[1];
+    const w2 = closestWeights[2];
+    const dx = positions[point * 4] - (w0 * positions[t0 * 4] + w1 * positions[t1 * 4] + w2 * positions[t2 * 4]);
+    const dy = positions[point * 4 + 1] - (w0 * positions[t0 * 4 + 1] + w1 * positions[t1 * 4 + 1] + w2 * positions[t2 * 4 + 1]);
+    const dz = positions[point * 4 + 2] - (w0 * positions[t0 * 4 + 2] + w1 * positions[t1 * 4 + 2] + w2 * positions[t2 * 4 + 2]);
+    pushConstraint(POINT_TRIANGLE, point, t0, t1, t2, 1, w0, w1, w2, dx, dy, dz, thickness);
+  }
+}
+
+function addEdgeEdgeConstraints(
+  positions: Float32Array, topology: ClothTopology, triA: number, triB: number, thickness: number
+): void {
+  const {triangles} = topology;
+  for (let cornerA = 0; cornerA < 3; cornerA++) {
+    const a0 = triangles[triA * 3 + cornerA];
+    const a1 = triangles[triA * 3 + (cornerA + 1) % 3];
+    for (let cornerB = 0; cornerB < 3; cornerB++) {
+      const b0 = triangles[triB * 3 + cornerB];
+      const b1 = triangles[triB * 3 + (cornerB + 1) % 3];
+      closestPointsOnSegments(positions, a0, a1, b0, b1, closestParameters);
+      const s = closestParameters[0];
+      const t = closestParameters[1];
+      const dx = lerp(positions[a0 * 4], positions[a1 * 4], s) - lerp(positions[b0 * 4], positions[b1 * 4], t);
+      const dy = lerp(positions[a0 * 4 + 1], positions[a1 * 4 + 1], s) - lerp(positions[b0 * 4 + 1], positions[b1 * 4 + 1], t);
+      const dz = lerp(positions[a0 * 4 + 2], positions[a1 * 4 + 2], s) - lerp(positions[b0 * 4 + 2], positions[b1 * 4 + 2], t);
+      pushConstraint(EDGE_EDGE, a0, a1, b0, b1, 1 - s, s, 1 - t, t, dx, dy, dz, thickness);
+    }
+  }
+}
+
+function pushConstraint(
+  type: number, id0: number, id1: number, id2: number, id3: number, w0: number, w1: number, w2: number, w3: number,
+  dx: number, dy: number, dz: number, thickness: number
+): void {
+  const distance = Math.hypot(dx, dy, dz);
+  if (distance >= thickness || distance < MIN_DISTANCE) {
+    return;
+  }
+  if (constraintLength + PROXIMITY_STRIDE > constraintBuffer.length) {
+    const grown = new Float32Array(constraintBuffer.length * 2);
+    grown.set(constraintBuffer);
+    constraintBuffer = grown;
+  }
+  const buffer = constraintBuffer;
+  const base = constraintLength;
+  buffer[base] = type;
+  buffer[base + 1] = id0;
+  buffer[base + 2] = id1;
+  buffer[base + 3] = id2;
+  buffer[base + 4] = id3;
+  buffer[base + 5] = w0;
+  buffer[base + 6] = w1;
+  buffer[base + 7] = w2;
+  buffer[base + 8] = w3;
+  buffer[base + 9] = dx / distance;
+  buffer[base + 10] = dy / distance;
+  buffer[base + 11] = dz / distance;
+  constraintLength += PROXIMITY_STRIDE;
+}
+
+function lerp(a: number, b: number, t: number): number {
+  return a + (b - a) * t;
+}
+
+/** XPBD with compliance 0, weighted by inverse mass (positions[i * 4 + 3]). Pushes only while the pair is closer than `thickness`. */
+export function solveProximityConstraints(positions: Float32Array, constraints: Float32Array, thickness: number): void {
+  for (let base = 0; base < constraints.length; base += PROXIMITY_STRIDE) {
+    fillSolveTerms(constraints, base);
+    const nx = constraints[base + 9];
+    const ny = constraints[base + 10];
+    const nz = constraints[base + 11];
+
+    let distance = 0;
+    let weightedInverseMass = 0;
+    for (let slot = 0; slot < 4; slot++) {
+      const vertex = solveIds[slot] * 4;
+      const coefficient = solveCoefficients[slot];
+      distance += coefficient * (positions[vertex] * nx + positions[vertex + 1] * ny + positions[vertex + 2] * nz);
+      weightedInverseMass += positions[vertex + 3] * coefficient * coefficient;
+    }
+    if (distance >= thickness || weightedInverseMass <= 0) {
+      continue;
+    }
+
+    const lambda = (thickness - distance) / weightedInverseMass;
+    for (let slot = 0; slot < 4; slot++) {
+      const vertex = solveIds[slot] * 4;
+      const move = positions[vertex + 3] * solveCoefficients[slot] * lambda;
+      positions[vertex] += nx * move;
+      positions[vertex + 1] += ny * move;
+      positions[vertex + 2] += nz * move;
+    }
+  }
+}
+
+/** Signed coefficients so that the sum of coefficient * position is (first primitive - second primitive). */
+function fillSolveTerms(constraints: Float32Array, base: number): void {
+  const isPointTriangle = constraints[base] === POINT_TRIANGLE;
+  for (let slot = 0; slot < 4; slot++) {
+    solveIds[slot] = constraints[base + 1 + slot];
+    const weight = constraints[base + 5 + slot];
+    const isSecondPrimitive = isPointTriangle ? slot > 0 : slot > 1;
+    solveCoefficients[slot] = isSecondPrimitive ? -weight : weight;
+  }
+}
+
+/** Ericson, Real-Time Collision Detection 5.1.5: barycentric weights of the closest point of the triangle to `point`. */
+function closestPointOnTriangle(positions: Float32Array, point: number, a: number, b: number, c: number, out: Float64Array): void {
+  const pa = a * 4;
+  const pb = b * 4;
+  const pc = c * 4;
+  const pp = point * 4;
+  const abx = positions[pb] - positions[pa];
+  const aby = positions[pb + 1] - positions[pa + 1];
+  const abz = positions[pb + 2] - positions[pa + 2];
+  const acx = positions[pc] - positions[pa];
+  const acy = positions[pc + 1] - positions[pa + 1];
+  const acz = positions[pc + 2] - positions[pa + 2];
+
+  const apx = positions[pp] - positions[pa];
+  const apy = positions[pp + 1] - positions[pa + 1];
+  const apz = positions[pp + 2] - positions[pa + 2];
+  const d1 = abx * apx + aby * apy + abz * apz;
+  const d2 = acx * apx + acy * apy + acz * apz;
+  if (d1 <= 0 && d2 <= 0) {
+    return setWeights(out, 1, 0, 0);
+  }
+
+  const bpx = positions[pp] - positions[pb];
+  const bpy = positions[pp + 1] - positions[pb + 1];
+  const bpz = positions[pp + 2] - positions[pb + 2];
+  const d3 = abx * bpx + aby * bpy + abz * bpz;
+  const d4 = acx * bpx + acy * bpy + acz * bpz;
+  if (d3 >= 0 && d4 <= d3) {
+    return setWeights(out, 0, 1, 0);
+  }
+
+  const vc = d1 * d4 - d3 * d2;
+  if (vc <= 0 && d1 >= 0 && d3 <= 0) {
+    const v = d1 / (d1 - d3);
+    return setWeights(out, 1 - v, v, 0);
+  }
+
+  const cpx = positions[pp] - positions[pc];
+  const cpy = positions[pp + 1] - positions[pc + 1];
+  const cpz = positions[pp + 2] - positions[pc + 2];
+  const d5 = abx * cpx + aby * cpy + abz * cpz;
+  const d6 = acx * cpx + acy * cpy + acz * cpz;
+  if (d6 >= 0 && d5 <= d6) {
+    return setWeights(out, 0, 0, 1);
+  }
+
+  const vb = d5 * d2 - d1 * d6;
+  if (vb <= 0 && d2 >= 0 && d6 <= 0) {
+    const w = d2 / (d2 - d6);
+    return setWeights(out, 1 - w, 0, w);
+  }
+
+  const va = d3 * d6 - d5 * d4;
+  if (va <= 0 && d4 - d3 >= 0 && d5 - d6 >= 0) {
+    const w = (d4 - d3) / (d4 - d3 + (d5 - d6));
+    return setWeights(out, 0, 1 - w, w);
+  }
+
+  const denominator = 1 / (va + vb + vc);
+  const v = vb * denominator;
+  const w = vc * denominator;
+  return setWeights(out, 1 - v - w, v, w);
+}
+
+function setWeights(out: Float64Array, a: number, b: number, c: number): void {
+  out[0] = a;
+  out[1] = b;
+  out[2] = c;
+}
+
+function clamp01(value: number): number {
+  return Math.min(1, Math.max(0, value));
+}
+
+/** Ericson 5.1.9: parameters s (segment a0-a1) and t (segment b0-b1) of the closest points. */
+function closestPointsOnSegments(positions: Float32Array, a0: number, a1: number, b0: number, b1: number, out: Float64Array): void {
+  const pa0 = a0 * 4;
+  const pa1 = a1 * 4;
+  const pb0 = b0 * 4;
+  const pb1 = b1 * 4;
+  const d1x = positions[pa1] - positions[pa0];
+  const d1y = positions[pa1 + 1] - positions[pa0 + 1];
+  const d1z = positions[pa1 + 2] - positions[pa0 + 2];
+  const d2x = positions[pb1] - positions[pb0];
+  const d2y = positions[pb1 + 1] - positions[pb0 + 1];
+  const d2z = positions[pb1 + 2] - positions[pb0 + 2];
+  const rx = positions[pa0] - positions[pb0];
+  const ry = positions[pa0 + 1] - positions[pb0 + 1];
+  const rz = positions[pa0 + 2] - positions[pb0 + 2];
+
+  const a = d1x * d1x + d1y * d1y + d1z * d1z;
+  const e = d2x * d2x + d2y * d2y + d2z * d2z;
+  const f = d2x * rx + d2y * ry + d2z * rz;
+  const epsilon = 1e-12;
+
+  let s = 0;
+  let t = 0;
+  if (a <= epsilon && e <= epsilon) {
+    s = 0;
+    t = 0;
+  } else if (a <= epsilon) {
+    t = clamp01(f / e);
+  } else {
+    const c = d1x * rx + d1y * ry + d1z * rz;
+    if (e <= epsilon) {
+      s = clamp01(-c / a);
+    } else {
+      const b = d1x * d2x + d1y * d2y + d1z * d2z;
+      const denominator = a * e - b * b;
+      s = denominator > epsilon ? clamp01((b * f - c * e) / denominator) : 0;
+      t = (b * s + f) / e;
+      if (t < 0) {
+        t = 0;
+        s = clamp01(-c / a);
+      } else if (t > 1) {
+        t = 1;
+        s = clamp01((b - c) / a);
+      }
+    }
+  }
+  out[0] = s;
+  out[1] = t;
+}

src/app/pages/algorithms/cloth-icm/response/icm-response.spec.ts

@@ -0,0 +1,382 @@
+import {initialPositions} from '../cloth-icm-scenes';
+import {ClothBvh} from '../collision/cloth-bvh';
+import {buildClothTopology, ClothTopology} from '../collision/cloth-topology';
+import {detectCollisions} from '../collision/collision-detection';
+import {Contour, traceContours} from '../collision/intersection-contours';
+import {IntersectionSet} from '../collision/intersection-set';
+import {
+  accumulateIcmDisplacement, applyAveragedDisplacement, contourDisplacement, contourSides, distributeDisplacement, IcmParams, pairGradient
+} from './icm-response';
+
+const PARAMS: IcmParams = {h0: 0.025, g0: 4};
+
+interface Analysis {
+  topology: ClothTopology;
+  positions: Float32Array;
+  set: IntersectionSet;
+  contours: Contour[];
+  normals: Float32Array;
+}
+
+function analyze(topology: ClothTopology, positions: Float32Array): Analysis {
+  const bvh = new ClothBvh(topology);
+  bvh.refit(positions);
+  const set = new IntersectionSet();
+  detectCollisions(bvh, positions, topology, 0, set, []);
+  return {topology, positions, set, contours: traceContours(set), normals: bvh.triangleNormals};
+}
+
+/**
+ * The tangled scene has a ridge (open contours). These specs need one closed contour, so the back half
+ * of the scene gets a round bump instead.
+ */
+function tangled(gridSize: number): Analysis {
+  const topology = buildClothTopology(gridSize, gridSize);
+  const positions = initialPositions('tangled', topology, 0.1, 1);
+  const frontRows = Math.floor(gridSize / 2);
+  const backStart = frontRows * gridSize;
+  let centerY = 0;
+  for (let vertex = backStart; vertex < topology.vertexCount; vertex++) {
+    centerY += positions[vertex * 4 + 1] / (topology.vertexCount - backStart);
+  }
+  const radius = 0.8;
+  for (let vertex = backStart; vertex < topology.vertexCount; vertex++) {
+    const distance = Math.hypot(positions[vertex * 4], positions[vertex * 4 + 1] - centerY);
+    const falloff = Math.max(0, 1 - (distance * distance) / (radius * radius));
+    positions[vertex * 4 + 2] = -0.1 + 0.2 * falloff ** 3;
+  }
+  return analyze(topology, positions);
+}
+
+const FLAP_COLUMNS = 12;
+const FLAP_SHEET_ROWS = 6;
+const FLAP_SHEET_VERTICES = FLAP_COLUMNS * FLAP_SHEET_ROWS;
+
+/**
+ * Two sheets of one grid. Sheet A (rows 0-5) lies in z = 0. Sheet B (rows 6-11) starts behind the lower border of A,
+ * is shifted to the left, and rises through A along a straight line. One open contour: it ends on the left border of A
+ * and on the right border of B. (Two ends on the same sheet would cancel in D: a straight contour has no length gradient
+ * against a relative translation of the sheets.)
+ */
+function flap(): Analysis {
+  const topology = buildClothTopology(FLAP_COLUMNS, FLAP_SHEET_ROWS * 2);
+  const positions = new Float32Array(topology.vertexCount * 4);
+  const angle = 0.9;
+  for (let row = 0; row < FLAP_SHEET_ROWS * 2; row++) {
+    for (let column = 0; column < FLAP_COLUMNS; column++) {
+      const base = (row * FLAP_COLUMNS + column) * 4;
+      if (row < FLAP_SHEET_ROWS) {
+        positions.set([(column - 5.5) * 0.1, -row * 0.1, 0, 1], base);
+      } else {
+        const along = (row - FLAP_SHEET_ROWS) * 0.1;
+        positions.set([(column - 7) * 0.1 + 0.0371, -0.5137 + along * Math.cos(angle), -0.1 + along * Math.sin(angle), 1], base);
+      }
+    }
+  }
+  return analyze(topology, positions);
+}
+
+function isSheetA(vertex: number): boolean {
+  return vertex < FLAP_SHEET_VERTICES;
+}
+
+function sub(a: ArrayLike<number>, b: ArrayLike<number>): number[] {
+  return [a[0] - b[0], a[1] - b[1], a[2] - b[2]];
+}
+
+function dot(a: ArrayLike<number>, b: ArrayLike<number>): number {
+  return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
+}
+
+function displacementOf(analysis: Analysis, contour: Contour, params = PARAMS): Float32Array {
+  const sides = contourSides(analysis.set, contour, analysis.topology);
+  const out = new Float32Array(3);
+  contourDisplacement(analysis.positions, analysis.topology, analysis.normals, analysis.set, contour, sides, params, out);
+  return out;
+}
+
+function accumulate(analysis: Analysis, params = PARAMS): {delta: Float32Array; count: Float32Array} {
+  const delta = new Float32Array(analysis.topology.vertexCount * 3);
+  const count = new Float32Array(analysis.topology.vertexCount);
+  accumulateIcmDisplacement(analysis.positions, analysis.topology, analysis.normals, analysis.set, analysis.contours, params, delta, count);
+  return {delta, count};
+}
+
+function totalSegmentLength(set: IntersectionSet): number {
+  let length = 0;
+  for (let segment = 0; segment < set.segmentCount; segment++) {
+    const a = set.segments[segment * 2];
+    const b = set.segments[segment * 2 + 1];
+    length += Math.hypot(set.points[a * 3] - set.points[b * 3], set.points[a * 3 + 1] - set.points[b * 3 + 1], set.points[a * 3 + 2] - set.points[b * 3 + 2]);
+  }
+  return length;
+}
+
+describe('pairGradient', () => {
+  // Grid 4x2: triangle 0 (vertices 0, 4, 5) lies in z = 0; triangle 4 (vertices 2, 6, 7) crosses it with two hits.
+  const topology = buildClothTopology(4, 2);
+
+  function crossing(swapEdgeEnds: boolean): Analysis {
+    const positions = new Float32Array([
+      0, 0, 0, 1, 2, 0, 0, 1, 0.6, 0.7, -1, 1, 1.0, 1.3, -1, 1,
+      0, 2, 0, 1, 2, 2, 0, 1, 0.6, 0.7, 1, 1, 1.0, 1.3, 1, 1
+    ]);
+    if (swapEdgeEnds) {
+      const first = positions.slice(2 * 4, 2 * 4 + 3);
+      positions.copyWithin(2 * 4, 6 * 4, 6 * 4 + 3);
+      positions.set(first, 6 * 4);
+    }
+    const set = new IntersectionSet();
+    set.addTrianglePair(positions, topology, 0, 4);
+    const normals = new Float32Array(topology.triangles.length);
+    const bvh = new ClothBvh(topology);
+    bvh.refit(positions);
+    normals.set(bvh.triangleNormals);
+    return {topology, positions, set, contours: [], normals};
+  }
+
+  function gradientOf(analysis: Analysis, hit: number, mate: number, normalSign = 1): Float32Array {
+    const normals = analysis.normals.map(value => value * normalSign);
+    const out = new Float32Array(3);
+    pairGradient(analysis.positions, analysis.topology, normals, analysis.set, hit, mate, out);
+    return out;
+  }
+
+  /** Where the contour end moves along the collision line when the edge is displaced by `move`. */
+  function movedHitPoint(analysis: Analysis, hit: number, move: Float32Array, normal: number[]): number[] {
+    const edge = analysis.set.edge[hit];
+    const a = analysis.topology.edges[edge * 2] * 4;
+    const b = analysis.topology.edges[edge * 2 + 1] * 4;
+    const e = sub(analysis.positions.subarray(b, b + 3), analysis.positions.subarray(a, a + 3));
+    const slide = dot(move, normal) / dot(e, normal);
+    return [move[0] - slide * e[0], move[1] - slide * e[1], move[2] - slide * e[2]];
+  }
+
+  it('moves the contour end towards its mate when the edge follows G', () => {
+    const analysis = crossing(false);
+    const g = gradientOf(analysis, 0, 1);
+    const toMate = sub(analysis.set.points.slice(3, 6), analysis.set.points.slice(0, 3));
+    const normal = Array.from(analysis.normals.subarray(0, 3));
+    const progress = movedHitPoint(analysis, 0, g, normal);
+    expect(dot(progress, toMate)).toBeGreaterThan(0);
+  });
+
+  it('gives the same G for a reversed triangle normal and a reversed edge', () => {
+    const base = crossing(false);
+    const flipped = crossing(true);
+    const reference = gradientOf(base, 0, 1);
+    const reversedNormal = gradientOf(base, 0, 1, -1);
+    const reversedEdge = gradientOf(flipped, 0, 1);
+    for (let axis = 0; axis < 3; axis++) {
+      expect(reversedNormal[axis]).toBeCloseTo(reference[axis], 5);
+      expect(reversedEdge[axis]).toBeCloseTo(reference[axis], 5);
+    }
+  });
+
+  it('gives opposite R parts for the two ends of one segment and stays finite', () => {
+    const analysis = crossing(false);
+    const g01 = gradientOf(analysis, 0, 1);
+    const g10 = gradientOf(analysis, 1, 0);
+    expect(g01.every(Number.isFinite)).toBe(true);
+    expect(g10.every(Number.isFinite)).toBe(true);
+    expect(Math.hypot(g01[0], g01[1], g01[2])).toBeGreaterThan(0);
+  });
+});
+
+describe('contourSides', () => {
+  it('splits a closed loop into the two sheets that cross', () => {
+    const analysis = tangled(30);
+    expect(analysis.contours.length).toBeGreaterThan(0);
+    const frontRows = 15;
+    for (const contour of analysis.contours) {
+      const sides = contourSides(analysis.set, contour, analysis.topology);
+      expect(sides.includes(0)).toBe(true);
+      expect(sides.includes(1)).toBe(true);
+      const sheetOfSide = new Map<number, boolean>();
+      contour.points.forEach((hit, index) => {
+        const firstVertex = analysis.topology.edges[analysis.set.edge[hit] * 2];
+        const isFront = Math.floor(firstVertex / 30) < frontRows;
+        const known = sheetOfSide.get(sides[index]);
+        expect(known === undefined || known === isFront).toBe(true);
+        sheetOfSide.set(sides[index], isFront);
+      });
+      expect(sheetOfSide.get(0)).not.toBe(sheetOfSide.get(1));
+    }
+  });
+
+  it('puts the edges of each sheet on their own side along an open contour', () => {
+    const analysis = flap();
+    expect(analysis.contours.length).toBe(1);
+    const [contour] = analysis.contours;
+    expect(contour.closed).toBe(false);
+    const sides = contourSides(analysis.set, contour, analysis.topology);
+    const sideOfSheetA = new Set<number>();
+    const sideOfSheetB = new Set<number>();
+    contour.points.forEach((hit, index) => {
+      const vertex = analysis.topology.edges[analysis.set.edge[hit] * 2];
+      (isSheetA(vertex) ? sideOfSheetA : sideOfSheetB).add(sides[index]);
+    });
+    expect(sideOfSheetA.size).toBe(1);
+    expect(sideOfSheetB.size).toBe(1);
+    expect([...sideOfSheetA][0]).not.toBe([...sideOfSheetB][0]);
+  });
+});
+
+describe('contourDisplacement', () => {
+  it('stays below h0 and reaches h0 when g0 is small against |D|', () => {
+    const analysis = flap();
+    const [contour] = analysis.contours;
+    const smoothed = displacementOf(analysis, contour);
+    const sharp = displacementOf(analysis, contour, {h0: PARAMS.h0, g0: 1e-6});
+    expect(Math.hypot(smoothed[0], smoothed[1], smoothed[2])).toBeLessThan(PARAMS.h0);
+    expect(Math.hypot(sharp[0], sharp[1], sharp[2])).toBeCloseTo(PARAMS.h0, 6);
+  });
+
+  it('gives a longer H for a larger |D| (shorter g0 effect)', () => {
+    const analysis = flap();
+    const [contour] = analysis.contours;
+    const strong = displacementOf(analysis, contour, {h0: PARAMS.h0, g0: 0.5});
+    const weak = displacementOf(analysis, contour, {h0: PARAMS.h0, g0: 50});
+    expect(Math.hypot(strong[0], strong[1], strong[2])).toBeGreaterThan(Math.hypot(weak[0], weak[1], weak[2]));
+  });
+});
+
+describe('accumulateIcmDisplacement', () => {
+  it('shrinks the closed contour of the tangled scene', () => {
+    const analysis = tangled(30);
+    expect(analysis.contours.some(contour => contour.closed)).toBe(true);
+    const before = totalSegmentLength(analysis.set);
+    const {delta, count} = accumulate(analysis);
+    applyAveragedDisplacement(analysis.positions, delta, count, 1);
+
+    const after = totalSegmentLength(analyze(analysis.topology, analysis.positions).set);
+    expect(after).toBeLessThan(before);
+  });
+
+  it('moves every touched particle along one common direction (the contour has one H)', () => {
+    const analysis = tangled(30);
+    const closed = analysis.contours.filter(contour => contour.closed);
+    expect(analysis.contours.length).toBe(closed.length);
+    const single = {...analysis, contours: [closed[0]]};
+    const {delta, count} = accumulate(single);
+    const h = displacementOf(single, closed[0]);
+
+    let touched = 0;
+    for (let vertex = 0; vertex < count.length; vertex++) {
+      if (count[vertex] === 0) {
+        continue;
+      }
+      touched++;
+      const move = [delta[vertex * 3], delta[vertex * 3 + 1], delta[vertex * 3 + 2]];
+      const cross = [
+        move[1] * h[2] - move[2] * h[1],
+        move[2] * h[0] - move[0] * h[2],
+        move[0] * h[1] - move[1] * h[0]
+      ];
+      expect(Math.hypot(cross[0], cross[1], cross[2])).toBeLessThan(1e-6);
+    }
+    expect(touched).toBeGreaterThan(0);
+  });
+
+  it('moves the two sheets of a straight open contour in opposite directions, everywhere along it', () => {
+    const analysis = flap();
+    const [contour] = analysis.contours;
+    const h = displacementOf(analysis, contour);
+    const {delta, count} = accumulate(analysis);
+
+    const alongH = (vertex: number): number => dot([delta[vertex * 3], delta[vertex * 3 + 1], delta[vertex * 3 + 2]], h);
+    const touchedA: number[] = [];
+    const touchedB: number[] = [];
+    for (let vertex = 0; vertex < count.length; vertex++) {
+      if (count[vertex] > 0) {
+        (isSheetA(vertex) ? touchedA : touchedB).push(vertex);
+      }
+    }
+    expect(touchedA.length).toBeGreaterThan(0);
+    expect(touchedB.length).toBeGreaterThan(0);
+    const signA = Math.sign(alongH(touchedA[0]));
+    expect(signA).not.toBe(0);
+    for (const vertex of touchedA) {
+      expect(Math.sign(alongH(vertex))).toBe(signA);
+    }
+    for (const vertex of touchedB) {
+      expect(Math.sign(alongH(vertex))).toBe(-signA);
+    }
+  });
+
+  it('gives edge and triangle of every hit opposite displacements that add up to H', () => {
+    const analysis = flap();
+    const [contour] = analysis.contours;
+    const h = displacementOf(analysis, contour);
+    const hit = contour.points[Math.floor(contour.points.length / 2)];
+    const delta = new Float32Array(analysis.topology.vertexCount * 3);
+    const count = new Float32Array(analysis.topology.vertexCount);
+    distributeDisplacement(analysis.positions, analysis.topology, analysis.set, hit, h, 1, delta, count);
+
+    const sum = (vertices: ArrayLike<number>): number[] => {
+      const total = [0, 0, 0];
+      for (const vertex of Array.from(vertices)) {
+        for (let axis = 0; axis < 3; axis++) {
+          total[axis] += delta[vertex * 3 + axis];
+        }
+      }
+      return total;
+    };
+    const edge = analysis.set.edge[hit];
+    const edgeMove = sum(analysis.topology.edges.subarray(edge * 2, edge * 2 + 2));
+    const triangle = analysis.set.triangle[hit];
+    const triangleMove = sum(analysis.topology.triangles.subarray(triangle * 3, triangle * 3 + 3));
+    for (let axis = 0; axis < 3; axis++) {
+      expect(edgeMove[axis] - triangleMove[axis]).toBeCloseTo(h[axis], 6);
+      expect(edgeMove[axis] * triangleMove[axis]).toBeLessThanOrEqual(1e-12);
+    }
+  });
+
+  it('never moves a pinned particle but still moves its partner', () => {
+    const analysis = flap();
+    const first = accumulate(analysis);
+    const pinned = first.count.findIndex(value => value > 0);
+    analysis.positions[pinned * 4 + 3] = 0;
+
+    const second = accumulate(analysis);
+    expect(second.count[pinned]).toBe(0);
+    expect(Math.hypot(second.delta[pinned * 3], second.delta[pinned * 3 + 1], second.delta[pinned * 3 + 2])).toBe(0);
+    const others = Array.from(second.count).filter(value => value > 0).length;
+    expect(others).toBeGreaterThan(0);
+  });
+
+  it('gives the whole H to the edge when the triangle is fully pinned', () => {
+    const analysis = flap();
+    const [contour] = analysis.contours;
+    const h = displacementOf(analysis, contour);
+    const hit = contour.points[1];
+    const triangle = analysis.set.triangle[hit];
+    const triangleVertices = Array.from(analysis.topology.triangles.subarray(triangle * 3, triangle * 3 + 3));
+    for (const vertex of triangleVertices) {
+      analysis.positions[vertex * 4 + 3] = 0;
+    }
+    const delta = new Float32Array(analysis.topology.vertexCount * 3);
+    const count = new Float32Array(analysis.topology.vertexCount);
+    distributeDisplacement(analysis.positions, analysis.topology, analysis.set, hit, h, 1, delta, count);
+
+    const edge = analysis.set.edge[hit];
+    for (let axis = 0; axis < 3; axis++) {
+      const edgeSum = delta[analysis.topology.edges[edge * 2] * 3 + axis] + delta[analysis.topology.edges[edge * 2 + 1] * 3 + axis];
+      expect(edgeSum).toBeCloseTo(h[axis], 6);
+    }
+    for (const vertex of triangleVertices) {
+      expect(count[vertex]).toBe(0);
+    }
+  });
+});
+
+describe('applyAveragedDisplacement', () => {
+  it('averages a particle that is in several pairs and skips particles without a count', () => {
+    const positions = new Float32Array([0, 0, 0, 1, 5, 5, 5, 1]);
+    const delta = new Float32Array([0.3, 0, 0, 9, 9, 9]);
+    const count = new Float32Array([3, 0]);
+    applyAveragedDisplacement(positions, delta, count, 0.5);
+    expect(positions[0]).toBeCloseTo(0.05);
+    expect(Array.from(positions.slice(4, 7))).toEqual([5, 5, 5]);
+  });
+});

src/app/pages/algorithms/cloth-icm/response/icm-response.ts

@@ -0,0 +1,268 @@
+import {Contour} from '../collision/intersection-contours';
+import {IntersectionSet} from '../collision/intersection-set';
+import {ClothTopology} from '../collision/cloth-topology';
+
+/**
+ * Intersection Contour Minimization (Volino and Magnenat-Thalmann 2006), global scheme.
+ * Formulas follow the diploma thesis of Dahm (2010), file Dahm2010-Diplomarbeit.pdf of the diploma project.
+ * Page numbers are the printed ones (PDF page = printed + 6).
+ *
+ * Per contour: a gradient G for every segment end, one common displacement H for the whole contour,
+ * and a spread of H over the particles of every pierced triangle and piercing edge.
+ */
+export interface IcmParams {
+  /** Maximum displacement of one contour per frame (h0). */
+  h0: number;
+  /** Smooths H for short contours (g0). */
+  g0: number;
+}
+
+const MIN_LENGTH = 1e-12;
+
+const gradientScratch = new Float32Array(3);
+const distanceSquaredScratch = new Float64Array(3);
+const distanceShareScratch = new Float64Array(3);
+const edgeVertexScratch = new Int32Array(2);
+const triangleVertexScratch = new Int32Array(3);
+
+/** Side 0 or 1 per contour point: which of the two polygon chains of a self-contour the hit's edge belongs to. */
+export function contourSides(set: IntersectionSet, contour: Contour, topology: ClothTopology): Int8Array {
+  // Diploma 5.1 p. 42, 4.1.4 p. 30-31. The two triangles of a pair never share a vertex, so two hits of one
+  // segment are on the same side exactly when their edges belong to the same triangle (share a vertex).
+  const sides = new Int8Array(contour.points.length);
+  for (let index = 1; index < contour.points.length; index++) {
+    const previousEdge = set.edge[contour.points[index - 1]];
+    const edge = set.edge[contour.points[index]];
+    const sameSide = edgesShareVertex(topology, previousEdge, edge);
+    sides[index] = sameSide ? sides[index - 1] : 1 - sides[index - 1];
+  }
+  return sides;
+}
+
+function edgesShareVertex(topology: ClothTopology, edgeA: number, edgeB: number): boolean {
+  const {edges} = topology;
+  const a0 = edges[edgeA * 2];
+  const a1 = edges[edgeA * 2 + 1];
+  const b0 = edges[edgeB * 2];
+  const b1 = edges[edgeB * 2 + 1];
+  return a0 === b0 || a0 === b1 || a1 === b0 || a1 === b1;
+}
+
+/**
+ * Gradient G of one segment end (diploma 3.5.3 p. 19, Volino's formula): the direction in which the edge E of the hit
+ * moves relative to the pierced triangle with unit normal N so that the contour gets shorter.
+ *   G = R - 2 * (E . R) / (E . N) * N
+ * R is the unit direction of the collision line from the hit point towards the other end of the segment (`mate`).
+ * The mate lies in the triangle that owns E, so R points from the edge into that triangle (the orientation check of the diploma).
+ * E and N may point either way: both flip together and G stays the same.
+ */
+export function pairGradient(
+  positions: Float32Array, topology: ClothTopology, triangleNormals: Float32Array,
+  set: IntersectionSet, hit: number, mate: number, out: Float32Array
+): void {
+  const edge = set.edge[hit];
+  const vertexA = topology.edges[edge * 2] * 4;
+  const vertexB = topology.edges[edge * 2 + 1] * 4;
+  const ex = positions[vertexB] - positions[vertexA];
+  const ey = positions[vertexB + 1] - positions[vertexA + 1];
+  const ez = positions[vertexB + 2] - positions[vertexA + 2];
+
+  let rx = set.points[mate * 3] - set.points[hit * 3];
+  let ry = set.points[mate * 3 + 1] - set.points[hit * 3 + 1];
+  let rz = set.points[mate * 3 + 2] - set.points[hit * 3 + 2];
+  const rLength = Math.hypot(rx, ry, rz);
+  if (rLength < MIN_LENGTH) {
+    out[0] = 0;
+    out[1] = 0;
+    out[2] = 0;
+    return;
+  }
+  rx /= rLength;
+  ry /= rLength;
+  rz /= rLength;
+
+  const normalBase = set.triangle[hit] * 3;
+  const nx = triangleNormals[normalBase];
+  const ny = triangleNormals[normalBase + 1];
+  const nz = triangleNormals[normalBase + 2];
+  const edgeDotNormal = ex * nx + ey * ny + ez * nz;
+  const edgeLength = Math.hypot(ex, ey, ez);
+  if (Math.abs(edgeDotNormal) < MIN_LENGTH * edgeLength) {
+    out[0] = rx;
+    out[1] = ry;
+    out[2] = rz;
+    return;
+  }
+
+  const factor = (2 * (ex * rx + ey * ry + ez * rz)) / edgeDotNormal;
+  out[0] = rx - factor * nx;
+  out[1] = ry - factor * ny;
+  out[2] = rz - factor * nz;
+}
+
+/**
+ * Displacement H of one contour (diploma 3.5.5 p. 21, 5.1 p. 42, 3.5.4 p. 21, 6.4.2 p. 61-62).
+ *   D = sum of s_i * G_i over all segment ends of the contour, s_i = +1 on side 0 and -1 on side 1
+ *   H = h0 * D / sqrt(|D|^2 + g0^2)
+ * |H| stays below h0 and grows towards h0 with |D|. All particles of one contour move along the same direction.
+ */
+export function contourDisplacement(
+  positions: Float32Array, topology: ClothTopology, triangleNormals: Float32Array,
+  set: IntersectionSet, contour: Contour, sides: Int8Array, params: IcmParams, out: Float32Array
+): void {
+  let dx = 0;
+  let dy = 0;
+  let dz = 0;
+  const {points} = contour;
+  const segmentCount = contour.closed ? points.length : points.length - 1;
+  for (let index = 0; index < segmentCount; index++) {
+    const nextIndex = (index + 1) % points.length;
+    const hitA = points[index];
+    const hitB = points[nextIndex];
+    const signA = sides[index] === 0 ? 1 : -1;
+    const signB = sides[nextIndex] === 0 ? 1 : -1;
+
+    pairGradient(positions, topology, triangleNormals, set, hitA, hitB, gradientScratch);
+    dx += signA * gradientScratch[0];
+    dy += signA * gradientScratch[1];
+    dz += signA * gradientScratch[2];
+
+    pairGradient(positions, topology, triangleNormals, set, hitB, hitA, gradientScratch);
+    dx += signB * gradientScratch[0];
+    dy += signB * gradientScratch[1];
+    dz += signB * gradientScratch[2];
+  }
+
+  const scale = params.h0 / Math.sqrt(dx * dx + dy * dy + dz * dz + params.g0 * params.g0);
+  out[0] = dx * scale;
+  out[1] = dy * scale;
+  out[2] = dz * scale;
+}
+
+/**
+ * Spreads the impulse I = sign * H of one hit over the particles of its edge and its triangle (diploma 3.6 p. 25, 4.1.5 p. 31-33).
+ * The edge gets +I and the triangle -I, so the two primitives move apart by H relative to each other.
+ *   1. Split between edge and triangle by the sum of inverse masses (p. 32): with equal masses the triangle gets 3/5 and the edge 2/5. A fully pinned primitive gets nothing.
+ *   2. Inside a primitive, particle x gets f_x = (sM_x + sD_x) / 2, the mean of its mass share and its distance share
+ *      (nearer to the hit point = more). Pinned particles get 0, the distance shares are renormalised over the free ones.
+ * `outDelta` (3 per vertex) and `outCount` collect the result; see applyAveragedDisplacement.
+ */
+export function distributeDisplacement(
+  positions: Float32Array, topology: ClothTopology, set: IntersectionSet, hit: number,
+  h: Float32Array, sign: number, outDelta: Float32Array, outCount: Float32Array
+): void {
+  const edge = set.edge[hit];
+  const triangle = set.triangle[hit];
+  const edgeVertices = edgeVertexScratch;
+  edgeVertices[0] = topology.edges[edge * 2];
+  edgeVertices[1] = topology.edges[edge * 2 + 1];
+  const triangleVertices = triangleVertexScratch;
+  triangleVertices.set(topology.triangles.subarray(triangle * 3, triangle * 3 + 3));
+  const edgeInverseMass = inverseMassOf(positions, edgeVertices);
+  const triangleInverseMass = inverseMassOf(positions, triangleVertices);
+  const totalInverseMass = edgeInverseMass + triangleInverseMass;
+  if (totalInverseMass <= 0) {
+    return;
+  }
+
+  const impulseX = sign * h[0];
+  const impulseY = sign * h[1];
+  const impulseZ = sign * h[2];
+  const edgeShare = edgeInverseMass / totalInverseMass;
+  const triangleShare = triangleInverseMass / totalInverseMass;
+  const hitBase = hit * 3;
+
+  addToPrimitive(positions, set.points, hitBase, edgeVertices, edgeShare, impulseX, impulseY, impulseZ, outDelta, outCount);
+  addToPrimitive(positions, set.points, hitBase, triangleVertices, -triangleShare, impulseX, impulseY, impulseZ, outDelta, outCount);
+}
+
+function inverseMassOf(positions: Float32Array, vertices: Int32Array): number {
+  let sum = 0;
+  for (const vertex of vertices) {
+    sum += positions[vertex * 4 + 3];
+  }
+  return sum;
+}
+
+function addToPrimitive(
+  positions: Float32Array, points: number[], hitBase: number, vertices: Int32Array, primitiveShare: number,
+  impulseX: number, impulseY: number, impulseZ: number, outDelta: Float32Array, outCount: Float32Array
+): void {
+  if (primitiveShare === 0) {
+    return;
+  }
+  const hitX = points[hitBase];
+  const hitY = points[hitBase + 1];
+  const hitZ = points[hitBase + 2];
+  const count = vertices.length;
+
+  let inverseMassSum = 0;
+  const distancesSquared = distanceSquaredScratch;
+  let distanceSum = 0;
+  for (let index = 0; index < count; index++) {
+    const base = vertices[index] * 4;
+    const dx = positions[base] - hitX;
+    const dy = positions[base + 1] - hitY;
+    const dz = positions[base + 2] - hitZ;
+    const distanceSquared = dx * dx + dy * dy + dz * dz;
+    distancesSquared[index] = distanceSquared;
+    distanceSum += distanceSquared;
+    inverseMassSum += positions[base + 3];
+  }
+
+  // Distance share: edge = the other particle's distance, triangle = half of the other two (p. 32-33).
+  const distanceShares = distanceShareScratch;
+  let freeDistanceShareSum = 0;
+  for (let index = 0; index < count; index++) {
+    const share = distanceSum > 0 ? (distanceSum - distancesSquared[index]) / ((count - 1) * distanceSum) : 1 / count;
+    distanceShares[index] = positions[vertices[index] * 4 + 3] > 0 ? share : 0;
+    freeDistanceShareSum += distanceShares[index];
+  }
+
+  for (let index = 0; index < count; index++) {
+    const vertex = vertices[index];
+    const inverseMass = positions[vertex * 4 + 3];
+    if (inverseMass <= 0) {
+      continue;
+    }
+    const massShare = inverseMass / inverseMassSum;
+    const distanceShare = freeDistanceShareSum > 0 ? distanceShares[index] / freeDistanceShareSum : massShare;
+    const factor = primitiveShare * (massShare + distanceShare) / 2;
+    outDelta[vertex * 3] += factor * impulseX;
+    outDelta[vertex * 3 + 1] += factor * impulseY;
+    outDelta[vertex * 3 + 2] += factor * impulseZ;
+    outCount[vertex] += 1;
+  }
+}
+
+/** Runs sides, displacement and distribution for every contour. `outDelta` and `outCount` must be zeroed by the caller. */
+export function accumulateIcmDisplacement(
+  positions: Float32Array, topology: ClothTopology, triangleNormals: Float32Array,
+  set: IntersectionSet, contours: Contour[], params: IcmParams, outDelta: Float32Array, outCount: Float32Array
+): void {
+  const displacement = new Float32Array(3);
+  for (const contour of contours) {
+    const sides = contourSides(set, contour, topology);
+    contourDisplacement(positions, topology, triangleNormals, set, contour, sides, params, displacement);
+    for (let index = 0; index < contour.points.length; index++) {
+      const sign = sides[index] === 0 ? 1 : -1;
+      distributeDisplacement(positions, topology, set, contour.points[index], displacement, sign, outDelta, outCount);
+    }
+  }
+}
+
+/**
+ * A particle in many pairs is averaged, not summed: positions += scale * delta / count (diploma p. 29, the result is stored
+ * per particle and applied once). Pinned particles never have a count.
+ */
+export function applyAveragedDisplacement(positions: Float32Array, delta: Float32Array, count: Float32Array, scale: number): void {
+  for (let vertex = 0; vertex < count.length; vertex++) {
+    if (count[vertex] <= 0) {
+      continue;
+    }
+    const factor = scale / count[vertex];
+    positions[vertex * 4] += delta[vertex * 3] * factor;
+    positions[vertex * 4 + 1] += delta[vertex * 3 + 1] * factor;
+    positions[vertex * 4 + 2] += delta[vertex * 3 + 2] * factor;
+  }
+}

src/app/pages/algorithms/fluid/fluid-simulation.spec.ts

@@ -1,4 +1,4 @@
-import {FluidParams, FluidSimulation} from './fluid-simulation';
+import type {FluidParams} from './fluid-simulation';
 import {DEFAULT_SETTINGS, MAX_PRESSURE_ITERATIONS, Splat} from './fluid.model';
 
 const renderCalls: string[] = [];
@@ -62,7 +62,11 @@ function params(pressureIterations: number, splats: Splat[] = []): FluidParams {
   return {settings: {...DEFAULT_SETTINGS, pressureIterations}, splats};
 }
 
-function start(readParams: () => FluidParams) {
+// gpu-state.ts is shared between specs and keeps the @babylonjs/core of whichever spec loaded it first (the real
+// one, or another spec's mock). A fresh module graph per test makes it use this spec's mock.
+async function start(readParams: () => FluidParams) {
+  vi.resetModules();
+  const {FluidSimulation} = await import('./fluid-simulation');
   renderCalls.length = 0;
   disposedTargets = 0;
   disposedWrappers = 0;
@@ -78,9 +82,9 @@ function count(name: string): number {
 }
 
 describe('FluidSimulation', () => {
-  it('drawsNothingUntilAllShadersReady', () => {
+  it('drawsNothingUntilAllShadersReady', async () => {
     const readParams = vi.fn(() => params(1, [splat]));
-    const {scene} = start(readParams);
+    const {scene} = await start(readParams);
     ready['fluidAdvect'] = false;
 
     scene.frame();
@@ -93,9 +97,9 @@ describe('FluidSimulation', () => {
     expect(renderCalls).toContain('fluidSplat');
   });
 
-  it('readsParamsAtFrameStart', () => {
+  it('readsParamsAtFrameStart', async () => {
     let iterations = 5;
-    const {scene} = start(() => params(iterations));
+    const {scene} = await start(() => params(iterations));
     scene.frame();
     expect(count('fluidPressure')).toBe(5);
 
@@ -105,14 +109,14 @@ describe('FluidSimulation', () => {
     expect(count('fluidPressure')).toBe(12);
   });
 
-  it('clampsPressureIterations', () => {
+  it('clampsPressureIterations', async () => {
-    const {scene} = start(() => params(1000));
+    const {scene} = await start(() => params(1000));
     scene.frame();
     expect(count('fluidPressure')).toBe(MAX_PRESSURE_ITERATIONS);
   });
 
-  it('passOrder', () => {
+  it('passOrder', async () => {
-    const {scene} = start(() => params(2, [splat]));
+    const {scene} = await start(() => params(2, [splat]));
     scene.frame();
     renderCalls.length = 0;
     scene.frame();
@@ -123,8 +127,8 @@ describe('FluidSimulation', () => {
     ]);
   });
 
-  it('resetClearsBeforeSplats', () => {
+  it('resetClearsBeforeSplats', async () => {
-    const {scene, simulation} = start(() => params(1, [splat]));
+    const {scene, simulation} = await start(() => params(1, [splat]));
     scene.frame();
     renderCalls.length = 0;
     scene.frame();
@@ -136,8 +140,8 @@ describe('FluidSimulation', () => {
     expect(renderCalls.slice(0, 4)).toEqual(['fluidClear', 'fluidClear', 'fluidClear', 'fluidSplat']);
   });
 
-  it('disposeFreesAllTargets', () => {
+  it('disposeFreesAllTargets', async () => {
-    const {simulation} = start(() => params(1));
+    const {simulation} = await start(() => params(1));
     simulation.dispose();
     expect(disposedTargets).toBe(8);
     expect(disposedWrappers).toBe(8);

src/app/pages/algorithms/gray-scott/gray-scott-simulation.spec.ts

@@ -1,4 +1,4 @@
-import {GrayScottParams, GrayScottSimulation} from './gray-scott-simulation';
+import type {GrayScottParams} from './gray-scott-simulation';
 import {DEFAULT_SETTINGS} from './gray-scott.model';
 
 const renderCalls: string[] = [];
@@ -55,9 +55,13 @@ function paramsWithSteps(stepsPerFrame: number): GrayScottParams {
 }
 
 describe('GrayScottSimulation', () => {
-  it('readsParamsAtFrameStart', () => {
+  it('readsParamsAtFrameStart', async () => {
     let stepsPerFrame = 1;
     const scene = fakeScene();
+    // gpu-state.ts is shared between specs and keeps the @babylonjs/core of whichever spec loaded it first; a fresh
+    // module graph makes it use this spec's mock.
+    vi.resetModules();
+    const {GrayScottSimulation} = await import('./gray-scott-simulation');
     const simulation = new GrayScottSimulation();
     simulation.init(scene as never, {} as never, 'float', () => paramsWithSteps(stepsPerFrame));
     scene.frame();

src/app/pages/algorithms/information/source-paths.ts

@@ -10,6 +10,7 @@ export const SOURCE_PATHS = {
   FRACTAL3D: ['src/app/pages/algorithms/fractal3d/fractal3d.shader.ts'],
   PENDULUM: ['src/app/pages/algorithms/pendulum/pendulum.shader.ts', 'src/app/pages/algorithms/pendulum/strategies/pendulum-gpu.strategy.ts', 'src/app/pages/algorithms/pendulum/strategies/pendulum-cpu-physics.ts'],
   CLOTH: ['src/app/pages/algorithms/cloth/cloth.shader.ts', 'src/app/pages/algorithms/cloth/strategies/cloth-gpu.strategy.ts', 'src/app/pages/algorithms/cloth/strategies/cloth-cpu-physics.ts'],
+  CLOTH_ICM: ['src/app/pages/algorithms/cloth-icm/response/icm-response.ts', 'src/app/pages/algorithms/cloth-icm/collision/cloth-bvh.ts', 'src/app/pages/algorithms/cloth-icm/collision/intersection-contours.ts'],
   BOIDS: ['src/app/pages/algorithms/boids/boids.shader.ts', 'src/app/pages/algorithms/boids/strategies/boids-gpu.strategy.ts', 'src/app/pages/algorithms/boids/strategies/boids-cpu-physics.ts'],
   PHYSARUM: ['src/app/pages/algorithms/physarum/physarum.shader.ts', 'src/app/pages/algorithms/physarum/strategies/physarum-gpu.strategy.ts', 'src/app/pages/algorithms/physarum/strategies/physarum-cpu-physics.ts'],
   GRAY_SCOTT: ['src/app/pages/algorithms/gray-scott/gray-scott.shader.ts', 'src/app/pages/algorithms/gray-scott/gray-scott-simulation.ts'],

src/app/pages/algorithms/lenia/lenia-simulation.spec.ts

@@ -1,4 +1,4 @@
-import {LeniaParams, LeniaSimulation} from './lenia-simulation';
+import type {LeniaParams} from './lenia-simulation';
 import {DEFAULT_SETTINGS} from './lenia.model';
 
 const renderCalls: {name: string; target: unknown}[] = [];
@@ -64,11 +64,15 @@ function paramsWithSteps(stepsPerFrame: number): LeniaParams {
 }
 
 describe('LeniaSimulation', () => {
-  it('keepsSeedWhileStepShaderCompiles', () => {
+  it('keepsSeedWhileStepShaderCompiles', async () => {
     renderCalls.length = 0;
     ready['leniaCopy'] = true;
     ready['leniaStep'] = false;
     const scene = fakeScene();
+    // gpu-state.ts is shared between specs and keeps the @babylonjs/core of whichever spec loaded it first; a fresh
+    // module graph makes it use this spec's mock.
+    vi.resetModules();
+    const {LeniaSimulation} = await import('./lenia-simulation');
     const simulation = new LeniaSimulation();
     simulation.init(scene as never, {} as never, 'float', () => paramsWithSteps(1));
     simulation.seed(new Uint8Array(4));
@@ -85,12 +89,16 @@ describe('LeniaSimulation', () => {
     expect(firstStep.target).not.toBe(seededTarget);
   });
 
-  it('readsParamsAtFrameStart', () => {
+  it('readsParamsAtFrameStart', async () => {
     renderCalls.length = 0;
     ready['leniaCopy'] = true;
     ready['leniaStep'] = true;
     let stepsPerFrame = 1;
     const scene = fakeScene();
+    // gpu-state.ts is shared between specs and keeps the @babylonjs/core of whichever spec loaded it first; a fresh
+    // module graph makes it use this spec's mock.
+    vi.resetModules();
+    const {LeniaSimulation} = await import('./lenia-simulation');
     const simulation = new LeniaSimulation();
     simulation.init(scene as never, {} as never, 'float', () => paramsWithSteps(stepsPerFrame));
     simulation.seed(new Uint8Array(4));

src/assets/i18n/de.json

@@ -712,6 +712,43 @@
       "ITEM_5": "Wind: feste Richtung (5, 0, 15) plus ein Sinus-Flattern aus Position und Zeit. Der Flatter-Anteil gibt zusätzlich leichten Auftrieb (y: Flattern · 2)."
     }
   },
+  "CLOTH_ICM": {
+    "TITLE": "Stoff-Selbstkollision (ICM)",
+    "SCENE": "Szene",
+    "SCENE_TANGLED": "Verheddigter Start",
+    "SCENE_DROP": "Fallender Stoff",
+    "ICM_ON": "ICM einschalten",
+    "ICM_OFF": "ICM ausschalten",
+    "BVH_ON": "BVH anzeigen",
+    "BVH_OFF": "BVH ausblenden",
+    "CULLING_ON": "Culling an",
+    "CULLING_OFF": "Culling aus",
+    "RESTART_SIMULATION": "Simulation neu starten",
+    "ICM_STEP": "ICM-Schritt (x Abstand)",
+    "BVH_DEPTH": "BVH-Ebene",
+    "HINT": "Rote Linien sind die Schnittkonturen. BVH-Boxen: grün = Bereich, den der Normalenkegel-Test überspringt, grau = Bereich, der geprüft wird. Schalte das Culling aus, um zu sehen, wie die Zahl der geprüften Paare steigt. Im verheddigten Start entwirrt sich der Stoff mit ICM und bleibt ohne ICM verheddert.",
+    "EXPLANATION": {
+      "TITLE": "Stoff-Selbstdurchdringungen auf der CPU auflösen",
+      "ICM_TITLE": "Intersection Contour Minimization",
+      "ICM": "Ein Stoff, der sich schon durchdringt, lässt sich nicht reparieren, indem man Punkte auseinanderschiebt: Niemand weiß, welche Seite „außen“ ist. Volino und Magnenat-Thalmann (2006) betrachten stattdessen die Schnittkontur, die Kurve, an der der Stoff sich selbst durchschneidet. Sie bewegen den Stoff so, dass diese Kurve kürzer wird, bis sie verschwindet. Die Methode braucht keine Historie früherer Frames.",
+      "BVH_TITLE": "BVH mit Normalenkegel-Culling",
+      "BVH": "Jedes Dreieck gegen jedes andere zu testen ist viel zu langsam. Ein Baum aus Hüllkörpern (hier 18-DOPs, eine Box mit abgeschnittenen Ecken) findet die Paare, die nahe beieinander liegen. Ein Bereich, dessen Normalen grob in dieselbe Richtung zeigen (ein Kegel unter 90 Grad) und dessen Rand sich in dieser Blickrichtung nicht selbst überlappt, kann sich nicht selbst durchdringen. Der Baum überspringt so einen Bereich komplett (Volino 1994).",
+      "PROXIMITY_TITLE": "Abstands-Constraints",
+      "PROXIMITY": "Stoff, der sich noch nicht durchdringt, wird durch eine kleine Dicke auseinandergehalten: Punkt-Dreieck- und Kante-Kante-Constraints in der XPBD-Schleife schieben nahe Teile voneinander weg. Neben einer Kontur wird ein Schub verworfen, der gegen die ICM-Verschiebung arbeitet (Diplomarbeit, Abschnitt 5.3), sonst würde er die Durchdringung vergrößern.",
+      "DISCLAIMER": "Grenzen dieser Demo:",
+      "DISCLAIMER_1": "Konturen können in einer Sackgasse enden: einem lokalen Minimum, in dem die Kontur ohne einen weiteren Zug nicht kürzer wird. ICM kann so einen Zustand nicht entwirren. Die beiden Szenen bauen absichtlich keinen auf.",
+      "DISCLAIMER_2": "Treffer genau durch einen Eckpunkt oder eine Kante werden für einen Frame verworfen. Eine Kontur kann sich kurz in zwei Teile spalten.",
+      "DISCLAIMER_3": "Kontaktreibung (Kugel und Boden) und ein Geschwindigkeitslimit halten den Stoff ruhig. Ohne sie rutscht der Stoff zu einem Haufen zusammen, den die Erkennung nicht mehr sortieren kann."
+    },
+    "HOW_IT_WORKS": {
+      "INTRO": "Der Stoff ist ein Gitter aus Partikeln (CPU, XPBD, 10 Solver-Iterationen pro Frame): 50 × 50 = 2.500 im verheddigten Start, 35 × 35 = 1.225 beim fallenden Stoff. Pro Frame: Positionen vorhersagen, BVH anpassen, alle Schnitte finden, Konturen verfolgen, ICM-Verschiebung berechnen, dann die Solver-Iterationen laufen lassen und die Verschiebung in kleinen Teilen anwenden. Beim fallenden Stoff ist der Abstand zwischen den Lagen 3,5-mal größer (0,035 statt 0,01) und die Partikel haben mehr Abstand (0,145 statt 0,1), damit der Stoff gleich groß bleibt. Die nahen Paare werden mitten in den Iterationen noch einmal gesucht (nach 5 von 10), weil sich der fallende Stoff in einem Frame stark bewegt.",
+      "ITEM_1": "BVH: ein Blatt pro Gitterquadrat (2 Dreiecke). Jeder Knoten speichert ein 18-DOP und einen Normalenkegel. Der Baum behält seine Form und wird jeden Frame von unten nach oben angepasst. Knotenpaare mit überlappenden Volumen werden bis zu Blattpaaren abgelaufen.",
+      "ITEM_2": "Narrow Phase: Die 3 Kanten jedes Dreiecks werden gegen das andere Dreieck getestet (Möller-Trumbore). Zwei Treffer eines Dreieckspaars bilden ein Segment der Schnittkontur. Dreiecke mit gemeinsamem Eckpunkt werden übersprungen.",
+      "ITEM_3": "Konturen: Die Segmente werden zu geschlossenen Schleifen oder zu offenen Pfaden verbunden, die am Stoffrand enden. Die roten Linien zeigen sie.",
+      "ITEM_4": "ICM: Pro Kontur bekommt jedes Segmentende einen Gradienten (die Richtung, in der sich seine Kante relativ zum Dreieck bewegen muss, damit die Kontur kürzer wird). Die Gradienten werden mit Vorzeichen für die beiden Stoffseiten summiert und zu H = h0 · D / sqrt(|D|² + g0²) skaliert, mit h0 = ICM-Schritt · Abstand und g0 = 4. Die Kantenseite bewegt sich um +H, die Dreiecksseite um −H. Die Verschiebung wird nach inverser Masse und Abstand zum Trefferpunkt auf die Partikel verteilt und pro Partikel gemittelt.",
+      "ITEM_5": "Verheddigter Start: keine Schwerkraft, die obere Reihe ist fixiert, die untere Hälfte ist zurückgefaltet und ein Wulst quer darüber drückt durch die vordere Hälfte. Fallender Stoff: Der Stoff fällt auf eine kleine Kugel auf dem Boden. Beide Szenen nutzen denselben Solver. Schalte ICM aus, um den Unterschied zu sehen: beim fallenden Stoff bleiben dann die Schnitte (etwa hundert statt fast keine)."
+    }
+  },
   "BOIDS": {
     "TITLE": "Boids & Particle Life",
     "MODE_BOIDS": "Boids",
@@ -944,6 +981,10 @@
       "TITLE": "Stoffsimulation",
       "DESCRIPTION": "Stoffsimulation mit XPBD auf der GPU (WebGPU, Fallback WebGL + CPU)."
     },
+    "CLOTH_ICM": {
+      "TITLE": "Stoff-Selbstkollision (ICM)",
+      "DESCRIPTION": "Ein Stoff, der sich selbst durchschneidet, wird entwirrt, indem die Schnittkontur verkürzt wird (Volino und Magnenat-Thalmann 2006). CPU, BVH mit Normalenkegel-Culling, XPBD."
+    },
     "BOIDS": {
       "TITLE": "Boids & Particle Life",
       "DESCRIPTION": "Schwarmverhalten und emergentes Verhalten aus lokalen Regeln, berechnet auf der GPU (WebGPU, Fallback WebGL + CPU)."
@@ -1206,6 +1247,11 @@
     "STEPS_PER_SECOND": "Schritte/s",
     "MASS": "Masse",
     "PRESSURE_ITERATIONS": "Druck-Iterationen",
-    "ITERATION": "Iteration"
+    "ITERATION": "Iteration",
+    "INTERSECTIONS": "Schnittpunkte",
+    "CONTOURS": "Konturen",
+    "TESTED_PAIRS": "Geprüfte Dreieckspaare",
+    "CULLED_NODES": "Übersprungene BVH-Knoten",
+    "COLLISION_MS": "Kollisionszeit pro Frame"
   }
 }

src/assets/i18n/en.json

@@ -712,6 +712,43 @@
       "ITEM_5": "Wind: fixed direction (5, 0, 15) plus a sine flutter from position and time. The flutter part also adds a slight lift (y: flutter · 2)."
     }
   },
+  "CLOTH_ICM": {
+    "TITLE": "Cloth self-collision (ICM)",
+    "SCENE": "Scene",
+    "SCENE_TANGLED": "Tangled start",
+    "SCENE_DROP": "Dropped cloth",
+    "ICM_ON": "ICM On",
+    "ICM_OFF": "ICM Off",
+    "BVH_ON": "Show BVH",
+    "BVH_OFF": "Hide BVH",
+    "CULLING_ON": "Culling On",
+    "CULLING_OFF": "Culling Off",
+    "RESTART_SIMULATION": "Restart Simulation",
+    "ICM_STEP": "ICM step (x spacing)",
+    "BVH_DEPTH": "BVH level",
+    "HINT": "Red lines are the intersection contours. BVH boxes: green = region skipped by the normal-cone test, grey = region that is tested. Switch the culling off to see the number of tested pairs rise. In the tangled start the cloth untangles with ICM on and stays tangled with ICM off.",
+    "EXPLANATION": {
+      "TITLE": "Resolving cloth self-intersections on the CPU",
+      "ICM_TITLE": "Intersection Contour Minimization",
+      "ICM": "Cloth that already intersects itself cannot be repaired by pushing points apart: nobody knows which side is 'outside'. Volino and Magnenat-Thalmann (2006) look at the intersection contour instead, the curve where the cloth cuts through itself. They move the cloth so that this curve gets shorter until it vanishes. The method needs no history of earlier frames.",
+      "BVH_TITLE": "BVH with normal-cone culling",
+      "BVH": "Testing every triangle against every other one is far too slow. A tree of bounding volumes (here 18-DOPs, a box with cut corners) finds the pairs that are close. A region whose normals all point roughly the same way (a cone narrower than 90 degrees) and whose border does not overlap itself when seen along that direction cannot intersect itself. The tree skips such a region completely (Volino 1994).",
+      "PROXIMITY_TITLE": "Proximity constraints",
+      "PROXIMITY": "Cloth that does not intersect yet is kept apart by a small thickness: point-triangle and edge-edge constraints inside the XPBD loop push close parts away from each other. Next to a contour a push that works against the ICM displacement is dropped (diploma, section 5.3), otherwise it would widen the intersection.",
+      "DISCLAIMER": "Limits of this demo:",
+      "DISCLAIMER_1": "Contours can end in a dead end: a local minimum where the contour cannot get shorter without another pull. ICM cannot untangle such a state. The two scenes do not build one on purpose.",
+      "DISCLAIMER_2": "Hits exactly through a vertex or an edge are dropped for one frame. A contour can split into two pieces for a moment.",
+      "DISCLAIMER_3": "Contact friction (sphere and ground) and a speed limit keep the cloth calm. Without them the cloth slides into a pile that the detection can no longer sort out."
+    },
+    "HOW_IT_WORKS": {
+      "INTRO": "The cloth is a grid of particles (CPU, XPBD, 10 solver iterations per frame): 50 × 50 = 2,500 in the tangled start, 35 × 35 = 1,225 in the dropped cloth. Per frame: predict positions, refit the BVH, find all intersections, trace the contours, compute the ICM displacement, then run the solver iterations and apply the displacement in small parts. In the dropped cloth the gap between the layers is 3.5 times larger (0.035 instead of 0.01) and the particles are spaced wider (0.145 instead of 0.1) so the cloth keeps its size. The close pairs are found again in the middle of the iterations (after 5 of 10), because the falling cloth moves a lot within one frame.",
+      "ITEM_1": "BVH: one leaf per grid quad (2 triangles). Every node stores an 18-DOP and a normal cone. The tree keeps its shape and is refit bottom-up every frame. Node pairs with overlapping volumes are walked down to leaf pairs.",
+      "ITEM_2": "Narrow phase: the 3 edges of each triangle are tested against the other triangle (Möller-Trumbore). Two hits of one triangle pair form one segment of the intersection contour. Triangles that share a vertex are skipped.",
+      "ITEM_3": "Contours: the segments are joined into closed loops or into open paths that end at the cloth border. The red lines show them.",
+      "ITEM_4": "ICM: per contour, every segment end gets a gradient (the direction in which its edge must move relative to the triangle to shorten the contour). The gradients are summed with a sign for the two sides of the cloth and scaled to H = h0 · D / sqrt(|D|² + g0²), with h0 = ICM step · spacing and g0 = 4. The edge side moves by +H, the triangle side by −H. The displacement is spread over the particles by inverse mass and by distance to the hit point, and averaged per particle.",
+      "ITEM_5": "Tangled start: no gravity, the top row is pinned, the lower half is folded back and a ridge across it pokes through the front half. Dropped cloth: the cloth falls on a small sphere on the ground. Both scenes share the same solver. Switch ICM off to see the difference: in the dropped cloth the intersections then stay (about a hundred instead of almost none)."
+    }
+  },
   "BOIDS": {
     "TITLE": "Boids & Particle Life",
     "MODE_BOIDS": "Boids",
@@ -944,6 +981,10 @@
       "TITLE": "Cloth simulation",
       "DESCRIPTION": "Cloth simulation with XPBD on the GPU (WebGPU, fallback WebGL + CPU)."
     },
+    "CLOTH_ICM": {
+      "TITLE": "Cloth self-collision (ICM)",
+      "DESCRIPTION": "Cloth that cuts through itself is untangled by shrinking the intersection contour (Volino and Magnenat-Thalmann 2006). CPU, BVH with normal-cone culling, XPBD."
+    },
     "BOIDS": {
       "TITLE": "Boids & Particle Life",
       "DESCRIPTION": "Flocking and emergent behavior from local rules, computed on the GPU (WebGPU, fallback WebGL + CPU)."
@@ -1206,6 +1247,11 @@
     "STEPS_PER_SECOND": "Steps/s",
     "MASS": "Mass",
     "PRESSURE_ITERATIONS": "Pressure iterations",
-    "ITERATION": "Iteration"
+    "ITERATION": "Iteration",
+    "INTERSECTIONS": "Intersections",
+    "CONTOURS": "Contours",
+    "TESTED_PAIRS": "Tested triangle pairs",
+    "CULLED_NODES": "Culled BVH nodes",
+    "COLLISION_MS": "Collision time per frame"
   }
 }