Resources Numerical accuracy Calculation policy

Cross sections and end-area volumes along an alignment

Cuts the ground and the design along sample lines square to an alignment, works out cut and fill areas at each section and the volumes between sections.

How the precision is kept

A section has a vertex wherever the sample line crosses a triangle edge, so it follows the triangulated surface without sampling error. Areas are trapezoids between the two section lines, split where they cross. All of it is 64-bit. Volumes use the average end area method, which is exact only when the area varies linearly between sections. Sample stations on the step are rounded to 6 decimals, and the CSV export rounds areas to 2 decimals and volumes to 1.

Open item. Scale the 1e-15 parallel test in segCross to the lengths of the two segments.

Calculations under this policy

10 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.

sections.js CivilCraft 0.8.1 64-bit 10 calculations Reviewed

Tests bin/test-civilcraft-sections.mjs, bin/test-civilcraft-masshaul.mjs

Sample stations

64-bit

The stations where sections are cut: every multiple of the step, both ends, the curve ends and any extra stations.

sta = round₆(k × step) for whole k; neighbors closer than 1e-6 merged

Each station is k × step, one multiplication, then rounded to 6 decimals with toFixed(6) before it is kept. The first k uses ceil(lo ÷ step − 1e-9), so a start that is a multiple of the step within rounding is not skipped.

  • sections.js
  • lines 48 to 58
  • CivilCraft 0.8.1
export function sampleStations(al, { step, extra = [], from, to } = {}) {
  if (!al?.ok) return [];
  const lo = Math.max(al.start, Number.isFinite(from) ? from : al.start), hi = Math.min(al.end, Number.isFinite(to) ? to : al.end);
  if (!(hi >= lo)) return [];
  const st = new Set([lo, hi]);
  if (step > 0) for (let k = Math.ceil(lo / step - 1e-9); k * step < hi - 1e-9 && st.size < MAX_SECTIONS; k++) st.add(Number((k * step).toFixed(6)));
  for (const k of al.keyStations?.() || []) if (k.sta > lo && k.sta < hi) st.add(k.sta);
  for (const s of extra) if (s >= lo && s <= hi) st.add(s);
  const out = [...st].sort((a, b) => a - b).filter((s, i, a) => !i || s - a[i - 1] > 1e-6);
  return out.slice(0, MAX_SECTIONS);
}

Sample line at a station

64-bit

The two ends of the line square to the alignment at a station.

right = (sin dir, −cos dir); a = c − right × left; b = c + right × right offset

64-bit sine and cosine of the alignment direction, each end one multiply and add from the center point.

  • sections.js
  • lines 60 to 64
  • CivilCraft 0.8.1
export function sampleLine(al, sta, left, right) {
  const p = al.pointAt(sta);
  const rx = Math.sin(p.dir), ry = -Math.cos(p.dir);   // right of the forward direction (y up)
  return { sta, c: { x: p.x, y: p.y }, a: { x: p.x - rx * left, y: p.y - ry * left }, b: { x: p.x + rx * right, y: p.y + ry * right }, dir: p.dir, left, right };
}

Surface along a sample line

64-bit

The ground or design line along a sample line, with a point at every triangle edge it crosses.

z at a + (b − a)·t for every edge-crossing parameter t; o = offA + (offB − offA)·t

Coordinates are shifted by the surface origin before the intersections, which keeps the subtractions small. Crossing parameters closer than 1e-12 are merged; elevations come from the surface's own lookup with tolerance 1e-7. Each offset is one interpolation from the parameter, not an accumulated distance.

  • sections.js
  • lines 66 to 110
  • CivilCraft 0.8.1
