Design: Node Canonicalization for the Geometry Model#

Status: Draft Branch: node_canon/design_document Author: Connor Ferster Date: 2026-07-17

1. Problem#

papermodels converts marked-up drawings into a directed “geometry graph” where lower-rank elements (joists, beams) transfer load down to higher-rank elements (walls, columns, footings). The load path is derived from where geometries intersect and which geometries correspond between adjacent planes.

The library works but is fragile under floating-point error. The fragility is not caused by GEOS/shapely predicates (intersects, contains) — those are robust. It is caused by construction followed by re-testing: every intersection(), nearest_points(), interpolate(), centroid, and affine rotation produces coordinates rounded to the nearest double that almost never lie exactly on the input geometry, and the result is then fed back into an exact predicate or an exact float comparison, where a branch flips or an assert fires.

The deeper structural cause is that the same physical crossing is computed many times, independently, and the results are never reconciled. Concretely, in element.py::get_geometry_intersections (the O(n²) build loop):

  • For each ordered pair with j_rank > i_rank, geom_ops.get_intersection is called twice — once to store the crossing on i as an intersection_below, once to store it on j as an intersection_above (element.py:1152-1164). These two evaluations of one physical point can differ in their low bits.

  • Each Intersection NamedTuple carries a copy of other_geometry (element.py:43). The same wall is copied into every element that touches it.

  • Downstream methods (get_transfer_extents, align_frames_to_centroids, _get_support_locations, get_collector_extents) recompute self.geometry.intersection(other_geometry) again, producing further independent evaluations of the same crossing.

A single physical node therefore exists as several slightly-different floats scattered across the model, with nothing forcing them to agree. The scattered tolerance workarounds in geom_ops.py (grid_size=1e-3, eps=1e-6, buffer(1e-3), rel_tol=5e-2, exact == 0.0 / == 1.0 branch tests) patch individual symptoms but do not remove the cause; the instability simply moves.

2. Goal#

Establish a single canonical identity for every node (intersection point) in the model, computed once, shared by every element incident to it, and accelerated by spatial indexing so the model can be built and queried efficiently.

“Intersecting geometries share the same node location” must become true by construction, not by hoping two independent .intersection() calls round the same way.

Non-goals#

  • Replacing shapely or moving to exact arithmetic (CGAL, rational coordinates). The domain (drawing markup, ~1 mm meaningful precision, mostly rectangles and straight lines) does not justify the dependency or performance cost.

  • Changing the structural/analysis semantics (rank, reaction types, load transfer). This work changes how geometry identity is managed, not what a beam or column means.

  • Changing the networkx DiGraph as the element→element structural graph. It stays; its edges become backed by shared node ids.

3. Key Insight#

The spatial index is not the fix — it is only acceleration. The fix is node canonicalization: a single registry into which every computed intersection point is inserted, snapping to an existing node within an absolute tolerance if one exists, otherwise registering a new node. Every element meeting at that crossing then references the same node id.

This is equivalent to building a noded arrangement (a planar graph with shared vertices) over the input geometry.

4. Data-structure choice (and one correction)#

A KD-tree indexes points and answers nearest-neighbor / radius queries. It does not index extended geometry (line segments, polygons) for intersection queries. The two halves of this problem need two different structures:

Need

Structure

Notes

“Which lines/polygons cross this one?” (broad phase)

shapely.STRtree (R-tree)

Already available in shapely ≥ 2.0, GEOS-backed. Replaces the O(n²) loop. No new dependency.

“Is there already a node within abs_tol of this point?”

Quantized-coordinate hash (preferred) or KD-tree

Nodes are points. A dict keyed by coordinates quantized to abs_tol gives O(1) incremental insert/lookup; probe the 9-cell neighborhood to catch points straddling a bucket boundary. A scipy.spatial.cKDTree is only worth adding for a one-shot batch clustering pass.

The “lines/polygons KD-tree” from the original idea is really an STRtree. The “nodes KD-tree” is a real point index, but a quantized hash is a better fit because node insertion is incremental.

5. Proposed Architecture#

5.1 GeometryModel — the single source of truth#

A new class owns everything geometry-related in one place:

GeometryModel
├── geometries : dict[GeomId, shapely geometry]     # single source of truth, NO copies
├── ranks      : dict[GeomId, int]
├── planes     : dict[GeomId, PlaneId]              # page/plane id
├── index      : shapely.STRtree                    # broad-phase over geometries
├── nodes      : NodeRegistry                        # canonical node identities
└── incidence  : dict[NodeId, list[Incidence]]      # which geoms meet at each node

  • Element stops owning copies of other elements’ geometry. It becomes a view: it holds its own GeomId and references NodeIds, and answers intersections_below() / corresponding_above() by querying the model.

  • The networkx DiGraph is retained as the structural graph; its edges are backed by shared NodeIds.

  • Walls are stored in the STRtree as centerlines, not polygons (see §5.4). The wall’s polygon is retained separately for correspondence (§5.5).

