hex_roof source

← hex_roof · API reference · version 0.1.9

3 file(s), 1266 lines. This is the source the published version ships, and it is loft — the library is written in the same language you call it from.

Public items

A worked example is a call site the author chose as the thing to read first. Call sites are every line naming the item, found mechanically over the source with comments excluded — so an item with no worked example is not an item nobody uses.

src/hex_roof.loft

// Copyright (c) 2026 Jurjen Stellingwerff
// SPDX-License-Identifier: LGPL-3.0-or-later
// hex_roof — ROOFS, VAULTS AND ARCHES: one mechanism, several profiles.
//
// Once a cell carries a height, every roof form anyone names is the same function of a
// DISTANCE, and only the distance SOURCE changes:
//
//     d = |p - c|            a point    ->  cone, or a pyramid on a polygonal base
//     d = dist to a line     a line     ->  gable
//     d = dist to a SEGMENT  a segment  ->  hip — the ends slope because past the ridge's
//                                           end the distance stops running across the
//                                           ridge and starts running to its tip
//     d = 0                             ->  flat
//
// So "the list of roof forms" is not a list of special cases; it is a list of distance
// sources the way layer already provides. A hip is a gable with a shorter ridge —
// literally the same call with a shorter Track.
//
// A roof and an arch are the same object seen from opposite sides: a roof is the upper
// surface of a solid, an arch the lower surface of an opening. That is why they share
// everything but the sign of the profile's curvature, and why a cell wants a floor height
// AND a soffit height rather than a single z.
use hex_field::*;
use hex_way::*;

fn cellx(q: integer, r: integer) -> float {
  ((lattice_k(q, r)) as float) * 0.8660254037844386
}
fn celly(q: integer, r: integer) -> float {
  ((lattice_m(q, r)) as float) * 0.5
}

// z = peak - slope*d about a point.  Clamping d at dmax is not a detail: it is what makes
// an eave LEVEL on a quantised footprint, where the boundary cells sit at a spread of true
// radii (`eave_spread` measures the difference).
//
// AND THE CLAMP HAS ITS OWN COST, worth knowing before reaching for it: the cells it clamps
// become a flat plateau, and an interior cell inside that plateau PONDS.  On a quantised
// footprint the boundary and the interior overlap in radius, so there is no `dmax` that
// levels the eave and still drains — measured below.  `roof_hip` has neither problem,
// because it never measures a radius: the boundary is ring 0 by definition.
//
// `cx`/`cy` are WORLD coordinates, not cell indices.
//
// Example: @HXR-002
// Example: @HXR-003
pub fn roof_cone(s: HexSet, f: Heights, cx: float, cy: float,
                 peak: float, slope: float, dmax: float) {
  for r in hexset_r0(s)..hexset_r0(s) + hexset_h(s) {
    for q in hexset_q0(s)..hexset_q0(s) + hexset_w(s) {
      if hexset_get(s, q, r) {
        dx = cellx(q, r) - cx;
        dy = celly(q, r) - cy;
        d = sqrt(dx * dx + dy * dy) ?? 0.0;
        if d > dmax { d = dmax; }
        height_set(f, q, r, peak - slope * d);
      }
    }
  }
}

// The same profile about a ridge.  A full-length ridge gives a gable; shorten it and the
// ends hip in, because track_distance clamps to the segment.
//
// Example: @HXR-001
pub fn roof_ridge(s: HexSet, f: Heights, t: Track,
                  peak: float, slope: float, dmax: float) {
  for r in hexset_r0(s)..hexset_r0(s) + hexset_h(s) {
    for q in hexset_q0(s)..hexset_q0(s) + hexset_w(s) {
      if hexset_get(s, q, r) {
        d = track_distance(t, cellx(q, r), celly(q, r));
        if d > dmax { d = dmax; }
        height_set(f, q, r, peak - slope * d);
      }
    }
  }
}

// The SOFFIT of a barrel vault or an arch: spring + sqrt(R^2 - d^2), d measured from the
// crown line.  Beyond the springing the opening has closed, so the soffit is the springing
// itself.  Pass a radius larger than the half-span for a segmental arch; equal to it for a
// semicircular one.
pub fn vault_arc(s: HexSet, f: Heights, t: Track, spring: float, radius: float) {
  for r in hexset_r0(s)..hexset_r0(s) + hexset_h(s) {
    for q in hexset_q0(s)..hexset_q0(s) + hexset_w(s) {
      if hexset_get(s, q, r) {
        d = track_distance(t, cellx(q, r), celly(q, r));
        if d >= radius {
          height_set(f, q, r, spring);
        } else {
          height_set(f, q, r, spring + (sqrt(radius * radius - d * d) ?? 0.0));
        }
      }
    }
  }
}

// THE DISTANCE-TO-BOUNDARY ROOF — the one that actually matters on a quantised footprint.
//
// Every cell rises at a constant slope from its NEAREST EAVE, which is the definition of a
// hip roof. On a circular footprint that is a cone; on an octagonal one a pyramid; on an
// L-shaped hall a properly mitred hip-and-valley. The FORM is entirely the footprint's
// shape — there is one function here, not a catalogue, and that is the whole reason "roof
// forms" stopped being a list.
//
// Driving it from the FOOTPRINT rather than from the ideal circle is also what makes the
// eave level. Boundary cells are at ring 0 by definition, whatever their true radius. Drive
// a cone from each cell's own radius instead and the eave inherits the footprint's
// roundness deviation, wobbling by slope x that spread — `eave_spread` measures it, because
// the failure is invisible until you look at the eave rather than the apex.
//
// Multi-source BFS inward from the boundary, on the pre-allocated array + index-write
// idiom (appending while a head pointer walks the same vector is the broken shape).
//
// Example: @HXR-002
pub fn roof_hip(s: HexSet, f: Heights, eave: float, per_ring: float) -> integer {
  cap = hexset_w(s) * hexset_h(s);
  bq: vector<integer> = [];
  br: vector<integer> = [];
  bd: vector<integer> = [];
  for _ in 0..cap { bq += [0]; br += [0]; bd += [0]; }
  vis = hexset_chunk(hexset_q0(s), hexset_r0(s), hexset_w(s), hexset_h(s));
  n = 0;
  for r in hexset_r0(s)..hexset_r0(s) + hexset_h(s) {
    for q in hexset_q0(s)..hexset_q0(s) + hexset_w(s) {
      if hexset_get(s, q, r) {
        edge = false;
        for i in 0..6 {
          if !hexset_get(s, nb_q(q, r, i), nb_r(q, r, i)) { edge = true; }
        }
        if edge && n < cap {
          height_set(f, q, r, eave);
          hexset_set(vis, q, r, true);
          bq[n] = q; br[n] = r; bd[n] = 0;
          n = n + 1;
        }
      }
    }
  }
  head = 0;
  top = 0;
  while head < n {
    cq = bq[head] ?? 0;
    cr = br[head] ?? 0;
    cd = bd[head] ?? 0;
    head = head + 1;
    for i in 0..6 {
      nq = nb_q(cq, cr, i);
      nr = nb_r(cq, cr, i);
      if hexset_get(s, nq, nr) && !hexset_get(vis, nq, nr) && n < cap {
        hexset_set(vis, nq, nr, true);
        height_set(f, nq, nr, eave + ((cd + 1) as float) * per_ring);
        bq[n] = nq; br[n] = nr; bd[n] = cd + 1;
        n = n + 1;
        if cd + 1 > top { top = cd + 1; }
      }
    }
  }
  top
}

// A DOME — a point source with a CIRCULAR profile, where a cone is a point source with a
// linear one.  That single swap is the whole difference, and it is also why a dome is the
// harder object: a cone's slope is constant, a dome's diverges at the springing
// (dz/dr = -r/sqrt(R^2-r^2)).  Past the 45-degree point r = R/sqrt(2) the surface is more
// wall than roof and a per-cell height stops being able to carry it — which is exactly why
// real domes sit on a vertical DRUM and only the cap is vaulted.
pub fn dome(s: HexSet, f: Heights, cx: float, cy: float, base: float, radius: float) {
  for r in hexset_r0(s)..hexset_r0(s) + hexset_h(s) {
    for q in hexset_q0(s)..hexset_q0(s) + hexset_w(s) {
      if hexset_get(s, q, r) {
        dx = cellx(q, r) - cx;
        dy = celly(q, r) - cy;
        d2 = dx * dx + dy * dy;
        if d2 >= radius * radius {
          height_set(f, q, r, base);
        } else {
          height_set(f, q, r, base + (sqrt(radius * radius - d2) ?? 0.0));
        }
      }
    }
  }
}