/** The surface along a→b, exact: params where the segment crosses triangle edges, z by the surface itself. */
export function surfaceSection(surface, a, b, offA, offB) {
  if (!surface) return [];
  const { X, Y, tris, ox, oy } = surface;
  const ax = a.x - ox, ay = a.y - oy, bx = b.x - ox, by = b.y - oy, dx = bx - ax, dy = by - ay;
  const L = Math.hypot(dx, dy);
  if (!(L > 0)) return [];
  const ts = [0, 1];
  const g = surface.grid();
  const clampI = (v, n) => Math.max(0, Math.min(n - 1, Math.floor(v)));
  const x0 = Math.min(ax, bx), x1 = Math.max(ax, bx), y0 = Math.min(ay, by), y1 = Math.max(ay, by);
  if (!(x1 < g.b.minX || x0 > g.b.maxX || y1 < g.b.minY || y0 > g.b.maxY)) {
    const i0 = clampI((x0 - g.b.minX) / g.cw, g.nx), i1 = clampI((x1 - g.b.minX) / g.cw, g.nx), j0 = clampI((y0 - g.b.minY) / g.ch, g.ny), j1 = clampI((y1 - g.b.minY) / g.ch, g.ny);
    const seen = new Set();
    // cells along the segment: walk the cell rows/columns the segment's box covers but test each cell against the line
    for (let j = j0; j <= j1; j++) for (let i = i0; i <= i1; i++) {
      const cx0 = g.b.minX + i * g.cw, cy0 = g.b.minY + j * g.ch;
      if (!segHitsBox(ax, ay, bx, by, cx0, cy0, cx0 + g.cw, cy0 + g.ch)) continue;
      const k = j * g.nx + i;
      for (let s = g.start[k]; s < g.start[k + 1]; s++) {
        const tri = g.items[s];
        if (seen.has(tri)) continue;
        seen.add(tri);
        for (let e = 0; e < 3; e++) {
          const p = tris[3 * tri + e], q = tris[3 * tri + ((e + 1) % 3)];
          const tt = segCross(ax, ay, dx, dy, X[p], Y[p], X[q] - X[p], Y[q] - Y[p]);
          if (tt !== null) ts.push(tt);
        }
      }
    }
  }
  ts.sort((u, v) => u - v);
  const runs = [];
  let run = null, last = -1;
  for (const tt of ts) {
    if (tt - last < 1e-12) continue;
    last = tt;
    const z = surface.elevationAt(a.x + (b.x - a.x) * tt, a.y + (b.y - a.y) * tt, 1e-7);
    const o = offA + (offB - offA) * tt;
    if (z === null) { run = null; continue; }
    if (!run) runs.push((run = []));
    run.push({ o, z });
  }
  return runs.filter((r) => r.length >= 2);
}

Segment crossing parameter

64-bit

Where one segment crosses another, as a fraction along the first.

den = r × s; t = ((q − p) × s) ÷ den; u = ((q − p) × r) ÷ den

Plain 64-bit cross products, not an exact predicate. Segments are treated as parallel when |den| < 1e-15, an absolute value that is not scaled to the segment lengths. Parameters within 1e-12 outside [0, 1] are accepted and clamped.

  • sections.js
  • lines 111 to 116
  • CivilCraft 0.8.1
function segCross(px, py, rx, ry, qx, qy, sx, sy) {   // segment p + t·r vs q + u·s → t in [0,1] or null
  const den = rx * sy - ry * sx;
  if (Math.abs(den) < 1e-15) return null;
  const t0 = ((qx - px) * sy - (qy - py) * sx) / den, u = ((qx - px) * ry - (qy - py) * rx) / den;
  return t0 >= -1e-12 && t0 <= 1 + 1e-12 && u >= -1e-12 && u <= 1 + 1e-12 ? Math.max(0, Math.min(1, t0)) : null;
}

Elevation of a section line at an offset

64-bit

The elevation of a section at any offset, between its points.

z = z₀ + (z₁ − z₀)·(o − o₀) ÷ (o₁ − o₀)

Linear interpolation in 64-bit; a zero-width interval returns the second point's elevation instead of dividing. Offsets are matched with an absolute 1e-9.

  • sections.js
  • lines 128 to 137
  • CivilCraft 0.8.1
/** z of a piecewise-linear section (runs) at offset o, or null. */
export function zAt(runs, o) {
  for (const r of runs) {
    if (o < r[0].o - 1e-9 || o > r[r.length - 1].o + 1e-9) continue;
    for (let i = 1; i < r.length; i++) {
      if (o <= r[i].o + 1e-9) { const p = r[i - 1], q = r[i], d = q.o - p.o; return d > 0 ? p.z + ((q.z - p.z) * (o - p.o)) / d : q.z; }
    }
  }
  return null;
}

Cut and fill area of a section

64-bit

The area of cut and of fill between the ground line and the design line of one section.

per interval: area = (e₀ + e₁) ÷ 2 × w; where the lines cross: o× = o₀ + w·e₀ ÷ (e₀ − e₁), two triangles

The interval ends are the union of both lines' offsets, so each piece is a true trapezoid or a pair of triangles. The crossing division e₀ ÷ (e₀ − e₁) is only reached when the two differences have opposite signs, so the divisor is not near zero by cancellation. Areas are added in offset order.

  • sections.js
  • lines 139 to 167
  • CivilCraft 0.8.1
