hex_draw source

← hex_draw · API reference · version 0.1.4

5 file(s), 1251 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_draw.loft

// Copyright (c) 2026 Jurjen Stellingwerff
// SPDX-License-Identifier: LGPL-3.0-or-later
// hex_draw — draw a model into the field, and read its surface back.
//
// The `draw` map: a `Plan` becomes floor cells, wall edges, openings and a roof height field.
// And the inverse reading — a wall's ANALYTIC SURFACE recovered as the exact average of its
// stored edges, so a wall renders as one flat quad rather than a strip of them.
//
// Entry point: re-exports `housedraw`, `hexsurf`.
pub use housedraw::*;
pub use hexsurf::*;

src/hexsurf.loft

// Copyright (c) 2026 Jurjen Stellingwerff
// SPDX-License-Identifier: LGPL-3.0-or-later
// hexsurf — THE WALL SURFACE, BY AVERAGING.  Nothing here is fitted, so nothing has a tolerance.
//
// `ROUNDTRIP` §6.1: the straight surface a wall run approximates is DERIVED by averaging the edges
// that are already there, never fitted to them.
//
//   direction   the SUM of the edge vectors — an exact INTEGER vector in doubled coordinates
//   position    the mean of the edge MIDPOINTS.  A corner is an integer in (k,m), so a midpoint is
//               a half-integer and the mean is an exact RATIONAL — kept here as numerator and
//               denominator, never as a float
//   band        the exact perpendicular extent of the strip about that mean.  NOT an error term:
//               it is the geometry of the strip, and an exact algebraic constant per family
//
// This is what makes domain-A straight-wall recovery **R1**.  Least squares would introduce a
// residual to threshold; an average of exact rationals IS the answer.  `P4` therefore holds all
// the way through the presentation layer — there is no `ε` anywhere between `𝕋` and the drawn wall.
use hex_field::*;
use hex_grid::*;
use hex_form::*;

// ── the exact band constants (`ROUNDTRIP` §6.2), in world units ────────────
//
// 1 u = one hex edge = the circumradius; 1 hex step = sqrt(3) u = 1.5 m.  These are the CORNER
// bands — the full extent of the strip — which is what §6.2's table names.
pub const BAND_TOPS  = 0.5;                    // east family, dm = 0
pub const BAND_SIDES = 0.8660254037844386;     // north family, dk = 0 — sqrt(3)/2
pub const BAND_RATIO = 1.7320508075688772;     // sqrt(3), exactly

// The widening rule: both walls are presented straight and equally thick, at the LARGER band.
// Only the tops widen, and the amount is a closed form in Q(sqrt3) — nothing measured or tuned.
pub const WIDEN_TOTAL = 0.36602540378443865;   // (sqrt3 - 1)/2  u
pub const WIDEN_FACE  = 0.18301270189221933;   // (sqrt3 - 1)/4  u, applied to each face

// The averaged surface of one wall run.  Direction is exact integer; position is an exact
// rational held as `(num_k, num_m) / den`.
pub struct WallSurface {
  pub ws_dk:    integer,   // summed edge vector, doubled coordinates — exact
  pub ws_dm:    integer,
  pub ws_num_k: integer, // mean midpoint = (num_k / den, num_m / den) in doubled coordinates
  pub ws_num_m: integer,
  pub ws_den:   integer, // 2 * n
  pub ws_n:     integer, // edges averaged
}

// The two corners bounding the boundary edge between cell (qa,ra) and its outside neighbour.
fn edge_corners_of(qa: integer, ra: integer, qb: integer, rb: integer) -> (integer, integer) {
  hex_edge_corners(hex_neighbor_dir(qa, ra, qb, rb))
}

// THE DERIVATION.  Sum the edge vectors; sum the doubled midpoints.  Both are integer sums, so the
// surface is exact before any float is involved.
//
// Example: @HXD-001
pub fn surface_of(p: Plan, cells: HexSet, side: integer) -> WallSurface {
  sr = side_edges(p, cells, side);
  n = side_run_len(sr);
  sdk = 0;
  sdm = 0;
  smk = 0;
  smm = 0;
  for si in 0..n {
    qa = sr.sr_qa[si] ?? 0;
    ra = sr.sr_ra[si] ?? 0;
    (c1, c2) = edge_corners_of(qa, ra, sr.sr_qb[si] ?? 0, sr.sr_rb[si] ?? 0);
    sdk = sdk + corner_k(c2) - corner_k(c1);
    sdm = sdm + corner_m(c2) - corner_m(c1);
    smk = smk + 2 * lattice_k(qa, ra) + corner_k(c1) + corner_k(c2);
    smm = smm + 2 * lattice_m(qa, ra) + corner_m(c1) + corner_m(c2);
  }
  WallSurface { ws_dk: sdk, ws_dm: sdm, ws_num_k: smk, ws_num_m: smm,
                ws_den: 2 * n, ws_n: n }
}

// Which of the 12 HEADINGS is the mean direction an exact positive multiple of?  -1 if none.
//
// This is the exactness claim: the averaged direction equals the nominal side direction EXACTLY,
// as an integer multiple of a heading step.  Any other value would mean the average had drifted
// off the lattice, which is what a fit would do and an average must not.
//
// PARALLEL, not an integer multiple, and not "one component is zero".  Both of those tests
// reject geometry that is exact:
//
//   "one component is zero"   only 4 of the 12 headings have one — a side running NE is `(1,3)`,
//                             both components non-zero and perfectly exact.  Failed 16 of 24.
//   "an integer multiple"     the summed vector is `(0,14)` where the heading step is `(0,6)`;
//                             the SUM over a run is not a whole number of steps.  Failed 12 of 24.
//
// What §6.1 actually claims is that the DIRECTION is exact — no residual perpendicular component —
// which is a zero CROSS PRODUCT against the heading step, in exact integers.  The magnitude is the
// run's length and carries no claim.
//
// Example: @HXD-002
pub fn surface_heading(w: WallSurface) -> integer {
  for sh in 0..HEAD_N {
    st = head_step(sh);
    if w.ws_dk * st.1 - w.ws_dm * st.0 == 0 {
      if w.ws_dk * st.0 + w.ws_dm * st.1 > 0 { return sh; }
    }
  }
  -1
}

// Does this surface run along one of the exact headings?  An inexact one is still drawable
// but cannot be reasoned about on the lattice, so ask before treating it as integer.
pub fn surface_is_exact(w: WallSurface) -> boolean { surface_heading(w) >= 0 }

// World position of the mean, and the unit perpendicular to the mean direction.
pub fn surface_mean_x(w: WallSurface) -> float {
  (((w.ws_num_k as float) / (w.ws_den as float)) ?? 0.0) * 0.8660254037844386
}

// The surface's mean point in world y — the companion of `surface_mean_x`, and where a
// label or a normal is anchored.
pub fn surface_mean_y(w: WallSurface) -> float {
  (((w.ws_num_m as float) / (w.ws_den as float)) ?? 0.0) * 0.5
}

fn perp_x(w: WallSurface) -> float {
  wx = (w.ws_dk as float) * 0.8660254037844386;
  wy = (w.ws_dm as float) * 0.5;
  ln = sqrt(wx * wx + wy * wy) ?? 1.0;
  (0.0 - wy / ln) ?? 0.0
}

fn perp_y(w: WallSurface) -> float {
  wx = (w.ws_dk as float) * 0.8660254037844386;
  wy = (w.ws_dm as float) * 0.5;
  ln = sqrt(wx * wx + wy * wy) ?? 1.0;
  (wx / ln) ?? 0.0
}

// The perpendicular extent of the run's CORNERS about the mean line — the band of §6.2.
//
// Example: @HXD-001
pub fn surface_band_corners(p: Plan, cells: HexSet, side: integer) -> float {
  w = surface_of(p, cells, side);
  sr = side_edges(p, cells, side);
  px = perp_x(w);
  py = perp_y(w);
  cx = surface_mean_x(w);
  cy = surface_mean_y(w);
  lo = 1.0e9;
  hi = -1.0e9;
  for bi in 0..side_run_len(sr) {
    qa = sr.sr_qa[bi] ?? 0;
    ra = sr.sr_ra[bi] ?? 0;
    (c1, c2) = edge_corners_of(qa, ra, sr.sr_qb[bi] ?? 0, sr.sr_rb[bi] ?? 0);
    for wi in 0..2 {
      cc = c1;
      if wi == 1 { cc = c2; }
      ex = ((lattice_k(qa, ra) + corner_k(cc)) as float) * 0.8660254037844386;
      ey = ((lattice_m(qa, ra) + corner_m(cc)) as float) * 0.5;
      ev = (ex - cx) * px + (ey - cy) * py;
      if ev < lo { lo = ev; }
      if ev > hi { hi = ev; }
    }
  }
  hi - lo
}