5.2 NodeRegistry — canonical node identities#

NodeRegistry
├── node_abs_tol : float                             # absolute snapping tolerance (real-world units); separate from noding_abs_tol (§7)
├── coord   : dict[NodeId, tuple[float, float]]      # canonical coordinate per node
└── _buckets: dict[tuple[int, int], list[NodeId]]    # quantized spatial hash

  get_or_create(point) -> NodeId
      q = quantize(point, node_abs_tol)
      for candidate in neighbor_buckets(q):          # 9-cell neighborhood
          if distance(point, coord[candidate]) <= node_abs_tol:
              return candidate                        # snap to existing node
      return register_new(point)                      # otherwise new node

Canonicalization can also be done as a batch clustering pass (collect all raw intersection points, cluster within abs_tol via union-find over radius-neighbor pairs, assign one canonical coordinate — e.g. the cluster centroid — per cluster). The batch form is where a cKDTree would earn its place; the incremental get_or_create form needs no new dependency.

5.3 Incidence — what meets at a node#

Incidence(geom_id: GeomId, rank: int, role: str)   # role ∈ {"endpoint", "interior", "boundary"}

intersections_above / intersections_below collapse into a single query: “elements incident to this node, filtered by rank.” One structure replaces the two independently-computed lists maintained today.

5.4 Wall representation: polygon in, centerline in the index#

Walls are drawn and stored as polygons. This is deliberate — a polygon gives a forgiving, area-based way to determine correspondence between levels (adjacent PDF pages), which is robust to the wall being drawn slightly differently from one plane to the next (§5.5).

But a polygon is the wrong thing to intersect for load-path purposes: what a joist or beam actually bears on is the wall’s centerline. The historical sticking point has been computing that centerline. The existing geom_ops.get_rectangle_centerline explodes the polygon’s own edges and assumes a clean, effectively axis-aligned 4-vertex rectangle; it degrades on rotated walls, walls with vertex noise, or walls with extra/clipped vertices, and the surrounding orientation logic then falls back to an “is it horizontal or vertical?” test that fails for anything off-axis.

Generalized approach (implemented): oriented bounding box. geom_ops.get_wall_centerline fits the minimum-area rotated bounding rectangle (OBB, shapely.minimum_rotated_rectangle) around the wall polygon, then connects the midpoints of the OBB’s two short (opposite) end faces:

  • The OBB is a true rectangle regardless of the input’s orientation, so the result is correct at any angle — no axis-alignment assumption.

  • Because the OBB is fitted rather than read off the raw edges, it absorbs vertex noise, extra vertices, and clipped corners; an approximately rectangular wall yields the same spine as the clean rectangle would.

  • Selecting an opposite edge pair (not merely “the two shortest edges”) means a near-square wall degenerates to a valid long-axis spine instead of a diagonal.

Validated on axis-aligned, 30°-rotated, noisy (extra vertex + jitter), clipped-corner, and near-square walls: all produce a correct long-axis centerline of the expected length.

In the build pipeline, each wall element contributes its centerline to the STRtree and the node registry (so crossings are computed against the spine), while its polygon is retained on the element for correspondence.

5.5 Correspondents stay ratio-based (decided)#

Cross-plane correspondence remains computed from polygon area/length overlap ratios, not from shared nodes. This is a deliberate decision, not a limitation:

  • Walls (and columns) keep their polygon precisely so correspondence can be forgiving to plane-to-plane drawing differences.

  • The overlap ratio is load-bearing information: when multiple elements converge at one location between levels, the ratio provides the basis for choosing a specific load path (see prioritize_correspondents in element.py). Collapsing correspondence to a shared node would discard the signal used to disambiguate.

Correspondents therefore benefit from the shared STRtree broad phase (fewer pairwise checks) but not from the node registry. The node registry governs intersections (same-plane crossings); correspondence governs between-plane transfer and stays ratio-based.

6. Build Pipeline#

1. Ingest + snap to grid            (existing: annotations.py:45, set_precision grid=1e-3)
1a. Wall polygons -> centerlines    (geom_ops.get_wall_centerline; polygon kept for correspondence — §5.4)
2. Tolerance noding (Phase 2)       (extend/trim endpoints within abs_tol — §7)
3. Build STRtree over geometries    (walls contribute their centerline)
4. Broad phase: per geom, query index for candidates; filter by rank + plane
5. Compute each crossing ONCE; register its point via NodeRegistry.get_or_create
6. Build incidence graph keyed by NodeId; derive above/below by incident ranks
7. Build networkx DiGraph edges from incidence (backed by shared NodeIds)