// Two barrels over one bay.  Which way you combine them is the entire distinction between
// the two classical vaults:
//
//   GROIN    = max — the higher surface wins, so both crown lines stay at full height and
//              every wall keeps its arch.  The groins (concave) run along the diagonals.
//   CLOISTER = min — the lower wins, so the vault springs continuously from all four walls
//              and there is no wall opening anywhere.  The arrises (convex) run diagonally.
//
// Same two inputs, opposite operator, opposite building.
fn barrel_at(t: Track, x: float, y: float, spring: float, radius: float) -> float {
  d = track_distance(t, x, y);
  if d >= radius { return spring; }
  spring + (sqrt(radius * radius - d * d) ?? 0.0)
}

//
// Example: @HXR-004
pub fn vault_groin(s: HexSet, f: Heights, t1: Track, t2: Track,
                   spring: float, radius: float) {
  for r in hexset_r0(s)..hexset_r0(s) + hexset_h(s) {
    for q in hexset_q0(s)..hexset_q0(s) + hexset_w(s) {
      if hexset_get(s, q, r) {
        x = cellx(q, r); y = celly(q, r);
        a = barrel_at(t1, x, y, spring, radius);
        b = barrel_at(t2, x, y, spring, radius);
        if b > a { a = b; }
        height_set(f, q, r, a);
      }
    }
  }
}

//
// Example: @HXR-004
pub fn vault_cloister(s: HexSet, f: Heights, t1: Track, t2: Track,
                      spring: float, radius: float) {
  for r in hexset_r0(s)..hexset_r0(s) + hexset_h(s) {
    for q in hexset_q0(s)..hexset_q0(s) + hexset_w(s) {
      if hexset_get(s, q, r) {
        x = cellx(q, r); y = celly(q, r);
        a = barrel_at(t1, x, y, spring, radius);
        b = barrel_at(t2, x, y, spring, radius);
        if b < a { a = b; }
        height_set(f, q, r, a);
      }
    }
  }
}

// What can pass here: the gap between floor and soffit.
pub fn clear_height(fl: Heights, so: Heights, q: integer, r: integer) -> float {
  height_get(so, q, r) - height_get(fl, q, r)
}

// PONDING — the roof's correctness invariant, and the analogue of a stair's monotonicity.
// An INTERIOR cell (all six neighbours in the footprint) every one of whose neighbours is
// strictly higher is a local minimum: water collects and never leaves. A roof that sheds
// has none. Boundary cells are exempt — they drain off the eave, which is their job.
//
// A zero here is not "no local minima": a hollow on the BOUNDARY is not counted.
//
// Example: @HXR-005
pub fn roof_ponds(s: HexSet, f: Heights) -> integer {
  n = 0;
  for r in hexset_r0(s)..hexset_r0(s) + hexset_h(s) {
    for q in hexset_q0(s)..hexset_q0(s) + hexset_w(s) {
      if hexset_get(s, q, r) {
        interior = true;
        lower = false;
        z = height_get(f, q, r);
        for i in 0..6 {
          nq = nb_q(q, r, i);
          nr = nb_r(q, r, i);
          if !hexset_get(s, nq, nr) { interior = false; }
          if hexset_get(s, nq, nr) && height_get(f, nq, nr) < z - 0.000001 { lower = true; }
        }
        if interior && !lower { n = n + 1; }
      }
    }
  }
  n
}

// The spread of eave heights around a footprint's boundary — 0 for a level eave.
//
// THE INSTRUMENT FOR THE FAILURE NOTHING ELSE SEES: a roof built from the wrong distance
// source still has an apex, still is not flat and still sheds water.  This is the number
// that moves.
//
// Example: @HXR-002
// Example: @HXR-003
pub fn eave_spread(s: HexSet, f: Heights) -> float {
  lo = 1000000.0;
  hi = -1000000.0;
  for r in hexset_r0(s)..hexset_r0(s) + hexset_h(s) {
    for q in hexset_q0(s)..hexset_q0(s) + hexset_w(s) {
      if hexset_get(s, q, r) {
        edge = false;
        for i in 0..6 {
          if !hexset_get(s, nb_q(q, r, i), nb_r(q, r, i)) { edge = true; }
        }
        if edge {
          z = height_get(f, q, r);
          if z < lo { lo = z; }
          if z > hi { hi = z; }
        }
      }
    }
  }
  if hi < lo { return 0.0; }
  hi - lo
}

// ── RECOVERING THE FORM: the inverse of everything above ────────────────────
//
// Cells are STORAGE; the analytic surface is the truth.  A height field is read back the
// way a traced outline is matched to straights and arcs: read the per-cell z, decide the
// profile, recover its exact
// parameters, and hand the renderer a real cone — not a staircase of hexes wearing one.
//
// Nothing downstream should ever draw the quantised heights directly.  Once the form is
// recovered, the drawn surface is exact at any resolution, and the grid has served its
// purpose as a compact way to STORE the roof rather than as the shape itself.
pub const ROOF_UNKNOWN = 0;
pub const ROOF_PLANE   = 1;
pub const ROOF_CONE    = 2;
pub const ROOF_DOME    = 3;   // rf_z = base, rf_slope = sphere radius
pub const ROOF_RIDGE   = 4;   // rf_x,rf_y = one end; rf_x2,rf_y2 = the other

// A roof surface fitted to a footprint: its kind (`ROOF_*`), that kind's parameters and the
// worst residual over the footprint.  `roof_plane_fit`, `roof_cone_fit`, `roof_dome_fit`
// and `roof_ridge_fit` answer one.
pub struct RoofFit {
  pub rf_kind:  integer,
  pub rf_x:     float,     // plane: a       cone: centre x
  pub rf_y:     float,     // plane: b       cone: centre y
  pub rf_z:     float,     // plane: c       cone: apex z
  pub rf_slope: float,     // cone and ridge: dz per unit of distance
  pub rf_res:   float,     // worst residual over the footprint
  // ⚠ **THE RIDGE'S SECOND END — a segment needs four numbers where a point needs two.**
  // Zero for every other kind, which is what the existing literals give them.
  pub rf_x2:    float = 0.0,
  pub rf_y2:    float = 0.0,
}

// Least squares z = a*x + b*y + c.  Covers flat and every shed/mono-pitch.
pub fn roof_plane_fit(s: HexSet, f: Heights) -> RoofFit {
  sxx = 0.0; sxy = 0.0; sx = 0.0; syy = 0.0; sy = 0.0; sn = 0.0;
  sxz = 0.0; syz = 0.0; sz = 0.0;
  for r in hexset_r0(s)..hexset_r0(s) + hexset_h(s) {
    for q in hexset_q0(s)..hexset_q0(s) + hexset_w(s) {
      if hexset_get(s, q, r) {
        x = cellx(q, r); y = celly(q, r); z = height_get(f, q, r);
        sxx = sxx + x * x; sxy = sxy + x * y; sx = sx + x;
        syy = syy + y * y; sy = sy + y; sn = sn + 1.0;
        sxz = sxz + x * z; syz = syz + y * z; sz = sz + z;
      }
    }
  }
  det = sxx * (syy * sn - sy * sy) - sxy * (sxy * sn - sy * sx)
      + sx * (sxy * sy - syy * sx);
  ad = det; if ad < 0.0 { ad = 0.0 - ad; }
  if ad < 0.000000001 {
    return RoofFit { rf_kind: ROOF_UNKNOWN, rf_x: 0.0, rf_y: 0.0, rf_z: 0.0,
                     rf_slope: 0.0, rf_res: 1000000.0 };
  }
  da = sxz * (syy * sn - sy * sy) - sxy * (syz * sn - sy * sz)
     + sx * (syz * sy - syy * sz);
  db = sxx * (syz * sn - sy * sz) - sxz * (sxy * sn - sy * sx)
     + sx * (sxy * sz - syz * sx);
  dc = sxx * (syy * sz - syz * sy) - sxy * (sxy * sz - syz * sx)
     + sxz * (sxy * sy - syy * sx);
  a = da / det ?? 0.0;
  b = db / det ?? 0.0;
  c = dc / det ?? 0.0;
  res = 0.0;
  for r in hexset_r0(s)..hexset_r0(s) + hexset_h(s) {
    for q in hexset_q0(s)..hexset_q0(s) + hexset_w(s) {
      if hexset_get(s, q, r) {
        e = a * cellx(q, r) + b * celly(q, r) + c - height_get(f, q, r);
        if e < 0.0 { e = 0.0 - e; }
        if e > res { res = e; }
      }
    }
  }
  RoofFit { rf_kind: ROOF_PLANE, rf_x: a, rf_y: b, rf_z: c, rf_slope: 0.0, rf_res: res }
}

