diff --git a/lib/solvers/RectDiffSolver.ts b/lib/solvers/RectDiffSolver.ts
index eaed0f5..c2fe4a0 100644
--- a/lib/solvers/RectDiffSolver.ts
+++ b/lib/solvers/RectDiffSolver.ts
@@ -13,31 +13,39 @@ import {
   computeProgress,
 } from "./rectdiff/engine"
 import { rectsToMeshNodes } from "./rectdiff/rectsToMeshNodes"
+import type { GapFillOptions } from "./rectdiff/gapfill/types"
+import {
+  findUncoveredPoints,
+  calculateCoverage,
+} from "./rectdiff/gapfill/engine"
+import { GapFillSubSolver } from "./rectdiff/subsolvers/GapFillSubSolver"
 
-// A streaming, one-step-per-iteration solver.
-// Tests that call `solver.solve()` still work because BaseSolver.solve()
-// loops until this.solved flips true.
-
+/**
+ * A streaming, one-step-per-iteration solver for capacity mesh generation.
+ */
 export class RectDiffSolver extends BaseSolver {
   private srj: SimpleRouteJson
-  private mode: "grid" | "exact"
   private gridOptions: Partial<GridFill3DOptions>
+  private gapFillOptions: Partial<GapFillOptions>
   private state!: RectDiffState
   private _meshNodes: CapacityMeshNode[] = []
 
+  /** Active subsolver for GAP_FILL phases. */
+  declare activeSubSolver: GapFillSubSolver | null
+
   constructor(opts: {
     simpleRouteJson: SimpleRouteJson
-    mode?: "grid" | "exact"
     gridOptions?: Partial<GridFill3DOptions>
+    gapFillOptions?: Partial<GapFillOptions>
   }) {
     super()
     this.srj = opts.simpleRouteJson
-    this.mode = opts.mode ?? "grid"
     this.gridOptions = opts.gridOptions ?? {}
+    this.gapFillOptions = opts.gapFillOptions ?? {}
+    this.activeSubSolver = null
   }
 
   override _setup() {
-    // For now "exact" mode falls back to grid; keep switch if you add exact later.
     this.state = initState(this.srj, this.gridOptions)
     this.stats = {
       phase: this.state.phase,
@@ -45,12 +53,51 @@ export class RectDiffSolver extends BaseSolver {
     }
   }
 
-  /** IMPORTANT: exactly ONE small step per call */
+  /** Exactly ONE small step per call. */
   override _step() {
     if (this.state.phase === "GRID") {
       stepGrid(this.state)
     } else if (this.state.phase === "EXPANSION") {
       stepExpansion(this.state)
+    } else if (this.state.phase === "GAP_FILL") {
+      // Initialize gap fill subsolver if needed
+      if (
+        !this.activeSubSolver ||
+        !(this.activeSubSolver instanceof GapFillSubSolver)
+      ) {
+        const minTrace = this.srj.minTraceWidth || 0.15
+        const minGapSize = Math.max(0.01, minTrace / 10)
+        const boundsSize = Math.min(
+          this.state.bounds.width,
+          this.state.bounds.height,
+        )
+        this.activeSubSolver = new GapFillSubSolver({
+          placed: this.state.placed,
+          options: {
+            minWidth: minGapSize,
+            minHeight: minGapSize,
+            scanResolution: Math.max(0.05, boundsSize / 100),
+            ...this.gapFillOptions,
+          },
+          layerCtx: {
+            bounds: this.state.bounds,
+            layerCount: this.state.layerCount,
+            obstaclesByLayer: this.state.obstaclesByLayer,
+            placedByLayer: this.state.placedByLayer,
+          },
+        })
+      }
+
+      this.activeSubSolver.step()
+
+      if (this.activeSubSolver.solved) {
+        // Transfer results back to main state
+        const output = this.activeSubSolver.getOutput()
+        this.state.placed = output.placed
+        this.state.placedByLayer = output.placedByLayer
+        this.activeSubSolver = null
+        this.state.phase = "DONE"
+      }
     } else if (this.state.phase === "DONE") {
       // Finalize once
       if (!this.solved) {
@@ -65,47 +112,101 @@ export class RectDiffSolver extends BaseSolver {
     this.stats.phase = this.state.phase
     this.stats.gridIndex = this.state.gridIndex
     this.stats.placed = this.state.placed.length
+    if (this.activeSubSolver instanceof GapFillSubSolver) {
+      const output = this.activeSubSolver.getOutput()
+      this.stats.gapsFilled = output.filledCount
+    }
   }
 
-  // Let BaseSolver update this.progress automatically if present.
+  /** Compute solver progress (0 to 1). */
   computeProgress(): number {
-    return computeProgress(this.state)
+    if (this.solved || this.state.phase === "DONE") {
+      return 1
+    }
+    if (
+      this.state.phase === "GAP_FILL" &&
+      this.activeSubSolver instanceof GapFillSubSolver
+    ) {
+      return 0.85 + 0.1 * this.activeSubSolver.computeProgress()
+    }
+    return computeProgress(this.state) * 0.85
   }
 
   override getOutput(): { meshNodes: CapacityMeshNode[] } {
     return { meshNodes: this._meshNodes }
   }
 
-  // Helper to get color based on z layer
+  /** Get coverage percentage (0-1). */
+  getCoverage(sampleResolution: number = 0.05): number {
+    return calculateCoverage(
+      { sampleResolution },
+      {
+        bounds: this.state.bounds,
+        layerCount: this.state.layerCount,
+        obstaclesByLayer: this.state.obstaclesByLayer,
+        placedByLayer: this.state.placedByLayer,
+      },
+    )
+  }
+
+  /** Find uncovered points for debugging gaps. */
+  getUncoveredPoints(
+    sampleResolution: number = 0.05,
+  ): Array<{ x: number; y: number; z: number }> {
+    return findUncoveredPoints(
+      { sampleResolution },
+      {
+        bounds: this.state.bounds,
+        layerCount: this.state.layerCount,
+        obstaclesByLayer: this.state.obstaclesByLayer,
+        placedByLayer: this.state.placedByLayer,
+      },
+    )
+  }
+
+  /** Get color based on z layer for visualization. */
   private getColorForZLayer(zLayers: number[]): {
     fill: string
     stroke: string
   } {
     const minZ = Math.min(...zLayers)
     const colors = [
-      { fill: "#dbeafe", stroke: "#3b82f6" }, // blue (z=0)
-      { fill: "#fef3c7", stroke: "#f59e0b" }, // amber (z=1)
-      { fill: "#d1fae5", stroke: "#10b981" }, // green (z=2)
-      { fill: "#e9d5ff", stroke: "#a855f7" }, // purple (z=3)
-      { fill: "#fed7aa", stroke: "#f97316" }, // orange (z=4)
-      { fill: "#fecaca", stroke: "#ef4444" }, // red (z=5)
+      { fill: "#dbeafe", stroke: "#3b82f6" },
+      { fill: "#fef3c7", stroke: "#f59e0b" },
+      { fill: "#d1fae5", stroke: "#10b981" },
+      { fill: "#e9d5ff", stroke: "#a855f7" },
+      { fill: "#fed7aa", stroke: "#f97316" },
+      { fill: "#fecaca", stroke: "#ef4444" },
     ]
     return colors[minZ % colors.length]!
   }
 
-  // Streaming visualization: board + obstacles + current placements.
+  /** Streaming visualization: board + obstacles + current placements. */
   override visualize(): GraphicsObject {
+    // If a subsolver is active, delegate to its visualization
+    if (this.activeSubSolver) {
+      return this.activeSubSolver.visualize()
+    }
+
     const rects: NonNullable<GraphicsObject["rects"]> = []
     const points: NonNullable<GraphicsObject["points"]> = []
 
+    // Board bounds - use srj bounds which is always available
+    const boardBounds = {
+      minX: this.srj.bounds.minX,
+      maxX: this.srj.bounds.maxX,
+      minY: this.srj.bounds.minY,
+      maxY: this.srj.bounds.maxY,
+    }
+
     // board
     rects.push({
       center: {
-        x: (this.srj.bounds.minX + this.srj.bounds.maxX) / 2,
-        y: (this.srj.bounds.minY + this.srj.bounds.maxY) / 2,
+        x: (boardBounds.minX + boardBounds.maxX) / 2,
+        y: (boardBounds.minY + boardBounds.maxY) / 2,
       },
-      width: this.srj.bounds.maxX - this.srj.bounds.minX,
-      height: this.srj.bounds.maxY - this.srj.bounds.minY,
+      width: boardBounds.maxX - boardBounds.minX,
+      height: boardBounds.maxY - boardBounds.minY,
       fill: "none",
       stroke: "#111827",
       label: "board",
@@ -158,7 +259,7 @@ export class RectDiffSolver extends BaseSolver {
     }
 
     return {
-      title: "RectDiff (incremental)",
+      title: `RectDiff (${this.state?.phase ?? "init"})`,
       coordinateSystem: "cartesian",
       rects,
       points,
diff --git a/lib/solvers/rectdiff/candidates.ts b/lib/solvers/rectdiff/candidates.ts
index c597bf2..4bc59e1 100644
--- a/lib/solvers/rectdiff/candidates.ts
+++ b/lib/solvers/rectdiff/candidates.ts
@@ -2,13 +2,17 @@
 import type { Candidate3D, XYRect } from "./types"
 import { EPS, clamp, containsPoint, distancePointToRectEdges } from "./geometry"
 
-function isFullyOccupiedAllLayers(
-  x: number,
-  y: number,
-  layerCount: number,
-  obstaclesByLayer: XYRect[][],
-  placedByLayer: XYRect[][],
-): boolean {
+/**
+ * Check if a point is occupied on all layers.
+ */
+function isFullyOccupiedAllLayers(params: {
+  x: number
+  y: number
+  layerCount: number
+  obstaclesByLayer: XYRect[][]
+  placedByLayer: XYRect[][]
+}): boolean {
+  const { x, y, layerCount, obstaclesByLayer, placedByLayer } = params
   for (let z = 0; z < layerCount; z++) {
     const obs = obstaclesByLayer[z] ?? []
     const placed = placedByLayer[z] ?? []
@@ -20,14 +24,25 @@ function isFullyOccupiedAllLayers(
   return true
 }
 
-export function computeCandidates3D(
-  bounds: XYRect,
-  gridSize: number,
-  layerCount: number,
-  obstaclesByLayer: XYRect[][],
-  placedByLayer: XYRect[][], // all current nodes (soft + hard)
-  hardPlacedByLayer: XYRect[][], // only full-stack nodes, treated as hard
-): Candidate3D[] {
+/**
+ * Compute candidate seed points for a given grid size.
+ */
+export function computeCandidates3D(params: {
+  bounds: XYRect
+  gridSize: number
+  layerCount: number
+  obstaclesByLayer: XYRect[][]
+  placedByLayer: XYRect[][]
+  hardPlacedByLayer: XYRect[][]
+}): Candidate3D[] {
+  const {
+    bounds,
+    gridSize,
+    layerCount,
+    obstaclesByLayer,
+    placedByLayer,
+    hardPlacedByLayer,
+  } = params
   const out = new Map<string, Candidate3D>() // key by (x,y)
 
   for (let x = bounds.x; x < bounds.x + bounds.width; x += gridSize) {
@@ -44,13 +59,13 @@ export function computeCandidates3D(
 
       // New rule: Only drop if EVERY layer is occupied (by obstacle or node)
       if (
-        isFullyOccupiedAllLayers(
+        isFullyOccupiedAllLayers({
           x,
           y,
           layerCount,
           obstaclesByLayer,
           placedByLayer,
-        )
+        })
       )
         continue
 
@@ -58,16 +73,16 @@ export function computeCandidates3D(
       let bestSpan: number[] = []
       let bestZ = 0
       for (let z = 0; z < layerCount; z++) {
-        const s = longestFreeSpanAroundZ(
+        const s = longestFreeSpanAroundZ({
           x,
           y,
           z,
           layerCount,
-          1,
-          undefined, // no cap here
+          minSpan: 1,
+          maxSpan: undefined,
           obstaclesByLayer,
-          hardPlacedByLayer, // IMPORTANT: ignore soft nodes
-        )
+          placedByLayer: hardPlacedByLayer,
+        })
         if (s.length > bestSpan.length) {
           bestSpan = s
           bestZ = z
@@ -113,17 +128,30 @@ export function computeCandidates3D(
   return arr
 }
 
-/** Longest contiguous free span around z (optionally capped) */
-export function longestFreeSpanAroundZ(
-  x: number,
-  y: number,
-  z: number,
-  layerCount: number,
-  minSpan: number,
-  maxSpan: number | undefined,
-  obstaclesByLayer: XYRect[][],
-  placedByLayer: XYRect[][],
-): number[] {
+/**
+ * Find the longest contiguous free span around z (optionally capped).
+ */
+export function longestFreeSpanAroundZ(params: {
+  x: number
+  y: number
+  z: number
+  layerCount: number
+  minSpan: number
+  maxSpan: number | undefined
+  obstaclesByLayer: XYRect[][]
+  placedByLayer: XYRect[][]
+}): number[] {
+  const {
+    x,
+    y,
+    z,
+    layerCount,
+    minSpan,
+    maxSpan,
+    obstaclesByLayer,
+    placedByLayer,
+  } = params
+
   const isFreeAt = (layer: number) => {
     const blockers = [
       ...(obstaclesByLayer[layer] ?? []),
@@ -150,18 +178,25 @@ export function longestFreeSpanAroundZ(
   return res.length >= minSpan ? res : []
 }
 
+/**
+ * Compute default grid sizes based on bounds.
+ */
 export function computeDefaultGridSizes(bounds: XYRect): number[] {
   const ref = Math.max(bounds.width, bounds.height)
   return [ref / 8, ref / 16, ref / 32]
 }
 
-/** Compute exact uncovered segments along a 1D line given a list of covering intervals */
-function computeUncoveredSegments(
-  lineStart: number,
-  lineEnd: number,
-  coveringIntervals: Array<{ start: number; end: number }>,
-  minSegmentLength: number,
-): Array<{ start: number; end: number; center: number }> {
+/**
+ * Compute exact uncovered segments along a 1D line.
+ */
+function computeUncoveredSegments(params: {
+  lineStart: number
+  lineEnd: number
+  coveringIntervals: Array<{ start: number; end: number }>
+  minSegmentLength: number
+}): Array<{ start: number; end: number; center: number }> {
+  const { lineStart, lineEnd, coveringIntervals, minSegmentLength } = params
+
   if (coveringIntervals.length === 0) {
     const center = (lineStart + lineEnd) / 2
     return [{ start: lineStart, end: lineEnd, center }]
@@ -220,15 +255,26 @@ function computeUncoveredSegments(
   return uncovered
 }
 
-/** Exact edge analysis: find uncovered segments along board edges and blocker edges */
-export function computeEdgeCandidates3D(
-  bounds: XYRect,
-  minSize: number,
-  layerCount: number,
-  obstaclesByLayer: XYRect[][],
-  placedByLayer: XYRect[][], // all nodes
-  hardPlacedByLayer: XYRect[][], // full-stack nodes
-): Candidate3D[] {
+/**
+ * Compute edge candidates using exact edge analysis.
+ */
+export function computeEdgeCandidates3D(params: {
+  bounds: XYRect
+  minSize: number
+  layerCount: number
+  obstaclesByLayer: XYRect[][]
+  placedByLayer: XYRect[][]
+  hardPlacedByLayer: XYRect[][]
+}): Candidate3D[] {
+  const {
+    bounds,
+    minSize,
+    layerCount,
+    obstaclesByLayer,
+    placedByLayer,
+    hardPlacedByLayer,
+  } = params
+
   const out: Candidate3D[] = []
   // Use small inset from edges for placement
   const δ = Math.max(minSize * 0.15, EPS * 3)
@@ -237,13 +283,13 @@ export function computeEdgeCandidates3D(
     `${z}|${x.toFixed(6)}|${y.toFixed(6)}`
 
   function fullyOcc(x: number, y: number) {
-    return isFullyOccupiedAllLayers(
+    return isFullyOccupiedAllLayers({
       x,
       y,
       layerCount,
       obstaclesByLayer,
       placedByLayer,
-    )
+    })
   }
 
   function pushIfFree(x: number, y: number, z: number) {
@@ -273,16 +319,16 @@ export function computeEdgeCandidates3D(
     dedup.add(k)
 
     // Approximate z-span strength at this z (ignoring soft nodes)
-    const span = longestFreeSpanAroundZ(
+    const span = longestFreeSpanAroundZ({
       x,
       y,
       z,
       layerCount,
-      1,
-      undefined,
+      minSpan: 1,
+      maxSpan: undefined,
       obstaclesByLayer,
-      hardPlacedByLayer,
-    )
+      placedByLayer: hardPlacedByLayer,
+    })
     out.push({ x, y, z, distance: d, zSpanLen: span.length, isEdgeSeed: true })
   }
 
@@ -314,12 +360,12 @@ export function computeEdgeCandidates3D(
         end: Math.min(bounds.x + bounds.width, b.x + b.width),
       }))
     // Find uncovered segments that are large enough to potentially fill
-    const topUncovered = computeUncoveredSegments(
-      bounds.x + δ,
-      bounds.x + bounds.width - δ,
-      topCovering,
-      minSize * 0.5,
-    )
+    const topUncovered = computeUncoveredSegments({
+      lineStart: bounds.x + δ,
+      lineEnd: bounds.x + bounds.width - δ,
+      coveringIntervals: topCovering,
+      minSegmentLength: minSize * 0.5,
+    })
     for (const seg of topUncovered) {
       const segLen = seg.end - seg.start
       if (segLen >= minSize) {
@@ -340,12 +386,12 @@ export function computeEdgeCandidates3D(
         start: Math.max(bounds.x, b.x),
         end: Math.min(bounds.x + bounds.width, b.x + b.width),
       }))
-    const bottomUncovered = computeUncoveredSegments(
-      bounds.x + δ,
-      bounds.x + bounds.width - δ,
-      bottomCovering,
-      minSize * 0.5,
-    )
+    const bottomUncovered = computeUncoveredSegments({
+      lineStart: bounds.x + δ,
+      lineEnd: bounds.x + bounds.width - δ,
+      coveringIntervals: bottomCovering,
+      minSegmentLength: minSize * 0.5,
+    })
     for (const seg of bottomUncovered) {
       const segLen = seg.end - seg.start
       if (segLen >= minSize) {
@@ -365,12 +411,12 @@ export function computeEdgeCandidates3D(
         start: Math.max(bounds.y, b.y),
         end: Math.min(bounds.y + bounds.height, b.y + b.height),
       }))
-    const leftUncovered = computeUncoveredSegments(
-      bounds.y + δ,
-      bounds.y + bounds.height - δ,
-      leftCovering,
-      minSize * 0.5,
-    )
+    const leftUncovered = computeUncoveredSegments({
+      lineStart: bounds.y + δ,
+      lineEnd: bounds.y + bounds.height - δ,
+      coveringIntervals: leftCovering,
+      minSegmentLength: minSize * 0.5,
+    })
     for (const seg of leftUncovered) {
       const segLen = seg.end - seg.start
       if (segLen >= minSize) {
@@ -390,12 +436,12 @@ export function computeEdgeCandidates3D(
         start: Math.max(bounds.y, b.y),
         end: Math.min(bounds.y + bounds.height, b.y + b.height),
       }))
-    const rightUncovered = computeUncoveredSegments(
-      bounds.y + δ,
-      bounds.y + bounds.height - δ,
-      rightCovering,
-      minSize * 0.5,
-    )
+    const rightUncovered = computeUncoveredSegments({
+      lineStart: bounds.y + δ,
+      lineEnd: bounds.y + bounds.height - δ,
+      coveringIntervals: rightCovering,
+      minSegmentLength: minSize * 0.5,
+    })
     for (const seg of rightUncovered) {
       const segLen = seg.end - seg.start
       if (segLen >= minSize) {
@@ -420,12 +466,12 @@ export function computeEdgeCandidates3D(
             start: Math.max(b.y, bl.y),
             end: Math.min(b.y + b.height, bl.y + bl.height),
           }))
-        const obLeftUncovered = computeUncoveredSegments(
-          b.y,
-          b.y + b.height,
-          obLeftCovering,
-          minSize * 0.5,
-        )
+        const obLeftUncovered = computeUncoveredSegments({
+          lineStart: b.y,
+          lineEnd: b.y + b.height,
+          coveringIntervals: obLeftCovering,
+          minSegmentLength: minSize * 0.5,
+        })
         for (const seg of obLeftUncovered) {
           pushIfFree(obLeftX, seg.center, z)
         }
@@ -445,12 +491,12 @@ export function computeEdgeCandidates3D(
             start: Math.max(b.y, bl.y),
             end: Math.min(b.y + b.height, bl.y + bl.height),
           }))
-        const obRightUncovered = computeUncoveredSegments(
-          b.y,
-          b.y + b.height,
-          obRightCovering,
-          minSize * 0.5,
-        )
+        const obRightUncovered = computeUncoveredSegments({
+          lineStart: b.y,
+          lineEnd: b.y + b.height,
+          coveringIntervals: obRightCovering,
+          minSegmentLength: minSize * 0.5,
+        })
         for (const seg of obRightUncovered) {
           pushIfFree(obRightX, seg.center, z)
         }
@@ -467,12 +513,12 @@ export function computeEdgeCandidates3D(
             start: Math.max(b.x, bl.x),
             end: Math.min(b.x + b.width, bl.x + bl.width),
           }))
-        const obTopUncovered = computeUncoveredSegments(
-          b.x,
-          b.x + b.width,
-          obTopCovering,
-          minSize * 0.5,
-        )
+        const obTopUncovered = computeUncoveredSegments({
+          lineStart: b.x,
+          lineEnd: b.x + b.width,
+          coveringIntervals: obTopCovering,
+          minSegmentLength: minSize * 0.5,
+        })
         for (const seg of obTopUncovered) {
           pushIfFree(seg.center, obTopY, z)
         }
@@ -493,12 +539,12 @@ export function computeEdgeCandidates3D(
             start: Math.max(b.x, bl.x),
             end: Math.min(b.x + b.width, bl.x + bl.width),
           }))
-        const obBottomUncovered = computeUncoveredSegments(
-          b.x,
-          b.x + b.width,
-          obBottomCovering,
-          minSize * 0.5,
-        )
+        const obBottomUncovered = computeUncoveredSegments({
+          lineStart: b.x,
+          lineEnd: b.x + b.width,
+          coveringIntervals: obBottomCovering,
+          minSegmentLength: minSize * 0.5,
+        })
         for (const seg of obBottomUncovered) {
           pushIfFree(seg.center, obBottomY, z)
         }
diff --git a/lib/solvers/rectdiff/engine.ts b/lib/solvers/rectdiff/engine.ts
index 68a92c2..c497409 100644
--- a/lib/solvers/rectdiff/engine.ts
+++ b/lib/solvers/rectdiff/engine.ts
@@ -22,6 +22,9 @@ import {
 } from "./geometry"
 import { buildZIndexMap, obstacleToXYRect, obstacleZs } from "./layers"
 
+/**
+ * Initialize the RectDiff solver state from SimpleRouteJson.
+ */
 export function initState(
   srj: SimpleRouteJson,
   opts: Partial<GridFill3DOptions>,
@@ -106,7 +109,9 @@ export function initState(
   }
 }
 
-// Build per-layer list of "hard" placed rects = nodes spanning all layers
+/**
+ * Build per-layer list of "hard" placed rects (nodes spanning all layers).
+ */
 function buildHardPlacedByLayer(state: RectDiffState): XYRect[][] {
   const out: XYRect[][] = Array.from({ length: state.layerCount }, () => [])
   for (const p of state.placed) {
@@ -117,23 +122,27 @@ function buildHardPlacedByLayer(state: RectDiffState): XYRect[][] {
   return out
 }
 
+/**
+ * Check if a point is occupied on ALL layers.
+ */
 function isFullyOccupiedAtPoint(
   state: RectDiffState,
-  x: number,
-  y: number,
+  point: { x: number; y: number },
 ): boolean {
   for (let z = 0; z < state.layerCount; z++) {
     const obs = state.obstaclesByLayer[z] ?? []
     const placed = state.placedByLayer[z] ?? []
     const occ =
-      obs.some((b) => containsPoint(b, x, y)) ||
-      placed.some((b) => containsPoint(b, x, y))
+      obs.some((b) => containsPoint(b, point.x, point.y)) ||
+      placed.some((b) => containsPoint(b, point.x, point.y))
     if (!occ) return false
   }
   return true
 }
 
-/** Shrink/split any *soft* (non-full-stack) nodes overlapped by `newIndex` */
+/**
+ * Shrink/split any soft (non-full-stack) nodes overlapped by the newcomer.
+ */
 function resizeSoftOverlaps(state: RectDiffState, newIndex: number) {
   const newcomer = state.placed[newIndex]!
   const { rect: newR, zLayers: newZs } = newcomer
@@ -198,7 +207,9 @@ function resizeSoftOverlaps(state: RectDiffState, newIndex: number) {
   }
 }
 
-/** One micro-step during the GRID phase: handle (or fetch) exactly one candidate */
+/**
+ * One micro-step during the GRID phase: handle exactly one candidate.
+ */
 export function stepGrid(state: RectDiffState): void {
   const {
     gridSizes,
@@ -216,14 +227,14 @@ export function stepGrid(state: RectDiffState): void {
 
   // Ensure candidates exist for this grid
   if (state.candidates.length === 0 && state.consumedSeedsThisGrid === 0) {
-    state.candidates = computeCandidates3D(
-      state.bounds,
-      grid,
-      state.layerCount,
-      state.obstaclesByLayer,
-      state.placedByLayer, // all nodes (soft + hard) for fully-occupied test
-      hardPlacedByLayer, // hard blockers for ranking/span
-    )
+    state.candidates = computeCandidates3D({
+      bounds: state.bounds,
+      gridSize: grid,
+      layerCount: state.layerCount,
+      obstaclesByLayer: state.obstaclesByLayer,
+      placedByLayer: state.placedByLayer,
+      hardPlacedByLayer,
+    })
     state.totalSeedsThisGrid = state.candidates.length
     state.consumedSeedsThisGrid = 0
   }
@@ -238,14 +249,14 @@ export function stepGrid(state: RectDiffState): void {
     } else {
       if (!state.edgeAnalysisDone) {
         const minSize = Math.min(minSingle.width, minSingle.height)
-        state.candidates = computeEdgeCandidates3D(
-          state.bounds,
+        state.candidates = computeEdgeCandidates3D({
+          bounds: state.bounds,
           minSize,
-          state.layerCount,
-          state.obstaclesByLayer,
-          state.placedByLayer, // for fully-occupied test
+          layerCount: state.layerCount,
+          obstaclesByLayer: state.obstaclesByLayer,
+          placedByLayer: state.placedByLayer,
           hardPlacedByLayer,
-        )
+        })
         state.edgeAnalysisDone = true
         state.totalSeedsThisGrid = state.candidates.length
         state.consumedSeedsThisGrid = 0
@@ -262,16 +273,16 @@ export function stepGrid(state: RectDiffState): void {
   state.consumedSeedsThisGrid += 1
 
   // Evaluate attempts — multi-layer span first (computed ignoring soft nodes)
-  const span = longestFreeSpanAroundZ(
-    cand.x,
-    cand.y,
-    cand.z,
-    state.layerCount,
-    minMulti.minLayers,
-    maxMultiLayerSpan,
-    state.obstaclesByLayer,
-    hardPlacedByLayer, // ignore soft nodes for span
-  )
+  const span = longestFreeSpanAroundZ({
+    x: cand.x,
+    y: cand.y,
+    z: cand.z,
+    layerCount: state.layerCount,
+    minSpan: minMulti.minLayers,
+    maxSpan: maxMultiLayerSpan,
+    obstaclesByLayer: state.obstaclesByLayer,
+    placedByLayer: hardPlacedByLayer,
+  })
 
   const attempts: Array<{
     kind: "multi" | "single"
@@ -303,16 +314,16 @@ export function stepGrid(state: RectDiffState): void {
       if (hardPlacedByLayer[z]) hardBlockers.push(...hardPlacedByLayer[z]!)
     }
 
-    const rect = expandRectFromSeed(
-      cand.x,
-      cand.y,
-      grid,
-      state.bounds,
-      hardBlockers, // soft nodes DO NOT block expansion
+    const rect = expandRectFromSeed({
+      startX: cand.x,
+      startY: cand.y,
+      gridSize: grid,
+      bounds: state.bounds,
+      blockers: hardBlockers,
       initialCellRatio,
       maxAspectRatio,
-      attempt.minReq,
-    )
+      minReq: attempt.minReq,
+    })
     if (!rect) continue
 
     // Place the new node
@@ -325,7 +336,7 @@ export function stepGrid(state: RectDiffState): void {
 
     // New: relax candidate culling — only drop seeds that became fully occupied
     state.candidates = state.candidates.filter(
-      (c) => !isFullyOccupiedAtPoint(state, c.x, c.y),
+      (c) => !isFullyOccupiedAtPoint(state, { x: c.x, y: c.y }),
     )
 
     return // processed one candidate
@@ -334,10 +345,13 @@ export function stepGrid(state: RectDiffState): void {
   // Neither attempt worked; drop this candidate for now.
 }
 
-/** One micro-step during the EXPANSION phase: expand exactly one placed rect */
+/**
+ * One micro-step during the EXPANSION phase: expand exactly one placed rect.
+ */
 export function stepExpansion(state: RectDiffState): void {
   if (state.expansionIndex >= state.placed.length) {
-    state.phase = "DONE"
+    // Transition to gap fill phase instead of done
+    state.phase = "GAP_FILL"
     return
   }
 
@@ -355,16 +369,16 @@ export function stepExpansion(state: RectDiffState): void {
   }
 
   const oldRect = p.rect
-  const expanded = expandRectFromSeed(
-    p.rect.x + p.rect.width / 2,
-    p.rect.y + p.rect.height / 2,
-    lastGrid,
-    state.bounds,
-    hardBlockers,
-    0, // seed bias off
-    null, // no aspect cap in expansion pass
-    { width: p.rect.width, height: p.rect.height },
-  )
+  const expanded = expandRectFromSeed({
+    startX: p.rect.x + p.rect.width / 2,
+    startY: p.rect.y + p.rect.height / 2,
+    gridSize: lastGrid,
+    bounds: state.bounds,
+    blockers: hardBlockers,
+    initialCellRatio: 0,
+    maxAspectRatio: null,
+    minReq: { width: p.rect.width, height: p.rect.height },
+  })
 
   if (expanded) {
     // Update placement + per-layer index (replace old rect object)
@@ -382,6 +396,9 @@ export function stepExpansion(state: RectDiffState): void {
   state.expansionIndex += 1
 }
 
+/**
+ * Finalize placed rectangles into output format.
+ */
 export function finalizeRects(state: RectDiffState): Rect3d[] {
   // Convert all placed (free space) nodes to output format
   const out: Rect3d[] = state.placed.map((p) => ({
@@ -423,7 +440,9 @@ export function finalizeRects(state: RectDiffState): Rect3d[] {
   return out
 }
 
-/** Optional: rough progress number for BaseSolver.progress */
+/**
+ * Calculate rough progress number for BaseSolver.progress.
+ */
 export function computeProgress(state: RectDiffState): number {
   const grids = state.options.gridSizes.length
   if (state.phase === "GRID") {
diff --git a/lib/solvers/rectdiff/gapfill/detection.ts b/lib/solvers/rectdiff/gapfill/detection.ts
new file mode 100644
index 0000000..7126d89
--- /dev/null
+++ b/lib/solvers/rectdiff/gapfill/detection.ts
@@ -0,0 +1,3 @@
+// lib/solvers/rectdiff/gapfill/detection.ts
+export * from "./detection/findAllGaps"
+// findGapsOnLayer is not exported as it's only used by findAllGaps
diff --git a/lib/solvers/rectdiff/gapfill/detection/deduplicateGaps.ts b/lib/solvers/rectdiff/gapfill/detection/deduplicateGaps.ts
new file mode 100644
index 0000000..0f45164
--- /dev/null
+++ b/lib/solvers/rectdiff/gapfill/detection/deduplicateGaps.ts
@@ -0,0 +1,28 @@
+// lib/solvers/rectdiff/gapfill/detection/deduplicateGaps.ts
+import { rectsEqual } from "../../../../../utils/rectsEqual"
+import { rectsOverlap } from "../../../../../utils/rectsOverlap"
+import type { GapRegion } from "../types"
+
+export function deduplicateGaps(gaps: GapRegion[]): GapRegion[] {
+  const result: GapRegion[] = []
+
+  for (const gap of gaps) {
+    // Check if we already have a gap at the same location with overlapping layers
+    const existing = result.find(
+      (g) =>
+        rectsEqual(g.rect, gap.rect) ||
+        (rectsOverlap(g.rect, gap.rect) &&
+          gap.zLayers.some((z) => g.zLayers.includes(z))),
+    )
+
+    if (!existing) {
+      result.push(gap)
+    } else if (gap.zLayers.length > existing.zLayers.length) {
+      // Replace with the one that has more layers
+      const idx = result.indexOf(existing)
+      result[idx] = gap
+    }
+  }
+
+  return result
+}
diff --git a/lib/solvers/rectdiff/gapfill/detection/findAllGaps.ts b/lib/solvers/rectdiff/gapfill/detection/findAllGaps.ts
new file mode 100644
index 0000000..d4cc9ec
--- /dev/null
+++ b/lib/solvers/rectdiff/gapfill/detection/findAllGaps.ts
@@ -0,0 +1,83 @@
+// lib/solvers/rectdiff/gapfill/detection/findAllGaps.ts
+import type { XYRect } from "../../types"
+import type { GapRegion, LayerContext } from "../types"
+import { EPS } from "../../geometry"
+import { findGapsOnLayer } from "./findGapsOnLayer"
+import { rectsOverlap } from "../../../../../utils/rectsOverlap"
+import { deduplicateGaps } from "./deduplicateGaps"
+
+/**
+ * Find gaps across all layers and return GapRegions with z-layer info.
+ */
+export function findAllGaps(
+  {
+    scanResolution,
+    minWidth,
+    minHeight,
+  }: {
+    scanResolution: number
+    minWidth: number
+    minHeight: number
+  },
+  ctx: LayerContext,
+): GapRegion[] {
+  const { bounds, layerCount, obstaclesByLayer, placedByLayer } = ctx
+
+  // Find gaps on each layer
+  const gapsByLayer: XYRect[][] = []
+  for (let z = 0; z < layerCount; z++) {
+    const obstacles = obstaclesByLayer[z] ?? []
+    const placed = placedByLayer[z] ?? []
+    const gaps = findGapsOnLayer({ bounds, obstacles, placed, scanResolution })
+    gapsByLayer.push(gaps)
+  }
+
+  // Convert to GapRegions with z-layer info
+  const allGaps: GapRegion[] = []
+
+  for (let z = 0; z < layerCount; z++) {
+    for (const gap of gapsByLayer[z]!) {
+      // Filter out gaps that are too small
+      if (gap.width < minWidth - EPS || gap.height < minHeight - EPS) continue
+
+      // Check if this gap exists on adjacent layers too
+      const zLayers = [z]
+
+      // Look up
+      for (let zu = z + 1; zu < layerCount; zu++) {
+        const hasOverlap = gapsByLayer[zu]!.some((g) => rectsOverlap(g, gap))
+        if (hasOverlap) zLayers.push(zu)
+        else break
+      }
+
+      // Look down (if z > 0 and not already counted)
+      for (let zd = z - 1; zd >= 0; zd--) {
+        const hasOverlap = gapsByLayer[zd]!.some((g) => rectsOverlap(g, gap))
+        if (hasOverlap && !zLayers.includes(zd)) zLayers.unshift(zd)
+        else break
+      }
+
+      allGaps.push({
+        rect: gap,
+        zLayers: zLayers.sort((a, b) => a - b),
+        centerX: gap.x + gap.width / 2,
+        centerY: gap.y + gap.height / 2,
+        area: gap.width * gap.height,
+      })
+    }
+  }
+
+  // Deduplicate gaps that are essentially the same across layers
+  const deduped = deduplicateGaps(allGaps)
+
+  // Sort by priority: prefer larger gaps and multi-layer gaps
+  deduped.sort((a, b) => {
+    // Prefer multi-layer gaps
+    const layerDiff = b.zLayers.length - a.zLayers.length
+    if (layerDiff !== 0) return layerDiff
+    // Then prefer larger area
+    return b.area - a.area
+  })
+
+  return deduped
+}
diff --git a/lib/solvers/rectdiff/gapfill/detection/findGapsOnLayer.ts b/lib/solvers/rectdiff/gapfill/detection/findGapsOnLayer.ts
new file mode 100644
index 0000000..895ed9d
--- /dev/null
+++ b/lib/solvers/rectdiff/gapfill/detection/findGapsOnLayer.ts
@@ -0,0 +1,100 @@
+// lib/solvers/rectdiff/gapfill/detection/findGapsOnLayer.ts
+import type { XYRect } from "../../types"
+import { EPS } from "../../geometry"
+
+import { mergeUncoveredCells } from "./mergeUncoveredCells"
+
+/**
+ * Sweep-line algorithm to find maximal uncovered rectangles on a single layer.
+ */
+export function findGapsOnLayer({
+  bounds,
+  obstacles,
+  placed,
+  scanResolution,
+}: {
+  bounds: XYRect
+  obstacles: XYRect[]
+  placed: XYRect[]
+  scanResolution: number
+}): XYRect[] {
+  const blockers = [...obstacles, ...placed]
+
+  // Collect all unique x-coordinates
+  const xCoords = new Set<number>()
+  xCoords.add(bounds.x)
+  xCoords.add(bounds.x + bounds.width)
+
+  for (const b of blockers) {
+    if (b.x > bounds.x && b.x < bounds.x + bounds.width) {
+      xCoords.add(b.x)
+    }
+    if (b.x + b.width > bounds.x && b.x + b.width < bounds.x + bounds.width) {
+      xCoords.add(b.x + b.width)
+    }
+  }
+
+  // Also add intermediate points based on scan resolution
+  for (let x = bounds.x; x <= bounds.x + bounds.width; x += scanResolution) {
+    xCoords.add(x)
+  }
+
+  const sortedX = Array.from(xCoords).sort((a, b) => a - b)
+
+  // Similarly for y-coordinates
+  const yCoords = new Set<number>()
+  yCoords.add(bounds.y)
+  yCoords.add(bounds.y + bounds.height)
+
+  for (const b of blockers) {
+    if (b.y > bounds.y && b.y < bounds.y + bounds.height) {
+      yCoords.add(b.y)
+    }
+    if (
+      b.y + b.height > bounds.y &&
+      b.y + b.height < bounds.y + bounds.height
+    ) {
+      yCoords.add(b.y + b.height)
+    }
+  }
+
+  for (let y = bounds.y; y <= bounds.y + bounds.height; y += scanResolution) {
+    yCoords.add(y)
+  }
+
+  const sortedY = Array.from(yCoords).sort((a, b) => a - b)
+
+  // Build a grid of cells and mark which are uncovered
+  const uncoveredCells: Array<{ x: number; y: number; w: number; h: number }> =
+    []
+
+  for (let i = 0; i < sortedX.length - 1; i++) {
+    for (let j = 0; j < sortedY.length - 1; j++) {
+      const cellX = sortedX[i]!
+      const cellY = sortedY[j]!
+      const cellW = sortedX[i + 1]! - cellX
+      const cellH = sortedY[j + 1]! - cellY
+
+      if (cellW <= EPS || cellH <= EPS) continue
+
+      // Check if this cell is covered by any blocker
+      const cellCenterX = cellX + cellW / 2
+      const cellCenterY = cellY + cellH / 2
+
+      const isCovered = blockers.some(
+        (b) =>
+          cellCenterX >= b.x - EPS &&
+          cellCenterX <= b.x + b.width + EPS &&
+          cellCenterY >= b.y - EPS &&
+          cellCenterY <= b.y + b.height + EPS,
+      )
+
+      if (!isCovered) {
+        uncoveredCells.push({ x: cellX, y: cellY, w: cellW, h: cellH })
+      }
+    }
+  }
+
+  // Merge adjacent uncovered cells into maximal rectangles
+  return mergeUncoveredCells(uncoveredCells)
+}
diff --git a/lib/solvers/rectdiff/gapfill/detection/mergeUncoveredCells.ts b/lib/solvers/rectdiff/gapfill/detection/mergeUncoveredCells.ts
new file mode 100644
index 0000000..940c75c
--- /dev/null
+++ b/lib/solvers/rectdiff/gapfill/detection/mergeUncoveredCells.ts
@@ -0,0 +1,75 @@
+// lib/solvers/rectdiff/gapfill/detection/mergeUncoveredCells.ts
+import type { XYRect } from "../../types"
+import { EPS } from "../../geometry"
+
+/**
+ * Merge adjacent uncovered cells into larger rectangles using a greedy approach.
+ */
+export function mergeUncoveredCells(
+  cells: Array<{ x: number; y: number; w: number; h: number }>,
+): XYRect[] {
+  if (cells.length === 0) return []
+
+  // Group cells by their left edge and width
+  const byXW = new Map<string, typeof cells>()
+  for (const c of cells) {
+    const key = `${c.x.toFixed(9)}|${c.w.toFixed(9)}`
+    const arr = byXW.get(key) ?? []
+    arr.push(c)
+    byXW.set(key, arr)
+  }
+
+  // Within each vertical strip, merge adjacent cells
+  const verticalStrips: XYRect[] = []
+  for (const stripCells of byXW.values()) {
+    // Sort by y
+    stripCells.sort((a, b) => a.y - b.y)
+
+    let current: XYRect | null = null
+    for (const c of stripCells) {
+      if (!current) {
+        current = { x: c.x, y: c.y, width: c.w, height: c.h }
+      } else if (Math.abs(current.y + current.height - c.y) < EPS) {
+        // Adjacent vertically, merge
+        current.height += c.h
+      } else {
+        // Gap, save current and start new
+        verticalStrips.push(current)
+        current = { x: c.x, y: c.y, width: c.w, height: c.h }
+      }
+    }
+    if (current) verticalStrips.push(current)
+  }
+
+  // Now try to merge horizontal strips with same y and height
+  const byYH = new Map<string, XYRect[]>()
+  for (const r of verticalStrips) {
+    const key = `${r.y.toFixed(9)}|${r.height.toFixed(9)}`
+    const arr = byYH.get(key) ?? []
+    arr.push(r)
+    byYH.set(key, arr)
+  }
+
+  const merged: XYRect[] = []
+  for (const rowRects of byYH.values()) {
+    // Sort by x
+    rowRects.sort((a, b) => a.x - b.x)
+
+    let current: XYRect | null = null
+    for (const r of rowRects) {
+      if (!current) {
+        current = { ...r }
+      } else if (Math.abs(current.x + current.width - r.x) < EPS) {
+        // Adjacent horizontally, merge
+        current.width += r.width
+      } else {
+        // Gap, save current and start new
+        merged.push(current)
+        current = { ...r }
+      }
+    }
+    if (current) merged.push(current)
+  }
+
+  return merged
+}
diff --git a/lib/solvers/rectdiff/gapfill/engine.ts b/lib/solvers/rectdiff/gapfill/engine.ts
new file mode 100644
index 0000000..a440663
--- /dev/null
+++ b/lib/solvers/rectdiff/gapfill/engine.ts
@@ -0,0 +1,7 @@
+// lib/solvers/rectdiff/gapfill/engine.ts
+export * from "./engine/calculateCoverage"
+export * from "./engine/findUncoveredPoints"
+export * from "./engine/getGapFillProgress"
+export * from "./engine/initGapFillState"
+
+export * from "./engine/stepGapFill"
diff --git a/lib/solvers/rectdiff/gapfill/engine/addPlacement.ts b/lib/solvers/rectdiff/gapfill/engine/addPlacement.ts
new file mode 100644
index 0000000..9c7bda7
--- /dev/null
+++ b/lib/solvers/rectdiff/gapfill/engine/addPlacement.ts
@@ -0,0 +1,27 @@
+// lib/solvers/rectdiff/gapfill/engine/addPlacement.ts
+import type { Placed3D, XYRect } from "../../types"
+import type { GapFillState } from "../types"
+
+/**
+ * Add a new placement to the state.
+ */
+export function addPlacement(
+  state: GapFillState,
+  {
+    rect,
+    zLayers,
+  }: {
+    rect: XYRect
+    zLayers: number[]
+  },
+): void {
+  const placed: Placed3D = { rect, zLayers: [...zLayers] }
+  state.placed.push(placed)
+
+  for (const z of zLayers) {
+    if (!state.placedByLayer[z]) {
+      state.placedByLayer[z] = []
+    }
+    state.placedByLayer[z]!.push(rect)
+  }
+}
diff --git a/lib/solvers/rectdiff/gapfill/engine/calculateCoverage.ts b/lib/solvers/rectdiff/gapfill/engine/calculateCoverage.ts
new file mode 100644
index 0000000..d9de4b3
--- /dev/null
+++ b/lib/solvers/rectdiff/gapfill/engine/calculateCoverage.ts
@@ -0,0 +1,44 @@
+// lib/solvers/rectdiff/gapfill/engine/calculateCoverage.ts
+import type { LayerContext } from "../types"
+
+/**
+ * Calculate coverage percentage (0-1).
+ */
+export function calculateCoverage(
+  { sampleResolution = 0.1 }: { sampleResolution?: number },
+  ctx: LayerContext,
+): number {
+  const { bounds, layerCount, obstaclesByLayer, placedByLayer } = ctx
+
+  let totalPoints = 0
+  let coveredPoints = 0
+
+  for (let z = 0; z < layerCount; z++) {
+    const obstacles = obstaclesByLayer[z] ?? []
+    const placed = placedByLayer[z] ?? []
+    const allRects = [...obstacles, ...placed]
+
+    for (
+      let x = bounds.x;
+      x <= bounds.x + bounds.width;
+      x += sampleResolution
+    ) {
+      for (
+        let y = bounds.y;
+        y <= bounds.y + bounds.height;
+        y += sampleResolution
+      ) {
+        totalPoints++
+
+        const isCovered = allRects.some(
+          (r) =>
+            x >= r.x && x <= r.x + r.width && y >= r.y && y <= r.y + r.height,
+        )
+
+        if (isCovered) coveredPoints++
+      }
+    }
+  }
+
+  return totalPoints > 0 ? coveredPoints / totalPoints : 1
+}
diff --git a/lib/solvers/rectdiff/gapfill/engine/findUncoveredPoints.ts b/lib/solvers/rectdiff/gapfill/engine/findUncoveredPoints.ts
new file mode 100644
index 0000000..c0a09b6
--- /dev/null
+++ b/lib/solvers/rectdiff/gapfill/engine/findUncoveredPoints.ts
@@ -0,0 +1,43 @@
+// lib/solvers/rectdiff/gapfill/engine/findUncoveredPoints.ts
+import type { LayerContext } from "../types"
+
+/**
+ * Find uncovered points for debugging gaps.
+ */
+export function findUncoveredPoints(
+  { sampleResolution = 0.05 }: { sampleResolution?: number },
+  ctx: LayerContext,
+): Array<{ x: number; y: number; z: number }> {
+  const { bounds, layerCount, obstaclesByLayer, placedByLayer } = ctx
+
+  const uncovered: Array<{ x: number; y: number; z: number }> = []
+
+  for (let z = 0; z < layerCount; z++) {
+    const obstacles = obstaclesByLayer[z] ?? []
+    const placed = placedByLayer[z] ?? []
+    const allRects = [...obstacles, ...placed]
+
+    for (
+      let x = bounds.x;
+      x <= bounds.x + bounds.width;
+      x += sampleResolution
+    ) {
+      for (
+        let y = bounds.y;
+        y <= bounds.y + bounds.height;
+        y += sampleResolution
+      ) {
+        const isCovered = allRects.some(
+          (r) =>
+            x >= r.x && x <= r.x + r.width && y >= r.y && y <= r.y + r.height,
+        )
+
+        if (!isCovered) {
+          uncovered.push({ x, y, z })
+        }
+      }
+    }
+  }
+
+  return uncovered
+}
diff --git a/lib/solvers/rectdiff/gapfill/engine/getGapFillProgress.ts b/lib/solvers/rectdiff/gapfill/engine/getGapFillProgress.ts
new file mode 100644
index 0000000..6d4d65a
--- /dev/null
+++ b/lib/solvers/rectdiff/gapfill/engine/getGapFillProgress.ts
@@ -0,0 +1,42 @@
+// lib/solvers/rectdiff/gapfill/engine/getGapFillProgress.ts
+import type { GapFillState } from "../types"
+
+/**
+ * Get progress as a number between 0 and 1.
+ * Accounts for four-stage processing (scan → select → expand → place for each gap).
+ */
+export function getGapFillProgress(state: GapFillState): number {
+  if (state.done) return 1
+
+  const iterationProgress = state.iteration / state.options.maxIterations
+  const gapProgress =
+    state.gapsFound.length > 0 ? state.gapIndex / state.gapsFound.length : 0
+
+  // Add sub-progress within current gap based on stage
+  let stageProgress = 0
+  switch (state.stage) {
+    case "scan":
+      stageProgress = 0
+      break
+    case "select":
+      stageProgress = 0.25
+      break
+    case "expand":
+      stageProgress = 0.5
+      break
+    case "place":
+      stageProgress = 0.75
+      break
+  }
+
+  const gapStageProgress =
+    state.gapsFound.length > 0
+      ? stageProgress / (state.gapsFound.length * 4) // 4 stages per gap
+      : 0
+
+  return Math.min(
+    0.999,
+    iterationProgress +
+      (gapProgress + gapStageProgress) / state.options.maxIterations,
+  )
+}
diff --git a/lib/solvers/rectdiff/gapfill/engine/initGapFillState.ts b/lib/solvers/rectdiff/gapfill/engine/initGapFillState.ts
new file mode 100644
index 0000000..8f5d886
--- /dev/null
+++ b/lib/solvers/rectdiff/gapfill/engine/initGapFillState.ts
@@ -0,0 +1,57 @@
+// lib/solvers/rectdiff/gapfill/engine/initGapFillState.ts
+import type { Placed3D } from "../../types"
+import type { GapFillState, GapFillOptions, LayerContext } from "../types"
+
+const DEFAULT_OPTIONS: GapFillOptions = {
+  minWidth: 0.1,
+  minHeight: 0.1,
+  maxIterations: 10,
+  targetCoverage: 0.999,
+  scanResolution: 0.5,
+}
+
+/**
+ * Initialize the gap fill state from existing rectdiff state.
+ */
+export function initGapFillState(
+  {
+    placed,
+    options,
+  }: {
+    placed: Placed3D[]
+    options?: Partial<GapFillOptions>
+  },
+  ctx: LayerContext,
+): GapFillState {
+  const opts = { ...DEFAULT_OPTIONS, ...options }
+
+  // Deep copy placed arrays to avoid mutation issues
+  const placedCopy = placed.map((p) => ({
+    rect: { ...p.rect },
+    zLayers: [...p.zLayers],
+  }))
+
+  const placedByLayerCopy = ctx.placedByLayer.map((layer) =>
+    layer.map((r) => ({ ...r })),
+  )
+
+  return {
+    bounds: { ...ctx.bounds },
+    layerCount: ctx.layerCount,
+    obstaclesByLayer: ctx.obstaclesByLayer,
+    placed: placedCopy,
+    placedByLayer: placedByLayerCopy,
+    options: opts,
+    iteration: 0,
+    gapsFound: [],
+    gapIndex: 0,
+    done: false,
+    initialGapCount: 0,
+    filledCount: 0,
+    // Four-stage visualization state
+    stage: "scan",
+    currentGap: null,
+    currentSeed: null,
+    expandedRect: null,
+  }
+}
diff --git a/lib/solvers/rectdiff/gapfill/engine/stepGapFill.ts b/lib/solvers/rectdiff/gapfill/engine/stepGapFill.ts
new file mode 100644
index 0000000..ad143c8
--- /dev/null
+++ b/lib/solvers/rectdiff/gapfill/engine/stepGapFill.ts
@@ -0,0 +1,128 @@
+// lib/solvers/rectdiff/gapfill/engine/stepGapFill.ts
+import type { GapFillState } from "../types"
+import { findAllGaps } from "../detection"
+import { tryExpandGap } from "./tryExpandGap"
+import { addPlacement } from "./addPlacement"
+
+/**
+ * Perform one step of gap filling with four-stage visualization.
+ * Stages: scan → select → expand → place
+ * Returns true if still working, false if done.
+ */
+export function stepGapFill(state: GapFillState): boolean {
+  if (state.done) return false
+
+  switch (state.stage) {
+    case "scan": {
+      // Stage 1: Gap detection/scanning
+
+      // Check if we need to find new gaps
+      if (
+        state.gapsFound.length === 0 ||
+        state.gapIndex >= state.gapsFound.length
+      ) {
+        // Check if we've hit max iterations
+        if (state.iteration >= state.options.maxIterations) {
+          state.done = true
+          return false
+        }
+
+        // Find new gaps
+        state.gapsFound = findAllGaps(
+          {
+            scanResolution: state.options.scanResolution,
+            minWidth: state.options.minWidth,
+            minHeight: state.options.minHeight,
+          },
+          {
+            bounds: state.bounds,
+            layerCount: state.layerCount,
+            obstaclesByLayer: state.obstaclesByLayer,
+            placedByLayer: state.placedByLayer,
+          },
+        )
+
+        if (state.iteration === 0) {
+          state.initialGapCount = state.gapsFound.length
+        }
+
+        state.gapIndex = 0
+        state.iteration++
+
+        // If no gaps found, we're done
+        if (state.gapsFound.length === 0) {
+          state.done = true
+          return false
+        }
+      }
+
+      // Move to select stage
+      state.stage = "select"
+      return true
+    }
+
+    case "select": {
+      // Stage 2: Show the gap being targeted
+      if (state.gapIndex >= state.gapsFound.length) {
+        // No more gaps in this iteration, go back to scan
+        state.stage = "scan"
+        return true
+      }
+
+      state.currentGap = state.gapsFound[state.gapIndex]!
+      state.currentSeed = {
+        x: state.currentGap.centerX,
+        y: state.currentGap.centerY,
+      }
+      state.expandedRect = null
+
+      // Move to expand stage
+      state.stage = "expand"
+      return true
+    }
+
+    case "expand": {
+      // Stage 3: Show expansion attempt
+      if (!state.currentGap) {
+        // Shouldn't happen, but handle gracefully
+        state.stage = "select"
+        return true
+      }
+
+      // Try to expand from the current seed
+      const expandedRect = tryExpandGap(state, {
+        gap: state.currentGap,
+        seed: state.currentSeed!,
+      })
+      state.expandedRect = expandedRect
+
+      // Move to place stage
+      state.stage = "place"
+      return true
+    }
+
+    case "place": {
+      // Stage 4: Show the placed result
+      if (state.expandedRect && state.currentGap) {
+        // Actually place the rectangle
+        addPlacement(state, {
+          rect: state.expandedRect,
+          zLayers: state.currentGap.zLayers,
+        })
+        state.filledCount++
+      }
+
+      // Move to next gap and reset to select stage
+      state.gapIndex++
+      state.currentGap = null
+      state.currentSeed = null
+      state.expandedRect = null
+      state.stage = "select"
+      return true
+    }
+
+    default:
+      state.stage = "scan"
+      return true
+  }
+}
diff --git a/lib/solvers/rectdiff/gapfill/engine/tryExpandGap.ts b/lib/solvers/rectdiff/gapfill/engine/tryExpandGap.ts
new file mode 100644
index 0000000..d3ae73c
--- /dev/null
+++ b/lib/solvers/rectdiff/gapfill/engine/tryExpandGap.ts
@@ -0,0 +1,78 @@
+// lib/solvers/rectdiff/gapfill/engine/tryExpandGap.ts
+import type { XYRect } from "../../types"
+import type { GapFillState, GapRegion } from "../types"
+import { expandRectFromSeed } from "../../geometry"
+
+/**
+ * Try to expand a rectangle from a seed point within the gap.
+ * Returns the expanded rectangle or null if expansion fails.
+ */
+export function tryExpandGap(
+  state: GapFillState,
+  {
+    gap,
+    seed,
+  }: {
+    gap: GapRegion
+    seed: { x: number; y: number }
+  },
+): XYRect | null {
+  // Build blockers for the gap's z-layers
+  const blockers: XYRect[] = []
+  for (const z of gap.zLayers) {
+    blockers.push(...(state.obstaclesByLayer[z] ?? []))
+    blockers.push(...(state.placedByLayer[z] ?? []))
+  }
+
+  // Try to expand from the seed point
+  const rect = expandRectFromSeed({
+    startX: seed.x,
+    startY: seed.y,
+    gridSize: Math.min(gap.rect.width, gap.rect.height),
+    bounds: state.bounds,
+    blockers,
+    initialCellRatio: 0,
+    maxAspectRatio: null,
+    minReq: { width: state.options.minWidth, height: state.options.minHeight },
+  })
+
+  if (!rect) {
+    // Try additional seed points within the gap
+    const seeds = [
+      { x: gap.rect.x + state.options.minWidth / 2, y: gap.centerY },
+      {
+        x: gap.rect.x + gap.rect.width - state.options.minWidth / 2,
+        y: gap.centerY,
+      },
+      { x: gap.centerX, y: gap.rect.y + state.options.minHeight / 2 },
+      {
+        x: gap.centerX,
+        y: gap.rect.y + gap.rect.height - state.options.minHeight / 2,
+      },
+    ]
+
+    for (const altSeed of seeds) {
+      const altRect = expandRectFromSeed({
+        startX: altSeed.x,
+        startY: altSeed.y,
+        gridSize: Math.min(gap.rect.width, gap.rect.height),
+        bounds: state.bounds,
+        blockers,
+        initialCellRatio: 0,
+        maxAspectRatio: null,
+        minReq: {
+          width: state.options.minWidth,
+          height: state.options.minHeight,
+        },
+      })
+
+      if (altRect) {
+        return altRect
+      }
+    }
+
+    return null
+  }
+
+  return rect
+}
diff --git a/lib/solvers/rectdiff/gapfill/types.ts b/lib/solvers/rectdiff/gapfill/types.ts
new file mode 100644
index 0000000..59bdbf6
--- /dev/null
+++ b/lib/solvers/rectdiff/gapfill/types.ts
@@ -0,0 +1,60 @@
+// lib/solvers/rectdiff/gapfill/types.ts
+import type { XYRect, Placed3D } from "../types"
+
+export interface GapFillOptions {
+  /** Minimum width for gap-fill rectangles (can be smaller than main solver) */
+  minWidth: number
+  /** Minimum height for gap-fill rectangles */
+  minHeight: number
+  /** Maximum iterations to prevent infinite loops */
+  maxIterations: number
+  /** Target coverage percentage (0-1) to stop early */
+  targetCoverage: number
+  /** Grid resolution for gap detection */
+  scanResolution: number
+}
+
+export interface GapRegion {
+  /** Bounding box of the gap */
+  rect: XYRect
+  /** Z-layers where this gap exists */
+  zLayers: number[]
+  /** Center point for seeding */
+  centerX: number
+  centerY: number
+  /** Approximate area of the gap */
+  area: number
+}
+
+export interface GapFillState {
+  bounds: XYRect
+  layerCount: number
+  obstaclesByLayer: XYRect[][]
+  placed: Placed3D[]
+  placedByLayer: XYRect[][]
+  options: GapFillOptions
+
+  // Progress tracking
+  iteration: number
+  gapsFound: GapRegion[]
+  gapIndex: number
+  done: boolean
+
+  // Stats
+  initialGapCount: number
+  filledCount: number
+
+  // Four-stage visualization state
+  stage: "scan" | "select" | "expand" | "place"
+  currentGap: GapRegion | null
+  currentSeed: { x: number; y: number } | null
+  expandedRect: XYRect | null
+}
+
+/** Context for layer-based operations shared across gap fill functions */
+export interface LayerContext {
+  bounds: XYRect
+  layerCount: number
+  obstaclesByLayer: XYRect[][]
+  placedByLayer: XYRect[][]
+}
diff --git a/lib/solvers/rectdiff/geometry.ts b/lib/solvers/rectdiff/geometry.ts
index 19d0e9f..f2b3b52 100644
--- a/lib/solvers/rectdiff/geometry.ts
+++ b/lib/solvers/rectdiff/geometry.ts
@@ -209,17 +209,28 @@ function maxExpandUp(
   return Math.max(0, e)
 }
 
-/** Grow a rect around (startX,startY), honoring bounds/blockers/aspect/min sizes */
-export function expandRectFromSeed(
-  startX: number,
-  startY: number,
-  gridSize: number,
-  bounds: XYRect,
-  blockers: XYRect[],
-  initialCellRatio: number,
-  maxAspectRatio: number | null | undefined,
-  minReq: { width: number; height: number },
-): XYRect | null {
+/** Grow a rect around a seed point, honoring bounds/blockers/aspect/min sizes. */
+export function expandRectFromSeed(params: {
+  startX: number
+  startY: number
+  gridSize: number
+  bounds: XYRect
+  blockers: XYRect[]
+  initialCellRatio: number
+  maxAspectRatio: number | null | undefined
+  minReq: { width: number; height: number }
+}): XYRect | null {
+  const {
+    startX,
+    startY,
+    gridSize,
+    bounds,
+    blockers,
+    initialCellRatio,
+    maxAspectRatio,
+    minReq,
+  } = params
+
   const minSide = Math.max(1e-9, gridSize * initialCellRatio)
   const initialW = Math.max(minSide, minReq.width)
   const initialH = Math.max(minSide, minReq.height)
@@ -297,6 +308,7 @@ export function expandRectFromSeed(
   return best
 }
 
+/** Find the intersection of two 1D intervals, or null if they don't overlap. */
 export function intersect1D(a0: number, a1: number, b0: number, b1: number) {
   const lo = Math.max(a0, b0)
   const hi = Math.min(a1, b1)
diff --git a/lib/solvers/rectdiff/subsolvers/GapFillSubSolver.ts b/lib/solvers/rectdiff/subsolvers/GapFillSubSolver.ts
new file mode 100644
index 0000000..ecc6cdc
--- /dev/null
+++ b/lib/solvers/rectdiff/subsolvers/GapFillSubSolver.ts
@@ -0,0 +1,253 @@
+// lib/solvers/rectdiff/subsolvers/GapFillSubSolver.ts
+import { BaseSolver } from "@tscircuit/solver-utils"
+import type { GraphicsObject } from "graphics-debug"
+import type { XYRect, Placed3D } from "../types"
+import type {
+  GapFillState,
+  GapFillOptions,
+  LayerContext,
+} from "../gapfill/types"
+import {
+  initGapFillState,
+  stepGapFill,
+  getGapFillProgress,
+} from "../gapfill/engine"
+
+/**
+ * A sub-solver that fills empty spaces (gaps) left by the main grid-based
+ * placement algorithm.
+ *
+ * The preceding grid-based placement is fast but can leave irregular un-placed
+ * areas. This solver maximizes board coverage by finding and filling these
+ * gaps, which is critical for producing a high-quality capacity mesh.
+ *
+ * The core of the algorithm is its gap-detection phase. It works by first
+ * collecting all unique x and y-coordinates from the edges of existing
+ * obstacles and placed rectangles. This set of coordinates is supplemented by a
+ * uniform grid based on the `scanResolution` parameter. Together, these form a
+ * non-uniform grid of cells. The solver then tests the center of each cell for
+ * coverage. Contiguous uncovered cells are merged into larger, maximal
+ * rectangles, which become the candidate gaps to be filled.
+ *
+ * Once a prioritized list of gaps is generated (favoring larger, multi-layer
+ * gaps), the solver iteratively attempts to fill each one by expanding a new
+ * rectangle from a seed point until it collides with an existing boundary.
+ *
+ * The time complexity is dominated by the gap detection, which is approximately
+ * O((N+1/R)^2 * B), where N is the number of objects, R is the scan
+ * resolution, and B is the number of blockers. The algorithm's performance is
+ * therefore highly dependent on the `scanResolution`. It is a heuristic
+ * designed to be "fast enough" by avoiding a brute-force search, instead
+ * relying on this grid-based cell checking to find significant gaps.
+ */
+export class GapFillSubSolver extends BaseSolver {
+  private state: GapFillState
+  private layerCtx: LayerContext
+
+  constructor(params: {
+    placed: Placed3D[]
+    options?: Partial<GapFillOptions>
+    layerCtx: LayerContext
+  }) {
+    super()
+    this.layerCtx = params.layerCtx
+    this.state = initGapFillState(
+      {
+        placed: params.placed,
+        options: params.options,
+      },
+      params.layerCtx,
+    )
+  }
+
+  /**
+   * Execute one step of the gap fill algorithm.
+   * Each gap goes through four stages: scan for gaps, select a target gap,
+   * expand a rectangle from seed point, then place the final result.
+   */
+  override _step() {
+    const stillWorking = stepGapFill(this.state)
+    if (!stillWorking) {
+      this.solved = true
+    }
+  }
+
+  /**
+   * Calculate progress as a value between 0 and 1.
+   * Accounts for iterations, gaps processed, and current stage within each gap.
+   */
+  computeProgress(): number {
+    return getGapFillProgress(this.state)
+  }
+
+  /**
+   * Get all placed rectangles including original ones plus newly created gap-fill rectangles.
+   */
+  getPlaced(): Placed3D[] {
+    return this.state.placed
+  }
+
+  /**
+   * Get placed rectangles organized by Z-layer for efficient layer-based operations.
+   */
+  getPlacedByLayer(): XYRect[][] {
+    return this.state.placedByLayer
+  }
+
+  override getOutput() {
+    return {
+      placed: this.state.placed,
+      placedByLayer: this.state.placedByLayer,
+      filledCount: this.state.filledCount,
+    }
+  }
+
+  /** Zen visualization: show four-stage gap filling process. */
+  override visualize(): GraphicsObject {
+    const rects: NonNullable<GraphicsObject["rects"]> = []
+    const points: NonNullable<GraphicsObject["points"]> = []
+
+    // Board bounds (subtle)
+    rects.push({
+      center: {
+        x: this.layerCtx.bounds.x + this.layerCtx.bounds.width / 2,
+        y: this.layerCtx.bounds.y + this.layerCtx.bounds.height / 2,
+      },
+      width: this.layerCtx.bounds.width,
+      height: this.layerCtx.bounds.height,
+      fill: "none",
+      stroke: "#e5e7eb",
+      label: "",
+    })
+
+    switch (this.state.stage) {
+      case "scan": {
+        // Stage 1: Show scanning/detection phase with light blue overlay
+        rects.push({
+          center: {
+            x: this.layerCtx.bounds.x + this.layerCtx.bounds.width / 2,
+            y: this.layerCtx.bounds.y + this.layerCtx.bounds.height / 2,
+          },
+          width: this.layerCtx.bounds.width,
+          height: this.layerCtx.bounds.height,
+          fill: "#dbeafe",
+          stroke: "#3b82f6",
+          label: "scanning",
+        })
+        break
+      }
+
+      case "select": {
+        // Stage 2: Show the gap being targeted (red outline)
+        if (this.state.currentGap) {
+          rects.push({
+            center: {
+              x:
+                this.state.currentGap.rect.x +
+                this.state.currentGap.rect.width / 2,
+              y:
+                this.state.currentGap.rect.y +
+                this.state.currentGap.rect.height / 2,
+            },
+            width: this.state.currentGap.rect.width,
+            height: this.state.currentGap.rect.height,
+            fill: "#fecaca",
+            stroke: "#ef4444",
+            label: "target gap",
+          })
+
+          // Show the seed point
+          if (this.state.currentSeed) {
+            points.push({
+              x: this.state.currentSeed.x,
+              y: this.state.currentSeed.y,
+              color: "#dc2626",
+              label: "seed",
+            })
+          }
+        }
+        break
+      }
+
+      case "expand": {
+        // Stage 3: Show expansion attempt (yellow growing rectangle + seed)
+        if (this.state.currentGap) {
+          // Show gap outline (faded)
+          rects.push({
+            center: {
+              x:
+                this.state.currentGap.rect.x +
+                this.state.currentGap.rect.width / 2,
+              y:
+                this.state.currentGap.rect.y +
+                this.state.currentGap.rect.height / 2,
+            },
+            width: this.state.currentGap.rect.width,
+            height: this.state.currentGap.rect.height,
+            fill: "none",
+            stroke: "#f87171",
+            label: "",
+          })
+        }
+
+        if (this.state.currentSeed) {
+          // Show seed point
+          points.push({
+            x: this.state.currentSeed.x,
+            y: this.state.currentSeed.y,
+            color: "#f59e0b",
+            label: "expanding",
+          })
+        }
+
+        if (this.state.expandedRect) {
+          // Show expanded rectangle
+          rects.push({
+            center: {
+              x: this.state.expandedRect.x + this.state.expandedRect.width / 2,
+              y: this.state.expandedRect.y + this.state.expandedRect.height / 2,
+            },
+            width: this.state.expandedRect.width,
+            height: this.state.expandedRect.height,
+            fill: "#fef3c7",
+            stroke: "#f59e0b",
+            label: "expanding",
+          })
+        }
+        break
+      }
+
+      case "place": {
+        // Stage 4: Show final placed rectangle (green)
+        if (this.state.expandedRect) {
+          rects.push({
+            center: {
+              x: this.state.expandedRect.x + this.state.expandedRect.width / 2,
+              y: this.state.expandedRect.y + this.state.expandedRect.height / 2,
+            },
+            width: this.state.expandedRect.width,
+            height: this.state.expandedRect.height,
+            fill: "#bbf7d0",
+            stroke: "#22c55e",
+            label: "placed",
+          })
+        }
+        break
+      }
+    }
+
+    const stageNames = {
+      scan: "scanning",
+      select: "selecting",
+      expand: "expanding",
+      place: "placing",
+    }
+
+    return {
+      title: `GapFill (${stageNames[this.state.stage]}): ${this.state.filledCount} filled`,
+      coordinateSystem: "cartesian",
+      rects,
+      points,
+    }
+  }
+}
diff --git a/lib/solvers/rectdiff/types.ts b/lib/solvers/rectdiff/types.ts
index d247c05..2ae50bd 100644
--- a/lib/solvers/rectdiff/types.ts
+++ b/lib/solvers/rectdiff/types.ts
@@ -34,7 +34,7 @@ export type Candidate3D = {
 }
 export type Placed3D = { rect: XYRect; zLayers: number[] }
 
-export type Phase = "GRID" | "EXPANSION" | "DONE"
+export type Phase = "GRID" | "EXPANSION" | "GAP_FILL" | "DONE"
 
 export type RectDiffState = {
   // static
diff --git a/tests/obstacle-extra-layers.test.ts b/tests/obstacle-extra-layers.test.ts
index e4e8602..fff701f 100644
--- a/tests/obstacle-extra-layers.test.ts
+++ b/tests/obstacle-extra-layers.test.ts
@@ -29,7 +29,7 @@ test("RectDiffSolver clamps extra layer names to available z indices", () => {
     ],
   }
 
-  const solver = new RectDiffSolver({ simpleRouteJson: srj, mode: "grid" })
+  const solver = new RectDiffSolver({ simpleRouteJson: srj })
   solver.setup()
 
   expect(srj.obstacles[0]?.zLayers).toEqual([1])
diff --git a/tests/obstacle-zlayers.test.ts b/tests/obstacle-zlayers.test.ts
index 5909da4..d4fa4c0 100644
--- a/tests/obstacle-zlayers.test.ts
+++ b/tests/obstacle-zlayers.test.ts
@@ -29,7 +29,7 @@ test("RectDiffSolver maps obstacle layers to numeric zLayers", () => {
     ],
   }
 
-  const solver = new RectDiffSolver({ simpleRouteJson: srj, mode: "grid" })
+  const solver = new RectDiffSolver({ simpleRouteJson: srj })
   solver.setup()
 
   expect(srj.obstacles[0]?.zLayers).toEqual([0])
diff --git a/tests/rect-diff-solver.test.ts b/tests/rect-diff-solver.test.ts
index 570ebf9..25bef1a 100644
--- a/tests/rect-diff-solver.test.ts
+++ b/tests/rect-diff-solver.test.ts
@@ -27,7 +27,6 @@ test("RectDiffSolver creates mesh nodes with grid-based approach", () => {
 
   const solver = new RectDiffSolver({
     simpleRouteJson,
-    mode: "grid",
   })
 
   solver.solve()
@@ -62,7 +61,6 @@ test("RectDiffSolver handles multi-layer spans", () => {
 
   const solver = new RectDiffSolver({
     simpleRouteJson,
-    mode: "grid",
     gridOptions: {
       minSingle: { width: 0.4, height: 0.4 },
       minMulti: { width: 1.0, height: 1.0, minLayers: 2 },
@@ -105,7 +103,6 @@ test("RectDiffSolver respects single-layer minimums", () => {
 
   const solver = new RectDiffSolver({
     simpleRouteJson,
-    mode: "grid",
     gridOptions: {
       minSingle: { width: minWidth, height: minHeight },
       minMulti: { width: 1.0, height: 1.0, minLayers: 2 },
@@ -131,7 +128,7 @@ test("disruptive placement resizes single-layer nodes", () => {
     layerCount: 3,
     minTraceWidth: 0.2,
   }
-  const solver = new RectDiffSolver({ simpleRouteJson: srj, mode: "grid" })
+  const solver = new RectDiffSolver({ simpleRouteJson: srj })
   solver.setup()
 
   // Manually seed a soft, single-layer node occupying center (simulate early placement)
diff --git a/utils/README.md b/utils/README.md
new file mode 100644
index 0000000..d7a66ac
--- /dev/null
+++ b/utils/README.md
@@ -0,0 +1,21 @@
+# Project Utilities
+
+This directory contains global utility functions that can be reused across different parts of the application.
+
+## Geometry
+
+### `rectsOverlap(a: XYRect, b: XYRect): boolean`
+
+- **Description:** Checks if two rectangles overlap.
+- **Parameters:**
+  - `a`: The first rectangle (`XYRect`).
+  - `b`: The second rectangle (`XYRect`).
+- **Returns:** `true` if the rectangles overlap, `false` otherwise.
+
+### `rectsEqual(a: XYRect, b: XYRect): boolean`
+
+- **Description:** Checks if two rectangles are equal within a small tolerance (EPS).
+- **Parameters:**
+  - `a`: The first rectangle (`XYRect`).
+  - `b`: The second rectangle (`XYRect`).
+- **Returns:** `true` if the rectangles are equal, `false` otherwise.
diff --git a/utils/rectsEqual.ts b/utils/rectsEqual.ts
new file mode 100644
index 0000000..93f43f2
--- /dev/null
+++ b/utils/rectsEqual.ts
@@ -0,0 +1,18 @@
+// utils/rectsEqual.ts
+import type { XYRect } from "../lib/solvers/rectdiff/types"
+import { EPS } from "../lib/solvers/rectdiff/geometry"
+
+/**
+ * Checks if two rectangles are equal within a small tolerance (EPS).
+ * @param a The first rectangle.
+ * @param b The second rectangle.
+ * @returns True if the rectangles are equal, false otherwise.
+ */
+export function rectsEqual(a: XYRect, b: XYRect): boolean {
+  return (
+    Math.abs(a.x - b.x) < EPS &&
+    Math.abs(a.y - b.y) < EPS &&
+    Math.abs(a.width - b.width) < EPS &&
+    Math.abs(a.height - b.height) < EPS
+  )
+}
diff --git a/utils/rectsOverlap.ts b/utils/rectsOverlap.ts
new file mode 100644
index 0000000..f986bdf
--- /dev/null
+++ b/utils/rectsOverlap.ts
@@ -0,0 +1,18 @@
+// utils/rectsOverlap.ts
+import type { XYRect } from "../lib/solvers/rectdiff/types"
+import { EPS } from "../lib/solvers/rectdiff/geometry"
+
+/**
+ * Checks if two rectangles overlap.
+ * @param a The first rectangle.
+ * @param b The second rectangle.
+ * @returns True if the rectangles overlap, false otherwise.
+ */
+export function rectsOverlap(a: XYRect, b: XYRect): boolean {
+  return !(
+    a.x + a.width <= b.x + EPS ||
+    b.x + b.width <= a.x + EPS ||
+    a.y + a.height <= b.y + EPS ||
+    b.y + b.height <= a.y + EPS
+  )
+}