// ── G2: THE WALL AS *ONE* QUAD ─────────────────────────────────────────────
//
// `G2` asks that every wall render as **one flat quad — its analytic surface — not the strip**.
// Today the shot draws the wall *as it is stored*: one strip per boundary edge, so a side that is
// nine edges long is nine quads that zigzag by the band `X47` measured.  The surface is already
// exact (direction exactly a heading, position an exact rational); what was missing is its EXTENT,
// so a renderer can emit the segment instead of the strip.
//
// The span is the run's corners **projected onto the mean line**, taking the extremes.  The two
// points returned lie ON the line by construction, so the quad the renderer builds from them has
// perpendicular spread exactly **0** — which is `G2`'s check — while the strip it replaces spreads
// by `surface_band_corners`.
//
// Example: @HXD-003
pub fn surface_span(p: Plan, cells: HexSet, side: integer) -> (float, float, float, float) {
  w = surface_of(p, cells, side);
  sr = side_edges(p, cells, side);
  dx = (w.ws_dk as float) * 0.8660254037844386;
  dy = (w.ws_dm as float) * 0.5;
  dl = sqrt(dx * dx + dy * dy) ?? 1.0;
  ux = (dx / dl) ?? 1.0;
  uy = (dy / dl) ?? 0.0;
  cx = surface_mean_x(w);
  cy = surface_mean_y(w);
  lo = 1.0e9;
  hi = -1.0e9;
  for si in 0..side_run_len(sr) {
    qa = sr.sr_qa[si] ?? 0;
    ra = sr.sr_ra[si] ?? 0;
    (c1, c2) = edge_corners_of(qa, ra, sr.sr_qb[si] ?? 0, sr.sr_rb[si] ?? 0);
    for sj in 0..2 {
      cc = c1;
      if sj == 1 { cc = c2; }
      ex = ((lattice_k(qa, ra) + corner_k(cc)) as float) * 0.8660254037844386;
      ey = ((lattice_m(qa, ra) + corner_m(cc)) as float) * 0.5;
      tv = (ex - cx) * ux + (ey - cy) * uy;
      if tv < lo { lo = tv; }
      if tv > hi { hi = tv; }
    }
  }
  (cx + lo * ux, cy + lo * uy, cx + hi * ux, cy + hi * uy)
}

// The perpendicular spread of the FITTED quad's own corners about its line — `G2`'s check.  Both
// endpoints come off the line, so this is the number `eave_spread(fitted) = 0` names.  Measured
// rather than asserted, so that if the span ever stops lying on the line the gate says so.
//
// ⚠ *EXACTLY* 0 IS A CLAIM ABOUT FLOATS, AND IT IS ORIENTATION-DEPENDENT.  The derivation above is
// exact — integer sums and rationals — but the span is projected through world coordinates, and on
// the orientations where that arithmetic does not cancel the round trip costs an ulp.  Measured
// over 12 orientations x 4 sides: 32 read exactly 0 and 16 read 2.22e-16.  So `== 0.0` passes only
// on the orientations it was tried on; what holds everywhere is the comparison that means
// something — flat against the STRIP's band, which is fifteen orders of magnitude larger.
//
// Example: @HXD-001
pub fn surface_fitted_spread(p: Plan, cells: HexSet, side: integer) -> float {
  w = surface_of(p, cells, side);
  (sx0, sy0, sx1, sy1) = surface_span(p, cells, side);
  px = perp_x(w);
  py = perp_y(w);
  cx = surface_mean_x(w);
  cy = surface_mean_y(w);
  e0 = (sx0 - cx) * px + (sy0 - cy) * py;
  e1 = (sx1 - cx) * px + (sy1 - cy) * py;
  lo = e0;
  hi = e0;
  if e1 < lo { lo = e1; }
  if e1 > hi { hi = e1; }
  hi - lo
}

// THE MITER POINT (`I-CORNER` part 4): where two adjacent fitted surfaces intersect.  This is the
// corner the RENDER should use — each quad running miter-to-miter, so consecutive sides share an
// endpoint exactly and the outline closes with no gap.  That closure is the invariant, and it is
// exact.
//
// ⚠ IT IS NOT THE `Plan`'S CORNER, and it never can be.  An earlier note here claimed the
// intersection *"must land on the exact model corner, or the surfaces have drifted off the shape
// they were fitted to"* — a check nobody ran.  Measured, it misses by 0.47 world units on the 5x4
// cottage, and swept over plan sizes the miss is neither constant (0.465 .. 0.827) nor even
// consistently outward.  The reason is structural rather than numerical: what these surfaces
// recover is the face of the wall AS DRAWN, and the drawn footprint is the plan quantised to
// cells, so the recovered corner is an exact rational of the lattice while the plan's half-depth
// is an irrational multiple of `sqrt(3)/2`.  Two number systems that meet only at the origin.
//
// So a caller must never diff a recovered outline against plan geometry, and must not try to
// correct for the offset either.  Use the recovered corner for anything that has to line up with
// what was drawn.
//
// `ok` is false when the two lines are parallel, which for adjacent sides of a closed form cannot
// happen — so a false here is a real defect, not a degenerate input to be shrugged at.
//
// Example: @HXD-003
pub fn surface_miter(p: Plan, cells: HexSet, sa: integer, sb: integer)
                     -> (float, float, boolean) {
  wa = surface_of(p, cells, sa);
  wb = surface_of(p, cells, sb);
  ax = surface_mean_x(wa);
  ay = surface_mean_y(wa);
  bx = surface_mean_x(wb);
  by = surface_mean_y(wb);
  aux = (wa.ws_dk as float) * 0.8660254037844386;
  auy = (wa.ws_dm as float) * 0.5;
  bux = (wb.ws_dk as float) * 0.8660254037844386;
  buy = (wb.ws_dm as float) * 0.5;
  den = aux * buy - auy * bux;
  if den < 0.000000001 && den > 0.0 - 0.000000001 { return (0.0, 0.0, false); }
  tt = (((bx - ax) * buy - (by - ay) * bux) / den) ?? 0.0;
  (ax + tt * aux, ay + tt * auy, true)
}

// THE MITERED QUAD (`G2` + `I-CORNER` part 4): the side's segment running MITER TO MITER, so
// consecutive sides share an endpoint exactly and the outline closes with no corner gap.  This is
// what a renderer should draw; `surface_span` (the projected extent) leaves the gaps the fitted
// render made visible.
//
// Example: @HXD-003
pub fn surface_quad(p: Plan, cells: HexSet, side: integer, nsides: integer)
                    -> (float, float, float, float) {
  // `%` is FALLIBLE in loft, so `nsides = 0` types `prev`/`next` as `integer?` — and a null index
  // arriving at `surface_miter` would read as a valid-looking side **0** rather than a refusal: a
  // silently wrong corner instead of a caught one.  A miter needs `prev != side != next`, so it
  // needs at least three sides; below that there is no corner to cut and the honest answer is the
  // un-mitered span, which is already this function's defined escape when a miter fails.
  //
  // Every caller passes 4, so this cannot fire today — but *"this cannot fire"* is a claim like
  // any other (`X57`: `wall_separates`' comment said exactly that, was wrong, and was
  // load-bearing under a green suite).  The `?? 0` below is unreachable PAST this guard; the
  // guard is what refuses, the `??` only discharges the type.
  if nsides < 3 { return surface_span(p, cells, side); }
  prev = ((side - 1) % nsides + nsides) % nsides ?? 0;
  next = (side + 1) % nsides ?? 0;
  (ax, ay, aok) = surface_miter(p, cells, prev, side);
  (bx, by, bok) = surface_miter(p, cells, side, next);
  if !aok || !bok { return surface_span(p, cells, side); }
  (ax, ay, bx, by)
}