// Worst |z - (apex - slope*r)| for one candidate centre, with apex and slope regressed.
fn cone_at(s: HexSet, f: Heights, cx: float, cy: float) -> (float, float, float) {
  sr = 0.0; szz = 0.0; srr = 0.0; srz = 0.0; sn = 0.0;
  for r in hexset_r0(s)..hexset_r0(s) + hexset_h(s) {
    for q in hexset_q0(s)..hexset_q0(s) + hexset_w(s) {
      if hexset_get(s, q, r) {
        dx = cellx(q, r) - cx; dy = celly(q, r) - cy;
        rad = sqrt(dx * dx + dy * dy) ?? 0.0;
        z = height_get(f, q, r);
        sr = sr + rad; szz = szz + z; srr = srr + rad * rad; srz = srz + rad * z;
        sn = sn + 1.0;
      }
    }
  }
  den = sn * srr - sr * sr;
  ad = den; if ad < 0.0 { ad = 0.0 - ad; }
  if ad < 0.000000001 { return (0.0, 0.0, 1000000.0); }
  grad = (sn * srz - sr * szz) / den ?? 0.0;     // dz/dr, negative for a cone
  apex = (szz - grad * sr) / sn ?? 0.0;
  res = 0.0;
  for r in hexset_r0(s)..hexset_r0(s) + hexset_h(s) {
    for q in hexset_q0(s)..hexset_q0(s) + hexset_w(s) {
      if hexset_get(s, q, r) {
        dx = cellx(q, r) - cx; dy = celly(q, r) - cy;
        rad = sqrt(dx * dx + dy * dy) ?? 0.0;
        e = apex + grad * rad - height_get(f, q, r);
        if e < 0.0 { e = 0.0 - e; }
        if e > res { res = e; }
      }
    }
  }
  (apex, 0.0 - grad, res)
}

// ── A RIDGE — the cone's own profile about a SEGMENT ────────────────────
//
// `roof_ridge` draws `peak - slope * d` where `d` is the distance to a LINE SEGMENT, exactly
// as `roof_cone` draws it about a point.  So the fit is `cone_at` with the radius replaced by
// `track_distance` — the same primitive the generator uses, rather than a second opinion about
// what distance-to-a-segment means.
fn ridge_at(s: HexSet, f: Heights, ax: float, ay: float, bx: float, by: float)
    -> (float, float, float) {
  t = track_new();
  track_straight(t, ax, ay, bx, by);
  sd = 0.0; szz = 0.0; sdd = 0.0; sdz = 0.0; sn = 0.0;
  for r in hexset_r0(s)..hexset_r0(s) + hexset_h(s) {
    for q in hexset_q0(s)..hexset_q0(s) + hexset_w(s) {
      if hexset_get(s, q, r) {
        d = track_distance(t, cellx(q, r), celly(q, r));
        z = height_get(f, q, r);
        sd = sd + d; szz = szz + z; sdd = sdd + d * d; sdz = sdz + d * z;
        sn = sn + 1.0;
      }
    }
  }
  den = sn * sdd - sd * sd;
  ad = den; if ad < 0.0 { ad = 0.0 - ad; }
  if ad < 0.000000001 { return (0.0, 0.0, 1000000.0); }
  grad = (sn * sdz - sd * szz) / den ?? 0.0;
  apex = (szz - grad * sd) / sn ?? 0.0;
  res = 0.0;
  for r in hexset_r0(s)..hexset_r0(s) + hexset_h(s) {
    for q in hexset_q0(s)..hexset_q0(s) + hexset_w(s) {
      if hexset_get(s, q, r) {
        d = track_distance(t, cellx(q, r), celly(q, r));
        e = apex + grad * d - height_get(f, q, r);
        if e < 0.0 { e = 0.0 - e; }
        if e > res { res = e; }
      }
    }
  }
  (apex, 0.0 - grad, res)
}

// Fit a ridge — the gable / hip profile `roof_ridge` draws — to the heights over these cells.
//
// FOUR SEEDS, AND THE BEST RESIDUAL WINS.  The one that works on a gable comes from the data:
// the ridge is the LOCUS OF THE MAXIMUM, so the two highest cells furthest apart give its
// direction, and the seed is that line extended across the whole footprint.  The climb then
// SHORTENS it for a hipped roof, which is the direction it handles well.  The other three —
// both ends at the highest cell, both at `roof_cone_fit`'s centre, one at each — are cheap
// insurance: a climb grown out of a POINT can lower its residual by running the segment
// ACROSS the ridge until the fit is "flat at the mean", and a cone fitted to a gable puts its
// apex far outside the footprint.  Taking the better of several starts is not a threshold;
// when they converge it is the same answer.
//
// Each end is climbed in the eight directions `roof_cone_fit` uses, on the same step schedule.
pub fn roof_ridge_fit(s: HexSet, f: Heights) -> RoofFit {
  // the two candidate anchors: the highest cell, which lies ON the ridge, and the cone fit's
  // centre, which is where a cone lands when asked to be a ridge
  hi = 0.0 - 1000000.0;
  lo = 1000000.0;
  hx = 0.0; hy = 0.0;
  for r in hexset_r0(s)..hexset_r0(s) + hexset_h(s) {
    for q in hexset_q0(s)..hexset_q0(s) + hexset_w(s) {
      if hexset_get(s, q, r) {
        z = height_get(f, q, r);
        if z > hi { hi = z; hx = cellx(q, r); hy = celly(q, r); }
        if z < lo { lo = z; }
      }
    }
  }
  // ⚠ **THE MAXIMUM IS ATTAINED TO WITHIN ROUNDING, NOT EXACTLY.** On a field `roof_ridge`
  // drew, `track_distance` is zero on the segment only up to rounding, so `z == hi` would let
  // one cell win, collapse the pair to a point and degenerate the fit.
  //
  // ⚠ **THE BAND IS THE FIELD'S OWN SPAN, NOT A TUNED EPSILON**: a millionth of `hi - lo` is
  // float equality made scale-aware, and it carries no units of its own.
  band = (hi - lo) * 0.000001;
  if band < 0.0 { band = 0.0 - band; }
  // THE RIDGE IS THE LOCUS OF THE MAXIMUM, so the cells attaining it ARE the ridge, and the
  // two furthest apart give its direction.
  ex1 = hx; ey1 = hy; ex2 = hx; ey2 = hy;
  far = 0.0 - 1.0;
  for r in hexset_r0(s)..hexset_r0(s) + hexset_h(s) {
    for q in hexset_q0(s)..hexset_q0(s) + hexset_w(s) {
      if hexset_get(s, q, r) && height_get(f, q, r) >= hi - band {
        for r2 in hexset_r0(s)..hexset_r0(s) + hexset_h(s) {
          for q2 in hexset_q0(s)..hexset_q0(s) + hexset_w(s) {
            if hexset_get(s, q2, r2) && height_get(f, q2, r2) >= hi - band {
              ddx = cellx(q2, r2) - cellx(q, r);
              ddy = celly(q2, r2) - celly(q, r);
              dd = ddx * ddx + ddy * ddy;
              if dd > far {
                far = dd;
                ex1 = cellx(q, r); ey1 = celly(q, r);
                ex2 = cellx(q2, r2); ey2 = celly(q2, r2);
              }
            }
          }
        }
      }
    }
  }
  // THE SEED IS EXTENDED TO SPAN THE FOOTPRINT, because growing is what the climb does badly:
  // the locus gives the right direction and too short a segment, and every cell past an end
  // would get a clamped distance.  A gable's ridge runs the length of its roof; a hipped roof
  // is the case the climb then SHORTENS.
  pdx = ex2 - ex1;
  pdy = ey2 - ey1;
  plen = sqrt(pdx * pdx + pdy * pdy) ?? 0.0;
  if plen > 0.000001 {
    ux = pdx / plen ?? 0.0;
    uy = pdy / plen ?? 0.0;
    tmin = 1000000.0;
    tmax = 0.0 - 1000000.0;
    for r in hexset_r0(s)..hexset_r0(s) + hexset_h(s) {
      for q in hexset_q0(s)..hexset_q0(s) + hexset_w(s) {
        if hexset_get(s, q, r) {
          tp = (cellx(q, r) - ex1) * ux + (celly(q, r) - ey1) * uy;
          if tp < tmin { tmin = tp; }
          if tp > tmax { tmax = tp; }
        }
      }
    }
    ex2 = ex1 + ux * tmax; ey2 = ey1 + uy * tmax;
    ex1 = ex1 + ux * tmin; ey1 = ey1 + uy * tmin;
  }
  c = roof_cone_fit(s, f);
  best = ridge_climb(s, f, ex1, ey1, ex2, ey2);
  t1 = ridge_climb(s, f, hx, hy, hx, hy);
  if t1.rf_res < best.rf_res { best = t1; }
  t2 = ridge_climb(s, f, c.rf_x, c.rf_y, c.rf_x, c.rf_y);
  if t2.rf_res < best.rf_res { best = t2; }
  t3 = ridge_climb(s, f, hx, hy, c.rf_x, c.rf_y);
  if t3.rf_res < best.rf_res { best = t3; }
  best
}

