Finds the pad level, and the tilt if allowed, at which cut fills the fill, and plans which cut goes to which fill and how far.
How the precision is kept
64-bit arithmetic on a grid of ground samples. Volumes are sums of depth times cell area, so they are known to the cell size, about 20,000 cells over the area's bounding box by default. The level comes from 80 bisections and the tilt from a stepped grid search, which finds a low point of the haul moment but does not prove it is the lowest.
Calculations under this policy
9 calculations in 1 source file. The code is read from the released app when this page loads. It is the calculation people run, not a copy kept for this page.
A grid of ground elevations inside the polygon, each standing for one cell of ground.
cell = √(bounding box area ÷ 20 000) unless given
Cell centers are a running sum of the cell size from the box corner. The inside test is the even-odd crossing rule in 64-bit; a center exactly on an edge can fall either way. The polygon's own area is not used: the area is the count of cells times the cell area.
earthwork.js
lines 23 to 37
CivilCraft 0.8.1
export function cellsInside(ground, polygon, { cell } = {}) {
let x0 = Infinity, y0 = Infinity, x1 = -Infinity, y1 = -Infinity;
for (const p of polygon) { x0 = Math.min(x0, p.x); y0 = Math.min(y0, p.y); x1 = Math.max(x1, p.x); y1 = Math.max(y1, p.y); }
if (!(x1 > x0 && y1 > y0)) return [];
const c = cell > 0 ? cell : Math.sqrt(((x1 - x0) * (y1 - y0)) / 20000);
const inside = (x, y) => { let k = false; for (let i = 0, j = polygon.length - 1; i < polygon.length; j = i++) { const p = polygon[i], q = polygon[j]; if ((p.y > y) !== (q.y > y) && x < ((q.x - p.x) * (y - p.y)) / (q.y - p.y) + p.x) k = !k; } return k; };
const out = [];
for (let y = y0 + c / 2; y < y1; y += c) for (let x = x0 + c / 2; x < x1; x += c) {
if (!inside(x, y)) continue;
const z = ground.elevationAt(x, y);
if (z !== null && z !== undefined && Number.isFinite(z)) out.push({ x, y, zg: z });
}
out.cell = c;
return out;
}
Cut and fill volumes
64-bit
Cut, fill, usable cut and the net at a design level.
fill = area × Σ max(design − ground, 0); cut = area × Σ max(ground − design, 0); net = cut × factor − fill
Depths are summed first and multiplied by the cell area once. Each sum is in cell order with all terms of one sign, so nothing cancels inside a sum; the net is a difference of the two totals.
earthwork.js
lines 39 to 44
CivilCraft 0.8.1
export function volumes(cells, designOf, z0, { area = 1, factor = 1 } = {}) {
let cut = 0, fill = 0;
for (const p of cells) { const dz = designOf(p.x, p.y, z0) - p.zg; if (dz > 0) fill += dz; else cut -= dz; }
cut *= area; fill *= area;
return { cut, fill, usable: cut * factor, net: cut * factor - fill };
}
Balanced level
64-bit
The design level at which usable cut equals fill.
Bisection on the sign of the net, always 80 halvings, starting 1 unit below and above the range of ground minus design. Eighty halvings take the bracket below 64-bit resolution, so the level is as fine as the sampled volumes allow. The net at the returned level is reported, not assumed zero.
earthwork.js
lines 46 to 58
CivilCraft 0.8.1
/** net = usable cut − fill falls as z0 rises (less cut, more fill): bisection on the sign. */
export function balanceLevel(cells, designOf, { area = 1, factor = 1, lo = null, hi = null } = {}) {
if (!cells.length) return { z0: 0, cut: 0, fill: 0, net: 0, ok: false };
let a = lo, b = hi;
if (a === null || b === null) {
let mn = Infinity, mx = -Infinity;
for (const p of cells) { const d = designOf(p.x, p.y, 0); mn = Math.min(mn, p.zg - d); mx = Math.max(mx, p.zg - d); }
a = mn - 1; b = mx + 1; // at a every cell is fill, at b every cell is cut
}
for (let k = 0; k < 80; k++) { const m = (a + b) / 2; if (volumes(cells, designOf, m, { area, factor }).net > 0) a = m; else b = m; }
const z0 = (a + b) / 2, v = volumes(cells, designOf, z0, { area, factor });
return { z0, cut: v.cut, fill: v.fill, net: v.net, ok: true };
}
Design plane
64-bit
The design elevation at a point for a level and two grades.
z = z₀ + gx·(x − cx) + gy·(y − cy)
Coordinates are reduced by the plane's center before they are multiplied by the grades, so the products are of short distances, not full-size coordinates.
The volume moved times the distance between the centroid of cut and the centroid of fill.
moment = min(cut × factor, fill) × |centroid(cut) − centroid(fill)|
Centroids are sums of depth × coordinate divided by the sum of depth, using full-size coordinates, so each centroid carries the rounding of a sum of large products. The distance between the two is then a difference of nearly equal numbers when cut and fill are close together.
The level and tilt, within the grade limit, that balance with the least haul moment.
A 9 by 9 grid of grades is tried around the current best, each one balanced by its own bisection; the step halves each round, for at most 6 rounds. This is a search on samples: it returns the best tilt tried, not a proven minimum. When a least fall is asked for, 16 directions at that grade are compared instead.
earthwork.js
lines 70 to 96
CivilCraft 0.8.1
export function bestPlane(cells, { area = 1, factor = 1, maxGrade = 0.05, steps = 9, cx = null, cy = null, minGrade = 0 } = {}) {
if (!cells.length) return null;
const mx = cx ?? cells.reduce((s, p) => s + p.x, 0) / cells.length, my = cy ?? cells.reduce((s, p) => s + p.y, 0) / cells.length;
const score = (gx, gy) => {
if (Math.hypot(gx, gy) > maxGrade + 1e-12) return null;
const d = planeOf({ cx: mx, cy: my, gx, gy }), b = balanceLevel(cells, d, { area, factor });
return { gx, gy, ...b, ...moment(cells, d, b.z0, area, factor), haulMoment: moment(cells, d, b.z0, area, factor).m };
};
let best = score(0, 0), step = maxGrade / ((steps - 1) / 2 || 1);
for (let round = 0; round < 6 && maxGrade > 0; round++) {
let moved = false;
const base = best;
for (let i = -Math.floor(steps / 2); i <= Math.floor(steps / 2); i++) for (let j = -Math.floor(steps / 2); j <= Math.floor(steps / 2); j++) {
const s = score(base.gx + i * step, base.gy + j * step);
if (s && s.haulMoment < best.haulMoment - 1e-12) { best = s; moved = true; }
}
step = moved ? step : step / 2;
if (!moved && step < maxGrade * 1e-4) break;
if (moved) step /= 2;
}
if (minGrade > 0 && Math.hypot(best.gx, best.gy) < minGrade) { // drainage wants a fall: take the best of the allowed tilts at the least grade
let b2 = null;
for (let k = 0; k < 16; k++) { const a = (k * Math.PI) / 8, s = score(Math.cos(a) * minGrade, Math.sin(a) * minGrade); if (s && (!b2 || s.haulMoment < b2.haulMoment)) b2 = s; }
if (b2) best = b2;
}
return { z0: best.z0, gx: best.gx, gy: best.gy, cx: mx, cy: my, cut: best.cut, fill: best.fill, net: best.net, haulMoment: best.haulMoment, grade: Math.hypot(best.gx, best.gy) };
}
Blocks of cut and fill
64-bit
Cells grouped into blocks, each with its volume and its volume-weighted center.
center = Σ (p × v) ÷ Σ v
Block centers are weighted sums of full-size coordinates divided by the block volume. The cell's block index is a floor of a scaled coordinate, clamped to the grid.
earthwork.js
lines 101 to 112
CivilCraft 0.8.1
function blocksOf(list, n) {
if (list.length <= n) return list.map((p) => ({ x: p.x, y: p.y, v: p.v }));
let x0 = Infinity, y0 = Infinity, x1 = -Infinity, y1 = -Infinity;
for (const p of list) { x0 = Math.min(x0, p.x); y0 = Math.min(y0, p.y); x1 = Math.max(x1, p.x); y1 = Math.max(y1, p.y); }
const k = Math.max(1, Math.ceil(Math.sqrt(n * 1.3))), sx = (x1 - x0 || 1) / k, sy = (y1 - y0 || 1) / k, map = new Map();
for (const p of list) {
const i = Math.min(k - 1, Math.floor((p.x - x0) / sx)), j = Math.min(k - 1, Math.floor((p.y - y0) / sy)), key = j * k + i;
let b = map.get(key); if (!b) map.set(key, (b = { x: 0, y: 0, v: 0 }));
b.x += p.x * p.v; b.y += p.y * p.v; b.v += p.v;
}
return [...map.values()].filter((b) => b.v > 0).map((b) => ({ x: b.x / b.v, y: b.y / b.v, v: b.v }));
}
Minimum-cost transport
64-bit
How much goes from each block of cut to each block of fill at the least total volume times distance.
Successive shortest paths with real-valued capacities. A capacity or a cost improvement below 1e-12 is ignored and flows below 1e-9 are dropped; both are absolute, in volume and distance units. The loop stops after 50,000 augmentations whether or not it has finished.
earthwork.js
lines 114 to 141
CivilCraft 0.8.1
/** Transportation problem by successive shortest paths (Bellman-Ford on the residual network). supply[i], demand[j], cost[i][j]. */
export function transport(supply, demand, cost) {
const n = supply.length, m = demand.length, N = n + m + 2, S = n + m, T = n + m + 1;
const to = [], cap = [], cst = [], adj = Array.from({ length: N }, () => []);
const add = (u, v, c, w) => { adj[u].push(to.length); to.push(v); cap.push(c); cst.push(w); adj[v].push(to.length); to.push(u); cap.push(0); cst.push(-w); };
for (let i = 0; i < n; i++) add(S, i, supply[i], 0);
for (let j = 0; j < m; j++) add(n + j, T, demand[j], 0);
const first = to.length;
for (let i = 0; i < n; i++) for (let j = 0; j < m; j++) add(i, n + j, Infinity, cost[i][j]);
let flow = 0, guard = 0;
for (;;) {
const dist = new Float64Array(N).fill(Infinity), prev = new Int32Array(N).fill(-1), inq = new Uint8Array(N), q = [S];
dist[S] = 0;
while (q.length) {
const u = q.shift(); inq[u] = 0;
for (const e of adj[u]) if (cap[e] > 1e-12 && dist[u] + cst[e] < dist[to[e]] - 1e-12) { dist[to[e]] = dist[u] + cst[e]; prev[to[e]] = e; if (!inq[to[e]]) { inq[to[e]] = 1; q.push(to[e]); } }
}
if (!Number.isFinite(dist[T])) break;
let f = Infinity;
for (let v = T; v !== S; v = to[prev[v] ^ 1]) f = Math.min(f, cap[prev[v]]);
for (let v = T; v !== S; v = to[prev[v] ^ 1]) { cap[prev[v]] -= f; cap[prev[v] ^ 1] += f; }
flow += f;
if (++guard > 50000) break;
}
const flows = [];
for (let i = 0; i < n; i++) for (let j = 0; j < m; j++) { const e = first + 2 * (i * m + j), f = cap[e ^ 1]; if (f > 1e-9) flows.push({ i, j, vol: f }); }
return { flows, flow };
}
Haul plan totals
64-bit
Volume moved on site, total haul, average haul distance, and what is left to take away or bring in.
haul = Σ volume × distance; average = haul ÷ moved
Distances are straight lines between block centers by hypot. Totals are running sums over the flows. Spare and short volumes below 1e-9 are treated as none.
earthwork.js
lines 143 to 170
CivilCraft 0.8.1
export function haulPlan(cells, designOf, z0, { area = 1, factor = 1, blocks = 60, offsite = null } = {}) {
const cuts = [], fills = [];
for (const p of cells) { const dz = designOf(p.x, p.y, z0) - p.zg; if (dz < 0) cuts.push({ x: p.x, y: p.y, v: -dz * area * factor }); else if (dz > 0) fills.push({ x: p.x, y: p.y, v: dz * area }); }
const C = blocksOf(cuts, blocks), F = blocksOf(fills, blocks);
const cutV = C.reduce((s, b) => s + b.v, 0), fillV = F.reduce((s, b) => s + b.v, 0);
let far = offsite;
if (far === null) { let x0 = Infinity, y0 = Infinity, x1 = -Infinity, y1 = -Infinity; for (const p of cells) { x0 = Math.min(x0, p.x); y0 = Math.min(y0, p.y); x1 = Math.max(x1, p.x); y1 = Math.max(y1, p.y); } far = 2 * Math.hypot(x1 - x0, y1 - y0); }
// dummy sink takes spare cut, dummy source brings missing fill; both at distance `far` so on-site moves win
const supply = C.map((b) => b.v), demand = F.map((b) => b.v);
const spareTotal = Math.max(0, cutV - fillV), shortTotal = Math.max(0, fillV - cutV);
if (spareTotal > 1e-9) demand.push(spareTotal);
if (shortTotal > 1e-9) supply.push(shortTotal);
const cost = supply.map((_, i) => demand.map((_, j) => {
const a = i < C.length ? C[i] : null, b = j < F.length ? F[j] : null;
if (a && b) return Math.hypot(a.x - b.x, a.y - b.y);
return a || b ? far : 1e9; // dummy to dummy never
}));
const r = transport(supply, demand, cost);
const flows = [];
let moved = 0, haul = 0, spare = 0, short = 0;
for (const f of r.flows) {
const a = f.i < C.length ? C[f.i] : null, b = f.j < F.length ? F[f.j] : null;
if (a && b) { const d = Math.hypot(a.x - b.x, a.y - b.y); flows.push({ from: { x: a.x, y: a.y }, to: { x: b.x, y: b.y }, vol: f.vol, dist: d }); moved += f.vol; haul += f.vol * d; }
else if (a) spare += f.vol; else if (b) short += f.vol;
}
flows.sort((p, q) => q.vol - p.vol);
return { flows, moved, haul, avg: moved > 0 ? haul / moved : 0, spare, short, totals: { cut: cutV, fill: fillV }, blocks: { cut: C.length, fill: F.length } };
}
In detail
Method
Cell-centered prisms for cut and fill. Bisection on the sign of usable cut minus fill. Tilt by repeated grid search on the two grades with halving steps. Haul by a minimum-cost transport between blocks of cut and fill, solved by successive shortest paths (Bellman-Ford on the residual network), as the file header states.
Number types
64-bit numbers; the shortest-path distances are a Float64Array and the path links an Int32Array of indexes.
Tolerances
Level: 80 bisections from 1 unit outside the range of ground minus design. Tilt: grade limit + 1e-12; a better tilt must lower the haul moment by more than 1e-12; search stops when the step is below a ten-thousandth of the grade limit or after 6 rounds. Transport: capacities and cost improvements below 1e-12 ignored, flows below 1e-9 dropped, at most 50,000 augmentations.
Limits
Cells are grouped into at most 60 blocks of cut and 60 of fill by default, and haul distances are straight lines between block centroids. Leftover cut or fill goes to a dummy site at twice the area's diagonal.
Approach: 64-bit, stated limit
The arithmetic is 64-bit, but the method or the data sets the accuracy: a solver tolerance, a mesh, a step size, a datum shift, a table. The limit is stated where the result is shown.