// One stored edge's INTERVAL along the fitted line — where edge `i` of the run projects onto the
// surface.  This is what makes a feature renderable as `FEATURES` describes it: *"a door is
// `[s0,s1]` along an exact straight or arc"*, an interval on the analytic surface rather than a
// stored edge.  The renderer draws the wall as one quad (`surface_span`) and each feature as a
// sub-interval of it, so the zigzag never reaches the picture.
//
// Example: @HXD-006
pub fn surface_edge_span(p: Plan, cells: HexSet, side: integer, i: integer)
                         -> (float, float, float, float) {
  w = surface_of(p, cells, side);
  sr = side_edges(p, cells, side);
  dx = (w.ws_dk as float) * 0.8660254037844386;
  dy = (w.ws_dm as float) * 0.5;
  dl = sqrt(dx * dx + dy * dy) ?? 1.0;
  ux = (dx / dl) ?? 1.0;
  uy = (dy / dl) ?? 0.0;
  cx = surface_mean_x(w);
  cy = surface_mean_y(w);
  qa = sr.sr_qa[i] ?? 0;
  ra = sr.sr_ra[i] ?? 0;
  (c1, c2) = edge_corners_of(qa, ra, sr.sr_qb[i] ?? 0, sr.sr_rb[i] ?? 0);
  lo = 1.0e9;
  hi = -1.0e9;
  for ej in 0..2 {
    cc = c1;
    if ej == 1 { cc = c2; }
    ex = ((lattice_k(qa, ra) + corner_k(cc)) as float) * 0.8660254037844386;
    ey = ((lattice_m(qa, ra) + corner_m(cc)) as float) * 0.5;
    tv = (ex - cx) * ux + (ey - cy) * uy;
    if tv < lo { lo = tv; }
    if tv > hi { hi = tv; }
  }
  (cx + lo * ux, cy + lo * uy, cx + hi * ux, cy + hi * uy)
}

// The material on stored edge `i` of a side — so the renderer can find the features on a run
// without re-deriving the run's order.
pub fn surface_edge_mat(p: Plan, cells: HexSet, edges: EdgeSet, side: integer,
                        i: integer) -> integer {
  sr = side_edges(p, cells, side);
  edge_mat(edges, sr.sr_qa[i] ?? 0, sr.sr_ra[i] ?? 0, sr.sr_qb[i] ?? 0, sr.sr_rb[i] ?? 0)
}

// The perpendicular extent of the MIDPOINTS about the mean line.  Not in §6.2, and it is not the
// same as the corner band: on the east family the midpoints are exactly COLINEAR (0), while on the
// north family they span the full band — the row stagger, which the corners hide.
//
// Example: @HXD-001
pub fn surface_band_mids(p: Plan, cells: HexSet, side: integer) -> float {
  w = surface_of(p, cells, side);
  sr = side_edges(p, cells, side);
  px = perp_x(w);
  py = perp_y(w);
  cx = surface_mean_x(w);
  cy = surface_mean_y(w);
  lo = 1.0e9;
  hi = -1.0e9;
  for mi in 0..side_run_len(sr) {
    qa = sr.sr_qa[mi] ?? 0;
    ra = sr.sr_ra[mi] ?? 0;
    (c1, c2) = edge_corners_of(qa, ra, sr.sr_qb[mi] ?? 0, sr.sr_rb[mi] ?? 0);
    mx = ((lattice_k(qa, ra) as float)
          + ((corner_k(c1) + corner_k(c2)) as float) / 2.0) * 0.8660254037844386;
    my = ((lattice_m(qa, ra) as float)
          + ((corner_m(c1) + corner_m(c2)) as float) / 2.0) * 0.5;
    mv = (mx - cx) * px + (my - cy) * py;
    if mv < lo { lo = mv; }
    if mv > hi { hi = mv; }
  }
  hi - lo
}

// The scatter a LEAST-SQUARES formulation would have to threshold: the sum of squared
// perpendicular deviations of the midpoints from the derived mean line.  Averaging names the same
// geometry as an exact band constant and has nothing to tune; this is the control that the
// distinction is real rather than rhetorical.
//
// Example: @HXD-001
pub fn surface_lsq_residual(p: Plan, cells: HexSet, side: integer) -> float {
  w = surface_of(p, cells, side);
  sr = side_edges(p, cells, side);
  px = perp_x(w);
  py = perp_y(w);
  cx = surface_mean_x(w);
  cy = surface_mean_y(w);
  acc = 0.0;
  for ri in 0..side_run_len(sr) {
    qa = sr.sr_qa[ri] ?? 0;
    ra = sr.sr_ra[ri] ?? 0;
    (c1, c2) = edge_corners_of(qa, ra, sr.sr_qb[ri] ?? 0, sr.sr_rb[ri] ?? 0);
    mx = ((lattice_k(qa, ra) as float)
          + ((corner_k(c1) + corner_k(c2)) as float) / 2.0) * 0.8660254037844386;
    my = ((lattice_m(qa, ra) as float)
          + ((corner_m(c1) + corner_m(c2)) as float) / 2.0) * 0.5;
    dv = (mx - cx) * px + (my - cy) * py;
    acc = acc + dv * dv;
  }
  acc
}

// ── A5 — FEATURES ON A STRAIGHT SIDE, AND WHAT OF THEM SURVIVES ────────────
//
// `OD-9` is closed: a door or window IS the material on the wall slot, not an annotation beside
// one (`ROUNDTRIP` §2.4.1).  The edge is never removed, so a run cannot fragment — the doored-tower
// defect is structurally impossible rather than merely avoided.  A5 asks the next question: does
// the feature survive the ROUND TRIP?
//
// THE MODEL STORES `(side, t)` (`I2`, affine-invariant).  THE FIELD STORES A MATERIAL ON EDGES.
// So recovery has to get `t` back out of a set of marked edges — and it cannot get an arbitrary
// `t` back, because every `t` in an interval of width `1/n` selects the same edge.  What IS exact:
//
//     the edge centres along a run of n edges sit at   t = (2i + 1) / 2n
//
// exact rationals, odd numerator over `2n`.  `place_opening` snaps to the nearest, so a feature
// authored at one of those `t` values round-trips exactly and one authored between them does not.
// That is the same shape as the line-endpoint rule (§10.10): the doorstep quantises the model to
// what the storage can distinguish, rather than letting a value in that cannot come back.
//
// The run is NOT stored in `t` order — `side_edges` fills it in cell-scan order — so a feature is
// indexed by its exact `t` numerator, never by its position in the run.

// The exact `t` of stored edge `i`, as the numerator over `2n`: an odd integer in `1 .. 2n-1`.
// `sr_t` is a float only because `SideRun` holds it that way; the value is an exact rational and
// this recovers the numerator without a tolerance being involved in the answer.
//
// Example: @HXD-006
pub fn feature_t_num(sr: SideRun, i: integer, n: integer) -> integer {
  (round((sr.sr_t[i] ?? 0.0) * 2.0 * (n as float)) ?? 0.0) as integer
}

// Is `num/2n` a legal feature position — an edge centre, i.e. an ODD numerator in range?
pub fn feature_t_admissible(num: integer, n: integer) -> boolean {
  num > 0 && num < 2 * n && ((num % 2) + 2) % 2 == 1
}

// Recover a feature from the field: the `t` numerators of every edge on `side` carrying `kind`,
// in ascending `t` order (C2's ordering, not the storage order).
//
// Example: @HXD-006
pub fn feature_find(p: Plan, cells: HexSet, edges: EdgeSet,
                    side: integer, kind: integer) -> vector<integer> {
  sr = side_edges(p, cells, side);
  n = side_run_len(sr);
  out: vector<integer> = [];
  for fi in 0..n {
    if edge_mat(edges, sr.sr_qa[fi] ?? 0, sr.sr_ra[fi] ?? 0,
                       sr.sr_qb[fi] ?? 0, sr.sr_rb[fi] ?? 0) == kind {
      out += [feature_t_num(sr, fi, n)];
    }
  }
  // ascending t — insertion sort, the run is tens of edges at most
  for si in 1..len(out) {
    ps = si;
    while ps > 0 && (out[ps - 1] ?? 0) > (out[ps] ?? 0) {
      tmp = out[ps - 1] ?? 0;
      out[ps - 1] = out[ps] ?? 0;
      out[ps] = tmp;
      ps = ps - 1;
    }
  }
  out
}