// The climb itself, from a given pair of ends.  Each end moves in the eight directions
// `roof_cone_fit` uses, on the same step schedule — a second search rule here would be a second
// opinion about how this family searches.
fn ridge_climb(s: HexSet, f: Heights, sax: float, say: float, sbx: float, sby: float)
    -> RoofFit {
  ax = sax; ay = say; bx = sbx; by = sby;
  (ba, bs, br) = ridge_at(s, f, ax, ay, bx, by);
  step = 0.4;
  for _ in 0..5 {
    moved = true;
    while moved {
      moved = false;
      for i in 0..16 {
        ox = 0.0; oy = 0.0;
        k = i % 8 ?? 0;
        if k == 0 { ox = step; }
        if k == 1 { ox = 0.0 - step; }
        if k == 2 { oy = step; }
        if k == 3 { oy = 0.0 - step; }
        if k == 4 { ox = step; oy = step; }
        if k == 5 { ox = step; oy = 0.0 - step; }
        if k == 6 { ox = 0.0 - step; oy = step; }
        if k == 7 { ox = 0.0 - step; oy = 0.0 - step; }
        nax = ax; nay = ay; nbx = bx; nby = by;
        if i < 8 { nax = ax + ox; nay = ay + oy; }
        else { nbx = bx + ox; nby = by + oy; }
        (ca, cs, cr) = ridge_at(s, f, nax, nay, nbx, nby);
        if cr < br - 0.0000000001 {
          br = cr; ba = ca; bs = cs;
          ax = nax; ay = nay; bx = nbx; by = nby;
          moved = true;
        }
      }
    }
    step = step * 0.25;
  }
  RoofFit { rf_kind: ROOF_RIDGE, rf_x: ax, rf_y: ay, rf_z: ba, rf_slope: bs, rf_res: br,
            rf_x2: bx, rf_y2: by }
}

// Fit a cone to the heights over these cells: apex, slope and how well it fits.  The centre
// is seeded at the highest cell and then refined — the apex cell is only within one
// circumradius of the true apex, so the refinement is what makes the recovery exact.
// `roof_match` asks this, the plane, the dome and the ridge, and keeps what fits.
pub fn roof_cone_fit(s: HexSet, f: Heights) -> RoofFit {
  hi = -1000000.0;
  bx = 0.0; by = 0.0;
  for r in hexset_r0(s)..hexset_r0(s) + hexset_h(s) {
    for q in hexset_q0(s)..hexset_q0(s) + hexset_w(s) {
      if hexset_get(s, q, r) {
        z = height_get(f, q, r);
        if z > hi { hi = z; bx = cellx(q, r); by = celly(q, r); }
      }
    }
  }
  (ba, bs, br) = cone_at(s, f, bx, by);
  step = 0.4;
  for _ in 0..4 {
    moved = true;
    while moved {
      moved = false;
      for i in 0..8 {
        ox = 0.0; oy = 0.0;
        if i == 0 { ox = step; }
        if i == 1 { ox = 0.0 - step; }
        if i == 2 { oy = step; }
        if i == 3 { oy = 0.0 - step; }
        if i == 4 { ox = step; oy = step; }
        if i == 5 { ox = step; oy = 0.0 - step; }
        if i == 6 { ox = 0.0 - step; oy = step; }
        if i == 7 { ox = 0.0 - step; oy = 0.0 - step; }
        (ca, cs, cr) = cone_at(s, f, bx + ox, by + oy);
        if cr < br - 0.0000000001 {
          br = cr; ba = ca; bs = cs; bx = bx + ox; by = by + oy;
          moved = true;
        }
      }
    }
    step = step * 0.25;
  }
  RoofFit { rf_kind: ROOF_CONE, rf_x: bx, rf_y: by, rf_z: ba, rf_slope: bs, rf_res: br }
}

// A DOME.  Given a centre, the fit is LINEAR after one substitution: from
// z = base + sqrt(R^2 - r^2) we get (z^2 + r^2) = 2*base*z + (R^2 - base^2), so regressing
// (z^2 + r^2) on z hands back base as half the gradient and R from the intercept.  The
// centre is then refined by the same local search the cone uses.
fn dome_at(s: HexSet, f: Heights, cx: float, cy: float) -> (float, float, float) {
  sv = 0.0; su = 0.0; svv = 0.0; suv = 0.0; sn = 0.0;
  for r in hexset_r0(s)..hexset_r0(s) + hexset_h(s) {
    for q in hexset_q0(s)..hexset_q0(s) + hexset_w(s) {
      if hexset_get(s, q, r) {
        dx = cellx(q, r) - cx; dy = celly(q, r) - cy;
        z = height_get(f, q, r);
        v = z;
        u = z * z + dx * dx + dy * dy;
        sv = sv + v; su = su + u; svv = svv + v * v; suv = suv + u * v; sn = sn + 1.0;
      }
    }
  }
  den = sn * svv - sv * sv;
  ad = den; if ad < 0.0 { ad = 0.0 - ad; }
  if ad < 0.000000001 { return (0.0, 0.0, 1000000.0); }
  grad = (sn * suv - sv * su) / den ?? 0.0;
  icept = (su - grad * sv) / sn ?? 0.0;
  base = grad * 0.5;
  rr2 = icept + base * base;
  if rr2 <= 0.0 { return (0.0, 0.0, 1000000.0); }
  rad = sqrt(rr2) ?? 0.0;
  res = 0.0;
  for r in hexset_r0(s)..hexset_r0(s) + hexset_h(s) {
    for q in hexset_q0(s)..hexset_q0(s) + hexset_w(s) {
      if hexset_get(s, q, r) {
        dx = cellx(q, r) - cx; dy = celly(q, r) - cy;
        d2 = dx * dx + dy * dy;
        pred = base;
        if d2 < rad * rad { pred = base + (sqrt(rad * rad - d2) ?? 0.0); }
        e = pred - height_get(f, q, r);
        if e < 0.0 { e = 0.0 - e; }
        if e > res { res = e; }
      }
    }
  }
  (base, rad, res)
}