Step 5 is the crux: each physical crossing is evaluated exactly once and its identity is shared thereafter.

7. Phase 2 — Tolerance Extend/Trim (Noding)#

A preprocessing step that repairs “approximately intersecting” geometry, governed by an absolute tolerance noding_abs_tol in real-world units. This is separate from the canonicalization tolerance node_abs_tol (§5.2), so endpoint repair can reach wider than crossing-merging (decided, §9):

  • For each line endpoint, query the STRtree for geometry within noding_abs_tol.

  • If the nearest point on a candidate is within noding_abs_tol:

    • extend the line to that point if it falls short, or

    • trim the line to that point if it overshoots.

Absolute (not relative) tolerance is deliberate: a relative tolerance makes endpoint behavior depend on line length, reintroducing the scale-dependence we are trying to remove.

Caveat — cascading snaps. Moving endpoint A onto line L can pull it off line M it previously touched. A single sweep over a sparse, mostly-orthogonal structural drawing is usually sufficient, but the operation is not guaranteed to converge in one pass. The implementation must iterate to a fixed point (or a capped number of passes, warning if the cap is hit).

Ordering. Phase 2 runs before the STRtree/node-registry build (steps 3-6), so the arrangement is already topologically clean when crossings are canonicalized.

8. Migration / Impact#

This is a rewrite of the intersection-building layer, not an add-on. Call sites that today “compute a crossing” must change to “look up a node”:

  • element.py::get_geometry_intersections — the O(n²) double loop → STRtree broad phase + single-evaluation canonicalization.

  • Intersection / Correspondent construction — stop copying other_geometry; reference GeomId / NodeId.

  • Downstream recomputation sites that call self.geometry.intersection(other) — get_transfer_extents, align_frames_to_centroids, _get_support_locations, get_collector_extents — read canonical nodes from the model instead.

  • Many geom_ops.py snapping workarounds become unnecessary once every coordinate is a canonical node, and can be removed in a follow-up cleanup.

  • Wall-handling call sites that use get_rectangle_centerline with the axis-aligned assumption (and the "horizontal"/"vertical" orientation fallbacks around it) migrate to get_wall_centerline. get_rectangle_centerline may be retained for genuinely axis-aligned internal box geometry, or removed once no longer relied upon.

The existing test fixtures and the geometry-graph tests are the migration safety net; behavior (graph topology: nodes, edges, intersection counts) must be preserved.

9. Risks & Open Questions#

  • Batch vs incremental canonicalization — decided: incremental now, batch later. Phase 1 ships the incremental get_or_create + quantized hash form (dependency-free; canonical coordinate is “first point wins”). Batch clustering (union-find + cluster centroid, requiring a scipy cKDTree) is kept as a future option to be introduced only if stability testing shows “first point wins” produces unstable canonical coordinates.

  • Bucket-boundary straddling. Two points within abs_tol can fall in adjacent buckets; the 9-cell neighborhood probe handles the 2D case and must be covered by tests.

  • Tolerances — decided: separate. Noding (Phase 2 extend/trim) and node canonicalization use distinct tolerances rather than a single shared one, so endpoint repair can reach wider than crossing-merging without forcing the two to move together. (E.g. node_abs_tol for canonicalization, noding_abs_tol for Phase 2.) Both remain absolute, real-world units.

  • Cascade convergence bound — decided: cap + warning. Phase 2 iterates to a fixed point up to a capped number of passes; if the cap is hit it emits a warning and proceeds with the best-so-far arrangement rather than raising. Pathological inputs degrade gracefully instead of aborting the build.

  • Correspondents (cross-plane overlap) stay ratio-based — decided; see §5.5. The overlap ratio is retained deliberately because it is the basis for choosing a load path when multiple elements converge. They benefit from the shared STRtree broad phase but not from the node registry.

  • Wall centerline endpoints from the OBB — decided: accept overhang. OBB endpoints extend to the bounding-box extent, which for a clipped/irregular wall can slightly overhang the true polygon end. This is accepted as desirable for a bearing spine; it should still be confirmed (not clamped) against the extent/bearing calculations that consume it.