// How many boundary edges carry a wall material at all — the run's length, which a re-materialled
// feature must leave UNCHANGED.  Deleting an edge instead is what fragments a run.
//
// Example: @HXD-005
pub fn wall_edge_count(cells: HexSet, edges: EdgeSet) -> integer {
  n = 0;
  for wr in hexset_r0(cells)..hexset_r0(cells) + hexset_h(cells) {
    for wq in hexset_q0(cells)..hexset_q0(cells) + hexset_w(cells) {
      if hexset_get(cells, wq, wr) {
        // ALL SIX directions, as `draw_walls` does.  A boundary edge is visited exactly once
        // because only the inside cell is in the set — iterating just the three OWNED directions
        // misses every edge the outside cell owns and undercounts by about half (38 read as 19).
        for wd in 0..6 {
          (nq, nr) = hex_neighbor(wq, wr, wd);
          if !hexset_get(cells, nq, nr) {
            if edge_mat(edges, wq, wr, nq, nr) != 0 { n = n + 1; }
          }
        }
      }
    }
  }
  n
}

src/housedraw.loft

// Copyright (c) 2026 Jurjen Stellingwerff
// SPDX-License-Identifier: LGPL-3.0-or-later
// housedraw — DRAW A BUILDING IN ANY OF THE 12 ORIENTATIONS.
//
// The field primitives are hex_field's (HexSet/Labels/EdgeSet/Heights) and the roof
// profiles hex_roof's.  This is the layer above them: a LOCAL FRAME (hex_form's `Plan`) so
// the same building can be drawn at 6 rotations x 2 chiralities, and openings placed on its
// walls by (side, t) rather than by cell coordinates.
//
// FOUR ROUTINES, and each writes exactly one layer:
//
//   draw_floor    the massing  -> HexSet + Labels     (the floor you stand on)
//   draw_walls    its boundary -> EdgeSet             (thin walls: no floor is lost)
//   place_opening n edges at t -> EdgeSet material    (doors and windows)
//   draw_roof     a ridge      -> Heights             (a thin caller of hexroof)
//
// THE TWO RULES THAT MAKE 12 ORIENTATIONS WORK:
//
//  1. THE MIRROR IS APPLIED ONCE, IN THE FRAME.  Local -> world negates local u when
//     mirrored; a side's endpoints then only EXCHANGE.  Applying it in both places turns
//     every outward normal INWARD (measured: -1.000 at all 12 orientations).
//
//  2. A WALL IS THE BOUNDARY OF A FILLED REGION, never a band of cells buffered along a
//     line.  The boundary of a region is a closed loop BY CONSTRUCTION at any size,
//     orientation or tie; the buffer rule leaks (measured: 2 components, 0 cells enclosed
//     — a wall with a hole — while still passing an equivariance check 12/12).
//
// A DOOR IS A MATERIAL ON ITS EDGES, AND NEVER DELETES THEM.  Deleting fragments the wall
// run, which is the defect DRAWING.md names as a doored tower fitting 3 arcs instead of 1.
// The foxel stores ONE material per wall slot (ROUNDTRIP.md §2.4), so an opening REPLACES
// the wall material rather than sitting in a channel beside it.  Composition therefore
// belongs in the material vocabulary — "a door in a cottage wall" is a material — which is
// what the consumer reads to decide passability.
use hex_field::*;
use hex_form::*;
use hex_way::*;
use hex_roof::*;
use hex_grid::*;

const EPS = 0.000000001;

// Opening materials — values in the SAME vocabulary as the wall materials, because they
// share the one slot.  A door overrides `solid`; a window does not (it overrides opacity,
// which is an elevation property and never touches the footprint).
//
// ⚠ `OPEN_NONE` IS NOT "no opening" — it is the wall-material value **0, which means NO WALL**.
// Writing it DELETES the edge and breaks the run in two (`X51`: 37 edges, 2 dangling ends, and on
// a tower 3 arcs instead of 1).  It is kept because 0 is the slot's empty value, but it is never
// what an author means by an opening.
//
// THE TAXONOMY.  A door is not a gap in the wall; it is something that can be a collider or
// not.  A real gap is allowed too — the door itself is missing, but the wall continues as normal.
// Three distinct things live on a wall edge, and **all three keep the wall one continuous run**:
//
//   OPEN_DOOR    an opening fitted with a door leaf — a collider WHEN CLOSED, not when open
//   OPEN_WINDOW  an opening at elevation — never a floor-level collider, wall continues
//   OPEN_GAP     a real gap: the leaf is missing entirely, never a collider, wall continues
//
// The difference between `OPEN_GAP` and `OPEN_NONE` is the whole point: a gap is an **opening in a
// wall that is still there**; `OPEN_NONE` is the **absence of the wall**.
//
// THE RATIONALE, and it is what makes this not a matter of taste: a gap or a door is never
// rendered as a missing wall — there is a door frame or a ragged stone opening there.  An
// opening **has geometry of its own** — a frame, a lintel, a ragged
// jamb.  Absence has none.  So an opening must name a material that carries that geometry, which
// is precisely a wall id into the palette; writing 0 would ask the renderer to draw nothing and the
// collision layer to find no wall, both wrong.
pub const OPEN_NONE   = 0;
pub const OPEN_DOOR   = 1;
pub const OPEN_WINDOW = 2;
pub const OPEN_GAP    = 3;


// ── 1. THE FLOOR ────────────────────────────────────────────────────────────
//
// The massing is filled FIRST and walled second.  Comparisons are inclusive (<= bound
// + eps) so the fill never loses a column of cells to the structural tie: along u every
// cell-to-side distance is an exact multiple of sqrt(3)/2, which is exactly the default
// halfwidth, so a strict comparison sits on the last bit of a float for a whole column.
pub fn draw_floor(p: Plan, cells: HexSet, types: Labels, floor_type: integer) -> integer {
  hw = plan_hw(p);
  hd = plan_hd(p);
  n = 0;
  for fr in hexset_r0(cells)..hexset_r0(cells) + hexset_h(cells) {
    for fq in hexset_q0(cells)..hexset_q0(cells) + hexset_w(cells) {
      (lu, lv) = plan_to_local(p, cell_x(fq, fr), cell_y(fq, fr));
      if fabs(lu) <= hw + EPS && fabs(lv) <= hd + EPS {
        hexset_set(cells, fq, fr, true);
        label_set(types, fq, fr, floor_type);
        n = n + 1;
      }
    }
  }
  n
}

// ── 2. THE WALLS ────────────────────────────────────────────────────────────
//
// Every edge between a cell of the region and one outside it.  Closed by construction,
// and it costs no floor at all: a 7.5 x 6.0 m cottage keeps all 27 of its cells, where
// a wall drawn as a band of cells would eat 16 of them — 59% of the house.
//
// Example: @HXD-005
pub fn draw_walls(cells: HexSet, edges: EdgeSet, wall_mat: integer) -> integer {
  n = 0;
  for wr in hexset_r0(cells)..hexset_r0(cells) + hexset_h(cells) {
    for wq in hexset_q0(cells)..hexset_q0(cells) + hexset_w(cells) {
      if hexset_get(cells, wq, wr) {
        for d in 0..6 {
          nq = nb_q(wq, wr, d);
          nr = nb_r(wq, wr, d);
          if !hexset_get(cells, nq, nr) {
            edge_set_mat(edges, wq, wr, nq, nr, wall_mat);
            n = n + 1;
          }
        }
      }
    }
  }
  n
}