// Fit a dome to the same heights, reported the same way, so the two are comparable.
pub fn roof_dome_fit(s: HexSet, f: Heights) -> RoofFit {
  hi = -1000000.0;
  bx = 0.0; by = 0.0;
  for r in hexset_r0(s)..hexset_r0(s) + hexset_h(s) {
    for q in hexset_q0(s)..hexset_q0(s) + hexset_w(s) {
      if hexset_get(s, q, r) {
        z = height_get(f, q, r);
        if z > hi { hi = z; bx = cellx(q, r); by = celly(q, r); }
      }
    }
  }
  (bb, brad, bres) = dome_at(s, f, bx, by);
  step = 0.4;
  for _ in 0..4 {
    moved = true;
    while moved {
      moved = false;
      for i in 0..8 {
        ox = 0.0; oy = 0.0;
        if i == 0 { ox = step; }
        if i == 1 { ox = 0.0 - step; }
        if i == 2 { oy = step; }
        if i == 3 { oy = 0.0 - step; }
        if i == 4 { ox = step; oy = step; }
        if i == 5 { ox = step; oy = 0.0 - step; }
        if i == 6 { ox = 0.0 - step; oy = step; }
        if i == 7 { ox = 0.0 - step; oy = 0.0 - step; }
        (cb, crad, cres) = dome_at(s, f, bx + ox, by + oy);
        if cres < bres - 0.0000000001 {
          bres = cres; bb = cb; brad = crad; bx = bx + ox; by = by + oy;
          moved = true;
        }
      }
    }
    step = step * 0.25;
  }
  RoofFit { rf_kind: ROOF_DOME, rf_x: bx, rf_y: by, rf_z: bb, rf_slope: brad, rf_res: bres }
}

// Decide what the field IS.  A plane wins ties: it is the simpler statement, and a cone
// fitted to a flat roof would put its apex somewhere arbitrary.
//
// Example: @HXR-006
//
// It tries plane, cone, dome, then ridge, and answers the first whose worst residual is
// within `tol`.  The ridge is LAST because a cone is a ridge whose segment has zero length:
// trying it first would answer `ROOF_RIDGE` for a true cone, the less specific statement.
//
// ⛔ **WHEN NOTHING FITS, IT SAYS SO — `ROOF_UNKNOWN`**, never the least-bad shape under a
// kind it is not; `roof_eval` draws nothing for it.  Every other field is the closest
// candidate's — a ridge's both ends included — because which shape came closest and by how
// much is what a caller diagnoses with.
pub fn roof_match(s: HexSet, f: Heights, tol: float) -> RoofFit {
  p = roof_plane_fit(s, f);
  if p.rf_res <= tol { return p; }
  c = roof_cone_fit(s, f);
  if c.rf_res <= tol { return c; }
  d = roof_dome_fit(s, f);
  if d.rf_res <= tol { return d; }
  g = roof_ridge_fit(s, f);
  if g.rf_res <= tol { return g; }
  best = p;
  if c.rf_res < best.rf_res { best = c; }
  if d.rf_res < best.rf_res { best = d; }
  if g.rf_res < best.rf_res { best = g; }
  RoofFit { rf_kind: ROOF_UNKNOWN, rf_x: best.rf_x, rf_y: best.rf_y, rf_z: best.rf_z,
            rf_slope: best.rf_slope, rf_res: best.rf_res, rf_x2: best.rf_x2, rf_y2: best.rf_y2 }
}

// The recovered ANALYTIC surface — what a renderer should draw, at any resolution, instead
// of the cells it came from.
//
// Example: @HXR-006
pub fn roof_eval(fit: RoofFit, x: float, y: float) -> float {
  if fit.rf_kind == ROOF_PLANE { return fit.rf_x * x + fit.rf_y * y + fit.rf_z; }
  if fit.rf_kind == ROOF_CONE {
    dx = x - fit.rf_x; dy = y - fit.rf_y;
    return fit.rf_z - fit.rf_slope * (sqrt(dx * dx + dy * dy) ?? 0.0);
  }
  if fit.rf_kind == ROOF_RIDGE {
    rt = track_new();
    track_straight(rt, fit.rf_x, fit.rf_y, fit.rf_x2, fit.rf_y2);
    return fit.rf_z - fit.rf_slope * track_distance(rt, x, y);
  }
  if fit.rf_kind == ROOF_DOME {
    dx = x - fit.rf_x; dy = y - fit.rf_y;
    d2 = dx * dx + dy * dy;
    if d2 >= fit.rf_slope * fit.rf_slope { return fit.rf_z; }
    return fit.rf_z + (sqrt(fit.rf_slope * fit.rf_slope - d2) ?? 0.0);
  }
  0.0
}

tests/01-hex-roof.loft

// Copyright (c) 2026 Jurjen Stellingwerff
// SPDX-License-Identifier: LGPL-3.0-or-later
// hex_roof — the conformance tests: a roof is a HEIGHT FIELD over a cell set, and it must drain.
//
// These assert the library's CONTRACT, not its depth.  The exhaustive measurement lives in
// hexbody's `tests/house.loft` (roofs at 12 orientations, `eave_spread`, `roof_ponds = 0`), which
// prints its evidence and carries a control that must fire.  This file is what CI runs under
// `LOFT_DENY_WARNINGS=1`.
//
// ⚠ `roof_match` takes a TOLERANCE and is an R2 fit by design — recovering a profile from an
// arbitrary height field, not from one we authored.  That is licensed; it is not an `ε` smuggled
// into an exact path.
use hex_roof::*;
use hex_way::*;
use hex_field::*;

const W = 9;
const H = 9;

fn disc() -> HexSet {
  s = hexset_chunk(0, 0, W, H);
  for r in 0..H {
    for q in 0..W {
      dq = q - 4;
      dr = r - 4;
      if dq * dq + dr * dr <= 9 { hexset_set(s, q, r, true); }
    }
  }
  s
}

fn main() {
  s = disc();
  assert(hexset_count(s) > 0, "the test disc has cells");
}

fn test_a_cone_rises_to_a_single_peak() {
  // ⚠ `roof_cone`'s (cx, cy) is in WORLD units, not cell indices — cell (4,4) is NOT at world
  // (4.0, 4.0) on a pointy-top odd-r lattice.  Asserting "cell (4,4) is the apex" fails with 9
  // cells above it, which is a wrong assumption rather than a defect.  So assert the PROPERTY:
  // the cone is not flat, and its maximum is positive.
  s = disc();
  f = heights_new(0, 0, W, H);
  roof_cone(s, f, 4.0, 4.0, 5.0, 1.0, 9.0);
  hi = 0.0;
  lo = 1000000.0;
  for r in 0..H {
    for q in 0..W {
      if hexset_get(s, q, r) {
        h = height_get(f, q, r);
        if h > hi { hi = h; }
        if h < lo { lo = h; }
      }
    }
  }
  assert(hi > 0.0, "the cone has a positive maximum, got {hi}");
  assert(hi > lo, "the cone is not flat: max {hi} must exceed min {lo}");
}

fn test_a_cone_has_no_ponds() {
  // A pond is a cell with no downhill neighbour — water that cannot leave. A roof that ponds is
  // a roof that leaks, so this is the load-bearing property of every profile in this library.
  s = disc();
  f = heights_new(0, 0, W, H);
  roof_cone(s, f, 4.0, 4.0, 5.0, 1.0, 9.0);
  p = roof_ponds(s, f);
  assert(p == 0, "a cone must not pond, found {p}");
}

fn test_a_flat_field_does_pond_which_keeps_the_pond_check_honest() {
  // The control: if `roof_ponds` returned 0 for everything the test above would prove nothing.
  s = disc();
  f = heights_new(0, 0, W, H);
  p = roof_ponds(s, f);
  assert(p > 0, "a dead-flat field must pond, but roof_ponds returned {p}");
}