export function cutFill(ground, design) {
  const os = new Set();
  for (const r of [...ground, ...design]) for (const p of r) os.add(p.o);
  const xs = [...os].sort((a, b) => a - b);
  let cut = 0, fill = 0;
  const parts = [];
  const push = (o0, o1, kind, area) => {
    if (!(area > 0) || o1 - o0 <= 0) return;
    const l = parts[parts.length - 1];
    if (l && l.kind === kind && Math.abs(l.o1 - o0) < 1e-9) l.o1 = o1; else parts.push({ o0, o1, kind });
    if (kind === 'cut') cut += area; else fill += area;
  };
  for (let i = 1; i < xs.length; i++) {
    const o0 = xs[i - 1], o1 = xs[i], om = (o0 + o1) / 2;
    if (zAt(ground, om) === null || zAt(design, om) === null) continue;
    const g0 = zAt(ground, o0), g1 = zAt(ground, o1), d0 = zAt(design, o0), d1 = zAt(design, o1);
    if ([g0, g1, d0, d1].some((v) => v === null)) continue;
    const e0 = d0 - g0, e1 = d1 - g1, w = o1 - o0;
    if (e0 >= 0 && e1 >= 0) push(o0, o1, 'fill', ((e0 + e1) / 2) * w);
    else if (e0 <= 0 && e1 <= 0) push(o0, o1, 'cut', ((-e0 - e1) / 2) * w);
    else {   // the lines cross inside the interval
      const om2 = o0 + (w * e0) / (e0 - e1);
      const a0 = (Math.abs(e0) / 2) * (om2 - o0), a1 = (Math.abs(e1) / 2) * (o1 - om2);
      push(o0, om2, e0 > 0 ? 'fill' : 'cut', a0);
      push(om2, o1, e1 > 0 ? 'fill' : 'cut', a1);
    }
  }
  return { cut, fill, parts };
}

Feature line crossings

64-bit

The offset and elevation where each corridor feature line crosses the sample line.

o = −left + (left + right)·t; z = z_p + (z_q − z_p)·u

Uses the same segment crossing as the surface section, called twice to get the fraction along each segment; the elevation is a linear interpolation along the feature line.

  • sections.js
  • lines 169 to 186
  • CivilCraft 0.8.1
export function featureMarks(lines, line) {
  const out = [];
  const { a, b, left, right } = line;
  const dx = b.x - a.x, dy = b.y - a.y;
  for (const [code, pts] of Object.entries(lines || {})) {
    if (!Array.isArray(pts)) continue;
    for (let i = 1; i < pts.length; i++) {
      const p = pts[i - 1], q = pts[i];
      const tt = segCross(a.x, a.y, dx, dy, p.x, p.y, q.x - p.x, q.y - p.y);
      if (tt === null) continue;
      const u = segCross(p.x, p.y, q.x - p.x, q.y - p.y, a.x, a.y, dx, dy);
      const z = Number.isFinite(p.z) && Number.isFinite(q.z) && u !== null ? p.z + (q.z - p.z) * u : null;
      out.push({ code, o: -left + (left + right) * tt, z });
      break;
    }
  }
  return out.sort((m, n) => m.o - n.o);
}

Spans inside the graded regions

64-bit

The parts of a sample line that lie inside the graded regions, which is where the design counts.

Crossing parameters with the region outlines, then an even-odd inside test at each span's middle in plain 64-bit. A point within 1e-6 × max(1, edge length) of an outline counts as inside. Spans shorter than 1e-12 in parameter are dropped.

  • sections.js
  • lines 188 to 208
  • CivilCraft 0.8.1
export function regionSpans(line, regions) {
  const { a, b, left, right } = line, dx = b.x - a.x, dy = b.y - a.y;
  const ts = new Set([0, 1]);
  for (const r of regions || []) for (let i = 0; i < r.length; i++) {
    const p = r[i], q = r[(i + 1) % r.length];
    const tt = segCross(a.x, a.y, dx, dy, p.x, p.y, q.x - p.x, q.y - p.y);
    if (tt !== null) ts.add(tt);
  }
  const xs = [...ts].sort((m, n) => m - n), out = [];
  // inside, or on the outline (a section along a corridor's end edge is the corridor's end section)
  const onEdge = (x, y, r) => r.some((p, i) => { const q = r[(i + 1) % r.length], ex = q.x - p.x, ey = q.y - p.y, L = ex * ex + ey * ey, u = L ? Math.max(0, Math.min(1, ((x - p.x) * ex + (y - p.y) * ey) / L)) : 0; return Math.hypot(x - p.x - u * ex, y - p.y - u * ey) <= 1e-6 * Math.max(1, Math.sqrt(L)); });
  const inside = (x, y) => (regions || []).some((r) => { let c = false; for (let i = 0, j = r.length - 1; i < r.length; j = i++) if ((r[i].y > y) !== (r[j].y > y) && x < ((r[j].x - r[i].x) * (y - r[i].y)) / (r[j].y - r[i].y) + r[i].x) c = !c; return c || onEdge(x, y, r); });
  for (let i = 1; i < xs.length; i++) {
    const tm = (xs[i - 1] + xs[i]) / 2;
    if (xs[i] - xs[i - 1] < 1e-12 || !inside(a.x + dx * tm, a.y + dy * tm)) continue;
    const o0 = -left + (left + right) * xs[i - 1], o1 = -left + (left + right) * xs[i];
    const l = out[out.length - 1];
    if (l && Math.abs(l[1] - o0) < 1e-9) l[1] = o1; else out.push([o0, o1]);
  }
  return out;
}