// ── 3. OPENINGS ─────────────────────────────────────────────────────────────
//
// A door is N EDGES centred at t, not a t-interval erased with a segment: erasing
// over-erases by ~halfwidth at EACH end (a 3.00 m interval opens 4.50 m of wall) and
// shrinking the segment over-corrects.  One hex edge is 1.0 world unit = 0.87 m of
// clear width — a real domestic door, which a cell-thick wall could not express at all.
//
// The opening material REPLACES the wall material on those edges (the slot holds one).
// With any kind but `OPEN_NONE` nothing is deleted, so the wall still fits as one run.
// Returns how many edges it wrote: fewer than `nedges` when the side is shorter.
//
// Example: @HXD-004
pub fn place_opening(p: Plan, cells: HexSet, edges: EdgeSet,
                     side: integer, t: float, nedges: integer, kind: integer) -> integer {
  sr = side_edges(p, cells, side);
  total = len(sr.sr_t);
  if total == 0 { return 0; }
  taken: vector<boolean> = [for _ in 0..total { false }];
  placed = 0;
  for _pick in 0..nedges {
    best = -1;
    bestd = 0.0;
    for i in 0..total {
      if !(taken[i] ?? true) {
        d = fabs((sr.sr_t[i] ?? 0.0) - t);
        if best < 0 || d < bestd { best = i; bestd = d; }
      }
    }
    if best >= 0 {
      taken[best] = true;
      edge_set_mat(edges, sr.sr_qa[best] ?? 0, sr.sr_ra[best] ?? 0,
                          sr.sr_qb[best] ?? 0, sr.sr_rb[best] ?? 0, kind);
      placed = placed + 1;
    }
  }
  placed
}

// ── 4. THE ROOF ─────────────────────────────────────────────────────────────
//
// One ring of halo around the massing, because that ring is what makes the EAVE
// interpolable: without it 38 of the 76 wall-edge samples face a cell with no stored
// height, and the wall top has nothing to meet.
//
// Example: @HXD-007
pub fn grow_ring(src: HexSet, dst: HexSet) {
  for gr in hexset_r0(src)..hexset_r0(src) + hexset_h(src) {
    for gq in hexset_q0(src)..hexset_q0(src) + hexset_w(src) {
      if hexset_get(src, gq, gr) {
        hexset_set(dst, gq, gr, true);
        for gd in 0..6 { hexset_set(dst, nb_q(gq, gr, gd), nb_r(gq, gr, gd), true); }
      }
    }
  }
}

// A gable: the ridge runs along local u THROUGH the anchor and PAST both gable ends.
// Stopping it at the wall makes the distance go radial past the tip, so the gable end
// rolls over like a cone quadrant — 0.43 units, 0.37 m of bulge.  A hip is the same
// call with the ridge shortened; then the conical ends are the feature, not a defect.
//
// Example: @HXD-007
pub fn draw_roof(p: Plan, cells: HexSet, hz: Heights,
                 eave: float, pitch: float, hip_steps: integer) {
  hw = plan_hw(p);
  hd = plan_hd(p);
  half = hw + hd + HEX_LEN;                       // gable: run well past the ends
  if hip_steps > 0 { half = hw - (hip_steps as float) * HEX_LEN * 0.5; }
  (x0, y0) = plan_to_world(p, -half, 0.0);
  (x1, y1) = plan_to_world(p, half, 0.0);
  tr = track_new();
  track_straight(tr, x0, y0, x1, y1);
  roof_ridge(cells, hz, tr, eave + pitch * hd, pitch, 1000000.0);
}

tests/01-hex-draw.loft

// Copyright (c) 2026 Jurjen Stellingwerff
// SPDX-License-Identifier: LGPL-3.0-or-later
// hex_draw — the conformance tests: a Plan draws into the field, and its wall surface reads back.
//
// Contract only; the exhaustive measurement is hexbody's `tests/house.loft` and `tests/surface.loft`.
use hex_draw::*;
use hex_form::*;
use hex_field::*;

fn plan5x4() -> Plan { plan_new(0, 0, 5, 4, 0, false) }

fn main() {
  p = plan5x4();
  cells = hexset_chunk(-20, -20, 41, 41);
  types = labels_new(-20, -20, 41, 41);
  n = draw_floor(p, cells, types, 1);
  assert(n == 27, "a 5x4 house has 27 floor cells, got {n}");
}

fn test_the_walls_are_one_connected_boundary() {
  p = plan5x4();
  cells = hexset_chunk(-20, -20, 41, 41);
  types = labels_new(-20, -20, 41, 41);
  edges = edgeset_new(-20, -20, 41, 41);
  draw_floor(p, cells, types, 1);
  ne = draw_walls(cells, edges, 7);
  assert(ne == 38, "the 5x4 house has 38 wall edges, got {ne}");
}

// A consumer reads a surface's exact rational and builds one as a control (hexbody's surface
// gate): every field of `WallSurface` is readable and buildable outside this library.
fn test_a_surface_is_its_exact_rational() {
  p = plan5x4();
  cells = hexset_chunk(-20, -20, 41, 41);
  types = labels_new(-20, -20, 41, 41);
  draw_floor(p, cells, types, 1);
  w = surface_of(p, cells, 0);
  assert(w.ws_n > 0 && w.ws_den == 2 * w.ws_n, "den is twice the edges averaged: {w.ws_den} vs {w.ws_n}");
  one = WallSurface { ws_dk: 2, ws_dm: 0, ws_num_k: 4, ws_num_m: 0, ws_den: 2, ws_n: 1 };
  assert(surface_mean_x(one) > 1.73 && surface_mean_x(one) < 1.74, "mean x of (4/2) is 2 * 0.866, got {surface_mean_x(one)}");
}

fn test_a_wall_recovers_as_one_exact_surface() {
  // `G2`: a wall is not a strip of edges, it is one analytic surface -- and its direction is
  // EXACTLY a heading, not a fitted approximation.  This is what a game emits as a collision proxy.
  p = plan5x4();
  cells = hexset_chunk(-20, -20, 41, 41);
  types = labels_new(-20, -20, 41, 41);
  draw_floor(p, cells, types, 1);
  w = surface_of(p, cells, 0);
  assert(surface_is_exact(w), "side 0's surface is exact");
  spread = surface_fitted_spread(p, cells, 0);
  // Exactly 0 at this orientation, and that is luck: the derivation is exact but the span is
  // presented in floats, so 16 of the 12 orientations x 4 sides read one ulp instead (the
  // worked examples measure the census).  What holds everywhere is flat against the STRIP's
  // own width, which is fifteen orders of magnitude larger.
  assert(spread < BAND_TOPS * 1e-14,
         "the fitted surface is flat where the strip it replaces is {BAND_TOPS} wide, got {spread}");
}

tests/02-worked-examples.loft

// Copyright (c) 2026 Jurjen Stellingwerff
// SPDX-License-Identifier: LGPL-3.0-or-later
//
// Worked examples for `hex_draw` (@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: @HXD-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-draw.loft`.
//
// WHAT THIS PACKAGE ASKS OF A CALLER.  It draws a model into a field of cells
// and edges, and it reads a wall back out of that field as an analytic surface.
// Both halves live in TWO number systems at once, and every trap here is one of
// them leaking into the other:
//
//   the field    exact integers — a cell is (q,r), a corner is (k,m), an edge is
//                a pair of cells.  Everything recovered from it is an exact
//                rational, and nothing needs a tolerance
//   the plan     a rectangle in world units, whose half-extents are irrational
//                multiples of sqrt(3)/2 and land wherever they land on the lattice
//
// So a wall's DIRECTION is exact and its POSITION is a rational the plan cannot
// name; a door is an integer material and never a hole; a feature's `t` is an
// exact rational and never a float you chose.  Each of those is worked below.
use hex_draw::*;
use hex_field::*;
use hex_form::*;
use hex_grid::*;
use hex_shape::*;
use hex_way::*;
use hex_roof::*;

const W  = 41;
const H  = 41;
const Q0 = 0 - 20;
const R0 = 0 - 20;
const WALLMAT = 7;

// The cottage every example is drawn from: 5 hex steps by 4, at the origin, unrotated.
// 27 floor cells and 38 wall edges — 7.5 m by 6.0 m of real building.
fn cottage() -> Plan { plan_new(0, 0, 5, 4, 0, false) }

fn patch() -> HexSet { hexset_chunk(Q0, R0, W, H) }

// The cottage's floor, filled.
fn cottage_cells() -> HexSet {
  cells = patch();
  lb = labels_new(Q0, R0, W, H);
  draw_floor(cottage(), cells, lb, 1);
  cells
}

// How many floor cells touch the outside — the ring a wall drawn as a BAND OF CELLS
// would have to consume.
fn ring_cells(cells: HexSet) -> integer {
  n = 0;
  for r in R0..R0 + H {
    for q in Q0..Q0 + W {
      if hexset_get(cells, q, r) {
        out = false;
        for d in 0..6 { if !hexset_get(cells, nb_q(q, r, d), nb_r(q, r, d)) { out = true; } }
        if out { n = n + 1; }
      }
    }
  }
  n
}