// ── ⛔ A FALLBACK IS NOT A RECOVERY ────────────────────────────────────────
//
// `roof_match` returned the least-bad of plane/cone/dome with `rf_kind` naming a shape the
// field is NOT, and the residual was the only signal. It says `ROOF_UNKNOWN` now.
//
// ⚠ THE FIXTURES ARE FIELDS THIS PACKAGE ITSELF DRAWS, which is what makes the gap a fact
// about the package rather than about a fixture somebody invented: it draws a ridge, a hip
// and three vaults, and fits a plane, a cone and a dome.
fn test_a_ridge_recovers_as_a_ridge() {
  // ⛔ **WHAT `roof_match` USED TO ANSWER HERE WAS `ROOF_CONE` AT RESIDUAL 1.209**, with an
  // apex of 6.318 where the truth is 6.0 and a slope of 0.409 where it is 0.5 — a shape the
  // field is not, named with confidence. `roof_ridge_fit` is the cone's own profile about a
  // SEGMENT, so the answer is exact.
  s = disc();
  f = heights_new(0, 0, W, H);
  t = track_new();
  track_straight(t, 1.0, 3.0, 6.0, 3.0);
  roof_ridge(s, f, t, 6.0, 0.5, 100.0);
  m = roof_match(s, f, 0.001);
  assert(m.rf_kind == ROOF_RIDGE, "a ridge recovers as a ridge, got kind {m.rf_kind}");
  assert(m.rf_res <= 0.001, "…exactly: residual {m.rf_res}");
  // ⚠ **THE PARAMETERS, NOT JUST THE LABEL** — a kind that is right about nothing else is a
  // label, and `roof_eval` draws from these.
  assert(abs(m.rf_z - 6.0) < 0.01, "the peak comes back 6.0, got {m.rf_z}");
  assert(abs(m.rf_slope - 0.5) < 0.01, "and the slope 0.5, got {m.rf_slope}");
  // …and the recovered surface agrees with the field it was recovered from, at a cell centre.
  wx = (lattice_k(4, 4) as float) * 0.8660254037844386;
  wy = (lattice_m(4, 4) as float) * 0.5;
  assert(abs(roof_eval(m, wx, wy) - height_get(f, 4, 4)) < 0.01,
         "the recovered ridge evaluates to the field: {roof_eval(m, wx, wy)} against {height_get(f, 4, 4)}");
}

fn test_a_hip_is_still_unknown_and_that_keeps_the_refusal_honest() {
  // ⚠ **THE REFUSAL PATH MUST STAY REACHABLE.** With the ridge added, three of the shapes this
  // package draws now recover; a hip does not, and saying so is what stops `ROOF_UNKNOWN`
  // becoming a branch nothing takes — an unreachable refusal is one nobody would notice
  // breaking.
  s = disc();
  f = heights_new(0, 0, W, H);
  roof_hip(s, f, 2.0, 0.5);
  m = roof_match(s, f, 0.001);
  assert(m.rf_kind == ROOF_UNKNOWN, "a hip is none of the four, got {m.rf_kind}");
  assert(roof_eval(m, 0.0, 0.0) == 0.0, "and an unknown roof draws nothing");
}

fn test_an_unknown_roof_reports_the_closest_ridge_whole() {
  // When nothing fits and the ridge came closest, the answer must carry the ridge's segment —
  // BOTH ends.  A gable over a 7 x 3 block with one cell dented 0.3 fits no shape within 0.01,
  // and the ridge is the closest: its end on the far side of the block is what a caller needs
  // to see which ridge nearly fitted.  An end of (0, 0) would read as a ridge to the origin.
  s = hexset_new(4);
  for r in -1..2 { for q in -3..4 { hexset_set(s, q, r, true); } }
  full = track_new();
  track_straight(full, -3.0 * 1.7320508075688772, 0.0, 3.0 * 1.7320508075688772, 0.0);
  f = heights_new(-4, -4, 9, 9);
  roof_ridge(s, f, full, 5.0, 1.0, 100.0);
  height_set(f, 1, 0, height_get(f, 1, 0) - 0.3);
  g = roof_ridge_fit(s, f);
  m = roof_match(s, f, 0.01);
  assert(m.rf_kind == ROOF_UNKNOWN && m.rf_res == g.rf_res, "the ridge is closest: {m.rf_res} vs {g.rf_res}");
  assert(g.rf_x < 0.0 && g.rf_x2 > 0.0, "the fitted ridge spans the block: {g.rf_x} to {g.rf_x2}");
  assert(m.rf_x == g.rf_x && m.rf_x2 == g.rf_x2 && m.rf_y2 == g.rf_y2,
         "the unknown answer keeps the ridge's far end: ({m.rf_x2}, {m.rf_y2}) vs ({g.rf_x2}, {g.rf_y2})");
}

fn test_the_three_kinds_it_does_know_are_unaffected() {
  // ⚠ THE CONTROL, AND IT IS HALF THE POINT: a doorstep that answered UNKNOWN to everything
  // would pass both rows above. A cone this package drew must still come back a CONE, exactly.
  s = disc();
  fc = heights_new(0, 0, W, H);
  roof_cone(s, fc, 3.0, 1.5, 6.0, 0.5, 1000.0);
  mc = roof_match(s, fc, 0.001);
  assert(mc.rf_kind == ROOF_CONE, "a cone still recovers as a cone, got {mc.rf_kind}");
  assert(mc.rf_res <= 0.001, "…within tolerance: {mc.rf_res}");

  // ⛔ **IN WORLD UNITS, NOT CELL INDICES** — the same trap this file already warns about at
  // `test_a_cone_rises_to_a_single_peak`, and I walked into it one function later: on a
  // pointy-top odd-r lattice `q` is not the world x, so `2 + 0.25q` is a STAIRCASE and the
  // plane fit correctly refuses it. `cellx` is `lattice_k * √3/2`, which is what a plane is
  // linear in.
  fp = heights_new(0, 0, W, H);
  for r in 0..H { for q in 0..W {
    if hexset_get(s, q, r) {
      wx = (lattice_k(q, r) as float) * 0.8660254037844386;
      wy = (lattice_m(q, r) as float) * 0.5;
      height_set(fp, q, r, 2.0 + 0.25 * wx + 0.125 * wy);
    }
  } }
  mp = roof_match(s, fp, 0.001);
  assert(mp.rf_kind == ROOF_PLANE, "and a plane as a plane, got {mp.rf_kind}");
  // ⛔ **AND A CONE MUST NOT BECOME A RIDGE.** A cone IS a ridge whose segment has zero length,
  // so a fitter tried too early would answer the less specific shape for both — which is why
  // `roof_match` asks the ridge last.
  assert(mc.rf_kind != ROOF_RIDGE, "a cone is not answered as a degenerate ridge");
}

// ── ⛔ THE SHAPE A REAL EDITOR DRAWS — a gable over a rectangle ─────────────
//
// ⛔ **THE DISC FIXTURE ABOVE IS NOT ENOUGH, AND A CONSUMER FOUND THAT.** moros's `roof_over`
// writes a 27-cell gable — five rows across, CONSTANT along the ridge — and 0.1.2 returned a
// degenerate segment eight world units clear of it with slope 0 and residual 1.185, while the
// disc fixture passed. A suite written by the author cannot ask the questions the author does
// not have; this is that question, in the library.
//
// The field is exact: `7.75 - (2/3)|y - 6|`, which is what a 0.25-unit store holds for that
// gable — 5.75 / 6.75 / 7.75 / 6.75 / 5.75 on rows at y = 3, 4.5, 6, 7.5, 9.
fn gable_rows() -> HexSet {
  s = hexset_chunk(0, 0, W, H);
  for r in 2..7 { for q in 1..6 { hexset_set(s, q, r, true); } }
  s
}

fn test_a_gable_over_a_rectangle_recovers() {
  s = gable_rows();
  f = heights_new(0, 0, W, H);
  for r in 0..H { for q in 0..W {
    if hexset_get(s, q, r) {
      wy = (lattice_m(q, r) as float) * 0.5;
      d = wy - 6.0;
      if d < 0.0 { d = 0.0 - d; }
      height_set(f, q, r, 7.75 - 0.6666666666666666 * d);
    }
  } }
  g = roof_ridge_fit(s, f);
  assert(g.rf_res < 0.01,
         "the gable does not recover: residual {g.rf_res}, segment {g.rf_x},{g.rf_y} to {g.rf_x2},{g.rf_y2}, slope {g.rf_slope}");
  assert(abs(g.rf_slope - 0.6666666666666666) < 0.01, "the slope comes back, got {g.rf_slope}");
  assert(abs(g.rf_z - 7.75) < 0.02, "and the peak, got {g.rf_z}");
  // ⚠ **AND THE RIDGE RUNS ALONG THE ROWS, NOT ACROSS THEM** — a segment of near-zero length
  // is a cone wearing a ridge's name, and it fits this field badly.
  dx = g.rf_x2 - g.rf_x; dy = g.rf_y2 - g.rf_y;
  // ⚠ `??` BINDS LOOSEST, so the parentheses are load-bearing: `sqrt(..) ?? 0.0 > 1.0` parses
  // as `sqrt(..) ?? (0.0 > 1.0)` and is a type error, not a comparison.
  seglen = sqrt(dx * dx + dy * dy) ?? 0.0;
  assert(seglen > 1.0, "the recovered ridge has length, got {seglen}");
  m = roof_match(s, f, 0.01);
  assert(m.rf_kind == ROOF_RIDGE, "and `roof_match` names it a ridge, got kind {m.rf_kind}");
}