10. Validation Strategy#

  • Preserve existing tests. All current geometry-graph and integration tests must pass unchanged.

  • Jitter test. Perturb every input coordinate in the fixtures by ±1e-6…±1e-9 and assert the resulting graph topology (node set, edge set, intersection counts) is invariant. This converts “works on my drawings” into a checkable robustness property and is the acceptance criterion for the whole effort.

  • Idempotence. Running canonicalization twice yields identical node ids and coordinates.

11. Proposed Phasing#

  1. Phase 0 (this doc). Design + agreement. ✅

  2. Phase 1 (done). GeometryModel + NodeRegistry prototype behind the existing API, validated against current fixtures; jitter test added. Implemented in src/papermodels/datatypes/geometry_model.py, tested in tests/test_geometry_model.py:

    • NodeRegistry — incremental get_or_create, quantized-hash buckets, 9-cell neighborhood probe, nearest-wins snapping, “first point wins” canonical coordinate (the “incremental now, batch later” decision, §9).

    • GeometryModel.from_elements — STRtree broad phase, walls indexed by centerline (polygon retained), single-evaluation crossings canonicalized through the registry, incidence keyed by NodeId, rank-derived above/below queries. Not yet wired into GeometryGraph (parallel prototype).

    • Validation: reproduces all documented crossings of the intersections.pdf fixture; jitter test (perturb every coordinate by 1e-6…1e-9) shows the intersection topology — node set, edge set, incidence multiset — is invariant; build is idempotent.

  3. Phase 2 (done). Tolerance extend/trim noding preprocessing (noding_abs_tol). Implemented in src/papermodels/geometry/noding.py, tested in tests/test_noding.py:

    • node_geometries — per line endpoint, query the STRtree (dwithin) for geometry within noding_abs_tol and relocate the endpoint onto the nearest point on the nearest candidate (extend if short, trim if overshooting). Endpoints already touching a candidate are anchored and left alone; polygons (columns) and fixed_ids (wall centerlines) are snap targets but never moved.

    • Cascade handling: each pass works from a frozen snapshot; passes iterate to a fixed point, capped at max_passes with a warning (not an error) if the cap is hit — the cap+warning decision (§9). Returns a NodingReport (passes, converged, moves).

    • Wired into GeometryModel.from_elements as an opt-in pre-index step (noding_abs_tol=…), separate from node_abs_tol (§9). Validated: bridges a sub-tolerance gap into a real crossing, preserves already-clean fixture crossings, holds wall spines fixed, and the jitter test stays invariant with noding on.

  4. Phase 3 (substantive migration done). Migrate the intersection-building layer to canonical identities and robust constructions.

    • 3a (done). get_geometry_intersections computes each crossing once and snaps its region through a shared NodeRegistry, so the below view (on the lower-rank element) and the above view (on the higher-rank element) carry the byte-identical coordinate — removing the “computed twice, never reconciled” root cause. Regression test: test_intersection_region_shared_below_and_above.

    • 3b (done). O(n²) all-pairs scan → shapely.STRtree broad phase (candidates queried by the extent polygon when present, iterated in annotation order so topology is identical). get_rectangle_centerline reimplemented on the oriented bounding box (the robust get_wall_centerline construction, on_long_edge preserved), migrating every wall/rectangle centerline call site at once (§5.4). Downstream overlap recomputes (get_transfer_extents, align_frames_to_centroids) now read the canonical stored other_overlap instead of re-evaluating self.geometry.intersection(other).

    • Acceptance met. The design’s jitter criterion (§10) now holds on the production GeometryGraph, not just the prototype: test_production_graph_topology_invariant_under_jitter perturbs every input coordinate (1e-6/1e-8) and asserts the node set and typed edge set are invariant across the intersections and sketch_to_scale fixtures. All previously-passing behavioural tests remain green.

    Deliberately deferred (with rationale):

    • Retire geom_ops.py tolerance workarounds (grid_size=1e-3, rel_tol, buffer(1e-3), exact == 0.0/== 1.0 tests). This doc always scoped these as “a follow-up cleanup” (§8). With the production jitter test now passing with them present, removing them offers no robustness upside and real regression risk to current numerics; it should be its own change, guarded by the jitter test now in place.

    • Stop copying other_geometry; reference GeomId/NodeId. Intersection already carries other_tag; replacing the geometry copy with a lookup requires threading the GeometryModel through every consumer (many Element methods have no graph handle today). That is a model-as-view re-architecture, high-churn and low correctness value now that node identity is canonical — best done as a focused follow-up, not folded into this migration.

    • test_wall_point_load_locations (pre-existing failure). Traced to JoistArrayModel laying the first joist of an array off-grid at the trib edge (all 12 WT4.0 point loads come from the one J4.0 array collector). It fails on main, lives in collector generation, and is orthogonal to node canonicalization.