// The wall-edge count written the tempting way: only the three directions a cell
// "owns", so every edge the OUTSIDE cell owns is missed.  Parameterised by which of
// the two neighbour conventions is used, because they disagree about which three.
fn count_three_dirs(cells: HexSet, edges: EdgeSet, hexgrid_order: boolean) -> integer {
  n = 0;
  for r in R0..R0 + H {
    for q in Q0..Q0 + W {
      if hexset_get(cells, q, r) {
        for d in 0..3 {
          nq = nb_q(q, r, d);
          nr = nb_r(q, r, d);
          if hexgrid_order { (nq, nr) = hex_neighbor(q, r, d); }
          if !hexset_get(cells, nq, nr) && edge_mat(edges, q, r, nq, nr) != 0 { n = n + 1; }
        }
      }
    }
  }
  n
}

// The distance between the end of side `a`'s drawn segment and the start of side `b`'s.
fn span_gap(p: Plan, cells: HexSet, a: integer, b: integer) -> float {
  (_x0, _y0, x1, y1) = surface_span(p, cells, a);
  (u0, v0, _u1, _v1) = surface_span(p, cells, b);
  sqrt((x1 - u0) * (x1 - u0) + (y1 - v0) * (y1 - v0)) ?? 0.0
}

fn quad_gap(p: Plan, cells: HexSet, a: integer, b: integer) -> float {
  (_x0, _y0, x1, y1) = surface_quad(p, cells, a, 4);
  (u0, v0, _u1, _v1) = surface_quad(p, cells, b, 4);
  sqrt((x1 - u0) * (x1 - u0) + (y1 - v0) * (y1 - v0)) ?? 0.0
}

// A ridge along the plan's local u, extending `half` either side of the anchor.
fn ridge_heights(p: Plan, cells: HexSet, half: float) -> Heights {
  hz = heights_new(Q0, R0, W, H);
  (x0, y0) = plan_to_world(p, 0.0 - half, 0.0);
  (x1, y1) = plan_to_world(p, half, 0.0);
  tr = track_new();
  track_straight(tr, x0, y0, x1, y1);
  roof_ridge(cells, hz, tr, 2.0 + 0.5 * plan_hd(p), 0.5, 1000000.0);
  hz
}

// The smallest call that does something real: draw the cottage, and read one of its
// walls back as a segment a renderer can draw.
fn main() {
  p = cottage();
  cells = cottage_cells();
  edges = edgeset_new(Q0, R0, W, H);
  ne = draw_walls(cells, edges, WALLMAT);
  assert(ne == 38, "the cottage has 38 wall edges, got {ne}");
  w = surface_of(p, cells, 0);
  assert(surface_is_exact(w), "and side 0 recovers as an exact surface");
  (x0, y0, x1, y1) = surface_quad(p, cells, 0, 4);
  len2 = sqrt((x1 - x0) * (x1 - x0) + (y1 - y0) * (y1 - y0)) ?? 0.0;
  assert(len2 > 0.0, "which a renderer draws as one segment, {len2} world units long");
}

// @HXD-001 — A WALL IS ONE FLAT QUAD, NOT THE STRIP IT IS STORED AS.  The field
// stores a wall as N separate boundary edges, and those edges ZIGZAG: they are
// the sides of hexagons, so consecutive ones step sideways by an exact amount.
// Draw them as they are stored and a nine-edge wall is nine quads with a visible
// comb; draw the analytic surface and it is one.
//
// The surface is DERIVED by averaging, never fitted — the direction is the sum of
// the edge vectors and the position the mean of their midpoints, both exact
// integer sums — so there is no residual anywhere to threshold.  What the two
// bands below measure is not error: it is the GEOMETRY of the strip, and it is an
// exact algebraic constant per family, `sqrt(3)` apart.
fn test_a_wall_is_one_quad_and_the_strip_it_replaces_has_a_width() {
  p = cottage();
  cells = cottage_cells();

  // The two families, and their exact bands.
  b0 = surface_band_corners(p, cells, 0);
  b1 = surface_band_corners(p, cells, 1);
  assert(abs(b0 - BAND_SIDES) < 1e-12, "a north-family side's strip is sqrt(3)/2 wide, got {b0}");
  assert(abs(b1 - BAND_TOPS) < 1e-12, "an east-family side's is 1/2, got {b1}");
  assert(abs((b0 / b1 ?? 0.0) - BAND_RATIO) < 1e-12, "and the ratio is exactly sqrt(3)");

  // The quad that replaces the strip has no width at all: both its endpoints come
  // off the derived line, so its perpendicular spread is zero.
  assert(surface_fitted_spread(p, cells, 0) == 0.0, "the fitted quad is flat");
  assert(surface_fitted_spread(p, cells, 1) == 0.0, "on both families");

  // …to within the float it is presented in, which is NOT the same as "exactly 0"
  // and is worth knowing before you write `== 0.0` into a gate.  The derivation is
  // exact (integers and rationals); the world coordinates are f64, and projecting
  // a point onto a line and back costs an ulp on the orientations where the
  // arithmetic does not cancel.  Over the 12 orientations x 4 sides, 32 read
  // exactly 0 and 16 read 2.22e-16 — fifteen orders of magnitude below the band
  // they are being distinguished from, which is the comparison that means something.
  worst = 0.0;
  nexact = 0;
  ninexact = 0;
  for rot in 0..6 {
    for mi in 0..2 {
      q = plan_new(0, 0, 5, 4, rot, mi == 1);
      c = patch();
      lb = labels_new(Q0, R0, W, H);
      draw_floor(q, c, lb, 1);
      for side in 0..4 {
        sp = surface_fitted_spread(q, c, side);
        if sp == 0.0 { nexact = nexact + 1; } else { ninexact = ninexact + 1; }
        if sp > worst { worst = sp; }
      }
    }
  }
  assert(nexact + ninexact == 48, "12 orientations x 4 sides");
  assert(ninexact > 0, "an `== 0.0` assertion passes only on the orientations you tried");
  assert(worst < 1e-15, "and the worst spread is still an ulp, got {worst}");
  assert((worst / BAND_TOPS ?? 1.0) < 1e-14, "which is the flat-vs-strip comparison that holds");

  // What the strip's zigzag actually is: on the east family the edge MIDPOINTS are
  // exactly colinear and the corners alternate; on the north family the midpoints
  // themselves stagger by the full band.  The corners hide that difference.
  assert(surface_band_mids(p, cells, 1) == 0.0, "east-family midpoints are colinear");
  assert(abs(surface_band_mids(p, cells, 0) - BAND_SIDES) < 1e-12,
         "north-family midpoints span the whole band — the row stagger");

  // And the control that "averaging, not fitting" is a real distinction: a least
  // squares formulation would have a residual here to threshold, and the average
  // has an exact answer instead.
  assert(surface_lsq_residual(p, cells, 1) == 0.0, "nothing to fit on the east family");
  assert(abs(surface_lsq_residual(p, cells, 0) - (11.0 / 12.0 ?? 0.0)) < 1e-12,
         "and exactly 11/12 to threshold on the north one, got {surface_lsq_residual(p, cells, 0)}");
}