fn test_the_package_recovers_what_it_draws() {
  // ⛔ **THE STRONGEST FIXTURE THERE IS, AND IT WAS MISSING.** Every ridge test above builds its
  // field by hand; this one hands `roof_ridge_fit` the output of `roof_ridge`. 0.1.3 failed it
  // at residual 1.132 — the locus was collected with `z == hi`, and `track_distance` returns
  // zero on the segment only up to ROUNDING, so one cell won the equality, the pair collapsed
  // to a point and the fit degenerated. A hand-built field attains its maximum exactly and
  // hid it.
  s = disc();
  f = heights_new(0, 0, W, H);
  t = track_new();
  track_straight(t, 0.0 - 12.0, 3.0, 12.0, 3.0);
  roof_ridge(s, f, t, 6.0, 0.5, 1000.0);
  g = roof_ridge_fit(s, f);
  assert(g.rf_res < 0.01,
         "the package cannot recover its own ridge: residual {g.rf_res}, segment {g.rf_x},{g.rf_y} to {g.rf_x2},{g.rf_y2}");
  assert(abs(g.rf_slope - 0.5) < 0.01, "the slope comes back, got {g.rf_slope}");
  assert(abs(g.rf_z - 6.0) < 0.02, "and the peak, got {g.rf_z}");
}

tests/02-worked-examples.loft

// Copyright (c) 2026 Jurjen Stellingwerff
// SPDX-License-Identifier: LGPL-3.0-or-later
//
// Worked examples for `hex_roof` (@PLN141).  Each `test_*` below is a REAL,
// run-every-CI demonstration of a contract the signature cannot state; the
// function documents itself by CITING the tag with `// Example: @HXR-00N`, and
// `check_doc_drift.sh examples` fails if a citation goes dangling.  These are
// teaching examples, distinct from the conformance suite in `01-hex-roof.loft`.
//
// THE ONE IDEA, AND WHERE IT LEAKS.  Every roof form here is the same function
// of a DISTANCE, and only the distance SOURCE changes: a point gives a cone, a
// line a gable, a segment a hip, the FOOTPRINT ITSELF a proper hip-and-valley.
// So this is not a catalogue of forms and there are no special cases to learn.
//
// What a caller has to get right instead is which source a profile reads, and
// that is where every example below lives.  The trouble is that a roof built
// from the wrong source still looks like a roof: it still has an apex, it still
// sheds water, and it still passes every property the conformance suite checks.
// The damage shows at the EAVE, which is the last place anyone looks.
use hex_roof::*;
use hex_way::*;
use hex_field::*;

const W = 9;
const H = 9;

// A quantised disc: 29 cells, and deliberately not a circle — that is the whole
// point of everything below.
fn disc() -> HexSet {
  s = hexset_chunk(0, 0, W, H);
  for r in 0..H {
    for q in 0..W {
      dq = q - 4;
      dr = r - 4;
      if dq * dq + dr * dr <= 9 { hexset_set(s, q, r, true); }
    }
  }
  s
}

// The lattice -> world map, which is what this module measures distances in.
fn wx(q: integer, r: integer) -> float { ((lattice_k(q, r)) as float) * 0.8660254037844386 }
fn wy(q: integer, r: integer) -> float { ((lattice_m(q, r)) as float) * 0.5 }

fn peak_of(s: HexSet, f: Heights) -> float {
  m = -1000000.0;
  for r in 0..H { for q in 0..W { if hexset_get(s, q, r) {
    z = height_get(f, q, r);
    if z > m { m = z; }
  } } }
  m
}

// @HXR-001 — A HIP IS A GABLE WITH A SHORTER RIDGE.  Not "like" one: the same
// call, with a shorter `Track`.  `track_distance` clamps to the segment, so past
// the ridge's end the distance stops running across the ridge and starts running
// to its tip — and the ends slope in on their own.
//
// This is what "the list of roof forms is a list of distance sources" means in
// practice.  If you are reaching for a `roof_gable` and a `roof_hip_roof`, stop:
// you are choosing a Track length.
fn test_a_hip_is_a_gable_with_a_shorter_ridge() {
  s = disc();

  // One ridge the full width of the footprint, and one covering its middle third.
  full = track_new();
  track_straight(full, wx(0, 4), wy(4, 4), wx(8, 4), wy(4, 4));
  short = track_new();
  track_straight(short, wx(3, 4), wy(4, 4), wx(5, 4), wy(4, 4));

  gable = heights_new(0, 0, W, H);
  hip = heights_new(0, 0, W, H);
  roof_ridge(s, gable, full, 5.0, 1.0, 100.0);
  roof_ridge(s, hip, short, 5.0, 1.0, 100.0);

  // Over the middle the two are identical — both cells sit on their ridge.
  assert(height_get(gable, 4, 4) == 5.0, "the gable is at full height on the ridge");
  assert(height_get(hip, 4, 4) == 5.0, "and so is the hip — same ridge, same height");

  // Near the end they are not: the gable still runs at full height to the wall,
  // the hip has sloped away because the ridge ran out.
  assert(height_get(gable, 1, 4) == 5.0, "the gable carries the ridge to the end wall");
  assert(height_get(hip, 1, 4) < 2.0, "the hip has sloped in, to {height_get(hip, 1, 4)}");

  // Both drain, because both are one profile of one distance.
  assert(roof_ponds(s, gable) == 0 && roof_ponds(s, hip) == 0, "neither ponds");
}

// @HXR-002 — ON A QUANTISED FOOTPRINT A POINT-SOURCE CONE CANNOT GET BOTH THE
// EAVE AND THE DRAINAGE RIGHT.  Drive the roof from each cell's own radius and
// the eave inherits the footprint's roundness deviation.  Clamp at `dmax` to
// level it — which is what the clamp is for — and the clamped cells become a
// flat plateau that ponds.  There is no `dmax` that escapes both, and the reason
// is measurable: on a quantised disc the boundary and the interior OVERLAP in
// radius, so no single radius separates them.
//
// `roof_hip` has neither problem, because it does not measure a radius at all:
// it walks inward from the boundary, and the boundary is ring 0 BY DEFINITION,
// whatever radius its cells happen to sit at.
//
// Reach for `roof_cone` when the footprint really is a circle you control, and
// for `roof_hip` when it is a footprint someone drew.
fn test_the_eave_is_where_a_cone_on_a_quantised_footprint_goes_wrong() {
  s = disc();
  cx = wx(4, 4);
  cy = wy(4, 4);

  // The measurement that explains everything below: the innermost BOUNDARY cell
  // is closer to the centre than the outermost INTERIOR one.
  blo = 1000.0;
  ihi = -1000.0;
  for r in 0..H { for q in 0..W { if hexset_get(s, q, r) {
    dx = wx(q, r) - cx;
    dy = wy(q, r) - cy;
    d = sqrt(dx * dx + dy * dy) ?? 0.0;
    edge = false;
    for i in 0..6 { if !hexset_get(s, nb_q(q, r, i), nb_r(q, r, i)) { edge = true; } }
    if edge { if d < blo { blo = d; } } else { if d > ihi { ihi = d; } }
  } } }
  assert(ihi > blo, "interior reaches {ihi}, the nearest eave cell is at {blo} — they overlap");

  // Unclamped: drains perfectly, and the eave wobbles by more than two units.
  loose = heights_new(0, 0, W, H);
  roof_cone(s, loose, cx, cy, 5.0, 1.0, 100.0);
  assert(roof_ponds(s, loose) == 0, "it sheds water");
  assert(peak_of(s, loose) == 5.0, "it reaches its stated peak");
  assert(eave_spread(s, loose) > 2.0, "but the eave is {eave_spread(s, loose)} out of level");

  // Clamped low enough to level the eave: now the plateau ponds.
  tight = heights_new(0, 0, W, H);
  roof_cone(s, tight, cx, cy, 5.0, 1.0, 3.0);
  assert(eave_spread(s, tight) < 0.000001, "the eave is level");
  assert(roof_ponds(s, tight) == 3, "and three interior cells now hold water");

  // In between, both are wrong at once — there is no setting that works.
  mid = heights_new(0, 0, W, H);
  roof_cone(s, mid, cx, cy, 5.0, 1.0, 3.25);
  assert(eave_spread(s, mid) > 0.0 && roof_ponds(s, mid) > 0, "the crossover is not a solution");

  // Driven from the footprint instead: level eave AND no ponding, by construction.
  hip = heights_new(0, 0, W, H);
  rings = roof_hip(s, hip, 0.0, 1.0);
  assert(rings == 2, "the disc is {rings} rings deep");
  assert(eave_spread(s, hip) == 0.0, "every boundary cell is ring 0, so the eave is exactly level");
  assert(roof_ponds(s, hip) == 0, "and it still drains");
}