Volumes between sections

64-bit

Cut and fill volume between neighboring sections and the running totals.

V = (A₁ + A₂) ÷ 2 × (sta₂ − sta₁); cumulative = Σ V

Average end area in 64-bit; totals are running sums in station order with no compensation. The method, not the arithmetic, sets the accuracy between sections.

  • sections.js
  • lines 268 to 276
  • CivilCraft 0.8.1
const rows = [];
let cumCut = 0, cumFill = 0;
sections.forEach((s, i) => {
  let cutVol = 0, fillVol = 0;
  if (i) { const p = sections[i - 1], L = s.sta - p.sta; cutVol = ((p.cut + s.cut) / 2) * L; fillVol = ((p.fill + s.fill) / 2) * L; }
  cumCut += cutVol; cumFill += fillVol;
  rows.push({ sta: s.sta, cut: s.cut, fill: s.fill, cutVol, fillVol, cumCut, cumFill });
});
return { sections, rows, totals: { cut: cumCut, fill: cumFill } };

Section table as CSV

64-bit

Writes station, areas, volumes and running totals to a CSV file.

areas to 2 decimals; volumes to 1 decimal; feet projects: volume ÷ 27 (cubic yards)

The exported values are rounded with toFixed: areas to 2 decimals, volumes and running totals to 1 decimal. The running totals are rounded from the unrounded sums, not added up from rounded rows. 27 cubic feet per cubic yard is exact.

  • sections.js
  • lines 280 to 289
  • CivilCraft 0.8.1
export function sectionsCsv(set, { units = 'ft', station = (s) => String(s) } = {}) {
  const ft = units !== 'm', volPer = ft ? 27 : 1;
  const head = ft
    ? [t('sec.csv.station'), t('sec.csv.cut.area.ft'), t('sec.csv.fill.area.ft'), t('sec.csv.cut.vol.ft'), t('sec.csv.fill.vol.ft'), t('sec.csv.cum.cut.ft'), t('sec.csv.cum.fill.ft')]
    : [t('sec.csv.station'), t('sec.csv.cut.area.m'), t('sec.csv.fill.area.m'), t('sec.csv.cut.vol.m'), t('sec.csv.fill.vol.m'), t('sec.csv.cum.cut.m'), t('sec.csv.cum.fill.m')];
  const f = (v, d) => v.toFixed(d);
  const lines = [head.map(csvCell).join(',')];
  for (const r of set.rows) lines.push([station(r.sta), f(r.cut, 2), f(r.fill, 2), f(r.cutVol / volPer, 1), f(r.fillVol / volPer, 1), f(r.cumCut / volPer, 1), f(r.cumFill / volPer, 1)].map(csvCell).join(','));
  return lines.join('\r\n') + '\r\n';
}

In detail

Method
Line against triangle edge intersection by parameter; trapezoid areas; V = (A₁ + A₂) ÷ 2 × (sta₂ − sta₁), as the file header states.
Number types
64-bit numbers.
Tolerances
Parallel test |denominator| < 1e-15 (absolute, in squared model units); crossing parameters accepted within 1e-12 of [0, 1] and merged when closer than 1e-12; offsets matched within 1e-9; elevation lookup tolerance 1e-7 passed to the surface; stations closer than 1e-6 are merged; a point counts as on a region outline within 1e-6 × max(1, edge length).
Limits
The average end area method overstates or understates volume where the section area does not vary linearly. The parallel test is an absolute number, so it behaves differently for very small and very large coordinates. Cubic yards in the CSV are cubic feet ÷ 27, which is exact.

Approach: 64-bit, direct

Every value is a 64-bit number from input to result. Formulas are closed form or a direct solve, and every comparison uses a stated tolerance.

Every policy that takes this approach