// @HXD-002 — THE EXACTNESS CHECK, AND THE TWO OBVIOUS WAYS TO WRITE IT WRONG.
// The claim `surface_heading` makes is that the AVERAGED direction is exactly one
// of the twelve headings — no residual perpendicular component — which in integers
// is a zero CROSS PRODUCT against that heading's step.
//
// Both simpler-looking tests fail on geometry nobody doubts, and they fail in a
// way that reads like a real defect:
//
//   "one component is zero"   only 4 of the 12 headings have one; a side running
//                             NE is (1,3), both components non-zero and exact
//   "an integer multiple"     the summed vector over a nine-edge run is (0,14)
//                             where the heading step is (0,6) — the SUM is not a
//                             whole number of steps, and never claimed to be
//
// The magnitude carries no claim at all: it is the run's length.
fn test_exactness_is_a_cross_product_not_a_multiple() {
  p = cottage();
  cells = cottage_cells();

  w = surface_of(p, cells, 0);
  assert(surface_heading(w) == 3, "side 0 runs on heading 3");
  assert(w.ws_dk == 0 && w.ws_dm == 14, "its summed vector is (0,14), got ({w.ws_dk},{w.ws_dm})");
  st = head_step(3);
  assert(st.0 == 0 && st.1 == 6, "and heading 3's step is (0,6)");
  assert(w.ws_dk * st.1 - w.ws_dm * st.0 == 0, "the cross product is exactly zero — that is the claim");
  assert((14 % 6 ?? 0) != 0, "while 14 is not a multiple of 6 — the sum is not a whole number of steps");

  // Both wrong tests, counted over the six rotations x four sides.
  nexact = 0;
  fail_component = 0;
  fail_multiple = 0;
  for rot in 0..6 {
    q = plan_new(0, 0, 5, 4, rot, false);
    c = patch();
    lb = labels_new(Q0, R0, W, H);
    draw_floor(q, c, lb, 1);
    for side in 0..4 {
      v = surface_of(q, c, side);
      hd = surface_heading(v);
      if hd >= 0 { nexact = nexact + 1; }
      if !(v.ws_dk == 0 || v.ws_dm == 0) { fail_component = fail_component + 1; }
      sp = head_step(if hd >= 0 { hd } else { 0 });
      ismult = false;
      if sp.0 != 0 { ismult = (v.ws_dk % sp.0 ?? 1) == 0 && v.ws_dm * sp.0 == v.ws_dk * sp.1; }
      if sp.0 == 0 { ismult = (v.ws_dm % sp.1 ?? 1) == 0 && v.ws_dk == 0; }
      if !ismult { fail_multiple = fail_multiple + 1; }
    }
  }
  assert(nexact == 24, "every one of the 24 surfaces is exact by the cross product, got {nexact}");
  assert(fail_component == 16, "'one component is zero' rejects 16 of them, got {fail_component}");
  assert(fail_multiple == 12, "'an integer multiple' rejects 12 of them, got {fail_multiple}");
}

// @HXD-003 — THE MITER CLOSES THE OUTLINE, AND ITS CORNER IS NOT THE PLAN'S.
// Two adjacent surfaces are two lines; where they cross is the corner a renderer
// should cut to, and `surface_quad` returns each side miter to miter.  Then
// consecutive quads share an endpoint EXACTLY and the outline closes.  The
// un-mitered `surface_span` — each side's own extent projected onto its own line —
// leaves a real gap at every corner, and it is not a small one.
//
// What that corner is NOT is the corner of the `Plan`.  The recovered outline is
// the face of the wall as DRAWN, and the drawn footprint is the plan quantised to
// cells.  The two live in different number systems: the recovered corner is an
// exact rational of the lattice (a multiple of 3/4 across the rows), while the
// plan's half-depth is an irrational multiple of `sqrt(3)/2`.  They agree only at
// the origin — so this is not an offset a caller can correct for, and swept over
// plan sizes the difference is neither constant nor even consistently outward.
fn test_the_miter_closes_the_outline_but_is_not_the_plans_corner() {
  p = cottage();
  cells = cottage_cells();

  for side in 0..4 {
    nx = (side + 1) % 4 ?? 0;
    assert(quad_gap(p, cells, side, nx) == 0.0, "mitered sides {side}->{nx} share their endpoint");
    assert(abs(span_gap(p, cells, side, nx) - 0.4589642122738423) < 1e-12,
           "the un-mitered spans leave a corner gap, got {span_gap(p, cells, side, nx)}");
  }

  // Swept over plan sizes: the recovered corner never lands on the plan's, it is
  // sometimes outside the plan rectangle and sometimes inside it, and the distance
  // wanders with how the fill quantises.
  noutside = 0;
  ninside = 0;
  nearest = 9.0;
  for wid in 2..9 {
    q = plan_new(0, 0, wid, wid - 1, 0, false);
    c = patch();
    lb = labels_new(Q0, R0, W, H);
    draw_floor(q, c, lb, 1);
    (mx, my, mok) = surface_miter(q, c, 3, 0);
    assert(mok, "adjacent sides of a closed form always cross");
    (cx, cy) = plan_side_point(q, 0, 0.0);
    d = sqrt((mx - cx) * (mx - cx) + (my - cy) * (my - cy)) ?? 0.0;
    if d < nearest { nearest = d; }
    (_lu, lv) = plan_to_local(q, mx, my);
    // the recovered corner sits on the row lattice: an exact multiple of 3/4
    rows = (abs(lv) / 0.75) ?? 0.0;
    assert(abs(rows - ((round(rows) ?? 0.0))) < 1e-12,
           "the recovered corner is a rational of the lattice, got {abs(lv)}");
    // the plan's half-depth is not — it is (dep * sqrt(3) / 2)
    steps = (plan_hd(q) / 0.75) ?? 0.0;
    assert(abs(steps - ((round(steps) ?? 0.0))) > 1e-6, "and the plan's half-depth is not");
    if abs(lv) > plan_hd(q) { noutside = noutside + 1; } else { ninside = ninside + 1; }
  }
  assert(nearest > 0.4, "no plan size brings the two corners closer than {nearest} world units");
  assert(noutside > 0 && ninside > 0,
         "and the recovered corner is outside the plan on some sizes and inside on others");
}

// @HXD-004 — A DOOR IS A MATERIAL, AND NEVER A HOLE.  The field stores ONE
// material per wall slot, so an opening REPLACES the wall material rather than
// sitting in a channel beside it.  That leaves exactly one way to get it wrong,
// and it is the way that looks simplest: write 0.
//
// `OPEN_NONE` is not "no opening" — it is the wall material 0, which means NO
// WALL.  Writing it deletes the edge, and the wall stops being one closed run:
// the boundary that had no loose ends grows two.  Everything downstream that
// treats a wall as one run — recovery, collision, the surface above — is reading
// a different shape from then on.
//
// `OPEN_GAP` is what "a hole in the wall" actually means here: the leaf is missing
// entirely, nothing collides, and THE WALL IS STILL THERE.  An opening has
// geometry of its own — a frame, a lintel, a ragged jamb — and absence has none.
fn test_a_door_is_a_material_and_never_a_hole() {
  p = cottage();
  cells = cottage_cells();

  for kind in [OPEN_DOOR, OPEN_WINDOW, OPEN_GAP, OPEN_NONE] {
    edges = edgeset_new(Q0, R0, W, H);
    ne = draw_walls(cells, edges, WALLMAT);
    assert(ne == 38 && wall_chain_ends(edges, Q0, R0, W, H) == 0,
           "the drawn wall is a closed loop before anything is placed");
    placed = place_opening(p, cells, edges, 0, 0.5, 1, kind);
    assert(placed == 1, "one edge takes the opening");
    left = wall_edge_count(cells, edges);
    ends = wall_chain_ends(edges, Q0, R0, W, H);
    if kind == OPEN_NONE {
      assert(left == 37, "material 0 DELETES the edge: 38 -> {left}");
      assert(ends == 2, "and the closed loop grows two loose ends, got {ends}");
    } else {
      assert(left == 38, "an opening keeps its edge, got {left}");
      assert(ends == 0, "so the wall is still one closed run, got {ends}");
    }
  }
}

// @HXD-005 — A WALL IS THE BOUNDARY OF A FILLED REGION, WHICH IS WHY IT COSTS NO
// FLOOR.  Walls are edges between a cell of the region and a cell outside it, so
// the boundary is a closed loop BY CONSTRUCTION at any size, orientation or tie —
// and every floor cell survives.  Drawn instead as a band of cells along the
// outline, this cottage would lose the 16 of its 27 cells that touch the outside:
// 59% of the house, and the remaining room is not the room that was designed.
//
// Counting those edges has one trap, and it comes from the field's own symmetry:
// a boundary edge belongs to a PAIR of cells, and the enumeration a caller writes
// first — the three directions each cell "owns" — misses every edge whose owner is
// the outside cell.  Both neighbour conventions in this family have such an
// enumeration, and they are not the same three, so the same mistake reads 19 in one
// and 17 in the other while the answer is 38.
fn test_the_wall_costs_no_floor_and_the_count_needs_all_six() {
  cells = cottage_cells();
  assert(hexset_count(cells) == 27, "the cottage keeps all 27 floor cells");
  edges = edgeset_new(Q0, R0, W, H);
  ne = draw_walls(cells, edges, WALLMAT);
  assert(ne == 38, "with 38 wall edges around them");
  assert(ring_cells(cells) == 16, "16 of the 27 touch the outside — what a band wall would eat");

  assert(wall_edge_count(cells, edges) == 38, "all six directions from every cell: 38");
  assert(count_three_dirs(cells, edges, true) == 19, "the three a cell owns, hex_grid order: 19");
  assert(count_three_dirs(cells, edges, false) == 17, "the same shortcut, hex_field order: 17");
}