// @HXR-003 — THE CENTRE AND THE TRACK ARE IN WORLD UNITS, NOT CELL INDICES.
// A pointy-top odd-r lattice is not the integer grid its indices look like:
// cell (4,4) sits at world (6.93, 6.0).  Pass `(4.0, 4.0)` meaning "the middle
// cell" and the cone is centred outside the middle of the building.
//
// What makes this expensive is that the result is still a valid roof.  It has a
// positive peak, it is not flat, and it does not pond — so a suite that checks
// the apex passes it.  What moved is the apex's LOCATION and the eave's level.
fn test_the_centre_is_a_world_point_not_a_cell_index() {
  s = disc();
  assert(wx(4, 4) > 6.9 && wx(4, 4) < 6.94, "cell (4,4) is at world x {wx(4, 4)}, not 4.0");
  assert(wy(4, 4) == 6.0, "and world y {wy(4, 4)}, not 4.0");

  right = heights_new(0, 0, W, H);
  wrong = heights_new(0, 0, W, H);
  roof_cone(s, right, wx(4, 4), wy(4, 4), 5.0, 1.0, 100.0);
  roof_cone(s, wrong, 4.0, 4.0, 5.0, 1.0, 100.0);          // cell indices — the mistake

  // The apex is on the centre cell when the centre is a world point.
  assert(height_get(right, 4, 4) == 5.0, "the peak lands on the middle cell");
  assert(height_get(wrong, 4, 4) < 2.0, "and does not: {height_get(wrong, 4, 4)}");

  // Every cheap check still passes on the wrong one, which is why it survives.
  assert(roof_ponds(s, wrong) == 0, "the wrong roof still sheds water");
  assert(peak_of(s, wrong) > 0.0, "it still has a positive peak");

  // The eave is the tell — nearly four times worse than the merely-quantised one.
  assert(eave_spread(s, wrong) > 7.0, "the eave gives it away: {eave_spread(s, wrong)}");
}

// @HXR-004 — GROIN AND CLOISTER ARE THE SAME TWO BARRELS AND THE OPPOSITE
// OPERATOR.  `vault_groin` takes the higher of the two surfaces, so both crown
// lines stay at full height and every wall keeps its arch.  `vault_cloister`
// takes the lower, so the vault springs continuously off all four walls and no
// wall has an opening at all.
//
// Identical argument lists, one word different at the call site, and a different
// building.  If you have picked the wrong one, nothing about the call will say
// so — look at whether the walls have openings.
fn test_groin_and_cloister_differ_only_in_which_surface_wins() {
  s = disc();
  ew = track_new();
  track_straight(ew, wx(0, 4), wy(4, 4), wx(8, 4), wy(4, 4));
  ns = track_new();
  track_straight(ns, wx(4, 1), wy(4, 1), wx(4, 7), wy(4, 7));

  groin = heights_new(0, 0, W, H);
  cloister = heights_new(0, 0, W, H);
  vault_groin(s, groin, ew, ns, 0.0, 3.0);
  vault_cloister(s, cloister, ew, ns, 0.0, 3.0);

  // The higher surface never drops below the lower one, anywhere.
  same = 0;
  higher = 0;
  for r in 0..H { for q in 0..W { if hexset_get(s, q, r) {
    g = height_get(groin, q, r);
    c = height_get(cloister, q, r);
    if g >= c { same = same + 1; }
    if g > c { higher = higher + 1; }
  } } }
  assert(same == 29, "the groin is never lower, on all {same} cells");
  assert(higher == 27, "and strictly higher on {higher} of them — this is a different vault");
}

// @HXR-005 — PONDING COUNTS INTERIOR CELLS ONLY, AND THAT IS DELIBERATE.  A
// pond is water that cannot leave.  A cell on the boundary always can — it
// drains off the eave, which is its job — so a boundary local minimum is not a
// defect and is not counted.
//
// So `roof_ponds(s, f) == 0` does NOT mean "the field has no local minima".  If
// you are looking for a dip at the edge of a roof, this is not the function that
// finds it.
fn test_a_low_spot_on_the_eave_is_not_a_pond() {
  s = disc();

  // Sink a BOUNDARY cell far below everything around it: still no ponds.
  edged = heights_new(0, 0, W, H);
  roof_hip(s, edged, 0.0, 1.0);
  height_set(edged, 4, 1, -50.0);
  assert(roof_ponds(s, edged) == 0, "a hollow on the eave drains off it");

  // Sink an INTERIOR cell by the same amount: that is a pond.
  inner = heights_new(0, 0, W, H);
  roof_hip(s, inner, 0.0, 1.0);
  height_set(inner, 4, 4, -50.0);
  assert(roof_ponds(s, inner) == 1, "the same hollow inside is water that cannot leave");

  // The control that keeps the number honest: a dead-flat field is all ponds.
  assert(roof_ponds(s, heights_new(0, 0, W, H)) == 10, "a flat roof ponds everywhere inside");
}

// @HXR-006 — RECOVER THE SURFACE AND DRAW THAT, NEVER THE CELLS.  The heights
// are STORAGE; the analytic surface is the truth.  `roof_match` reads a height
// field back and names the profile, and `roof_eval` is then exact at any
// resolution — where drawing the cells directly gives a staircase of hexes
// wearing a cone.
//
// Two things to know about the answer.  A PLANE WINS TIES, because it is the
// simpler statement and a cone fitted to a flat roof puts its apex somewhere
// arbitrary.  And the fit is an R2 recovery with a TOLERANCE, by design — it
// recovers a profile from an arbitrary field, not from one you authored, so the
// recovered apex is not bit-identical to the one you asked for.
fn test_the_recovered_surface_is_the_thing_to_draw() {
  s = disc();

  // A field that IS a plane comes back as a plane.
  plane = heights_new(0, 0, W, H);
  for r in 0..H { for q in 0..W { if hexset_get(s, q, r) {
    height_set(plane, q, r, 2.0 + 0.5 * wx(q, r) + 0.25 * wy(q, r));
  } } }
  pf = roof_match(s, plane, 0.001);
  assert(pf.rf_kind == ROOF_PLANE, "a plane is recognised as a plane");
  assert(pf.rf_res < 0.000001, "with a residual of {pf.rf_res}");

  // A field that is a cone comes back as a cone.
  cone = heights_new(0, 0, W, H);
  roof_cone(s, cone, wx(4, 4), wy(4, 4), 5.0, 1.0, 100.0);
  cf = roof_match(s, cone, 0.001);
  assert(cf.rf_kind == ROOF_CONE, "and a cone as a cone");

  // `roof_eval` reads the ANALYTIC surface — close to the authored apex, and
  // deliberately not identical to it: this is a fit, not a round trip.
  z = roof_eval(cf, wx(4, 4), wy(4, 4));
  assert(z > 4.999 && z < 5.001, "the recovered apex is {z}, next to the authored 5.0");
}