// @HXD-006 — A FEATURE IS INDEXED BY ITS EXACT `t`, NEVER BY ITS POSITION IN THE
// RUN.  The model puts a door at `(side, t)`; the field marks edges.  Along a run
// of `n` edges the recoverable positions are the edge centres `(2i+1)/2n` — exact
// rationals, odd numerator over `2n` — and `place_opening` snaps to the nearest,
// which is what makes a feature survive the round trip (`hex_fit`'s @HXI-005 is
// the doorstep that refuses the positions this cannot get back).
//
// The trap is the run's ORDER.  `side_edges` fills it in cell-scan order, which is
// not `t` order: the numerators of this side come out `1 5 3 7 9 13 11 15 17`.  So
// two edges that are ADJACENT along the wall can sit two apart in the run with a
// plain wall edge between them — which is exactly what a two-edge window does here.
// Read features with `feature_find`, which returns numerators in ascending `t`.
fn test_a_feature_is_indexed_by_its_t_not_its_position() {
  p = cottage();
  cells = cottage_cells();
  edges = edgeset_new(Q0, R0, W, H);
  draw_walls(cells, edges, WALLMAT);

  sr = side_edges(p, cells, 0);
  n = side_run_len(sr);
  assert(n == 9, "side 0 carries 9 edges, got {n}");
  order = "";
  for i in 0..n { order += "{feature_t_num(sr, i, n)} "; }
  assert(order == "1 5 3 7 9 13 11 15 17 ", "the STORAGE order of the t numerators is [{order}]");

  // A door in the middle, and a two-edge window near the far end.
  assert(place_opening(p, cells, edges, 0, 0.5, 1, OPEN_DOOR) == 1, "one door edge");
  assert(place_opening(p, cells, edges, 0, 0.1, 2, OPEN_WINDOW) == 2, "two window edges");

  door = feature_find(p, cells, edges, 0, OPEN_DOOR);
  window = feature_find(p, cells, edges, 0, OPEN_WINDOW);
  assert(len(door) == 1 && (door[0] ?? 0) == 9, "the door is at t = 9/18, exactly the middle");
  assert(len(window) == 2 && (window[0] ?? 0) == 1 && (window[1] ?? 0) == 3,
         "the window took the two nearest centres to 0.1, t = 1/18 and 3/18");
  assert(feature_t_admissible(9, 9) && !feature_t_admissible(8, 9),
         "a legal t has an ODD numerator over 2n — an even one is a corner, not a centre");

  // …and those two adjacent window edges are NOT adjacent in the run.
  assert(surface_edge_mat(p, cells, edges, 0, 0) == OPEN_WINDOW, "run position 0 is window");
  assert(surface_edge_mat(p, cells, edges, 0, 1) == WALLMAT, "position 1 is plain wall");
  assert(surface_edge_mat(p, cells, edges, 0, 2) == OPEN_WINDOW, "position 2 is the window again");

  // What a renderer draws is not the marked EDGE but its interval ON the surface:
  // the wall is one quad and each feature is a sub-interval of it, so the zigzag
  // never reaches the picture.  One edge is exactly 1.0 world unit of clear
  // width — 0.87 m, a real domestic door, which a cell-thick wall could not express.
  (ex0, ey0, ex1, ey1) = surface_edge_span(p, cells, 0, 0);
  wide = sqrt((ex1 - ex0) * (ex1 - ex0) + (ey1 - ey0) * (ey1 - ey0)) ?? 0.0;
  assert(abs(wide - 1.0) < 1e-12, "one edge opens 1.0 world unit of wall, got {wide}");
  assert(abs(wide * 0.8660254037844386 - 0.8660254037844386) < 1e-12, "which is 0.87 m");
  // and both its ends lie on the wall's own line, not on the zigzag
  (sx0, sy0, sx1, sy1) = surface_span(p, cells, 0);
  along = sqrt((sx1 - sx0) * (sx1 - sx0) + (sy1 - sy0) * (sy1 - sy0)) ?? 0.0;
  assert(abs(ex0 - sx0) < 1e-12 && abs(ex1 - sx1) < 1e-12,
         "the interval runs along the same line as the wall quad");
  assert(along > wide, "and is a sub-interval of it: {wide} of {along}");
}

// @HXD-007 — THE RIDGE HAS TO RUN PAST THE GABLE ENDS.  A roof height is the
// distance from a ridge LINE, so where a cell's nearest point on that line is an
// ENDPOINT rather than a perpendicular foot, the distance goes radial and the roof
// rolls over like a quarter cone.  Stopping the ridge at the wall — the obvious
// thing to do — does exactly that, and the gable end drops by sqrt(3)/4 world
// units, 0.375 m of bulge on a real building.
//
// It is also a lesson about the measurement.  Compare the two rules over the
// MASSING alone and they are identical to the last bit, because no cell of the
// massing is past the ridge's end.  The difference lives entirely in the halo ring
// `grow_ring` adds — which is there because the eave has to interpolate against
// cells outside the wall — so an instrument that skips the halo reports a perfect
// result over the empty set.
fn test_the_ridge_runs_past_the_gable_ends() {
  p = cottage();
  massing = cottage_cells();
  hw = plan_hw(p);
  hd = plan_hd(p);

  // On the massing alone the two rules agree — the comparison has nowhere to bite.
  past_m = ridge_heights(p, massing, hw + hd + HEX_LEN);
  stop_m = ridge_heights(p, massing, hw);
  worst_m = 0.0;
  for r in R0..R0 + H {
    for q in Q0..Q0 + W {
      if hexset_get(massing, q, r) {
        d = abs(height_get(past_m, q, r) - height_get(stop_m, q, r));
        if d > worst_m { worst_m = d; }
      }
    }
  }
  assert(worst_m < 1e-12, "over the massing the two ridges agree — an empty measurement");

  // With the halo the difference is the bulge.  `draw_roof` is the shipped rule —
  // it runs the ridge past both ends — and the hand-built ridge below is the
  // tempting one that stops at the wall.
  grown = patch();
  grow_ring(massing, grown);
  assert(hexset_count(grown) == 49, "the halo takes 27 cells to 49, got {hexset_count(grown)}");
  past = heights_new(Q0, R0, W, H);
  draw_roof(p, grown, past, 2.0, 0.5, 0);
  stop = ridge_heights(p, grown, hw);
  worst = 0.0;
  ndiff = 0;
  for r in R0..R0 + H {
    for q in Q0..Q0 + W {
      if hexset_get(grown, q, r) {
        d = abs(height_get(past, q, r) - height_get(stop, q, r));
        if d > 1e-9 { ndiff = ndiff + 1; }
        if d > worst { worst = d; }
      }
    }
  }
  assert(ndiff == 10, "10 of the 49 cells sit past the shortened ridge, got {ndiff}");
  assert(abs(worst - 0.4330127018922193) < 1e-12,
         "and the worst of them rolls over by sqrt(3)/4 units, got {worst}");
  assert(abs(worst * 0.8660254037844386 - 0.375) < 1e-12, "which is 0.375 m of bulge");

  // The rolling over is not a defect in `roof_ridge` — it is what a SHORT ridge
  // means, and on a hip roof it is the shape you are asking for.  Same call, same
  // mechanism, and the only difference is whether the ridge was meant to end there.
  hip = heights_new(Q0, R0, W, H);
  draw_roof(p, grown, hip, 2.0, 0.5, 2);
  nlower = 0;
  for r in R0..R0 + H {
    for q in Q0..Q0 + W {
      if hexset_get(grown, q, r) && height_get(hip, q, r) < height_get(past, q, r) - 1e-9 {
        nlower = nlower + 1;
      }
    }
  }
  assert(nlower > 0, "a hip's shortened ridge lowers the ends deliberately, {nlower} cells");
}