Building a CAD kernel

This book is a lecture series on how to build a modern CAD kernel. It starts from a linear algebra course and a little programming, and it ends with you able to read, debug and extend a real boundary-representation kernel: its intervals, its NURBS, its booleans, its sketch and assembly solvers, its fillets and shells, and the incremental build that runs them all.

The kernel it teaches is Geop, an open-source kernel in Rust. Every statement this book makes about the kernel is meant to hold against the code in this repository, and paths and function names are cited so that you can go and read them: core/geop-core-geometry/src/aabb.rs, compute_aabb. Where something is not implemented, the book says so. The kernel changes faster than any text, so if a citation has gone stale, the code is right. Please tell us on the Discord.

Who it is for

You should know:

• Linear algebra: vectors, matrices, solving systems, rank, the idea of a null space.
• Calculus: derivatives, Taylor expansion, the chain rule, partial derivatives. Newton’s method helps.
• Some programming: you can read code in a C-like language. Rust is used in the excerpts, but none of its finer points are needed for the first eight parts. Reading it will come with practice.

You do not need to know anything about CAD, computational geometry or numerical analysis. Everything is built up from there: what a CAD kernel is, what floating point does, what a Bézier curve is. The book is meant to be read in order, since later lectures use what earlier ones established. A reader with a CAD background might skim Part 0 and start with Lecture 3, as the interval discipline shapes everything that follows.

How the lectures are organised

Each lecture opens with what you will be able to do after it, then develops the idea from first principles with worked examples. Boxes marked In the kernel point at where it lives in the code. A recap and exercises close it, with solutions or hints folded away. Try the exercise before you open the solution: much of the understanding in this subject is in the attempt.

The appendices hold references, a glossary and a map of the crates, the last with a reference page for many of them. For a generated API reference, see docs.geop-cad.dev.

If you are short of time, here is a thread through the book that carries the essential idea: Lectures 4, 5, 6, 12, 14, 17, 21 and 42.

The demos

Many pages have interactive figures. They run in your browser with no network: drag a control point, slide a knot, press Step to see one iteration of an algorithm. Each runs the mathematics it shows, in a few hundred lines of plain JavaScript in book/src/demos/, and where it is simpler than what the kernel does, its caption says how. The pure mathematics of each demo is tested against its formula by book/scripts/test_demo_math.mjs. On a phone, a touch on a handle drags it, and a touch elsewhere scrolls the page. On a keyboard, focus a figure and press n or p to choose a point, and the arrow keys to move it. Animations do not start by themselves if your system asks for reduced motion.

Here is one of the figures that shows a real part as the real kernel meshed it. The other figures like it link to the example in the web app, where the kernel itself runs in your browser as WebAssembly.

Running the kernel

You can do almost everything in this book without building anything:

• In the browser: app.geop-cad.dev, the kernel as WebAssembly. Open an example with a link, such as ?example=rounded_box.

• On the command line: curl -fsSL https://geop-cad.dev/install.sh | sh installs geop (on Windows, irm https://geop-cad.dev/install.ps1 | iex). Then geop examples --out-dir examples writes every built-in example as a program and a mesh, geop compile part.geop builds a program, and geop open my-project edits a folder in the browser, on the kernel running on your machine.

• For an agent: geop mcp runs the editor as a Model Context Protocol server, described in docs/mcp.md.

• From source: this is a Cargo workspace. Whenever you run something that searches geometry, cap it, as this repository’s rule demands:

  timeout 600 prlimit --as=16000000000 cargo test -p geop-core-geometry
  

  A search that runs away does not stop on its own; it grows until the machine runs out of memory.

Notation

We use the following throughout.

The kernel’s three-valued comparisons appear often: could_be_equal, and definitely_greater. Lecture 4 defines them. In Rust excerpts, ? passes an error up to the caller, and GeopResult<T> is either a T or an error with its chain of context.

A word on how this was made

Much of this book’s most useful material comes from AGENTS.md in the repository: a running account of what went wrong, written by those who chased it. When a lecture says that a number was “measured”, it comes from there, or from a capped run made while writing the lecture, and the lecture says which. The kernel’s mistakes are taught here on purpose. They are the fastest way to understand why it is built the way it is.

Questions are welcome on the project’s Discord.

Lecture 1: What a CAD kernel is

This book teaches you how to build a CAD kernel. You will need linear algebra, calculus, and enough programming to read Rust. You do not need to know anything about CAD or computational geometry; we build all of that up.

We teach through one real kernel, Geop, whose code is in the same repository as this book. Whenever we say “the kernel does X”, there is a file and a function you can read, and it is named in the text. The kernel changes, so if the text and the code disagree, believe the code.

A solid in a computer

Imagine a steel block with a hole drilled through it and every edge rounded. You want a computer to hold this block in such a way that it can answer questions about it:

• What is its volume, and where is its centre of mass?
• What is the distance between the hole’s axis and the block’s side, exactly?
• Is the point $(3, 4, 5)$ inside the steel?
• What does the block become if I also cut a slot through it?
• What curve should a milling machine follow along that rounded edge?
• If I make the block 2 mm longer, what happens to everything built on it?

A geometric modelling kernel is the library that holds such a description and answers such questions. It has no buttons and no window. It takes in operations (“extrude this outline”, “subtract that solid”, “round this edge by 3”) and gives out models, measurements, meshes for display and files for other programs.

Commercial CAD systems are built over a small number of such kernels. Parasolid sits under Siemens NX and SolidWorks; ACIS is licensed by Spatial; Open CASCADE is open source and sits under FreeCAD. SMLib is another, and docs/smlib-comparison/README.md compares Geop against it feature by feature. What separates these kernels from each other is rarely the mathematics of a fillet. It is how often the fillet works.

The kernel is distinct from the rest of a CAD program in a useful way:

Geop’s repository contains both, since a kernel is hard to test or to use with nothing around it. But the separation is strict, and we will see below that a script checks it.

Four ways to represent a solid

Let us compare the serious candidates by asking each the same questions.

Triangle meshes

A mesh is a list of triangles. Games, film and 3-D printers run on them. Meshes are simple, uniform, and very fast to draw, and the algorithms are short.

But a mesh only approximates a curved surface. A circle of radius $r$ replaced by a regular $n$-gon misses by the sagitta $r(1-\cos(\pi/n))\approx r\pi^2/(2n^2)$. A hole of radius 5 drawn with 16 segments is $0.096$ off, which in a machine shop is an enormous error. Getting to $10^{-6}$ on a radius of 50 takes about 15 700 segments per circle (exercise 2), and a surface needs the square of that. There is no finite mesh of a cylinder from which the radius can be read back exactly, nor one that is tangent to anything exactly. And where two facetted surfaces meet, the intersection is a jagged polyline with thousands of tiny segments, whose next boolean will be worse.

Voxels and distance fields

A voxel grid stores inside or outside per cube of space. A signed distance field stores the distance to the surface at each point, or a formula for it. Booleans become trivial: the union of two solids given by distance functions $d_1,d_2$ is $\min(d_1,d_2)$. Fillets can even be had by smoothing that minimum. The trouble is that the field’s sharp edges are lost at the resolution of the grid, nothing about the result is exact, and a question like “what is the axis of this hole?” has no answer in the representation, since a hole is not a thing in it. These representations are excellent for some purposes and they are growing, but mechanical parts are defined by exact dimensions and exact relationships.

Constructive solid geometry

CSG stores the solid as a tree: (box - cylinder) ∪ other_box. It is compact and exact, and point classification is easy: evaluate the tree at the point. But a tree does not know its own surface. To draw it, measure it, or pick “this edge” in it you must evaluate the tree, which means computing where those solids’ surfaces meet. That is exactly the hard part, deferred. There is no such thing as “round this edge” in a tree, because a tree has no edges.

Boundary representation

A boundary representation (B-rep) describes a solid by its skin. It has two layers:

• Geometry: exact mathematical surfaces and curves. “This face lies on the cylinder of radius 5 about this axis.” “This edge lies on this circle.”
• Topology: how the pieces connect. “These twelve edges bound these six faces. This edge is shared by exactly those two faces. These faces enclose a solid.”

A cube is 8 vertices, 12 edges and 6 faces. A face’s surface extends without limit in principle; the loop of edges around the face trims it to the part that is actually the solid’s skin. So the top of our block with a drilled hole is one plane, an outer loop of four edges, and an inner loop of one or more edges around the hole.

Topology is what makes questions about the part discrete. You can point to a face, name it, ask for its neighbours. “All edges of this face” is an exact finite list, not something that has to be recovered from a cloud of triangles. Since the geometry is exact, measurements and tangency are exact too. B-rep is the representation of essentially all mechanical CAD. Its one defect is that its operations are hard, and this book is about that.

Look at a real one. This is a box with every edge rounded, as Geop builds it. The picture is a mesh, because a screen draws triangles, but the mesh is only a rendering. The model behind it has faces lying on planes, cylinders and spheres.

The mesh you see has 1328 triangles (-O quality=8). With the kernel’s default quality it has 17 328 (we ran both, see Lecture 2), from the same model. That default mesh’s volume, computed by summing tetrahedra, is $2.85566$. The exact volume of this shape is $2.85620$, which you can compute in closed form (exercise 3). Our mesh’s flat facets lie inside the curved surface, and that is the difference. A B-rep’s volume is an integral over exact surfaces, not a sum over an approximation: Model::mass_properties in core/geop-core-topology/src/mass.rs turns it into line integrals around the faces’ trim curves by Green’s theorem.

What makes it hard

All the algorithms we will meet are, at their mathematical core, not very complicated. A boolean difference can be explained on a blackboard in ten minutes: find where the boundaries of $A$ and $B$ cross, cut them there, keep the pieces of $\partial A$ outside $B$ and the pieces of $\partial B$ inside $A$. So why are industrial kernels the work of decades? There are three hard problems, and they feed each other.

Problem 1: decisions from inexact numbers

The computer holds coordinates in finite precision (Lecture 3). But the output of a kernel is mostly not a number. It is a decision: this curve crosses this face, this point is on this boundary, these two vertices are the same. Decisions are discrete, and they have to be consistent with one another.

Suppose three operations each ask whether point $P$ lies on plane $E$. If the answer comes from $|\text{distance}|<10^{-9}$, then a slight change in how $P$ was computed changes the answer. A boolean that believes in one place that an edge crosses a face, and in another place that its end point lies on the face, has built something that isn’t a solid. The result of inconsistency is rarely a slightly inaccurate model. It is a crash, a hole in the shell, an infinite loop, or a model that looks fine and fails three operations later. Accuracy can degrade gracefully; topology cannot.

Problem 2: intersection

The central geometric question is “where do these two things meet?”. For two lines or two planes, a formula answers it. For two free-form surfaces, the answer is a space curve with no closed form, of degree up to hundreds if you write it algebraically. It may have several disconnected branches, closed loops that touch no boundary, cusps, or isolated tangent points, and the two surfaces may coincide along a stretch. And whatever algorithm you pick has to find all of it, because a missed branch means a wrong solid, and do so with the inexact numbers of problem 1. Part III of this book is about this.

Problem 3: keeping topology and geometry consistent

An operation changes a B-rep through many small steps: split this edge, imprint this curve in that face, drop these faces. After each one the model has to remain a meaningful thing. Every edge lies on both its faces’ surfaces. Every loop closes. Each vertex is on all its edges’ ends. Yet the geometric facts in that list are never exactly true in finite precision: an edge computed as the intersection of two surfaces does not lie exactly on either. So what does “lies on” mean? Traditional kernels answer with a tolerance stored per entity. As we will see in Lecture 5, Geop answers differently. Parts IV and V are about this.

There is also a fourth problem, which belongs to parametric CAD and not to kernels as such: when the block gets longer and the model is rebuilt, “this edge” in a later operation must still mean the same edge. This goes by the name of the topological naming problem, and Lecture 2 shows how Geop handles it.

An example of how it goes wrong

Take a box and cut from it a slightly rotated box such that one of the tool’s edges passes very nearly through an edge of the first. The two edges cross at a shallow angle. Computing the crossing’s position is ill-conditioned, so the answer is uncertain along the edges by far more than across them. Suppose the piercing search reports two crossings where the truth has one, separated by $10^{-13}$. Now an edge of length $10^{-13}$ exists. Nothing can be located along it. A face bounded by it may have no interior point at all. The boolean’s classifier asks for a point inside each face, finds none, and returns an error about a face far from any near-coincidence. The failure report names the distant symptom rather than the cause.

Nearly all the rules this book teaches are about never getting into such a state, or at least about finding out at the source if we do.

Geop’s position

Geop takes distinct positions on each hard problem. They are its reason for existing, and each gets a part of the book.

Every number is an interval. A coordinate is not a value $x$ but an enclosure $[x_{lo},x_{hi}]$ that provably contains the true value, and a comparison answers “definitely”, “definitely not” or “could be”. No epsilon ever decides anything. This is Lecture 4 and the rules of Lecture 5, and it is the answer to problem 1. It is what a tolerance would be if it were computed rather than chosen.

Every curve and surface is a NURBS. No line class, plane class, cylinder class or torus class. A plane is a degree-1 NURBS surface. A cylinder is a surface that is a rational quadratic one way and linear the other. No algorithm branches on what kind of surface it was given. From AGENTS.md, the repository’s rulebook:

Every curve and surface stays a NURBS, and no algorithm branches on what kind of surface it was given: no plane × cylinder intersection formula, no offset of a revolution through its profile, no torus fillet because the faces happen to be coaxial.

This is an unusual choice. A traditional kernel has one analytic surface type per common shape, with intersection formulas for many pairs, and falls back on a general NURBS code path for the rest. Those formulas are fast and exact, and every one of them is a separate code path with its own corner cases: the plane × cylinder formula has separate branches for a plane parallel to the axis, tangent, perpendicular and oblique. The mature kernels are, as AGENTS.md puts it, full of such branches, and Geop means to “port the configurations they handle (tangency, vanishing faces, seams), never the type dispatch”. When the general path is wrong or too wide for a case, it is the general path that gets fixed. Since intersection and tracing are fast for linear cases anyway, a plane special case buys little. The kernel is meant to stay small: “a compressed B-rep kernel, simple rather than special-cased”.

Where something does need to know a shape, as a user’s sketch does when it wants to project a circular edge as a circle rather than a spline, the shape is recognized from the control points, locally, at the place it is used (as_line, as_arc, as_plane in core/geop-core-geometry/src/shape.rs), and never stored as a kind in the B-rep.

Subdivide to isolate, Newton to refine. Intersection is not a matter of solving formulas, but of a global search that uses the NURBS control points’ hull to discard regions with proof, and a local polish with Newton’s method. Parts II and III.

A part is a program. The stored object is not a model but a list of steps that build it, each step with arguments that refer to what earlier steps built by stable names. Lecture 2.

Refuse by name. When an operation can’t handle a configuration it says which entities and which condition, instead of returning a subtly invalid solid. For example, the shell, asked to do what it does not support, says that “face … is taken away where it meets face … at an inward corner, along edge …, and the walls grow outward: not supported” (lib/geop-lib-shell/src/shell.rs).

How we got here

DEVELOPERS.md opens with a candid list of what did not work. It is worth reading because each attempt rejected something you may be tempted by.

1. Geop 1 used lines, arcs, spheres, tori and spirals as primitives, with explicit intersection formulas. It “didn’t scale as equations became too complex”. Interval arithmetic was already there, “and 2d booleans worked reliably”.
2. Geop 2 went into algebra: symbolic math, ideals, Gröbner bases. It was “too slow and too much overhead”, and the decision was made to use only numerical methods from then on.
3. Geop 3 focused on numerical methods, with first proofs of concept for intersecting NURBS and solving nonlinear equations.
4. Geop 4 “got stuck in numerical solvers and topology. Booleans looked promising but tracing was too difficult, too many edge cases.”
5. Geop 5 went back to subdivision only, with fat lines and fat surfaces in place of convex hulls, intersection curves never traced explicitly, and no topology, only a bag of bounded patches.

The kernel you are reading about differs in turn from Geop 5: it keeps the fat-line clipping (Lecture 13), but it does have full topology and it does trace intersection curves (Lecture 16). The lesson is not that any of these ideas was foolish, but rather how long it takes to learn which complexity is essential.

A tour of the repository

Geop is a Rust workspace in four tiers, with a web front end on top. The rule of the layout is that a crate’s place follows from what it knows.

As of this writing, that is 7 core crates, 7 lib crates, 21 op crates and 3 cad crates, and about 250 000 lines of Rust, tests included (a line count of core, lib, ops and cad).

Dependencies point down only. An operation is a leaf: it may use core/, lib/ and ops/geop-ops-model, never another op crate. Thus a change to the fillet cannot break the shell. scripts/check-layering.sh checks these rules in continuous integration, and its allow-list is empty. Composition happens in one place, cad/geop-cad-base/src/registry.rs, which calls one register per op crate, in one line each. Lecture 40 shows what such a crate looks like.

A map of every crate is in Appendix C. For the bottom of the stack, here is how the layers of knowledge stack up, with the part of this book that teaches each.

a part as a program           Part IX     core/geop-core-program, core/geop-core-editor
operations                    Part VII    ops/*, lib/geop-lib-sweep, lib/geop-lib-shell
booleans                      Part V      lib/geop-lib-booleans
sketches                      Part VI     core/geop-core-sketch
B-rep topology                Part IV     core/geop-core-topology
intersection                  Part III    core/geop-core-geometry/src/intersection
NURBS curves and surfaces     Part II     core/geop-core-geometry
numbers                       Part I      core/geop-core-math

What Geop does not do

Honest accounts include what is missing. docs/cad-gap.md compares Geop with Fusion 360, Creo, Onshape, NX and SolidWorks, and marks each capability done, partial or missing. ROADMAP.md states its yardstick: the features an engineer reaches for when designing a real machine such as a robot, leaving out “CAM, CAE, rendering, PDM”. Individual operations refuse configurations they do not support, and we will mention those as we come to them. Geop is also younger than the kernels it is compared with: docs/smlib-comparison/README.md measures itself against about 870 000 lines of mature C++ and lists where it is behind as well as ahead.

The plan of this book

• Part I, numbers you can trust (Lectures 3 to 5): floating point, intervals, and the rules that keep interval answers honest.
• Part II, curves and surfaces (6 to 11): Bézier curves, B-splines, NURBS and exact circles, subdivision, surfaces, and evaluating them accurately far from the origin.
• Part III, intersection (12 to 16): the global search, fat-line clipping, Newton, tangency and overlap, tracing curves.
• Part IV, topology (17 to 20): the B-rep’s data structure, poles and seams, point classification and validation.
• Part V, booleans (21 and 22): how all of the above turns into the operation everyone has heard of.
• Part VI, sketches (23 to 25): solving constraints, and proving the solution.
• Part VII, operations (26 to 35): sweep, loft, fillet, shell and the rest.
• Part VIII, assemblies (36 and 37): mates and the solver behind them.
• Part IX, the program model (38 to 40): runs, incremental builds, the editor and how to write an operation.
• Part X, the craft (41 and 42): performance and debugging, from real case studies.

Read it in order. The lectures have a goal box at the start, and a recap and exercises at the end; do the exercises, since some of them prove things the text then uses. The demos are real computations: drag things, and see whether what you see is what the text said.

Exercises

1. Count the vertices, edges and faces of a box with a square hole through it, top to bottom. Check the Euler–Poincaré formula $V-E+F-(L-F)=2(S-G)$, where $L$ is the number of loops, $S$ the number of shells and $G$ the genus (the number of through holes). Why does the formula need $L$, when $V-E+F=2$ works for a plain box?

2. A circle of radius $r$ is replaced by the inscribed regular $n$-gon. Show that the maximum distance between them is $r(1-\cos(\pi/n))$, and how many segments it takes to stay within $10^{-6}$ for $r=50$. How many triangles does that imply for a sphere of that radius, roughly?

3. The rounded box above is the set of points within distance $r=0.2$ of a box of sides $x=1.6,y=1.1,z=0.6$. Derive its volume by adding the pieces: the inner box, six slabs, twelve quarter cylinders, eight sphere octants. Check that you get $2.8562$. Why is the mesh’s $2.85566$ less than this, always?

4. In a CSG system, how would you decide whether a point lies inside $(A\cup B)\setminus C$, given a point-inside test for each primitive? Now list what you would need to know to draw the edges of the same solid. Which of those can the tree give you without intersection computations?

5. Predict from the tier rule which of these a crate in ops/ may depend on: core/geop-core-topology, lib/geop-lib-booleans, ops/geop-ops-sketch, cad/geop-cad-base. Then read scripts/check-layering.sh and check. Where does a function go that both the fillet and the shell need?

6. A traditional kernel stores a tolerance on each edge: “this edge lies within $10^{-5}$ of its faces”. Think of one way this could be wrong to compare with an interval enclosure. (Think about who chooses the number, and what happens when edges with different tolerances meet at a vertex.)

Lecture 2: A part is a program of steps

Before we descend to numbers and geometry for the next many lectures, let us see the whole system working once. We will take one small part, read the whole file that defines it, and run the real kernel on it. Everything we say here will be taken apart later.

Why store a program?

An engineer designs a bracket: a plate, four holes, a rounded corner. Next month the plate has to be 2 mm longer, and the holes should stay 10 mm from the edge. If all you have saved is the resulting solid, you have 40 faces, and you must find and shift a dozen of them by hand. And the fillet’s surface has to change along with them.

So mechanical CAD stores the history: the steps that made the part, with their dimensions. Change a dimension and run again. This is called parametric modelling, and every serious system does it. Geop takes this idea further than usual and makes it literal: the document is a program, and the solid is what that program outputs. There is no other model file. Well, there is a cache next to it, but that can be deleted at will.

A program being plain data brings consequences, all good:

• A .geop file is JSON text. It diffs in version control and can be written by a script or a language model. There is an MCP server for the latter (geop mcp, docs/mcp.md).
• Undo is keeping old copies of the list.
• Rebuilding after an edit can reuse everything not affected (Lecture 38).
• “What does this part mean” has one answer, the same in the browser, on the command line and in a test.

From core/geop-core-program/src/program.rs, a program is little more than this:

pub struct ProgramStep {
    pub id: String,
    pub operation: String,
    pub args: serde_json::Value,
    pub active: Option<String>,
}

pub struct Program {
    pub steps: Vec<ProgramStep>,
    pub state: StateValues,
    pub display: Display,
}

A step has an id you choose, the kind of operation it runs, and its arguments as they are saved. The optional active is a formula deciding whether the step runs at all, so a program can have configurations. state holds where placed parts are and how joints are turned, which matters in assemblies (Lecture 36). display is how the user chose to show what the program builds. Both are left out of the file when empty.

A whole program, read line by line

This is the example rounded_box. We wrote it to disk with geop examples, and you can see the source of it, which is a Rust function writing the same three steps, as rounded_box in ops/geop-ops-fillet/src/examples.rs. Below is its whole .geop text, elided only where it repeats.

Step 1: a sketch

{
  "steps": [
    {
      "id": "outline",
      "operation": "add_sketch",
      "args": {
        "plane": { "datum": "origin:xy" },
        "sketch": {
          "points": {
            "0": { "x": 9.857072060165695e-19, "y": 7.191076644733477e-19 },
            "1": { "x": 1.9999999999999998,    "y": 8.217957318638042e-19 },
            "2": { "x": 1.9999999999999998,    "y": 1.5 },
            "3": { "x": 1.1259343236666495e-19, "y": 1.5 }
          },
          "curves": {
            "4": { "type": "Line", "start": 0, "end": 1, "construction": false },
            "5": { "type": "Line", "start": 1, "end": 2, "construction": false },
            "6": { "type": "Line", "start": 2, "end": 3, "construction": false },
            "7": { "type": "Line", "start": 3, "end": 0, "construction": false }
          },
          "constraints": {
            "8":  { "type": "Fix", "point": 0, "x": 0.0, "y": 0.0 },
            "9":  { "type": "Horizontal", "line": 4 },
            "10": { "type": "Horizontal", "line": 6 },
            "11": { "type": "Vertical", "line": 5 },
            "12": { "type": "Vertical", "line": 7 },
            "13": { "type": "Length", "curve": 4, "value": 2.0 },
            "14": { "type": "Length", "curve": 5, "value": 1.5 }
          },
          "next_id": 15
        }
      }
    },

(We have reordered the keys and joined lines for reading; the real file has them sorted.)

The first step, with id outline, adds a sketch on the plane origin:xy. A part starts with a datum called origin, whose planes one can sketch on. A sketch has points, curves between them, and constraints among them. Each has an id, and they share one counter (next_id). Here are four lines around four corners. The constraints say what the sketch means: corner 0 is fixed at the origin, two lines are horizontal, two vertical, and two lengths are 2.0 and 1.5. That is exactly enough to determine all four corners, with no freedom left over (Lecture 23).

Now look at the numbers in points. Corner 0 is at $x=9.857\times10^{-19}$. Corner 1 is at $x=1.9999999999999998$, not 2. These are the values where the constraint solver left them when the example was made (solved in ops/geop-ops-model/src/sketch/figures.rs calls Sketch::solve). Nothing was wrong: the solver stopped when the constraints held to its accuracy. This small untidiness in the very first file we opened is a preview of Part I, which explains why you must never test whether x is 2.0, and of what to do instead.

Step 2: an extrusion

    {
      "id": "box",
      "operation": "extrude",
      "args": {
        "sketch": { "target": "/:sketch:outline" },
        "extent": {
          "side1": { "blind": 1.0 },
          "reversed": false,
          "symmetric": false
        },
        "output": { "solid": { "mode": "new_body" } }
      }
    },

The second step, box, extrudes the sketch. "sketch": {"target": "/:sketch:outline"} is a reference to what an earlier step built, and it is worth parsing. sketch:outline names a target, a named product of a step, here the sketch made by the step outline. The prefix /: says “of the part itself”. An element of a part placed in this one would read /bolt:… instead, so a picker never needs to ask where an element is from (Lecture 36).

extent says how far: blind, a distance of 1.0 on one side. output says what to do with the result: make a new body rather than add it to or subtract it from an existing solid. If mode were a boolean, this is where the extrusion would call the machinery of Lecture 22.

Step 3: a fillet

    {
      "id": "round",
      "operation": "fillet",
      "args": {
        "edges": [
          "extrude(box,outline,c4,start)",
          "extrude(box,outline,c4,end)",
          "extrude(box,outline,c5,start)",
          "extrude(box,outline,c5,end)",
          "extrude(box,outline,c6,start)",
          "extrude(box,outline,c6,end)",
          "extrude(box,outline,c7,start)",
          "extrude(box,outline,c7,end)",
          "extrude(box,outline,p0)",
          "extrude(box,outline,p1)",
          "extrude(box,outline,p2)",
          "extrude(box,outline,p3)"
        ],
        "radius": 0.2
      }
    }
  ]
}

The third step rounds twelve edges by 0.2. Read the names of those edges, because they are the main idea of this lecture.

Names say how a thing came to be

Every vertex, edge, face and solid of a Geop model has a stable text name, of the form

kind(operation, arg, ...)

where kind is the operation that created it, operation the id of the program step that ran it, and the arguments identify the entity within that step. Our extrusion of the sketch outline in the step box makes:

So the twelve edges in the fillet step are the four edges around the bottom, the four around the top, and the four uprights. No one had to label them. They are built from the program’s own vocabulary: c4 is curve 4 of the sketch, p0 point 0. And outline in the name is a scope (Namer::scoped in ops/geop-ops-model/src/names.rs).

The rule, from the docs of ops/geop-ops-model/src/lib.rs:

A name says how an entity came to be, never when: it is built only from inputs that stay the same when a part is rebuilt after an edit upstream — operation ids, sketch element ids, the names of the entities an operation consumed, and positions counted along those. Never from an internal id or the order in which an algorithm happened to create things.

The topological naming problem

Why this much fuss? Inside the model, entities do have ids: FaceId(42) and so on. Consider a step whose arguments said “round edge number 17”. Now the user adds a hole in an earlier step. The hole’s edges get made first, ids shift, and edge 17 is now something else. The fillet moves to a nonsensical place, or silently changes the design. Or say a boolean rearranges a model so that the same face is now split in two. Which piece is the “same” face?

This is the topological naming problem. Every parametric system runs into it, and it is notorious among users of systems that handle it badly. Ids are positions in a history, and an upstream edit changes the history. Geometry cannot serve as an address either (“the edge nearest $(2,1.5,0.5)$”), because the whole point of the edit is to move the edges.

What remains is a description of derivation. “The edge swept by point 0 of the outline in the step box” remains true whatever the height of the box. Each operation has to define how the names of its outputs derive from the names of its inputs. For a boolean, for example (Lecture 22 has the full table), with N standing for the step’s name:

• a piece of edge E starting at vertex V is named N(E,V), and the first piece keeps the name E;
• an edge traced along faces F1 < F2 from P to Q is N(F1,F2,P,Q);
• a vertex where edges E1 < E2 cross is N(E1,E2,i,n), the i-th of n crossings along E1.

There are some traps in making it work. Three of them, from AGENTS.md:

• Order must not reach a name. A model’s entities live in HashMaps, whose iteration order differs from one run to the next. A corner planner took a corner’s edges in that order, and the ball’s patch began at a different contact on every build. The test caught it by building a program, saving it, loading it and building it again. Wherever order matters, sort by something stable.
• Names may grow. A piece of a split edge is named after the edge and the vertex it starts at, and that vertex after the edge again, so a name doubled per cut. Twenty cuts made names of tens of megabytes. An argument longer than LONGEST_ARGUMENT is now spelled by its digest.
• A trimmed face is still that face. A boolean splitting a face gives one piece the name N(F,E) and the other keeps F, whichever one the splice happened to leave it on. If the boolean then keeps only one piece, that piece is the face, trimmed, and takes the name back (name_lone_pieces, in lib/geop-lib-booleans/src/naming.rs). Otherwise a later draft “about the bottom face” could no longer find it because of an irrelevant detail of the construction.

Typed references

A reference in Rust is not a bare string. When a step refers to a target’s element, its type says which target and which kind of element: ElementReference<NamedModel, Edge> is an edge in the model target. Here is the head of the fillet’s arguments, from ops/geop-ops-fillet/src/operation.rs:

pub struct FilletArgs {
    /// The edges to round, all of one solid.
    pub edges: Vec<ElementReference<NamedModel, Edge>>,
    /// The radius: a number, or a formula of the part's parameters.
    pub radius: NumericInput,
    // ... end_radius, vertex_radii, setback
}

A reference of the wrong kind doesn’t compile. In the file it is saved in the shortest unambiguous form. There is only one model in a part, so a reference into it is saved as the element’s name alone; a sketch is saved as {"target": "/:sketch:outline"}. Note too that radius is a NumericInput, “a number, or a formula”: the 0.2 in the file could have been "thickness / 2", had the part defined parameters.

The types live in core/geop-core-program/src/element.rs. The crate defining a target (say the sketch vocabulary in ops/geop-ops-model) defines the kinds of element it holds, through the trait HasElements. The core crate knows no specifics: core/geop-core-program knows nothing about models and sketches at all.

Running it

The geop command line is the shortest way to run a program. It is built from cad/geop-cad-cli. Install it with the script on the home page or from a checkout by cargo install --path cad/geop-cad-cli, or just open the example in the web app, in which everything we do here can be done by clicking.

Here is the session we ran for this lecture, with a release build of the kernel at version 0.12.0. You have to save the program above as rb.geop first; geop examples --only rounded_box writes it for you.

$ geop compile rb.geop -o rb.stl
3 steps, 866484 bytes -> rb.stl

Time: 0.7 s, for 17 328 triangles. The bytes are $84+50\cdot17\,328$, the format of a binary STL file. Other extensions select other exporters: step, for instance. AGENTS.md asks that every run be capped in memory and time, so in the repository this looks like timeout 120 prlimit --as=6000000000 target/release/geop compile rb.geop. We will say why in Lecture 41.

Changing a dimension

Now the point of the whole exercise. Change the extrusion from "blind": 1.0 to "blind": 2.0, and nothing else:

$ geop compile tall.geop -o tall.stl
3 steps, 866484 bytes -> tall.stl

The twelve edge names have kept their meaning: the fillet found the same edges of a taller box and rounded them. We measured both meshes:

The difference between the mesh volumes, $2.9655$, is the mesh’s estimate of the rounded rectangle’s area, whose exact value is $3-(4-\pi)(0.2)^2=2.9657$. Both meshes happen to have 17 328 triangles, which is why the byte counts agree. The exact values come from exercise 3 of Lecture 1.

Renaming a step

A step’s id is part of every name it creates. If you rename the step box to body, you must rename the references too:

{ "id": "body", "operation": "extrude", ... },
{ "id": "round", ... "edges": ["extrude(body,outline,c4,start)", ...] }

This builds the identical solid: same bounding box, same volume in our measurement. There is no magic behind a name. It is text, chosen to say something true.

When a step fails

What about a design the kernel can’t build? Set the radius to 2.0, more than half the box’s size:

$ geop compile big.geop -o big.stl
error: GeopError
Backtrace: disabled backtrace
RootError: the blend is too large for face FaceId(42): it would meet the face's surface beyond its end
Context:
where the corner's ball touches face FaceId(42), at Vector { data: [Scalar { lo: -7.99360577730113e-14, hi: 7.99360577730113e-14 }, Scalar { lo: 1.9999999999999398, hi: 2.0000000000000613 }, Scalar { lo: -1.0169642905566435e-13, hi: 1.0191847366058938e-13 }] }
Context:
rounding the corner at vertex "extrude(box,outline,p1,start)"
Context:
fillet(round, FilletArgs { edges: [edge("/:model", "extrude(box,outline,c4,start)"), ...
Context:
making target `model` in step `round`
Context:
step "round"

Read this carefully, because it is the kernel’s main way of talking to you. It is a chain, from the root cause up through the stack of what was being done:

• The root says what went wrong and why (“the blend is too large for face … it would meet the face’s surface beyond its end”).
• The first context gives the numbers. Note that the point is not printed as (0, 2, 0) but as an enclosure of it: $x\in[-8\times10^{-14},8\times10^{-14}]$. That is what Part I is about. Those are not rounding-error garbage; each coordinate is a proven bound on the true value.
• The next ones say which corner, by stable name, then which step and which of its targets.

We will explain this style of error, and why it is permanent code and not a debugging aid, in Lecture 5.

If instead a reference names something that isn’t there, say extrude(box,outline,p9):

RootError: no entity is named "extrude(box,outline,p9)" (names of "box": ["extrude(box)", "extrude(box,end)", "extrude(box,outline,c4)", ...

The error lists names a step of that id did make. So a refusal comes with enough information to fix the program, which is especially useful when the one reading it is a script or an agent. In an editor, a step that fails does not stop the others: later steps that do not depend on it still build, and the run lists all failures in program order (RunOutput in core/geop-core-program/src/runner.rs).

Also notice that next to each program the run left a file ending in .geopc. That is the cache of what the program built, read back on the next run when everything it depends on is the same (Lecture 38). It can be deleted at any time and is never committed to the repository.

Other ways to run a program

All of these drive one and the same editor (core/geop-core-editor) over one registry of operations:

• geop examples writes every built-in example, as program and mesh. --only NAME selects some.
• The web app at app.geop-cad.dev runs the kernel as WebAssembly in a browser. A link https://app.geop-cad.dev/?example=rounded_box opens that example. Lecture 39 shows how.
• geop open DIR edits the files of a folder in the browser, with the kernel running natively on your machine.
• geop mcp runs the editor as a Model Context Protocol server, so an agent can add steps, read which ones failed, look up names and export the result. Tools and protocol are in docs/mcp.md. Setting it up is one command: claude mcp add geop -- geop mcp.
• From Rust, Program::push appends a step, which is how the examples are written. See rounded_box in ops/geop-ops-fillet/src/examples.rs.

You will find it helpful to keep one of these open while reading. Each lecture from here on can be tried, by changing a number and seeing what happens.

Exercises

1. Without running anything: what are the names of the top cap, the wall along the sketch’s line 6, and the edge between those two in the program above? What name would the side wall have had if the extrude’s step id were slab?

2. Open the example in the web app, or use geop compile on a modified file. Change "value": 2.0 in constraint 13 to 3.0. Which of the names in the fillet step change? Why does the fillet still work, though every wall moved?

3. Suppose a kernel names edges by the index they got on creation, and you insert a hole before the fillet in rounded_box. Describe what happens. Now suppose it names them by their midpoint’s coordinates and you change the box’s width. Describe what happens. What property does a good naming scheme need that both fail?

4. Can stable names be fully stable? Suppose you delete sketch curve 4 and the point it shares with curve 5, and connect things differently. What happens to extrude(box,outline,c4,start)? What should the kernel do: guess, or refuse? What does Geop do?

5. Read validate_operation_id in core/geop-core-program/src/program.rs. Why can a step id not contain a comma or a parenthesis? Give an example of two different entities whose names would collide if it could.

6. The sketch’s stored x of 1.9999999999999998 is one unit in the last place below 2 (you will define that in the next lecture). Suppose that a later step, written by a careless programmer, tested if x == 2.0. What would happen on this file? And what if the solver had happened to stop exactly on 2.0 in the previous version of the kernel?

Lecture 3: Floating point

Geop stores coordinates in f64, IEEE 754 double precision. Before we can build the rest of Part I on top of it, we need to know exactly what those numbers are and what they do. Most programmers carry around a vague picture (“floats are slightly inaccurate”). That picture leads directly to the wrong fix, namely comparing with an epsilon. The right picture is much more precise, and much more useful.

What a double is

A f64 is 64 bits: one sign bit $s$, eleven exponent bits, fifty-two fraction bits. For all but the extremes, the bits stand for the number

$$x = (-1)^s\left(1+\frac{f}{2^{52}}\right)2^{e},\qquad 0\le f<2^{52},\quad -1022\le e\le 1023.$$

So a double is a 53-bit binary integer scaled by a power of two. It is floating in that where the binary point sits depends on the exponent. And it is exact: every bit pattern (besides NaN and infinities) denotes one precise rational number. There is nothing fuzzy about the number 1.9999999999999998 that you met in Lecture 2. It is exactly $2-2^{-52}$.

All numbers with one exponent $e$ are equally spaced, $2^{e-52}$ apart. This spacing is the unit in the last place, or ulp. When the exponent goes up by one, the spacing doubles:

• near 1, numbers are $2^{-52}\approx2.2\times10^{-16}$ apart;
• near 1000, they are $1.1\times10^{-13}$ apart;
• near $10^5$, $1.5\times10^{-11}$ apart;
• near $10^{16}$, 2 apart. The number $10^{16}+1$ does not exist.

So floats are a ruler whose ticks spread out as you go away from zero. Relative to the number’s size the spacing is nearly constant: $\mathrm{ulp}(x)/\lvert x\rvert$ is between $2^{-53}$ and $2^{-52}$. Absolute spacing is not constant at all, and this will decide our question about epsilons. Below is the ruler itself, computed from the bits of the numbers in your browser.

At the two ends: the largest exponent gives about $1.8\times10^{308}$, beyond which an operation yields infinity. At the small end, subnormal numbers fill in the gap between $2^{-1022}\approx2.2\times10^{-308}$ and zero, with fixed spacing $2^{-1074}\approx5\times10^{-324}$. Zero has a sign, and -0.0 == 0.0.

Most decimals are not representable

Binary fractions terminate only for denominators that are powers of two. $\tfrac1{10}$ in binary is $0.000110011001100\ldots_2$, repeating, so the source text 0.1 becomes the nearest double. To all its digits, these are the actual values:

0.1 is   0.1000000000000000055511151231257827021181583404541015625
0.2 is   0.200000000000000011102230246251565404236316680908203125
0.3 is   0.299999999999999988897769753748434595763683319091796875

Now add the first two. The exact sum of those two doubles lies between two representable numbers, and rounds to the nearest one:

0.1 + 0.2 is   0.3000000000000000444089209850062616169452667236328125

That is not the double called 0.3, which is one ulp lower. This is why 0.1 + 0.2 == 0.3 is false, in every language. Nothing was wrong with the addition; it was as accurate as it can be. The inaccuracies were both there already: 0.1, 0.2 and 0.3 are all slightly off, in different directions.

The standard model of rounding

IEEE 754 demands that +, -, *, / and sqrt be correctly rounded: the result is what you would get by computing the exact real result, and then rounding it to the nearest float. This makes float arithmetic amenable to analysis. If $\circ$ is one of those operations and $\hat\circ$ its float version, then (barring overflow and underflow)

$$x\mathbin{\hat\circ}y = (x\circ y)(1+\delta),\qquad \lvert\delta\rvert\le u=2^{-53}\approx1.1\times10^{-16}.$$

The number $u$ is called the unit roundoff. Each operation, taken alone, is nearly perfect, with a relative error of at most $u$. So where do bad answers come from? From what errors do when they meet later operations. Functions such as sin and cos, by contrast, come from the math library and are not guaranteed correctly rounded at all. This matters for later: the kernel’s interval sin has to allow for it.

Accumulation and order

Consider $s_n=s_{n-1}+x_n$. Each addition contributes a relative error up to $u$ of the running sum. Summing $n$ positive numbers therefore has an error bound growing like $nu$ times the sum. The bound suggests the order matters: put the largest numbers last, so that small additions do not suffer rounding against a large running total. In the extreme, adding 1 to $10^{16}$ changes nothing at all, because the result rounds back:

$$(10^{16}+1)-10^{16}=0,\qquad(10^{16}-10^{16})+1=1.$$

So float addition is not associative. Even 0.1 + 0.2 + 0.3 depends on the order: we computed (0.1+0.2)+0.3 = 0.6000000000000001 and 0.1+(0.2+0.3) = 0.6.

At $N=10^6$, the true value is $14.39272672286\ldots$. In f32, largest first gives $14.3574$, wrong in the fourth digit. Smallest first gives $14.39265$, right to five digits. The two runs have the same numbers and the same operations, and differ only in order. (This demo is a calculation you could repeat in any language; book/scripts/demo_tests/float-lab.test.mjs compares it with an independent one.)

For a kernel, the lesson goes beyond summing. Two ways of computing the same value give different floats. So an evaluation by one route is not == to an evaluation by another, even if there has been no mistake at all, and even if the formula says they are equal.

Cancellation

The worst thing that happens to an honest computation is that it subtracts two nearly equal numbers. Each number has an error of perhaps one ulp. Subtraction itself is exact in this situation (a theorem by Sterbenz), but the leading digits vanish and what remains is mostly those errors, now promoted to the leading position. The relative error is no longer $u$, but $u$ times how much the magnitude has shrunk.

As an example, take

$$f(x)=\frac{1-\cos x}{x^2}.$$

Its limit as $x\to0$ is $\tfrac12$. Let us try it: at $x=10^{-8}$, $\cos x$ is $1-5\times10^{-17}$, closer to 1 than half an ulp, so it rounds to 1. The numerator is zero: $f(10^{-8})=0$. At $x=2\times10^{-8}$, the result is $0.5551$. No digit is right.

By the half-angle identity, $1-\cos x=2\sin^2(x/2)$, so the same function is $2\sin^2(x/2)/x^2$. That formula has no subtraction, and it gives $0.5$, accurate to the last bit, at those same arguments.

Rewriting a formula to avoid subtracting nearby numbers is a major technique. We will see it again in Lecture 11: a rational patch evaluated in world coordinates loses accuracy with its distance from the origin, and the fix is to subtract a nearby point first, where it can be done exactly, and add it back after.

Why no epsilon is correct

Now to the central claim. Since computed values are inexact, a natural thought is: to test whether $a$ equals $b$, test

(a - b).abs() < EPSILON

for some small constant. This looks like a reasonable engineering solution, and it is the source of much of the misery in geometric software. Let us see why, in three steps.

There is no right size. An epsilon has units, but the rounding error is relative to what was being computed and how. Take EPSILON = 1e-9.

• At coordinates near $10^8$, the ulp is $1.5\times10^{-8}$, greater than epsilon. The test is the same as a == b, and any rounding error at all makes it fail.
• A value that accumulated a thousand operations at coordinates near $10^5$, where the ulp is $1.5\times10^{-11}$ and errors can reach $10^{-8}$, fails the test though the true values are equal.
• For a part measured in nanometres, $10^{-9}$ could be most of a feature.

Use a relative epsilon, $\lvert a-b\rvert<\varepsilon\max(\lvert a\rvert,\lvert b\rvert)$? It is then useless near zero, where $10^{-20}$ and $-10^{-20}$ are “different” for any $\varepsilon<2$, though both are obviously rounding noise. In the file of Lecture 2, a point’s x was stored as $9.857\times10^{-19}$, standing for 0. And the magnitude of $a$ tells nothing about the size of the error anyway: $a$ could be the small difference of two huge numbers.

The kernel’s own experience makes the point. As AGENTS.md records, a STEP edge lay $2\times10^{-31}$ above its plane face, and that mattered. At the other end, a quarter arc cut forty times once came out $0.08$ wide. No single constant is both smaller than the first and larger than the second.

It is not transitive. Equality has a property you rely upon without saying so: $a=b$ and $b=c$ imply $a=c$. With an epsilon they do not. With EPSILON = 1e-9, take $a=0$, $b=0.9\times10^{-9}$, $c=1.8\times10^{-9}$. Then a “equals” b, b “equals” c, and a does not “equal” c.

So suppose three vertices with those coordinates lie along an edge. A boolean merges coincident vertices. If it merges a with b, then b’s neighbour c is no longer within epsilon of the merged vertex, or is within epsilon, depending on which one is kept. Thus the topology depends on the order of the merges, and the order came out of a hash map. This failure cannot be fixed with a better number.

It conflates two different questions. “Are these the same?” has three honest answers: definitely yes (exact), definitely no, and I cannot tell at the accuracy I have. An epsilon test collapses the third answer into one of the first two, without telling you. Any later code therefore acts on a belief it has no grounds for. That loss is the real harm, and Lecture 4 is about not losing it.

Geometry on floats: the orientation test

Here is the classic demonstration that this is not just about equality. The most basic predicate of computational geometry asks on which side of the line through $q$ and $r$ the point $p$ lies. The answer is the sign of a determinant:

$$\operatorname{orient}(p,q,r)=\operatorname{sign}\bigl((q_x-p_x)(r_y-p_y)-(q_y-p_y)(r_x-p_x)\bigr).$$

It takes two multiplications and five subtractions. If $p$, $q$ and $r$ are nearly collinear, the two products are nearly equal, and cancellation strikes. Kettner, Mehlhorn, Pion, Schirra and Yap collected examples of what follows in Classroom examples of robustness problems in geometric computations, from which the configuration below is taken.

Fix $q=(12,12)$ and $r=(24,24)$, and let $p$ take each of the $256\times256$ values $(0.5+i\,\mathrm{ulp},\,0.5+j\,\mathrm{ulp})$, with $\mathrm{ulp}=2^{-53}$. The true line is $y=x$. Thus the true sign is the sign of $j-i$: the diagonal is collinear, and each side has its own sign. To be careful, we compute the exact sign too, in rational arithmetic, using BigInt.

What you are looking at is the geometry of a line, as seen by floating point: a jagged mosaic, with blocks of wrong signs far from the true diagonal. We counted the failures: on the $16\times16$ grid, 114 of 256 points get the wrong sign, or a false “collinear”. On the $256\times256$ grid it is 11 972 of 65 536. Note also that it can’t be fixed by a threshold on the determinant. Any threshold just turns some of the wrong signs into “zero”, and it also turns right signs into “zero”.

Why does this break algorithms and not merely give small errors? Because a geometric algorithm builds on such predicates as axioms. A convex hull algorithm relies on the fact that if $p$ is left of $qr$ and $r$ is left of $pq$, certain other things follow. When the predicates are inconsistent, the algorithm’s invariants fail. The cited paper shows incremental convex hull algorithms that then lose points or produce a polygon that isn’t convex. The output isn’t a slightly wrong hull; it is not a hull.

What can be done about it

There are three broad strategies, and each CAD kernel is a compromise among them.

1. Exact arithmetic. Rational numbers or algebraic numbers keep each comparison exact. Shewchuk’s adaptive predicates make this efficient for determinants of points. But coordinates of exact intersection points of NURBS are algebraic numbers of very high degree. DEVELOPERS.md records that Geop 2 tried algebra, “ideals, Groebner bases”, and found it “too slow and too much overhead”.

2. Tolerances. Decide that things within $\varepsilon$ are the same. A kernel stores a tolerance with each entity and tries to keep the decisions consistent. Most classical kernels do so, and years of engineering have made them work remarkably often, at the cost of the problems shown above. docs/smlib-comparison/README.md describes SMLib as “a mature NURBS B-rep kernel built on tolerances”.

3. Enclosures. Keep the inexact numbers, but compute along with every value a rigorous bound on how inexact it may be, and let comparisons say when they can’t tell. This is Geop’s approach, and the next lecture.

Exercises

1. Show that 0.1 in binary is $0.0\overline{0011}_2$. What is $\mathrm{ulp}(0.1)$? Check the actual double’s error from the text’s expansion against half of it.

2. What are the following in f64, and why? (a) 1e16 + 1.0 - 1e16 (b) 1e16 - 1e16 + 1.0 (c) 1.0 + 1e-16 == 1.0 (d) 1.0 - 1e-16 == 1.0.

3. The smaller root of $x^2-10^8x+1=0$ is about $10^{-8}$. Compute it by the textbook formula $\frac{b-\sqrt{b^2-4c}}{2}$ with $b=10^8$, $c=1$. Then derive a formula with no cancellation. How many digits does each give?

4. A program merges vertices closer than $10^{-6}$, and processes vertices in the order given. Give four vertices on a line such that the result depends on the order in which they are processed. How would you make it independent of order? What would your fix do on a long chain of vertices each $0.9\times10^{-6}$ from the next?

5. Change the grid in the orientation demo to q=(3,3), r=(100,100) and to p=(0.1,0.1), q=(1.7,1.7), r=(2.9,2.9). Is the mosaic still there? Is the true line still the diagonal, in the second case? (Careful: the double called 1.7 is not exactly 1.7.)

6. A surface is evaluated at a parameter, in two ways, giving $z_1$ and $z_2$, both within a few ulps of the same value. Suppose your code selects the larger for one purpose and the smaller for another. Could the topology you build depend on the platform’s math library? On the compiler’s optimization level, on hardware with fused multiply-add? What does that imply for a test suite comparing files between machines?

Lecture 4: Interval arithmetic and enclosures

The previous lecture ended with a problem: a computed float is wrong by an amount it does not tell you, and comparisons based on it are guesses. The remedy is older than CAD. Moore’s Interval Analysis dates from 1966. Rather than one number, carry two: a lower and an upper bound on the true value. The bounds are themselves floats, and all arithmetic is arranged so that the truth always stays between them.

Enclosures

Let $x^*$ be a real number we would like to know: the height of some surface at some parameter. An enclosure of $x^*$ is an interval $[x_{lo},x_{hi}]$ with

$$x_{lo}\le x^*\le x_{hi}.$$

A narrow enclosure is more useful than a wide one, but a wide enclosure is still true. This reverses the usual relationship between numerics and correctness. Ordinarily an answer of $3.14159$ is more or less accurate, and you don’t know how much. An enclosure is $[3.14,3.15]$, and it is either valid or a bug. All the imprecision is out in the open, and it is measured, as its width $x_{hi}-x_{lo}$.

The arithmetic

To add two enclosures: if $x\in[a,b]$ and $y\in[c,d]$ then $x+y\in[a+c,b+d]$. Each operation asks for the least and greatest of $x\circ y$ over the box.

$$\begin{aligned} [a,b]+[c,d]&=[a+c,\;b+d]\\ [a,b]-[c,d]&=[a-d,\;b-c]\\ [a,b]\cdot[c,d]&=[\min(ac,ad,bc,bd),\;\max(ac,ad,bc,bd)]\\ [a,b]/[c,d]&=[a,b]\cdot[1/d,\;1/c]\qquad\text{if }0\notin[c,d] \end{aligned}$$

Subtraction crosses over: the lowest difference pairs the lowest $x$ with the highest $y$. Division by an interval containing zero is undefined, since the quotient could be anything; we return an error rather than something that merely looks like an answer. Monotone functions like $\sqrt{\cdot}$ take their end points to the result’s end points. Non-monotone ones need more care, $\sin$ for instance has to include 1 if the interval passes through $\pi/2$.

A few by hand:

The inclusion property

Why trust this? Because of one simple theorem, which is why interval analysis works at all.

Theorem (inclusion). Let $f$ be a function written as an expression in $+,-,\cdot,/$ and elementary functions, each variable replaced by an interval. Let $F(X_1,\dots,X_n)$ be the result of evaluating that expression with the interval operations above. Then for all $x_i\in X_i$, $$f(x_1,\dots,x_n)\in F(X_1,\dots,X_n).$$

Proof. By induction over the expression tree. It holds for a variable. For a node $A\circ B$, the induction hypothesis says that the true values of $A$ and $B$ lie in the computed intervals, and by construction of $\circ$, every $a\circ b$ for $a,b$ in those intervals lies in the output. $\square$

The theorem’s quiet power: it holds for every point in the box at once. Evaluate $f$ on $[0,1]$ and you have a claim about uncountably many inputs, using a finite number of operations. If the result excludes zero, then it is proved that $f$ has no root in $[0,1]$. No sample can do that. This is what makes subdivision searches in Part III possible: they evaluate over a whole region, and may discard the whole region with certainty.

Outward rounding

In a computer, $a+c$ must itself be rounded to a float. If it rounds downward for the upper bound, then the truth might fall outside. So lower bounds must round down and upper bounds up: outward rounding. Inclusion then still holds, whatever the rounding errors were.

IEEE 754 hardware offers rounding modes for this, but switching modes is awkward and slow, and WebAssembly, where Geop also runs, only has round-to-nearest. Geop does something simpler. Since round-to-nearest is within half an ulp of the exact result, stepping the result by one whole float outward is guaranteed to enclose it:

pub fn add(self, other: Self) -> Self {
    Scalar::new(next_down(self.lo + other.lo), next_up(self.hi + other.hi))
}

This is Scalar::add, from core/geop-core-math/src/scalars/mod.rs; Scalar is simply a struct of two public fields lo and hi, both f64. The other operations are the same with the formulas above. So:

Each row deserves a remark.

Row one. The exact sum lies at $0.3000000000000000166$, between the two bounds, so the enclosure is true. It isn’t the tightest interval possible, since rounding by a full ulp each side covers more than needed, but it’s cheap and it is what the kernel uses. And it does tell you what 0.1 + 0.2 did not: that the result is not to be compared with equality.

Note that from_f64(0.1) means exactly the double nearest 0.1, a sharp interval of width zero: the kernel takes that number at its word as the value meant. If you mean the true rational $1/10$, Scalar::from_ratio(1, 10) encloses it: one ulp to each side of the rounded quotient.

Row two. Multiplying by exactly 1 widened the value by an ulp each side. Rounding outward cannot tell that the operation was exact. This seems harmless, but as AGENTS.md records, rewriting extrude and revolve as one sweep builder multiplied every station row by a weight of exactly 1, and that alone broke six “inscribed cylinder” boolean tests, whose geometry was identical. The rule from there: never multiply or add by exact constants needlessly in construction code. The kernel has Scalar::is_exactly partly for this: “A weight that is exactly one is skipped rather than multiplied by, since the product would only be widened by rounding.”

Row three. $0.7-0.7$ is exactly 0, yet the result has width $10^{-323}$. In next_up, the step from zero goes to the smallest subnormal. So that result is not is_zero() (which needs both bounds to be 0 exactly), but it could_be_equal zero.

Functions from the math library get wider treatment, because they are not correctly rounded. Scalar::sin and Scalar::cos evaluate at both ends, step outward by an ulp, and recognize the extrema inside the interval structurally rather than numerically, through contains_periodic, then they clamp to $[-1,1]$. Scalar::atan2 is more interesting: it takes the bounds in plain f64 at the corners of the box, and proves them afterwards in outward-rounded arithmetic, moving a bound outward if the proof fails. A good pattern to remember: a float is a free choice of where to put the bound, as long as you prove it holds.

Scalar::div fails with “divisor interval contains zero” instead of returning an unbounded interval. A computation that loses its denominator in rounding errors therefore comes back as an error, with context, at the place it happened.

Three-valued comparison

Given two enclosures, what can we say about the true values? These are the kernel’s comparison operators, all one line each:

definite_sign sums it up for one number: Some(true), Some(false), or None if it could be zero.

Try them on some cases:

The operators give a logic with three outcomes in effect: yes, no, undecided. Two different kinds of question get asked of it, and the difference is crucial:

• “Could this curve cross this face?” If the answer is “could be”, the search must keep looking: discarding a region without proof could lose a real intersection.
• “Is this definitely inside?” If the answer isn’t “definitely”, you may not build on it, and need to subdivide, try another point, or stop and report.

Undecided is a real answer that must be handled. There is no if a < b { } else { }: you always have to think about the third case. This is more work to write, and exactly what makes the code correct in the bad cases. In the orientation demo of Lecture 3, the intervals give no wrong sign at any of the 65 536 positions of the largest grid. Where f64 gave a confident wrong answer, they say “cannot tell”.

You will also notice operations which are not comparisons:

• Scalar::width returns $hi-lo$ as a sharp value, so that a comparison of widths is always decided;
• Scalar::lower and Scalar::upper return the bounds as sharp values, for cutting a domain at the outer bound of an enclosure, which loses nothing;
• Scalar::sharpen collapses to the midpoint. It is dangerous, and has the next lecture to itself;
• Scalar::ENTIRE, $(-\infty,\infty)$, is the interval that could be anything. It stands in for a value not yet known;
• Scalar::max, Scalar::min and Scalar::clamped are honest versions of their namesakes: clamped moves a value to a bound only where it is definitely beyond it.

The dependency problem

Interval arithmetic has one great weakness, and much of the craft is about managing it. What is $x-x$ for $x\in[1,2]$? It is plainly 0. But the formula gives

$$[1,2]-[1,2]=[1-2,\;2-1]=[-1,1].$$

Interval subtraction treats the two operands as independent. The formula’s two appearances of $x$ might as well be different numbers, each of which may independently be anywhere in the interval. All correlation between them is forgotten. The result is true but needlessly wide. This is the dependency problem.

A second example: let $f(x)=x(1-x)$, for $x\in[0,1]$. The true range is $[0,\tfrac14]$. Interval evaluation gives $[0,1]\cdot[0,1]=[0,1]$, four times too wide. Now try it yourself.

There are three standard ways out, and the kernel uses all of them.

Rewrite the formula so each variable occurs once. If each variable appears only once in an expression, then interval evaluation gives the exact range (rounding aside). For $x(1-x)=\tfrac14-(x-\tfrac12)^2$, with a proper square that knows it’s non-negative, you get $[0,\tfrac14]$ on $[0,1]$. When you can, this is the cheapest fix. A related lesson is in Scalar::interpolate, with the two forms of linear interpolation:

$$a+\alpha(b-a)\qquad\text{and}\qquad(1-\alpha)a+\alpha b.$$

Mathematically the same; numerically opposites. The first mentions $a$ twice and counts its width twice, but $b-a$ is computed first. So if $a$ and $b$ are the same value, as with a weight of 1 that has to stay 1 through many subdivisions, the difference vanishes, whatever the width of $\alpha$. The second mentions each of $a,b$ once, so it doesn’t double-count wide control points, but splits $\alpha$ into two uncorrelated factors, so it can’t see the cancellation. Neither dominates. So the kernel computes both, and returns their intersection (see below). AGENTS.md recounts this as the fix to the first root cause of a real curve_curve_intersect blowup.

Mean value form. By the mean value theorem, for any $x,c\in X$, there is a $\xi$ between them with $f(x)=f(c)+f'(\xi)(x-c)$. Thus $$f(X)\subseteq F_{mv}(X)=f(c)+F'(X)\,(X-c).$$ Here the variable occurs in only two places, and $X-c$ is small and centred about zero. Worked example, for $f=x(1-x)$, $X=[0.4,0.6]$, $c=0.5$:

• $f(0.5)=0.25$
• $F'(X)=1-2X=[-0.2,0.2]$
• $X-c=[-0.1,0.1]$
• So $F_{mv}=0.25+[-0.02,0.02]=[0.23,0.27]$, width $0.04$

The true range is $[0.24,0.25]$, width $0.01$, and naive evaluation gives $[0.4,0.6]\cdot[0.4,0.6]=[0.16,0.36]$, width $0.2$.

How do they converge? Halve $X$ to $[0.45,0.55]$. Mean value gives $0.25\pm0.005$: width $0.01$, true width $0.0025$, excess $0.0075$, a quarter of what it was ($0.03$). Naive gives $[0.2025,0.3025]$: width $0.1$, excess $0.0975$, about half of what it was ($0.19$). In general:

• naive interval evaluation has an excess width linear in the width of the input;
• the mean value form has an excess quadratic in the width of the input.

This is why NurbSurface3D::evaluate in core/geop-core-geometry/src/nurb_surface/evaluate.rs intersects its direct enclosure with the mean value form (“the middle plus the partials over the box times the half widths”), as told in the section “A bound is as local as what it bounds” of AGENTS.md.

Subdivide. Split $X$ into pieces, evaluate on each, and take the union. For $X=[0.4,0.6]$ in two pieces: $[0.4,0.5]\cdot[0.5,0.6]=[0.2,0.3]$ and $[0.5,0.6]\cdot[0.4,0.5]=[0.2,0.3]$, so the union is $[0.2,0.3]$, half as wide as naive. Doubling the pieces halves the excess. It is expensive when used alone, but subdivision has another virtue: where the pieces are far from a solution, they can be thrown away entirely. That is Part III.

Rationals show the dependency problem in an especially costly form. Take $x/(1+x)$ on $[1000,1001]$ (the demo has this function, “x / (1 + x)”, over smaller values). Numerator and denominator move together, but the naive formula pairs the smallest numerator with the largest denominator: $[1000/1002,\;1001/1001]=[0.99800,1.0]$, when the true range is $[0.999001,0.999002]$. A NURBS curve evaluates to exactly such a ratio, $A/W$, and its width then scales with the distance from the origin. Lecture 11 is entirely about this one effect.

Union and intersection

An interval is a set, so it has lattice operations. Scalar::union is the smallest interval containing both. Scalar::intersect is their overlap. Be aware of what the kernel’s intersect does when they do not overlap: it returns the narrower input rather than an inverted interval. Its documentation says that callers must only intersect two enclosures of the same exact value. They then must overlap unless something has gone wrong.

Think about when each is right:

• Two enclosures of the same value, coming from different formulas: the truth is in both. So it is in the intersection, which is tighter than either. That is what Scalar::interpolate does with its two forms.
• A quantity that may be either of two things, or one value with two estimates that must both be covered: the union.
• Averaging: never. Lecture 5 explains.

The wrapping effect

The dependency problem has a geometric sibling. A point in space is a vector of three intervals, that is a box aligned with the axes. Rotate that box by 45 degrees. The true set is the same box, turned. But the only thing a vector of intervals can represent is an axis-aligned box, and the smallest one around a turned unit square is $\sqrt2$ times as wide. Repeat this, and the width multiplies by $\sqrt2$ per step, even though the true set never grows.

This effect is why interval methods work best with computations that aren’t long chains of transformations, and why construction code in a kernel thinks carefully about which frame to compute in. When growth of a width looks mysterious, check for wrapping.

Exercises

1. Evaluate by hand: $[-2,3]\cdot[-1,4]$, $[1,2]/[2,4]$, $[1,2]^2$ as a product, and $[-1,2]^2$ as a product. What is the true range of the last, and how would you write a proper “square” operation?

2. Prove inclusion for $f(x,y)=(x+y)x$, by following the theorem’s induction for $X=[0,1]$, $Y=[1,2]$. Compare with the true range. Then give an expression for the same function with less overestimation on that box.

3. Use the mean value form to enclose $f(x)=x^3-2x$ on $[0.9,1.1]$, with $c=1$, and compare with naive evaluation and the true range. Use the demo to check.

4. Why does Scalar::interpolate return a directly if a.is_exactly(b) and alpha lies in $[0,1]$? Why isn’t computing a + alpha*(b - a) enough? (Read the comment in the code.) What’s the risk if this were not done across forty repeated subdivisions?

5. Enclose $g(x)=(x+1)/(x+2)$ on $[100,101]$ naively. Then rewrite it so that $x$ occurs once and compare. What general trick for ratios is this?

6. A box $[-1,1]^2$ is rotated by $\theta=10°$ in an interval computation, 36 times, by the formula $x'=x\cos\theta-y\sin\theta$, $y'=x\sin\theta+y\cos\theta$, for a full turn. How wide is the result, ignoring rounding? How could you rotate 36 times by $10°$ without that?

Lecture 5: The rules of honest numerics

Intervals give us a number that does not lie. But a program can still make it lie: narrow an enclosure past what was established, compare it with a fudge factor, or ask it a question whose answer doesn’t apply to it. This lecture is a list of rules, each with the bug that taught it. They come from the section headings of AGENTS.md, and the stories, and numbers, are from there. You’ll meet most of the code involved again later, in its proper place. What matters now is the reasoning.

The rules boil down to one sentence:

We are always seeking the mathematical truth: bounds should be as tight as they can legitimately be made, but a bound must never be narrowed past what the computation actually established.

Keep that in mind. Dropping uncertainty does not make a result more accurate. It makes it wrong in a way nothing downstream can detect, and in a geometry kernel the concrete consequence is a missed intersection.

Rule 1: never use an epsilon

No (a - b).abs() < EPSILON, and no new constant to paper over a comparison. Lecture 3 showed why it cannot be right, and Lecture 4 what to write instead: could_be_equal, definitely_greater, and so on.

There are tunable numbers in the kernel, though. A search has a node budget (max_nodes) and a size below which it stops subdividing (min_subdivision_size). The distinction is stated like this: these numbers are allowed to affect how hard the search tries. More budget or a tighter tolerance changes the running time, and how quickly a search converges. It does not change what a converged answer means. A search that hasn’t converged returns an error, None or a wide interval, never a value that is merely “close enough”.

RemeshParams in lib/geop-lib-booleans/src/remesh/mod.rs collects all the numbers of a boolean. Read its comments: “Node budget for every subdivision search”, “An effort budget only: a trace that reaches its vertex does so in the same steps whatever the budget.”

So how do you tell a budget from an epsilon? Ask: if I changed this number, could a correct result become incorrect? If yes, it’s an epsilon in disguise. And running out of budget has to be an error. The geometry crate’s page says it plainly: “Running out of max_nodes is always an error, never read as ‘not contained’: the search is incomplete, so neither answer is justified.”

There is one number that looks like a tolerance and isn’t: MAX_WIDTH in core/geop-core-topology/src/validation/numerical_accuracy.rs, $10^{-4}$, the widest a coordinate of any stored entity may be. No comparison is made against it to decide geometry. It is a bound on how much uncertainty an entity may carry, because beyond it the searches can’t locate anything. Past that, validation reports a bug, typically “a missing refinement”.

If a result comes out wider than expected, the fix is to understand why, not to shrink a number until the symptom disappears. That is Rule 8.

Rule 2: union to cover, intersect to constrain, never average

Suppose a point on an intersection curve of two surfaces is computed twice, once by projecting onto each surface. You get two enclosures of the same point:

$$A=[1.000,1.002],\qquad B=[1.001,1.005].$$

What should you do with them? There are three candidates.

• Average the midpoints: $1.002$. This is a sharp value that claims more precision than either input had, and discards the most useful signal in the pair: how far apart the two answers were. That gap is the honest measure of how much the computations disagree, and later convergence tests need to see it.
• Union: $[1.000,1.005]$. Every value either computation allows. It keeps the gap as width.
• Intersection: $[1.001,1.002]$. Every value both computations allow.

Which of the last two is correct depends on what you know:

Union when you must cover both possibilities. A point on a face × face curve, projected once onto each surface, is combined this way: the gap between the two projections is what later “have I arrived?” tests need to see. merge_vertex (core/geop-core-topology/src/edit/merge_vertex.rs) combines vertex points this way too.

Intersect when both constraints must hold at once. If each is a valid enclosure of the very same number, then the truth lies in both, so in the intersection, and we rightly keep less width. Example: splice_edge_into_face pins each end of a new loop to the arriving end’s $(u,v)$ intersected with the vertex’s projection onto the surface. Using the end alone let the width of interpolate_enclosing spread around the loop until a face split came out degenerate.

Average never. It is not an enclosure of anything.

Rule 3: sharpen only where the value is a free choice

Scalar::sharpen collapses an interval to its midpoint. Its documentation warns that it is for “when you are free to pick any value within self”. It is the sharpest tool for violating honesty, and has a rule of its own. It is legitimate only where any value inside the interval would serve equally well, so that collapsing to one loses no mathematical accuracy. There are three such places:

1. Subdivision parameters. A search that splits an interval may split it anywhere; the split point is the algorithm’s choice, not a claim about geometry. Sharpening here is necessary, too: repeatedly re-deriving a split point from a wide domain compounds width through knot insertion, with $\alpha=(t-a)/(b-a)$, a span divided by a shrinking span, until the search is useless.
2. The predictor, before the corrector. A predicted point is a seed: the corrector converges to the same fixed point from any seed nearby. See the comment “Sharpened exactly once, here” in predictor_corrector_step (core/geop-core-geometry/src/intersection/branch.rs).
3. Where to cut an edge at a vertex, once refined. Every parameter of the refined enclosure ends a piece on the curve and within the vertex. Cut at a sharp one and check it. Cutting at all of them, by contrast, makes each piece wider than its parent. AGENTS.md measured a quarter arc cut forty times going from $3\times10^{-15}$ wide to $0.08$, where it now adds only one subdivision’s rounding per cut (NurbCurve::split_parameter_at_point).

Everywhere else, sharpening a value that is an answer is a bug. The first big example: NurbSurface::project, which finds the point of a surface nearest a target, sharpened every Newton iterate, including the last. For all but the last that is correct, since they are only seeds. But the final step’s du, dv inherit the target’s width through the residual, and that width honestly states how precisely an uncertain target pins down a foot point. Sharpening returned one sharp $(u,v)$ for a target that only determined a range, and every later comparison against anything else derived from the same uncertain point failed. Fixing that one thing took validate_fast failures across the scene sweep from 21 of 175 to 1 of 175.

A sharpened answer lies to every interval test downstream

The worst of it is that the damage shows up three modules away. predictor_corrector_step sharpened its iterate at every Newton step, including the one it returned. A traced intersection curve running exactly along a face’s trim boundary therefore came back with a sharp $(u,v)$ sitting $6\times10^{-17}$ off the boundary. So face_contains classified it Inside, since could_be_equal had no width to work with. The direction check whose job is to reject exactly such a curve waved it through. The curve was spliced in as a duplicate of the boundary, carving a zero-area sliver. Much later, this surfaced as a face that a boolean could not classify, because it had no interior point.

Notice what the sharpening did not do: it caused no geometric error. The point was where it should be, to within $6\times10^{-17}$. What it destroyed was the statement of uncertainty. A value that is honestly $[0,1.2\times10^{-16}]$ answers “could this be on the boundary?” correctly. The same value sharpened to $6\times10^{-17}$ answers it wrongly, with total confidence.

Keep two copies

The fix generalizes. Whenever a loop needs sharp seeds but produces an answer, keep two copies:

• the seed, sharpened, which starts the next step. Sharpening here is a free choice, and without it, 20 Newton steps compound width until evaluate fails outright (measured at 36 of 175 scenes);
• the answer, the same update with its width as computed, which is what you return.

The returned width is that of a single step from a sharp seed, so it is bounded per step rather than compounded or thrown away. Here’s the shape of it, from the loop in NurbSurface::project (core/geop-core-geometry/src/nurb_surface/project.rs):

|&(u, v)| (clamp((u.sharpen(), v.sharpen())), clamp((u, v))),

The pair is (seed, answer). core/geop-core-math/src/solvers/iterate.rs holds the loop that all of the kernel’s Newton-type methods share, and its module comment summarizes the idea: “Which of the two a caller returns says whether its value is a free choice … or an answer that later interval tests compare against.”

Now see it for yourself, in a real computation, with the kernel’s own rounding. We solve $(x+s)^3=c$ by Newton’s method, in interval arithmetic, where $c$ is computed as $(b+s)^3$ for a given $b$. So the true solution is exactly $b$. Then we ask the three-valued question of face_contains: is the answer inside, outside, or on the boundary $b$?

We ran this for $b=0.7$, $s=0.3$. The answer sharpened is 0.7000000000000001, one ulp above $b$, and says “outside”: wrong. Two copies give an interval $1.3\times10^{-15}$ wide around it, which says “on boundary”: right. Across a thousand values of $b$ in $[0.5,2)$, sharpening everything gives a wrong classification for 87 of them, and two copies for none, in every one of which the interval contains $b$. If instead the iterate is used as its own seed, width doubles at every step, and at step 50, for $b=0.7$, the divisor could be zero and the iteration fails with the iterate about $10^4$ wide.

Note what the “two copies” answer is, and is not. A single Newton step from a sharp seed encloses the Newton image of that seed, which differs from the root by second order. Here the rounding width covers that difference once the seed has converged, as the survey’s “contains” count confirms. But it’s not a proof that the root lies inside. For a proof you need a test like Krawczyk’s (Lecture 24). The kernel takes the same stance: the width states how well the data pins the answer down, and the point where existence must be established uses a proof.

Rule 4: validate the value you are about to use, not a wider one

Model::split_edge_at_vertex once checked curve.evaluate(edge_t).could_be_equal(&vertex_point) on the wide edge_t, a fat box that genuinely contains the vertex, and then split at the sharpened edge_t. The check answers “is the vertex somewhere along this stretch of curve?”, while the split asks “cut exactly here”. Those are different questions, and the gap between them is real geometry: the true crossing $t^*$ lies somewhere in edge_t, the midpoint generally isn’t $t^*$, and $C(t_{sharp})$ is a different point from $C(t^*)=$ vertex, off by $\lvert t_{sharp}-t^*\rvert\,\lvert C'(t)\rvert$.

Note what is not wrong there: the two halves are exact sub-arcs. What’s wrong is asserting a topological identity (“this split point is that vertex”) that the numbers never established. So when you sharpen a parameter, re-validate the sharpened one, or better, refine it until the identity is true, by Newton on $C(t)=$ vertex. That is NurbCurve::split_parameter_at_point today.

Do not patch a mismatch by moving geometry to match the assertion. Pinning a split curve’s end point onto the vertex was tried, and looked as if it worked for the 3-D curve. Applied to the pcurve, whose $(u,v)$ came from a search min_subdivision_size wide, it pinned a control point to a box so wide that the pcurve had no usable tangent (“Cannot normalize zero”), and took the sweep from 0 failing scenes of 175 to 12. Fix the parameter, not the geometry.

Rule 5: ask a question about a box only where its answer holds for a box

The converse trap. face_interior_point_where (core/geop-core-topology/src/contains/face.rs) tests a box around a candidate, so that the accepted point keeps clear of the boundary. Vertex and pcurve checks answer that honestly for a box. But the ray casting that came after built its ray from the box, which turned the ray into a strip. A strip merely grazing a pcurve overlaps it just as a crossing does, so that gives one hit, and the parity flips. A point $0.44$ from a hole’s axis, with the rim at $0.33$, was found “inside” the disc cut from a plate’s top. The boolean kept that disc, and the result was left open.

Parity counting only means something for a ray from one point. Once the box is known to touch no boundary, all of it lies on one side, and one sharp point of it decides which. So before you run an interval computation on a wide input, check that its answer still means something for every point in the input.

The inclusion theorem of the previous lecture guarantees an enclosure of $f$ over the box. It doesn’t guarantee that what you did with that enclosure, here counting, makes sense over a box.

Rule 6: proximity is not membership

Something is selected because it is near, and then used as if it belonged. A traced intersection curve ends when it reaches a known vertex, and candidate_within chose the nearest vertex within one marching stride. But a traced curve lies on both surfaces it is traced between, so its end vertex must lie on both too, and an oversized solid puts many unrelated vertices within a stride of the curve’s end. Adopting one of those ended the edge at a point on neither face. fit_pcurve then had nothing to project onto, and silently returned the nearest foot point, about $8\times10^{-3}$ away, which finally surfaced as 32 validate_fast errors about pcurve end points.

The fix requires what the geometry actually asserts: surface_could_contain on both surfaces. When you select by proximity, ask what property the code is about to assume, and test that directly.

A diagnostic that pays for itself: when such a check fails, report whether the point lies on the surface at all. “The pcurve names the wrong $(u,v)$ for a point that is on the surface” and “the edge was spliced onto a face it does not belong to” produce identical error text, but need opposite fixes. One extra surface_could_contain in the error path separated them: of those 32 failures, every one had exactly one end point on the surface and one off.

Rule 7: check an iterative answer against everything it was asked

An iteration that stops is not an iteration that converged. The shell moves a vertex onto the offsets of its faces, by Gauss-Newton in meet (core/geop-core-geometry/src/meet.rs). It returned a hemisphere’s pole at $(0,0.9,0.1)$, on the offset plane but off the offset sphere, and the failure surfaced two steps later as “not on the circle”. There were two silent causes, a dropped condition and a projection that never moved. Checking the returned point against every surface it was asked to lie on turned this into an error at the source, naming the surface it missed, with the numbers. It costs one projection per condition. It is what makes the next bug of the kind take minutes rather than an afternoon.

Rule 8: trace numerical bugs to their root cause

When a search misbehaves, there’s a strong pull towards a local workaround: trim the input away from the trouble spot, pad a result “to be safe”, seed the search with a known-good answer, add a special case for the one failing pair. These can make the symptom vanish while leaving the defect in place. Some (padding, seeding) smuggle an epsilon back in by a side door.

Instead, keep descending until the answer is a specific line of arithmetic, not a vague “it’s imprecise near shared vertices”. From AGENTS.md, a real example. A curve_curve_intersect blowup was first “fixed” by trimming curves near a shared vertex. That worked, but only by avoiding the region where the bug lived. The chain as eventually found, one level at a time:

1. Repeated NurbCurve::split calls compounded interval width in a control point’s weight. Fixed by a tighter interpolate formula.
2. The same with knot values divided by an ever-shrinking span. Fixed by sharpening a free-choice value immediately after the operation that introduced the rounding noise, not by padding the result.
3. Only with both fixed did it become clear that a third, independent thing was needed: a domain-width convergence signal. The first two had buried it under noise.

Name the entities before you name the cause

Explanations such as “near-tangent crossing” and “imprecise near a shared vertex” sound plausible, cannot be checked and are often wrong. Before you state why something failed, identify the exact entities at the failing vertex or edge by name, and check the explanation against them.

A revolve cut “up to next” failed with a vertex too wide to validate, and it was explained to the user as a near-tangency. The user looked at the scene: there was none. Naming the vertex showed it was where the drilled hole’s rim crossed the revolve’s end face, and that face was there because the operation’s own fallback had stopped the tool flat at the first point of contact, an existing corner. The bug was a logic error, not numerics.

So first ask whether something our code chose (a stop position, a cap, a split point) put the entities there. Blame numerics only after that is ruled out. Lecture 42 has more cases of this.

Rule 9: debug with context in the real call stack

When you need to know why something fails, do not write eprintln! or a scratch harness. Add context to the real path:

let ctx = |e: GeopError| {
    e.with_context(format!(
        "Model::split_edge_at_vertex(edge_id={edge_id}, edge_t={edge_t}, vertex_id={vertex_id}, max_nodes={max_nodes}, min_subdivision_size={min_subdivision_size})"
    ))
};

let edge = self.get_edge(edge_id).with_context(&ctx)?.clone();

That is real code from core/geop-core-topology/src/edit/split_edge_at_vertex.rs. The context is permanent: it stays and helps the next failure of the kind report itself. You saw what it gives in Lecture 2: the root cause, then a chain of what was being done, with stable names and numbers.

When a cause depends on actual numbers, print them with {:?}. Display on a Scalar truncates to three decimals of the midpoint, silently hiding exactly the width you’re trying to diagnose. GeopError in core/geop-core-math/src/geop_error.rs is the error type.

The sanctioned alternative is a regression test. Every bug fix begins by pinning the failure into its own standalone #[test], red for exactly that reason. Then fix, then require both the reproducer and the broader suite to be green. A scene that fails as one line of a 175-scene log tells you almost nothing, while the same scene as its own test can be run in a second. Fixes that green the reproducer but regress the sweep are not fixes, and it has happened here. Over time, the suite gets one test per real defect.

Rule 10: refuse what you do not support, early and by name

An operation handling only some configurations must recognize the others before it starts building. Otherwise the unsupported case appears deep in a boolean, as a degenerate splice, or worse, as an invalid model that nobody notices. The shell assumed a removed face’s inner copy lies inside the face, which is false where that face meets a kept face at an inward corner. It failed in five different ways before inward_corner (lib/geop-lib-shell/src/shell.rs) recognized that up front:

face … is taken away where it meets face … at an inward corner, along edge …, and the walls grow outward: not supported

A refusal names the entities and the condition, so the user can change the model and the next session can extend the code. This costs a little code. It saves a lot of despair.

Exercises

1. Classify each as a legitimate or illegitimate use of sharpen, and say why: (a) choosing the parameter at which to bisect a domain in a containment search; (b) storing the vertex where a curve pierces a face; (c) choosing the point at which to evaluate a surface’s normal in order to fix a plane for a corrector step; (d) returning the arc length of a curve; (e) evaluating a sample point to describe a problem in an error message.

2. Two computations of a parameter give $[0.30,0.34]$ and $[0.33,0.40]$. Give what you’d return, if (a) both are rigorous enclosures of one number; (b) the same rigorous enclosures, but of two different branches of a curve, and you want a box around both; (c) rigorous enclosures of one number, and the second is really $[0.35,0.40]$. What does (c) tell you, and what does Scalar::intersect do with it?

3. Write the two-copies loop as pseudocode for a one-dimensional Newton iteration with an interval residual f and derivative df. Which of seed and answer does each of these get: the evaluation of f? The value returned? The check of whether the seed has repeated? The intersection with the box that subdivision proved?

seed = start.sharpen()
answer = start
for k in 0..N:
    next = seed - f(seed) / df(seed)
    answer = next
    seed = next.sharpen()
    if seed was seen before: break
return answer.intersect(box)

4. In the sharpen demo, drag $b$ until “sharpen everything” says something other than “on boundary”, then look for a $b$ where it is right. The readout also counts over 1000 boundaries. Why can the answer be off by one ulp, though the problem is well-conditioned?

5. A search returns None when it exhausts max_nodes, which the caller reads as “no intersection”. Give a scenario for a wrong solid, and fix the interface. Then find what the kernel actually does: read the doc comment of ValidationParameters::max_nodes in core/geop-core-topology/src/validation/parameters.rs.

6. Which rule does each of these break? (a) “Project the vertex onto the surface, and if the foot point is within 1e-6, it is on the surface.” (b) “This vertex is the closest one, so it must be the end of the edge.” (c) “We ran Newton for 20 iterations, so (u, v) is the solution.” (d) “The ray from this box crossed the loop three times, so the box is inside.” (e) “The boolean fails near the corner, so we trim the tool 1e-5 away from it.”

Lecture 6: Bézier curves

Everything so far was about numbers. From here on we are concerned with shape. A CAD kernel has to represent a free-form curve exactly, change it cheaply, and above all answer questions about it with proof: does this curve cross that one? Where? Could it reach this point? The whole of the next ten lectures hangs off one remarkable representation that makes those questions tractable, and it starts with nothing more than repeated linear interpolation.

From a line to a curve

Take two points $P_0, P_1$ in the plane. The line segment between them is the set of all

$$(1-t)\,P_0 + t\,P_1,\qquad t\in[0,1].$$

We call this linear interpolation, and the only thing it asks of a space is that you can add vectors and scale them. Now take three points. Interpolating on each leg of the polygon gives two points moving along the two legs, and interpolating between those with the same $t$ gives a third:

$$\begin{aligned} P_0^1(t) &= (1-t)P_0 + tP_1,\\ P_1^1(t) &= (1-t)P_1 + tP_2,\\ P_0^2(t) &= (1-t)P_0^1(t) + tP_1^1(t). \end{aligned}$$

As $t$ runs from 0 to 1, $P_0^2(t)$ traces a parabola. Nothing in the construction is special to three points. Start from $n+1$ control points and repeat, each level one point shorter than the one before:

$$P_i^0 = P_i,\qquad P_i^{k}(t) = (1-t)\,P_i^{k-1}(t) + t\,P_{i+1}^{k-1}(t),\quad k=1,\dots,n,\ i=0,\dots,n-k.$$

The last point $C(t)=P_0^n(t)$ is the Bézier curve of degree $n$, and this procedure is de Casteljau’s algorithm. Play with it before reading on.

Things to notice while you drag:

1. The curve starts at $P_0$ and ends at $P_n$, but does not in general pass through the others.
2. The curve never leaves the shaded hull of the control points.
3. At $t=0$ the curve leaves along the first leg $P_0P_1$, and it arrives along the last leg.
4. Moving one control point moves the whole curve, though least where $t$ is far from $i/n$.

Each of these can be proved, and each of them will be used.

Bernstein polynomials

Unroll the recursion by hand for $n=2$:

$$C(t) = (1-t)^2 P_0 + 2t(1-t)P_1 + t^2 P_2.$$

The coefficients are the terms of $\bigl((1-t)+t\bigr)^2$. This is general: a path through the triangle of the construction from $P_i$ to the top takes $n$ steps, $i$ of them to the right (each with a factor $t$) and $n-i$ to the left (each with a factor $1-t$), and there are $\binom{n}{i}$ such paths. So

$$C(t) = \sum_{i=0}^{n} B_i^n(t)\,P_i,\qquad B_i^n(t)=\binom{n}{i}t^i(1-t)^{n-i}.$$

The $B_i^n$ are the Bernstein polynomials of degree $n$. The demo prints both evaluations under the plot: the two numbers agree, as book/scripts/demo_tests/de-casteljau.test.mjs also checks. Three facts about them carry the rest of this lecture.

They are a partition of unity. By the binomial theorem, $$\sum_i B_i^n(t) = \bigl((1-t)+t\bigr)^n = 1.$$

They are non-negative on the interval. For $t\in[0,1]$, every factor of $B_i^n(t)$ is at least zero.

They are a basis. The $n+1$ polynomials $B_0^n,\dots,B_n^n$ span the polynomials of degree at most $n$, so every polynomial curve of degree $n$ has exactly one set of Bézier control points. A Bézier curve is therefore not a new kind of curve, only a polynomial written in a clever basis.

Why is that basis clever? Because its coefficients are points that have geometric meaning. In the monomial basis $a_0+a_1t+a_2t^2+\dots$ the coefficient $a_3$ is not a location anywhere, and a small change in it can throw the curve far away. Worse, the monomial form is badly conditioned for evaluation and root finding on $[0,1]$; Farouki and Goodman proved that the Bernstein basis is the best conditioned among polynomial bases that are non-negative on the interval. A kernel living on rounding-error margins cares.

What the representation gives you

Convex hull. Since the Bernstein values are non-negative and sum to one, $C(t)$ is a convex combination of the control points, for every $t\in[0,1]$. So the curve lies in the convex hull of its control points.

This single property is the foundation of every intersection algorithm in this book. If the hulls of two curves are disjoint, then the curves do not meet, and that is a proof, not an estimate: it holds for all infinitely many points of the curve, using finitely many numbers. Compare it with sampling, which proves nothing about what lies between two samples (we return to that in Lecture 12).

Affine invariance. The weights sum to one, so $A\bigl(\sum B_iP_i\bigr)=\sum B_iA(P_i)$ for an affine map $A(x)=Mx+b$. Translating, rotating or scaling a curve means transforming its control points, not sampling and refitting. Note the weights must sum to one for the translation part to survive: this will be what the “homogeneous coordinates” of Lecture 8 are about.

End points and end tangents. $B_0^n(0)=1$ and all others vanish, so $C(0)=P_0$; likewise $C(1)=P_n$.

Derivative. Differentiating the Bernstein form and collecting terms gives $$C'(t) = n\sum_{i=0}^{n-1} B_i^{n-1}(t)\,(P_{i+1}-P_i).$$ The derivative of a Bézier curve is a Bézier curve of one degree less, the hodograph, whose control points are the differences of the original ones. At $t=0$ it is $n(P_1-P_0)$, which is observation 3 above. The convex hull applies to the hodograph too: if all differences $P_{i+1}-P_i$ lie in some open half plane, then no tangent of the curve can vanish, or turn back through that direction. In one line of arithmetic we have bounded every tangent of the curve.

Variation diminishing. A straight line crosses the curve at most as many times as it crosses the control polygon. We state this without proof (see Farin in Appendix A). It means that a curve whose control polygon has no wiggle cannot have a wiggle either.

Splitting

Press split at t in the demo. The two polygons drawn along the outside of de Casteljau’s triangle are themselves control polygons, and of the two pieces of the curve exactly:

$$\text{left: } P_0^0, P_0^1, \dots, P_0^n,\qquad \text{right: } P_0^n, P_1^{n-1}, \dots, P_n^0.$$

The right-hand piece, reparametrized to $[0,1]$, is the same point set as the original curve restricted to $[c,1]$. Why? Because a polynomial is determined by its values on any interval, and one can check by induction that the left polygon’s Bézier curve evaluates to $C(ct)$.

Splitting is what makes subdivision work. It costs $O(n^2)$ additions, it is numerically benign (every step is a convex combination, so errors are not amplified), and each piece’s control polygon is closer to the curve than the original’s. Indeed, for pieces of parameter width $h$, the distance between a piece and its control polygon’s chord shrinks as $h^2$. So repeated splitting converges quadratically to line segments, each boxed in by its hull. Everything in Part III is built from this.

A worked example

Let $P_0=(0,0)$, $P_1=(1,2)$, $P_2=(3,2)$, $P_3=(4,0)$ and $t=\tfrac12$.

So $C(\tfrac12)=(2,\tfrac32)$, and by the Bernstein sum $\tfrac18(0,0)+\tfrac38(1,2)+\tfrac38(3,2)+\tfrac18(4,0)=(2,\tfrac32)$. The left piece has control points $(0,0),(\tfrac12,1),(\tfrac54,\tfrac32),(2,\tfrac32)$, the right piece $(2,\tfrac32),(\tfrac{11}{4},\tfrac32),(\tfrac72,1),(4,0)$. The tangent at $t=\tfrac12$ is $3\bigl((\tfrac{11}{4},\tfrac32)-(\tfrac54,\tfrac32)\bigr)=(\tfrac92,0)$, horizontal, as the symmetry demands.

Where one polynomial is not enough

Bézier curves have two limitations that drive the next two lectures. First, moving any control point moves the whole curve: there is no local control. A long curve needs a high degree, and at high degree a control point’s influence becomes weak and global. The fix is to chain many low-degree pieces with an agreed continuity at the joints, which leads to B-splines (Lecture 7). Second, a polynomial cannot draw a circle: $x(t)^2+y(t)^2=1$ has no non-constant polynomial solution (exercise 5). The fix is a rational weight, which leads to NURBS (Lecture 8).

Recap

Exercises

1. Compute $C(\tfrac14)$ for the worked example by de Casteljau, and verify with the Bernstein sum.

2. Prove the derivative formula from the Bernstein form. Start from $\frac{d}{dt}B_i^n=n(B_{i-1}^{n-1}-B_i^{n-1})$.

3. Degree elevation. Show that the same curve has degree $n+1$ control points $Q_0=P_0$, $Q_{n+1}=P_n$, and $Q_i=\frac{i}{n+1}P_{i-1}+\bigl(1-\frac{i}{n+1}\bigr)P_i$. Check it with the demo’s “Add point” button, which uses this formula: the curve does not move.

4. Prove $B_i^mB_k^n=\binom{m}{i}\binom{n}{k}\big/\binom{m+n}{i+k}\,B_{i+k}^{m+n}$. Use it to write down the control points of the product of two scalar Bézier polynomials. This is what Bernstein in core/geop-core-geometry/src/bernstein.rs implements. Why can it do this with no evaluation and no sampling?

5. Show that no polynomial curve traces a piece of the unit circle: there are no non-constant polynomials $x,y$ with $x^2+y^2=1$.

6. Suppose a cubic’s control points are $(0,0),(1,1),(2,-1),(3,0)$. Without computing the curve, bound its $y$ and its slope $dy/dx$ for all $t$. How wide are those bounds, and how would you tighten them?

Lecture 7: B-splines

Lecture 6 ended with a defect. A Bézier curve is one polynomial, so a change to any control point changes all of it, and a long intricate curve needs a high degree. Neither is acceptable in design. An engineer who adjusts the bulge near one end of a wing section does not want the other end to move. And a degree-40 polynomial is a numerically fragile object in any basis.

The fix is obvious in outline: use many low-degree polynomial pieces, joined with an agreed smoothness at their joints. What is not obvious is that this can be done with a single formula, keeping every good property of the Bernstein form, in such a way that the joint constraints hold automatically for any choice of control points. That formula is the B-spline.

Why not simply chain Bézier curves?

You could. Two cubics $P_0,\dots,P_3$ and $Q_0,\dots,Q_3$ meet with a common tangent if $Q_0=P_3$ and $Q_1-Q_0=P_3-P_2$. Asking for a continuous second derivative as well adds a further relation among $P_1,P_2,P_3,Q_0,Q_1,Q_2$. So out of the eight control points of two pieces, only a few are actually free. If a user drags one of them, you must move others to preserve the joints, and the correction propagates along the chain.

What we would like is to work with only the free quantities: a set of numbers such that any values give a curve with the desired smoothness. For a chain of $s$ cubic pieces joined $C^2$, a count gives $4s$ coefficients per coordinate, minus $3(s-1)$ constraints, which is $s+3$ free numbers. So we seek $s+3$ basis functions, each piecewise cubic and $C^2$, such that any weighted sum has those joints. That is exactly what B-spline basis functions are.

The knot vector

Fix a degree $p$ and a non-decreasing sequence of numbers called the knot vector,

$$U = (u_0, u_1, \dots, u_{m}),\qquad u_i\le u_{i+1}.$$

The intervals $[u_k,u_{k+1})$ are the knot spans. They are where the curve is a single polynomial. A value may repeat; the number of times it does is its multiplicity. Repeating a knot will turn out to be how we ask for less smoothness at that joint.

The Cox-de Boor recursion

The basis functions are defined from the knots alone, by recursion on the degree:

$$N_{i,0}(t) = \begin{cases}1 & u_i\le t<u_{i+1}\\ 0 & \text{otherwise}\end{cases}$$

$$N_{i,p}(t) = \frac{t-u_i}{u_{i+p}-u_i}\,N_{i,p-1}(t) + \frac{u_{i+p+1}-t}{u_{i+p+1}-u_{i+1}}\,N_{i+1,p-1}(t),$$

with the convention that a fraction whose denominator is zero counts as zero ($0/0=0$): its partner $N$ is zero too, on an empty span.

Each level is a blend of two neighbours of the level below, with a weight that runs linearly from 0 to 1 across the support. Compare it with $P_i^k=(1-t)P_i^{k-1}+tP_{i+1}^{k-1}$: this is de Casteljau’s idea, with $t$ replaced by a local fraction depending on the knots.

Given $n+1$ control points $P_0,\dots,P_n$, the B-spline curve is

$$C(t) = \sum_{i=0}^{n} N_{i,p}(t)\,P_i.$$

The counting that makes this consistent: $n+1$ basis functions of degree $p$ need $n+p+2$ knots, $u_0,\dots,u_{n+p+1}$. The curve is defined on the domain $[u_p,u_{n+1}]$. Outside it fewer than $p+1$ basis functions are in play and they no longer sum to one.

Try these while reading on. Move a control point far from the red marker and notice which colours of the curve stay put. Drag two knots together in degree 3 and watch the curve acquire a visible crease in curvature. Press “Triple knot”: the curve now passes through a control point with a corner. Set the degree to 6: with seven control points, this is one span, and the basis functions are the Bernstein polynomials.

A small worked example

Take $p=2$ and knots $U=(0,0,0,1,2,2,2)$, so $n+1=4$ control points and domain $[0,2]$. On the first span $[0,1)$, at $t=\tfrac12$, work upward. Only $N_{2,0}(t)=1$. At degree 1:

$$N_{1,1}=\frac{u_3-t}{u_3-u_2}N_{2,0}=\tfrac12,\qquad N_{2,1}=\frac{t-u_2}{u_3-u_2}N_{2,0}=\tfrac12.$$

At degree 2:

$$N_{0,2}=\frac{u_3-t}{u_3-u_1}N_{1,1}=\tfrac12\cdot\tfrac12=\tfrac14,$$ $$N_{1,2}=\frac{t-u_1}{u_3-u_1}N_{1,1}+\frac{u_4-t}{u_4-u_2}N_{2,1}=\tfrac12\cdot\tfrac12+\tfrac34\cdot\tfrac12=\tfrac58,$$ $$N_{2,2}=\frac{t-u_2}{u_4-u_2}N_{2,1}=\tfrac14\cdot\tfrac12=\tfrac18.$$

They sum to $\tfrac14+\tfrac58+\tfrac18=1$, and $N_{3,2}(\tfrac12)=0$. With $P=(0,0),(1,2),(3,2),(4,0)$ we get $C(\tfrac12)=\tfrac58(1,2)+\tfrac18(3,2)=(1,\tfrac32)$.

Properties

Local support. By induction, $N_{i,p}$ is zero outside $[u_i,u_{i+p+1})$. So on a span $[u_k,u_{k+1})$, only the $p+1$ functions $N_{k-p,p},\dots,N_{k,p}$ can be nonzero. Moving $P_i$ affects at most $p+1$ spans of the curve. This is the local control we were after, and it is why evaluation costs $O(p^2)$ however many control points there are.

Non-negative. On its support every factor in the recursion is between 0 and 1, so $N_{i,p}\ge0$.

Partition of unity. On the domain, $\sum_i N_{i,p}(t)=1$. By induction on $p$: expanding the recursion and collecting the two terms in which $N_{i,p-1}$ appears gives the coefficient $\frac{t-u_i}{u_{i+p}-u_i}+\frac{u_{i+p}-t}{u_{i+p}-u_{i}}=1$, so $\sum_i N_{i,p}=\sum_i N_{i,p-1}$, and at degree zero exactly one function is 1.

The strong convex hull. Together the last three give something stronger than Bézier’s: on span $k$, the curve lies in the convex hull of just its $p+1$ active control points, $P_{k-p},\dots,P_k$. A long B-spline is enclosed in a chain of small hulls rather than one big one.

Affine invariance. As before: partition of unity means that transforming the control points transforms the curve.

It contains Bézier curves. Take $U=(\underbrace{0,\dots,0}_{p+1},\underbrace{1,\dots,1}_{p+1})$. The recursion then reproduces the Bernstein polynomials (exercise 2). The kernel uses precisely this: bezier_knots in core/geop-core-geometry/src/nurb_curve/construct.rs is that vector, and a straight line is NurbCurve::line over bezier_knots(1).

Multiplicity and continuity

Each $N_{i,p}$ is a polynomial of degree $p$ on each span. Across a knot of multiplicity $m$, it is $C^{p-m}$. Therefore so is the curve:

This table is a statement about the parameter, not the shape: a $C^1$ joint can still look perfectly smooth, and a curve can even be $C^2$ and have a cusp if its derivative vanishes there. Geometric continuity ($G^1$: the unit tangent is continuous; $G^2$: the curvature is too) is the shape-level notion, and we meet it in Lecture 29.

Multiplicity gives us the vocabulary to reconcile Bézier and B-splines. A knot of multiplicity $p$ at a joint makes the curve pass through a control point there, with each piece on its own. A knot of multiplicity $p+1$ cuts the control points into two independent groups. Raising every interior knot to full multiplicity turns any B-spline into separate Bézier curves, which is how NurbCurve::bezier_pieces (core/geop-core-geometry/src/nurb_curve/split.rs) gets them, splitting at each knot. Knot insertion is the subject of Lecture 9.

Clamped knot vectors

If the first $p+1$ knots are equal, and likewise the last $p+1$, the knot vector is clamped. Then, exactly as for Bézier curves, $C(u_p)=P_0$ and $C(u_{n+1})=P_n$, and the end tangents lie along the first and last legs of the polygon. An unclamped curve starts somewhere inside the hull of its first few control points, a place you cannot read off from any single number.

For a kernel this matters a lot. An edge’s curve has to start exactly at a vertex. If the first control point is the start of the curve, then the vertex is a stored number and no evaluation (and no rounding) is needed. The kernel relies on this in several places: NurbCurve3D::start and NurbCurve3D::end simply divide the first and last control points by their weights, for “the clamped knot vectors every curve here has”, and evaluate_into (core/geop-core-geometry/src/nurb_curve/evaluate.rs) returns the end control point when a parameter pins a clamped end sharply. is_clamped_at in core/geop-core-geometry/src/knot_insertion.rs is the test, and NurbCurve::extended refuses to carry on a curve past an end that is not clamped.

De Boor’s algorithm

The sum of basis functions is the definition, but not how to compute. Just as de Casteljau beats the Bernstein sum, de Boor’s algorithm evaluates by repeated linear interpolation of control points directly.

Let $t$ lie in span $k$, $u_k\le t<u_{k+1}$. Start with the active points, $d_j^0=P_j$ for $j=k-p,\dots,k$, and for $r=1,\dots,p$:

$$d_j^{r} = (1-\alpha_j^{r})\,d_{j-1}^{r-1}+\alpha_j^{r}\,d_j^{r-1},\qquad \alpha_j^{r}=\frac{t-u_j}{u_{j+p-r+1}-u_j},\quad j=k-p+r,\dots,k.$$

The answer is $C(t)=d_k^{p}$. Note that the denominators stretch over $p-r+1$ spans around $t$, and for $t$ in span $k$ every $\alpha_j^r$ lies in $[0,1]$, so each step is a convex combination. On one Bézier span, $\alpha_j^r$ reduces to $(t-a)/(b-a)$ for every $j,r$, and this is de Casteljau’s algorithm.

Run the worked example: $k=2$, $t=\tfrac12$. Level one: $d_1^1=\tfrac12P_0+\tfrac12P_1=(\tfrac12,1)$ (the denominator $u_3-u_1=1$), and $d_2^1=\tfrac34P_1+\tfrac14P_2=(\tfrac32,2)$ (the denominator $u_4-u_2=2$). Level two: $\alpha_2^2=(t-u_2)/(u_3-u_2)=\tfrac12$, so $d_2^2=\tfrac12(\tfrac12,1)+\tfrac12(\tfrac32,2)=(1,\tfrac32)$, as before. The demo’s grey polygons are the levels of this triangle.

Derivatives

Differentiating the recursion and simplifying gives a formula as clean as the Bézier one (Piegl and Tiller derive it in full; see Appendix A):

$$C'(t)=\sum_{i=0}^{n-1} N_{i+1,p-1}(t)\,Q_i,\qquad Q_i=\frac{p}{u_{i+p+1}-u_{i+1}}\,(P_{i+1}-P_i).$$

The derivative is a B-spline of degree $p-1$, on the same knots with the first and last dropped, whose control points are scaled differences. It has the same properties, so the hull of the $Q_i$ bounds every tangent on a span: a finite proof about all the tangents of a curve. For repeated knots the denominator vanishes and the term is dropped.

Why does a knot of multiplicity $p$ allow a corner? The derivative spline has degree $p-1$ and carries the same knot, now of multiplicity equal to its degree plus one: that is the row of the table where a spline may jump, here the tangent.

A parameter that is an interval

So far $t$ has been a number. In this kernel it is often an interval, the output of a search over a box. What does it mean to evaluate a spline over an interval that reaches across several knots?

The natural mistake is the one the kernel used to make. find_span picks one span. Evaluating that span’s polynomial over a $t$ that runs past the span’s end is not the spline any more, it is the polynomial continued, and a polynomial piece of a fitted curve continued across other spans goes wild. It is no enclosure of the spline there. It holds the right value only by luck.

AGENTS.md records the consequence in Count the work, not the seconds: a spline evaluated over an interval reaching past one knot span used that span’s polynomial for all of it, and came out 25 times too wide over a fitted curve, which a search then subdivided away at great cost. The comment on find_spans gives the same figure for a cubic fitted through 48 points and evaluated over its whole domain. And a test (a_wide_parameter_encloses_the_curve_tightly in core/geop-core-geometry/src/nurb_curve/evaluate.rs) records the extreme: a quarter circle of radius 2.25, interpolated through 49 points, came out 57 wide over its whole domain.

The fix is find_spans: when $t$ is an interval reaching past one span, it returns every span it reaches, with the piece of $t$ in that span; the curve’s evaluation is the union of the pieces. Union, not intersection, since these are different points of the curve and the answer must cover them all.

Notably this is not yet complete, and the source says so. Curves unite their spans (evaluate_into, cartesian_derivatives), and so do a surface’s partials, but NurbSurface::evaluate still extrapolates from find_span’s one span. The comment on find_spans records that uniting the pieces there left a boolean in chamfer_cylinder_flat_side with a face straddling a cylinder’s cap: something downstream leans on the extrapolated width, and that has not been tracked down. AGENTS.md adds that the tangent-branch search on blends pays for this. This is honest engineering: a known limit, written next to the code that has it.

Fitting a spline through points

A kernel makes splines two ways. A designer gives control points, or an algorithm gives points on the curve, typically the samples of a traced intersection (Lecture 16). Then we need the control points that pass through them.

Given data points $Q_0,\dots,Q_m$, the method (Piegl and Tiller, chapter 9) is:

1. Choose parameters $t_k$, one for each data point.
2. Choose knots so the system is well posed. The averaging technique sets $u_{j+p}=\frac1p\sum_{i=j}^{j+p-1}t_i$, which guarantees each span contains a parameter value; by the Schoenberg-Whitney condition, the matrix is then nonsingular.
3. Solve $\sum_i N_{i,p}(t_k)P_i=Q_k$ for the control points. Because of local support, this system is banded.

Step 1 looks like a free choice, and in the mathematical sense it is: any increasing sequence gives a curve through the points. But AGENTS.md asks that a free choice also be a good one, and here is a measured example. A marching trace takes even strides, then one short final leg onto the vertex where it terminates. With uniform parameters the spline must cover the long legs and the short one in equal parameter time, and overshoots. According to the comment on chord_length_params in core/geop-core-geometry/src/nurb_curve/interpolate.rs, the resulting edge left the true intersection by 7.1e-4, and the same arc traced from its two ends produced two curves 1.3e-4 apart, which no containment search can recognise as one curve. Chord-length parameters, each point getting the fraction of the polyline’s length before it, fixed this:

$$t_0=0,\qquad t_k=t_{k-1}+\frac{\lvert Q_k-Q_{k-1}\rvert}{\sum_j\lvert Q_j-Q_{j-1}\rvert}.$$

An interpolant lies between its points

A curve through samples of a true curve is not that curve. Between samples it drifts, and in a kernel where everything is an enclosure, that drift must be stated, or the kernel will conclude that a point on the true curve is not on its own interpolant.

So there are two functions. NurbCurve::interpolate passes exactly through its points. NurbCurve::interpolate_enclosing is also given true points between the samples; it fits through the samples’ sharp centres (a free choice, and solving through wide samples would amplify their width through the elimination), then measures how far each true point sticks out of the curve’s enclosure, and widens every control point by the largest such amount in each coordinate. Since the basis is non-negative and sums to one, this widens the curve by exactly that amount everywhere (widen_to_enclose, with the measured pad from enclosing_pad).

Note the honest caveat in that comment: the enclosure is exact at the checked points and “as good as that measurement in between”. This is a place the kernel measures rather than proves. Where the true points are placed is a matter of care: true_point_fractions uses one midpoint for an interior interval, where a cubic interpolant’s drift peaks, and seven for the two end intervals, where a clamped end makes the drift lopsided and several times larger.

Where does this leave us?

A B-spline does everything we asked of a piecewise curve, but it still cannot draw a circle: its pieces are polynomial, and exercise 5 of Lecture 6 still holds. And its derivative test says something about the parameter, not the shape. Lecture 8 fixes the first by adding one more coordinate.

Recap

Exercises

1. Show that a B-spline with $n+1$ control points and degree $p$ has $n-p+1$ knot spans in its domain (counting empty ones). How many free numbers per coordinate does a chain of $s$ cubic pieces joined $C^2$ have? Check that a clamped cubic B-spline with $s$ non-empty spans has that many control points.

2. Show that with $U=(0,\dots,0,1,\dots,1)$, $p+1$ of each, $N_{i,p}=B_i^p$. Start with $p=1$ and induct.

3. For the quadratic knots $U=(0,0,0,1,2,3,3,3)$ compute by hand all $N_{i,2}(1.5)$ and check that they sum to one. Which control points are active at $t=1.5$? Then use the demo with degree 2 to see the same active-point pattern.

4. Using the derivative formula, show that a clamped curve leaves $P_0$ with $C'(u_p)=\frac{p}{u_{p+1}-u_1}(P_1-P_0)$. Then explain in one sentence why this denominator is a single span’s width for a clamped curve, whatever the rest of the knots are.

5. Show by example that uniform parameters can overshoot. Take four collinear points at $x=0,1,2,2.01$ on the $x$-axis and fit a single cubic with $t_k=k/3$. How far past $x=2.01$ does the curve go? Then repeat with chord length. Why is this the same defect as the traced curve’s last leg?

6. Reading the code. Open find_spans in core/geop-core-geometry/src/spline.rs. Under what two conditions does it return a single span with $t$ unchanged? Why is the first condition, t.is_sharp(), safe even if the sharp $t$ lies exactly on a knot?

Lecture 8: NURBS and exact conics

Lecture 7 gave us piecewise polynomials, and left one gap: a polynomial curve cannot be a circle (exercise 5 of Lecture 6). This is no academic detail. A hole, a shaft, a fillet round a straight edge, a thread: machined parts are made largely of cylinders and cones, and their cross sections are circles. Approximating a circle by a cubic to $10^{-5}$ is fine for a picture, and worthless for a kernel, where two cylinders that are in truth tangent must be recognised as tangent to within a few ulps.

Older kernels answer this by storing a circle as a separate kind of thing, with its own intersection formulas against every other kind. The history in DEVELOPERS.md (“Geop 1”) records that this was tried and “didn’t scale as equations became too complex”. The other answer is to widen the representation slightly, until the circle is merely a particular set of control points. That is NURBS: non-uniform rational B-splines. “Non-uniform” we have already met, the knots may be anywhere. The new word is “rational”.

A curve seen from the origin

Take a polynomial curve in $\mathbb{R}^3$, and write its coordinates $(X(t),Y(t),W(t))$. Now look at it from the origin. Each point on it lies on a ray through the origin, and the ray crosses the plane $W=1$ in one point:

$$C(t)=\left(\frac{X(t)}{W(t)},\ \frac{Y(t)}{W(t)}\right).$$

So the shadow of a polynomial curve, projected through the origin onto a plane, is a curve whose coordinates are ratios of polynomials. That is what a rational curve is: a polynomial curve in a space of one more dimension, projected. The three-dimensional space is called homogeneous space.

Everything we know about polynomial curves applies upstairs, and only at the end do we divide. The kernel takes this literally: control points are stored upstairs. NurbCurve<4> has $(wx,wy,wz,w)$, NurbCurve<3> has $(wu,wv,w)$ (see the doc comment in core/geop-core-geometry/src/nurb_curve/mod.rs), and de_boor runs on them as a plain polynomial computation “before any division by the weight”, as the header of core/geop-core-geometry/src/spline.rs says.

From control points to weights

A control point upstairs is a point $(X_i,Y_i,W_i)$. Dividing by the last coordinate gives its shadow $P_i=(X_i/W_i,Y_i/W_i)$, and the value $w_i=W_i$ is its weight. So $(X_i,Y_i,W_i)=(w_ix_i,w_iy_i,w_i)$. Summing upstairs and dividing:

$$C(t)=\frac{\sum_i N_{i,p}(t)\,w_i\,P_i}{\sum_i N_{i,p}(t)\,w_i}.$$

Grouping the terms gives a rational basis:

$$C(t)=\sum_i R_i(t)\,P_i,\qquad R_i(t)=\frac{N_{i,p}(t)\,w_i}{\sum_j N_{j,p}(t)\,w_j}.$$

If all weights are equal, they cancel, and we have the polynomial B-spline back. If all weights are positive, the $R_i$ are non-negative and sum to one. So with positive weights all our properties survive: end-point interpolation for clamped knots, local support, and the strong convex hull: a span of the curve lies in the hull of the shadows of its active points.

Weights act like a pull. Raising $w_i$ draws the curve towards $P_i$ without changing $P_i$.

One property is even better than before. A polynomial curve is preserved by an affine map, but a rational one is preserved by any projective map: apply a matrix upstairs, and the curve follows. So a perspective view of a NURBS curve is again a NURBS curve, with no approximation.

Things to try: set all three weights to 2 and see nothing move in the upper plot, while upstairs the curve is scaled along its rays. Then try the preset “same curve” with weights $(2,\sqrt2,1)$. The curve in the plane is the same set of points as with weights $(1,1,1)$, but the red point has moved along it: the parametrization changed. Only one combination of the weights matters for the shape. We come to it now.

Rational quadratics are conics

Consider a single quadratic span with weights $w_0,w_1,w_2$:

$$C(t)=\frac{(1-t)^2w_0P_0+2t(1-t)w_1P_1+t^2w_2P_2}{(1-t)^2w_0+2t(1-t)w_1+t^2w_2}.$$

Assume $P_0,P_1,P_2$ are not collinear, and write any point $Q$ of the plane in barycentric coordinates relative to them: $Q=aP_0+bP_1+cP_2$, $a+b+c=1$. Reading the coefficients off the formula, a point on the curve has $a:b:c=(1-t)^2w_0:2t(1-t)w_1:t^2w_2$, so

$$b^2\,w_0w_2=4\,a\,c\,w_1^2.$$

This is a homogeneous quadratic equation in $(a,b,c)$, and those are affine functions of position. So every rational quadratic lies on a conic section, with $P_0P_1$ and $P_2P_1$ as tangents at its ends. Conversely, any conic arc whose end tangents meet in a finite point is one such span.

The equation has one parameter. Dividing through, only the ratio

$$k=\frac{w_1^2}{w_0w_2}$$

appears, so $k$ determines the shape alone. Changing the three weights while keeping $k$ is the reparametrization you tried in the demo, $t\mapsto\frac{\rho t}{1-t+\rho t}$. We can therefore normalize to $w_0=w_2=1$, leaving $w=w_1=\sqrt{k}$:

The parabola case is what you should expect: equal weights give a polynomial quadratic, and a polynomial quadratic is a parabola. As for the rest, look at the middle of the span. By the formula $C(\tfrac12)$ is on the segment from the chord’s midpoint $M=\tfrac12(P_0+P_2)$ to $P_1$, at the fraction $\frac{w}{1+w}$ of the way. So $w$ runs from the chord ($w=0$) towards $P_1$ ($w\to\infty$) with the parabola at the half.

The exact circle

Now the construction every circle in this kernel is built from. We want an arc of a circle of radius $r$ centred at the origin, turning through an angle $\theta<\pi$, symmetric about the $x$-axis. Take

$$P_0=r(\cos\tfrac\theta2,-\sin\tfrac\theta2),\qquad P_2=r(\cos\tfrac\theta2,\sin\tfrac\theta2).$$

A conic arc’s end tangents meet at $P_1$, so $P_1$ has to be where the circle’s tangents at $P_0$ and $P_2$ meet:

$$P_1=\left(\frac{r}{\cos\frac\theta2},\ 0\right).$$

By symmetry, the curve’s midpoint is on the $x$-axis, and it must be at $(r,0)$. The midpoint is $M+\frac{w}{1+w}(P_1-M)$, with $M=(r\cos\frac\theta2,0)$. Write $c=\cos\frac\theta2$:

$$r=rc+\frac{w}{1+w}\left(\frac rc-rc\right)\ \Longrightarrow\ 1-c=\frac{w}{1+w}\cdot\frac{1-c^2}{c}\ \Longrightarrow\ \frac{w}{1+w}=\frac{c}{1+c}\ \Longrightarrow\ w=c.$$

So the middle weight is

$$w=\cos\frac\theta2.$$

For a quarter turn, $\theta=90^\circ$, this is $\cos45^\circ=\frac{\sqrt2}{2}$. Since $w<1$ the conic is an ellipse, and we still have to check that this ellipse is a circle. It is: a conic is fixed by five conditions, we imposed two end points, two tangents, and one more point; the circle meets them all. For a direct check, take the quarter arc $P_0=(1,0)$, $P_1=(1,1)$, $P_2=(0,1)$, $w=\frac{\sqrt2}{2}$:

$$x=\frac{(1-t)^2+\sqrt2\,t(1-t)}{D},\quad y=\frac{\sqrt2\,t(1-t)+t^2}{D},\quad D=(1-t)^2+\sqrt2\,t(1-t)+t^2,$$

and a short expansion (exercise 2) gives $x^2+y^2=1$ for all $t$.

Try the circle presets in the demo: the largest deviation from the radius, over 2000 samples of every span, is of the order of $10^{-15}$. That is rounding in the double-precision sampling of $w$ and the arithmetic. It is not an approximation error, and doesn’t shrink or grow as you add spans.

A whole circle

A full circle needs more than one span, as $\theta<\pi$. It is usually four quarters. Each span is its own Bézier piece, so each interior knot has multiplicity two, and the curve passes through a control point at every joint:

$$U=(0,0,0,\ \tfrac14,\tfrac14,\ \tfrac12,\tfrac12,\ \tfrac34,\tfrac34,\ 1,1,1),$$

with nine control points (the first and last coincide) and weights $1,\frac{\sqrt2}2,1,\frac{\sqrt2}2,\dots$. This is what NurbSurface::revolve writes in $u$: $2n+1$ rows, knots at multiplicity two, in core/geop-core-geometry/src/nurb_surface/revolve.rs.

The price of the circle

Two things are lost.

The parameter does not run at constant speed. No rational parametrization of a circle has $\lvert C'\rvert$ constant. Over a quarter arc with $w=\frac{\sqrt2}2$ the speed in the middle is about 1.17 times the speed at the ends (exercise 3). So the parameter is not the angle, and finding “the point at 30 degrees” takes a little computation. The kernel’s helices run into this directly: the header of core/geop-core-geometry/src/nurb_curve/helix.rs explains that a helix cannot be a rational curve, because its height grows linearly with the angle, and “no rational parametrization of a circle turns at a constant rate”. So it keeps the circle exact, and approximates the height: quarter-turn spans stay within 0.53% of the pitch, “5.3e-3 times the pitch, measured” in the comment. The curve lies exactly on its cylinder; only its rise is approximate.

The joints are only $C^1$, even where the shape is perfect. Across a knot of multiplicity 2 in a quadratic, a rational curve is only guaranteed $C^0$, by the table in Lecture 7. With equal spans of a circle, the two speeds at a joint agree and so does the direction, so the first derivative is continuous. The second is not: the speed is rising on one side of the joint and falling on the other. The circle has constant curvature, so the shape is as smooth as can be; only its parametrization has kinks in its acceleration.

That detail has practical consequences. AGENTS.md’s section on hash maps ends with a passage on exactly this: the kernel’s arcs and helices are rational quadratics, “only C1 where their spans meet”. fit_pcurve used to fit one $C^2$ cubic through 48 samples of the whole curve, and across those joints the drift shrank only as $h^2$. AGENTS.md says a helix of pitch 1 came within 15% of the 1e-4 an entity may carry. The joints are known: they are the knots of the source curve. So it now fits each smooth piece separately and joins them. The lesson: use the structure you know before you sample blindly. We return to this in Lecture 16.

Positive weights, always

Negative or zero weights are legal in some systems. Geop refuses them. Look at the check in NurbCurve::try_new, in core/geop-core-geometry/src/nurb_curve/mod.rs:

#![allow(unused)]
fn main() {
if !p[D - 1].definitely_greater(Scalar::ZERO) { /* error */ }
}

The reason is everything built on the denominator. If all weights are positive, then $W(t)=\sum_iN_i(t)w_i$ is positive on the whole domain, because it is a convex combination of positive numbers. Then:

• The convex hull property holds, and every search’s first test is valid.
• $C(t)$ never divides by zero, and in fact the division can be avoided entirely. The containment search asks where $X_k(t)-p_kW(t)=0$ (contains::curve), a polynomial equation with the same zeros and no quotient. We will use this in Lecture 12.
• Splitting forms convex combinations of existing homogeneous points, so it preserves positivity.

Note the three-valued test: a weight must be definitely positive, not merely $>0$ at its midpoint. An interval weight of $[-10^{-20},1]$ is rejected.

What about a semicircle in one span? $\theta=\pi$ gives $w=0$, with $P_1$ at infinity. A homogeneous control point $(X,Y,0)$ is well defined upstairs; downstairs it is a direction. This is the classical trick, and it is refused here for the reasons above. Revolves and helices are therefore cut into spans of at most a quarter turn (NurbSurface::revolve, helix_rows).

One more thing worth noting: extrapolating a rational curve is no convex combination, so it can produce a weight that is not positive. NurbCurve::extended’s comment says it fails “where the continued weights stop being positive”, and since it ends in try_new, that is where such a failure would surface.

Derivatives of a rational curve

The quotient rule applies, but it pays to organize it upstairs, where everything is polynomial. Write $A(t)=W(t)\,C(t)$, with $A$ the numerator. Then by Leibniz

$$A^{(k)}=\sum_{i=0}^{k}\binom{k}{i}W^{(i)}C^{(k-i)},$$

and solving for the last term:

$$C^{(k)}=\frac{1}{W}\left(A^{(k)}-\sum_{i=1}^{k}\binom{k}{i}W^{(i)}C^{(k-i)}\right).$$

For $k=1$: $C'=(A'-W'C)/W$. This is algorithm A4.2 of Piegl and Tiller, and it is rational_derivatives in core/geop-core-geometry/src/spline.rs, line for line. $A^{(k)}$ and $W^{(k)}$ are derivatives of polynomials, computed by Lecture 7’s hodograph.

As an example, take the quarter arc at $t=0$. $A'(0)=2(w_1P_1-w_0P_0)$ with $w_0=1$, $w_1=\frac{\sqrt2}2$, $P_0=(1,0)$, $P_1=(1,1)$, so $A'(0)=2(\frac{\sqrt2}{2}-1,\frac{\sqrt2}{2})$, and $W'(0)=2(\frac{\sqrt2}2-1)$. Then $C'(0)=A'(0)-W'(0)P_0=(0,\sqrt2)$. So $\lvert C'(0)\rvert=\sqrt2$ for a unit circle. In general $C'(0)=2\frac{w_1}{w_0}(P_1-P_0)$, which is a useful fact: the end tangent points along the first leg, as in the polynomial case, scaled by the ratio of the weights.

The middle of the arc goes faster: by symmetry $W'(\tfrac12)=0$, so $\lvert C'(\tfrac12)\rvert=\lvert A'(\tfrac12)\rvert/W(\tfrac12)$, which exercise 3 works out to about 1.657, against 1.414 at the ends.

Look again at the formula for $C'$, and think about what happens for an interval $t$. $A'$, $W'$, $C$ and $W$ are all enclosures, each computed separately, and interval arithmetic cannot know that they are correlated. Far from the origin, $\lvert A\rvert$ is large, and this formula is in trouble. That is the subject of Lecture 11.

Recap

Exercises

1. Show that multiplying every weight of a NURBS curve by the same constant leaves the curve unchanged, parametrization included. What does this correspond to upstairs?

2. Check $x^2+y^2=1$ for the quarter arc given in the lecture.

3. Compute $\lvert C'(t)\rvert$ at $t=0$ and $t=\tfrac12$ for the unit quarter arc. What is their ratio? Do the same for a span of $120^\circ$.

4. Derive the formula $C(\tfrac12)=M+\frac{w}{1+w}(P_1-M)$ for a rational quadratic with $w_0=w_2=1$. What does the formula say about $w\to0$ and $w\to\infty$?

5. Write the nine homogeneous control points of a unit circle in the plane as $(wx,wy,w)$. Then verify by evaluating with the formula at $t=\tfrac18$ in the knot vector above.

6. Reading the code. NurbSurface::revolve writes a control row with homogeneous(middle, w.mul(middle_weight)). Why does a profile control point’s own weight w multiply the arc’s weight? What surface would you get by revolving the unit quarter arc, itself rational, about an axis through its centre?

Lecture 9: Knot insertion and subdivision

Every algorithm in Part III comes down to the same question: can we zoom in on a piece of a curve? We can evaluate a NURBS at a point, but that says nothing about what lies between points. What the searches need is another representation of the same curve, valid over a piece of it, whose control points give fresh, tighter hulls. That operation is knot insertion, and this lecture is about what it does when the numbers are enclosures. It is the place in this book where algebra and rounding interact most instructively.

The problem

Given a B-spline of degree $p$ with knots $U$ and control points $P_i$, we want to add a knot $\bar t$ to $U$ and find new control points $Q_i$ giving exactly the same curve. This must be possible: adding a knot only enlarges the spline space, since it allows the curve to have less smoothness at $\bar t$ but does not require it.

Since the curve must not change, we need

$$\sum_i N_{i,p}(t)P_i=\sum_i\bar N_{i,p}(t)Q_i,$$

where $\bar N$ are the basis functions on the new knots. There is one control point more than before. By local support, nothing far from $\bar t$ can change. The answer turns out to be pleasingly simple.

Boehm’s formula

Let $\bar t\in[u_k,u_{k+1})$. Then

$$Q_i=\begin{cases}P_i & i\le k-p,\\[2pt] \alpha_iP_i+(1-\alpha_i)P_{i-1} & k-p<i\le k,\\[2pt] P_{i-1} & i>k,\end{cases}\qquad \alpha_i=\frac{\bar t-u_i}{u_{i+p}-u_i}.$$

The new knot vector is $U$ with $\bar t$ inserted after $u_k$. So exactly $p$ old points are replaced, by points on the legs of the old polygon, and the rest are kept. No system of equations is solved, and because each $\alpha_i$ is in $[0,1]$, every new point is a convex combination of old ones. So the new polygon is the old polygon with its corners cut.

The derivation is the Cox-de Boor recursion run the other way: one shows that each old basis function is a combination of two new ones,

$$N_{i,p}=\frac{\bar t-\bar u_i}{\bar u_{i+p+1}-\bar u_i}\bar N_{i,p}+\frac{\bar u_{i+p+2}-\bar t}{\bar u_{i+p+2}-\bar u_{i+1}}\bar N_{i+1,p},$$

with the convention that a term whose denominator vanishes is dropped, and this by induction on $p$. Substituting into the left-hand sum and collecting the coefficient of each $\bar N_{i,p}$ gives the formula above (Piegl and Tiller, section 5.2, have it in full; Boehm’s paper of 1980 is the original). The Oslo algorithm generalizes this to several knots at once; the kernel does not use it, and inserts knots one at a time.

A worked example

Take the quadratic from Lecture 7, $U=(0,0,0,1,2,2,2)$, $P=(0,0),(1,2),(3,2),(4,0)$, and insert $\bar t=\tfrac12$. Here $k=2$ and $p=2$, so $i=1,2$ change:

$$\alpha_1=\frac{\frac12-u_1}{u_3-u_1}=\frac12,\qquad \alpha_2=\frac{\frac12-u_2}{u_4-u_2}=\frac14.$$

$$Q_1=\tfrac12P_1+\tfrac12P_0=(\tfrac12,1),\quad Q_2=\tfrac14P_2+\tfrac34P_1=(\tfrac32,2).$$

So $Q=(0,0),(\tfrac12,1),(\tfrac32,2),(3,2),(4,0)$, with knots $(0,0,0,\tfrac12,1,2,2,2)$. Those two new points are $d_1^1$ and $d_2^1$ of the de Boor triangle we computed for the same $t$. That is no accident: each knot insertion is one level of de Boor, and if we insert $\bar t$ $p$ times in all, we run through the whole triangle, and one control point is $C(\bar t)$ itself.

Splitting is inserting to full multiplicity

Insert $\bar t$ until its multiplicity reaches $p+1$. Then, by Lecture 7’s table, the curve is allowed to jump at $\bar t$, so the control points separate into two groups with nothing connecting them: the curve has been split into two curves. Each piece has a clamped knot vector.

For one Bézier span on $[a,b]$ all those $\alpha$ have the same value: every level of the triangle uses $\alpha=(\bar t-a)/(b-a)$. It is de Casteljau, and the two polygons along the sides of the triangle in Lecture 6 are the pieces.

Convergence

Why does subdividing help? Two facts.

First, with each insertion the polygon’s corners are cut, so the polygon moves toward the curve. Second, the rate is quadratic. For a curve with parameter span $h$, the distance between the curve and the chord of a piece is bounded by

$$\max_{t}\lvert C(t)-L(t)\rvert\le\frac18h^2\max\lvert C''\rvert,$$

the classical bound for linear interpolation. Applied to a piece’s control polygon, it shows the distance shrinking by about a factor of four when $h$ halves. Press “Halve every span” a few times and read the figures in the demo. With the starting polygon, we measured the farthest control point at $1.10$, $0.36$, $0.11$, $0.028$, $0.007$ over the rounds. The summed area of the boxes falls more slowly, $28.6$, $22.7$, $15.9$, $10.0$, $5.7$, tending to a halving per round: there are twice as many boxes each round, each a quarter of the area.

A hull of $p+1$ points around a nearly straight piece is therefore a thin sliver. An axis-aligned box of those points is not thin if the curve runs diagonally: it is $h\times h$. That is why a box test alone makes pairs of pieces multiply near a crossing, and why the kernel measures hulls across a direction when it can. See Lecture 13 and the branch and bound in AGENTS.md’s “Do not decide from samples”: “Measured across the direction between two pieces, the hull bound is off by the square of the pieces’ size… A box bound is off by the size itself, and with it the survivors multiply.”

Now the numbers are intervals

So far all this is classical. Here is the part specific to a kernel of enclosures. Let $t$ be an interval (because some search produced it), and let the knots $a,b$ be intervals (because they came from earlier cuts). Then

$$\alpha=\frac{t-a}{b-a}$$

is an interval, and its width has two sources:

$$\text{width}(\alpha)\approx\frac{\text{width}(t)+\text{width}(a)}{\lvert b-a\rvert}+\lvert\alpha\rvert\,\frac{\text{width}(a)+\text{width}(b)}{\lvert b-a\rvert}.$$

Every input width is divided by the span’s size. Then in a control point $(1-\alpha)P+\alpha P'$, it comes out multiplied by $\lvert P'-P\rvert$, which is the speed of the curve times the span’s size. Those cancel, and the control point widens by roughly the width of $t$ times the speed of the curve: a width in position, which is at least what you asked for by giving a wide $t$.

So one split at a wide $t$ widens control points by $\lvert C'\rvert\,\text{width}(t)$. No surprise. The trouble comes when the width returns. The next parameter is found on the widened piece, so its enclosure is at least the piece’s width divided by the speed, and then the new knot $a=t$ is wide as well and enters $\alpha$ twice. Each cut multiplies the width by a constant bigger than one, giving exponential growth in the number of cuts.

This is not a thought experiment. AGENTS.md records it under “Many booleans on one body”, in the section Width along a split edge:

An edge was cut at the whole refined enclosure of where it reaches the vertex. That kept every piece it could be cut into there in one: the parameter’s width went through Boehm’s α = (t - a) / (b - a) into the pieces’ control points and end knot, the next cut’s parameter, found on the wider piece, came out wider, and the width grew threefold per cut. A quarter arc cut forty times went from 3e-15 to 0.08; a disc’s rim cut by thirty-two slots ended 4e-10 below its cap, and gears stopped at 80 teeth, their rims drawn in twelve arcs to spread the cuts.

Before the fix, then, a width of 0.08 on an arc that started at 3e-15, with nothing wrong geometrically. Run it yourself:

With the default settings you can read the growth off the plot. In this simplified setting we measured a factor between 3 and 4 per cut, reaching 0.05 by the 25th cut, and two cuts later the division fails. The blue line stays below 1e-13 all the way to forty cuts. So the demo shows the same phenomenon as the AGENTS.md record, and it reproduces the rate, though not the same numbers. A real test in the kernel, a_rim_cut_many_times_stays_as_narrow_as_cut_once in core/geop-core-geometry/src/nurb_curve/split.rs, cuts an arc forty times and asserts that the widest control point stays under 1e-13.

Where to cut is a free choice

The fix is the rule from Lecture 5, and this is its canonical example. The vertex is an enclosure, and every parameter of the refined enclosure of where the curve reaches it ends a piece on the curve and within the vertex. Where to cut within it is therefore a free choice: any parameter there would serve. So choose one, sharply, and check that it does the job. The vertex keeps its width. The pieces carry only the rounding of the cut.

NurbCurve::split_parameter_at_point in core/geop-core-geometry/src/nurb_curve/refine.rs is that function: it cuts at the sharp middle of the refined enclosure and checks the curve reaches the point there. Compare what NurbCurve::split says about itself:

t is used exactly as given — it is not sharpened here. Whether t’s width may be discarded belongs to whoever produced it, since only they know what it means, so this function neither assumes nor imposes an answer.

The abstraction is careful. split cannot know. A subdivision search may split at a sharp midpoint. Code cutting an edge at a vertex goes through split_parameter_at_point. A caller that wanted the width to stay could give a wide $t$.

And the rule’s limit: AGENTS.md’s section “Validate the value you are about to use, not a wider one” tells the story of checking curve.evaluate(edge_t).could_be_equal(&vertex_point) on the wide parameter, then splitting at the sharpened one. The check answers “is the vertex somewhere along this stretch?”, while the split asks “cut exactly here”. Those are different questions. So when you sharpen, re-validate the sharp value: split_parameter_at_point does.

Pinning does not help

An idea you will surely have: after splitting, set the piece’s end control point to the vertex exactly. AGENTS.md records this as tried: “Pinning the split curve’s endpoint onto the vertex was tried and looked like it worked for the 3-D curve; applied to the pcurve — whose (u, v) came from a min_subdivision_size-bounded containment search — it pinned a control point to a box so wide the pcurve had no usable tangent left (“Cannot normalize zero”), and took the sweep from 0/175 failing scenes to 12/175.“ Fix the parameter, not the geometry.

Which formula for a lerp?

Look closer at the heart of it all: a point between two enclosures. In exact arithmetic these are the same:

$$a+\alpha(b-a)\qquad\text{and}\qquad(1-\alpha)a+\alpha b.$$

In interval arithmetic they are not, and neither is always better. Let us compute with $a=b=[1,1]$, as with a weight that should stay 1 through many subdivisions, and $\alpha=[0.3,0.4]$:

• $a+\alpha(b-a)$: $b-a=[0,0]$ (up to outward rounding), and so the result is $[1,1]$ (up to the same). Tight.
• $(1-\alpha)a+\alpha b$: $[0.6,0.7]+[0.3,0.4]=[0.9,1.1]$, of width $0.2$. The two alphas are treated as independent.

So for equal values the first wins clearly. Now take wide $a$ and $b$: $a=[0,1]$, $b=[0.1,1.1]$, with a sharp $\alpha=\tfrac12$:

• $a+\alpha(b-a)$: $[0,1]+\tfrac12[-0.9,1.1]=[-0.45,1.55]$, width 2. Since $a$ is mentioned twice, its width counts twice.
• $(1-\alpha)a+\alpha b$: $\tfrac12[0,1]+\tfrac12[0.1,1.1]=[0.05,1.05]$, width 1. Tight.

The first case is what a repeatedly split curve’s weights need; the second is what a patch with truly wide control points needs. Since both are valid enclosures of the same number, intersect them. This is Scalar::interpolate in core/geop-core-math/src/scalars/mod.rs:

#![allow(unused)]
fn main() {
if a.is_exactly(b) && alpha.lo >= 0.0 && alpha.hi <= 1.0 {
    return a;
}
let via_delta = a.add(alpha.mul(b.sub(a)));
let via_weights = Self::ONE.sub(alpha).mul(a).add(alpha.mul(b));
via_delta.intersect(via_weights)
}

The opening line is an extra observation: every point between two equal enclosures lies in that enclosure, whatever the two values stand for. Without it, $b-a$ rounds outward even where it is exactly zero, and $a$ plus that is $a$ widened by an ulp. The comment notes the benefit: “A polynomial patch cut in pieces keeps its weights exactly one, and stays a polynomial.” That has a practical consequence: core/geop-core-topology/src/mass.rs can then integrate with exact rules (Lecture 41).

The story of how the kernel got here is in AGENTS.md’s “Trace numerical bugs to their root cause”: a blowup in curve_curve_intersect was first “fixed” by trimming curves near a shared vertex, which only moved away from the region where the bug lived. Tracing down one level at a time found the real chain. Repeated NurbCurve::split calls compounded width in a control point’s weight (root cause one, fixed by the tighter formula above), and in knot values divided by an ever-shrinking span (root cause two, fixed by sharpening a free-choice value right after the operation that introduced the rounding noise). Use the second control of the demo to compare: with only one of the formulas, how do the two curves change?

Carrying a curve on

There is one more operation in this family, and it has a lesson. Booleans and shells often want a curve or surface longer than it is: a face that moves needs room on its neighbours’ surfaces. Mathematically the polynomial just continues. De Casteljau at $s>1$ extrapolates as exactly as it interpolates, in that the left polygon of a “split” at $s$ is the control polygon over $[0,s]$.

Note that this isn’t a convex combination any more. $\alpha>1$ amplifies, and a rational curve’s weights can fall.

The first implementation replaced the end piece by the piece over the longer span, recomputing all its control points by extrapolation. AGENTS.md’s “Carrying a patch on must leave it as it was” gives the result: “a sphere’s octant carried on evaluated its old region four orders of magnitude wider, and a shell’s sphere vertex went from 1e-15 to 7e-10 wide.”

The fix “was not to give curved surfaces less room but to keep the curve as it was, control point for control point, and join the continuation on as a piece of its own at a knot of full multiplicity. Room is then free: where the face stays, nothing changed.” NurbCurve::extended in core/geop-core-geometry/src/nurb_curve/extend.rs does that: it isolates the end piece by knot insertion, finds the continuation by continued and beyond, then appends the continuation’s points to the untouched original ones, with the shared end point the intersect of its two enclosures. It fails on a knot vector not clamped at that end, and where the continued weights stop being positive.

The general rule: any operation that only needs more of a curve should not touch the part it already had.

Recap

Exercises

1. Insert $\bar t=1$ into the worked example’s knots $U=(0,0,0,1,2,2,2)$, and write the new control points. What is the multiplicity of the knot now, and what does the curve do at $t=1$ relative to the new control polygon?

2. Prove that the number of changed control points in one insertion is $p$, and that the $\alpha_i$ lie in $[0,1]$. What happens to the formula when $\bar t$ is already a knot of multiplicity $s$? How many points really change?

3. In the width demo, take $\alpha$ with the two single formulas, and observe the two lines. Which formula is worse for the wide cuts, and which for the sharp ones? Give a reason using the two examples in the lecture.

4. Suppose a search returns $t$ of width $10^{-10}$ on a curve with $\lvert C'\rvert\approx1$. The width of the resulting split point in space is what? Compute $\text{width}(\alpha)$ for spans of $10^{-2}$ and of $10^{-8}$. Why, at a span of $10^{-8}$, does it hardly matter any more what the control points were?

5. Why does split_bezier_span make all levels use one $\alpha$, while Boehm’s insertion would compute $(t-t)/(b-t)$ at later levels? Compute $[t-\varepsilon,t+\varepsilon]-[t-\varepsilon,t+\varepsilon]$.

6. Show that the left polygon of de Casteljau at $s>1$ is the control polygon of the same polynomial over $[0,s]$. Compute the weights of that polygon for the unit quarter arc at $s=2$. For which $s$ does its middle weight reach zero? Then compute the weights of the polygon over $[1,s]$ only, which beyond in core/geop-core-geometry/src/nurb_curve/extend.rs keeps. What does that tell you about adding a piece rather than recomputing?

Lecture 10: Surfaces

A solid is bounded by surfaces. We now have everything needed to define one, with almost no new mathematics. The trick is called the tensor product: a surface is a curve whose control points move along curves.

From a moving curve to a net

Take a curve $C(v)=\sum_jN_{j,q}(v)Q_j$. Now let each control point move along its own curve, with a common parameter $u$ and a common basis:

$$Q_j(u)=\sum_iN_{i,p}(u)P_{ij}.$$

Substituting gives

$$S(u,v)=\sum_{i}\sum_{j}N_{i,p}(u)\,N_{j,q}(v)\,P_{ij}.$$

That is a tensor-product B-spline surface: a rectangular control net $P_{ij}$, degree $p$ and knots $U$ in $u$, degree $q$ and knots $V$ in $v$, over the rectangle $[u_p,u_{n+1}]\times[v_q,v_{m+1}]$. With weights it is a NURBS surface, the homogeneous points $(w_{ij}P_{ij},w_{ij})$ obeying exactly this polynomial formula upstairs:

$$S(u,v)=\frac{\sum_{i,j}N_{i,p}(u)N_{j,q}(v)w_{ij}P_{ij}}{\sum_{i,j}N_{i,p}(u)N_{j,q}(v)w_{ij}}.$$

What makes this powerful: fix one parameter, and the rest is a curve algorithm. Fix $u=u_0$, and $S(u_0,v)$ is a B-spline in $v$ with control points $Q_j(u_0)$. This is an iso-curve (or iso-parametric curve). So:

• Evaluate: de Boor along $u$ on each column of the net gives one row of points; de Boor along $v$ on that row gives the point. The order is irrelevant.
• Split: insert the same knot in every row. This is what “works on rows” meant in Lecture 9.
• Convex hull: with non-negative weights, the products $N_{i,p}N_{j,q}$ are non-negative and sum to one. So $S(u,v)$ lies in the hull of the active $(p+1)(q+1)$ net points.
• Boundary: for clamped knots, the boundary of the patch is four curves whose control points are the outer rows of the net, and the corners are the corner control points.

The demo shows the two-stage evaluation for a bicubic Bézier patch: sixteen control points, cubic in each direction.

Things to check in the demo: move $v$ only, with $u$ still: the orange curve does not change. Choose the corner $P_{00}$ and raise it: the surface’s corner follows exactly. Raise an inner point $P_{11}$: the surface follows only part of the way.

Partial derivatives

Differentiate the double sum with respect to $u$, holding $v$. By Lecture 7’s derivative formula, applied to every column, $S_u$ is again a tensor-product surface, whose net is differences along $u$:

$$S_u(u,v)=\sum_{i,j}N_{i+1,p-1}(u)\,N_{j,q}(v)\,\frac{p\,(P_{i+1,j}-P_{ij})}{u_{i+p+1}-u_{i+1}}.$$

For a Bézier patch on $[0,1]^2$, this is $p\sum B_i^{p-1}(u)B_j^q(v)(P_{i+1,j}-P_{ij})$. The same goes for $S_v$. Second derivatives come from differencing twice: $S_{uu}$ along $u$ twice, $S_{vv}$ along $v$ twice, and $S_{uv}$ once each way. A key fact is that all these nets are finite data, so their hulls bound the derivative over a whole box of $(u,v)$. Tangent directions and normals of a piece of surface can be bounded by proof, not by sampling.

For a rational surface one differentiates upstairs and applies Leibniz in two variables, as for curves. With $A=wS$:

$$S_u=\frac{A_u-w_uS}{w},\qquad S_{uv}=\frac{A_{uv}-w_{uv}S-w_uS_v-w_vS_u}{w}.$$

The normal

$S_u$ and $S_v$ span the tangent plane at a point, if they are independent. The unit normal is

$$n=\frac{S_u\times S_v}{\lvert S_u\times S_v\rvert}.$$

Two things to keep in mind.

Its direction has meaning, and the parametrization decides it. Swapping $u$ and $v$, or reversing $u$, flips the normal. In a solid, the normal must point out of the material. The kernel has no orientation flag on a face: “Orientation has no flag. It lives in the surface’s parametrization” (from the topology crate page). reverse_face mirrors the surface in $u$ and its pcurves with it.

It can fail. Where $S_u\times S_v=0$ the formula is $0/0$. It may be that the surface itself has a crease or a tip. It may equally well be that a perfectly smooth surface has a bad coordinate system there. We will now look at the latter, the most important kind.

Fundamental forms and curvature

Everything a kernel needs to know about the local shape of a surface comes from two quadratic forms. They are worth deriving rather than memorizing.

First fundamental form. Move in the domain at a velocity $(\dot u,\dot v)$. The velocity in space is $S_u\dot u+S_v\dot v$, of squared length

$$\lvert\dot S\rvert^2=E\dot u^2+2F\dot u\dot v+G\dot v^2,\quad E=S_u\cdot S_u,\ F=S_u\cdot S_v,\ G=S_v\cdot S_v.$$

So $E,F,G$ measure lengths and angles on the surface, using only the parameters. A patch of surface of domain size $\mathrm{d}u\,\mathrm{d}v$ has area $\sqrt{EG-F^2}\,\mathrm{d}u\,\mathrm{d}v=\lvert S_u\times S_v\rvert\,\mathrm{d}u\,\mathrm{d}v$.

Second fundamental form. Now take a curve on the surface, $\gamma(s)=S(u(s),v(s))$, with arc length $s$. Its acceleration $\ddot\gamma$ has a component in the tangent plane, which says how the curve turns within the surface, and one along $n$, which says how much the surface forces it to bend. Since $\dot\gamma\cdot n=0$ everywhere, differentiating gives $\ddot\gamma\cdot n=-\dot\gamma\cdot\dot n$, and expanding $\ddot\gamma$:

$$\ddot\gamma\cdot n=L\dot u^2+2M\dot u\dot v+N\dot v^2,\quad L=S_{uu}\cdot n,\ M=S_{uv}\cdot n,\ N=S_{vv}\cdot n.$$

This $\ddot\gamma\cdot n$ is the normal curvature $\kappa_n$ in the direction $(\dot u,\dot v)$. It depends only on the direction, not on the curve. For a general (non-unit) velocity:

$$\kappa_n=\frac{L\dot u^2+2M\dot u\dot v+N\dot v^2}{E\dot u^2+2F\dot u\dot v+G\dot v^2}.$$

As a direction turns, $\kappa_n$ varies between a maximum and a minimum, the principal curvatures $\kappa_1,\kappa_2$. They are the eigenvalues of $\mathrm{I}^{-1}\mathrm{II}$, where $\mathrm{I}=\begin{pmatrix}E&F\\F&G\end{pmatrix}$ and $\mathrm{II}=\begin{pmatrix}L&M\\M&N\end{pmatrix}$. Their determinant and half trace are the Gaussian and mean curvature:

$$K=\kappa_1\kappa_2=\frac{LN-M^2}{EG-F^2},\qquad H=\frac{\kappa_1+\kappa_2}{2}=\frac{LG-2MF+NE}{2(EG-F^2)}.$$

$K>0$ means the surface locally stays on one side of the tangent plane, like a dome. $K<0$ means a saddle. A surface with $K=0$ at every point, like a plane or a cylinder, is called developable: it can be unrolled to a flat sheet without stretching.

A worked example

Take the bicubic form of $z=xy$: $S(u,v)=(u,v,uv)$ on $[0,1]^2$ (with the control heights $\frac i3\frac j3$; try the “saddle” preset, which has this shape up to scale). Then

$$S_u=(1,0,v),\quad S_v=(0,1,u),\quad S_{uu}=S_{vv}=0,\quad S_{uv}=(0,0,1),$$

$$S_u\times S_v=(-v,-u,1),\quad n=\frac{(-v,-u,1)}{d},\quad d=\sqrt{1+u^2+v^2}.$$

$$E=1+v^2,\ F=uv,\ G=1+u^2;\quad EG-F^2=d^2;\quad L=N=0,\ M=\frac1d.$$

So $K=\frac{-M^2}{d^2}=-\frac{1}{d^4}$ and $H=\frac{-2MF}{2d^2}=-\frac{uv}{d^3}$. At $(0,0)$: $K=-1$, a true saddle. These two formulas are what book/scripts/demo_tests/surface-patch.test.mjs checks the demo against.

The surfaces of a CAD system, as NURBS

Everything a mechanical part is made of is exactly representable.

A plane or planar rectangle is a bilinear patch: degree (1,1), four control points, weights 1 (NurbSurface3D::bilinear). It is the simplest tensor product, and its normal is $(P_{10}-P_{00})\times(P_{01}-P_{00})$.

A ruled surface joins two curves by straight lines: $S(u,v)=(1-v)A(u)+vB(u)$. If $A$ and $B$ have the same degree and knots, the net is simply their two control polygons side by side: degree $(p,1)$. Making two arbitrary curves compatible (the same degree and knots) is by degree elevation and knot insertion, NurbCurve::compatible in nurb_curve/compatible.rs, and then NurbSurface3D::ruled. An extrusion is a ruled surface between a curve and its translate, NurbCurve::sweep in nurb_curve/sweep.rs. So a cylinder is a circle extruded: degree (2,1).

A surface of revolution turns a profile $P(v)=\sum_jN_j(v)w_jP_j$ about an axis. Each profile control point describes a circle, which Lecture 8 tells us how to write exactly:

$$\text{row } j\text{ of the net: the arc through the turn of } P_j,\text{ with weights } w_j\cdot(1,\ \cos\tfrac\theta2,\ 1,\ \cos\tfrac\theta2,\ \dots).$$

So the net is a product, and the weights multiply. NurbSurface::revolve in nurb_surface/revolve.rs does this: $u$ is the turn over $[0,1]$, with rational quadratic pieces of at most a quarter turn, and $v$ is the profile’s own parameter, so the degree is $(2,q)$. A cone revolves a straight line that meets the axis. A sphere revolves a half circle. A torus revolves a circle clear of the axis. They are not kinds in the kernel, only control points.

Poles: where the coordinates fail

Take a revolved surface whose profile touches the axis: a sphere’s pole or a cone’s tip. The profile’s control point on the axis does not move as it turns. So a whole row of the net collapses to one point. The surface at that boundary is the same point no matter what the turn $u$ is, so $S_u=0$ there, and $S_u\times S_v=0$: the formula for the normal fails.

The surface may be perfectly fine. For a sphere, the pole is a point like any other, with a good tangent plane; a different coordinate system would not notice it. This is the same situation as the Earth’s north pole in latitude and longitude: every longitude names one point, and no step in longitude moves you, but the ground is flat there as anywhere else.

AGENTS.md’s “Degenerate parametrizations: the surface is fine, the coordinates are not” gives three real bugs, and each comes from confusing the coordinates with the geometry:

• A seed on the pole never leaves it. NurbSurface::project skipped any Newton step whose Jacobian could be singular. From a seed on the pole every step is singular, so it returned the pole for every target. The fix is to step along the parameter that still moves the point, down a meridian: which meridian is a free choice, so take the one that heads most towards the target (off_pole in nurb_surface/project.rs).
• Decide tangency where the geometry is, not at moving foot points. Two quarters of a sphere meet at its pole, as one smooth surface. Comparing normals at foot points clamped to each patch’s domain made them two nearly equal conditions, and a singular system.
• Straight in $(u,v)$ is not preserved by offsetting. A sphere’s meridian cut by a plane comes out as a small circle on the inner copy.

So how does the kernel get a normal at a pole? From the net, directly. NurbSurface::normal first tries $S_u\times S_v$. If that cannot be normalized, it looks for a pole on a boundary row the point touches (rows_touching, pole, pole_normal). Leaving a pole, the derivative across the row is a positive combination of the spokes, the vectors from the pole to the next row’s control points. So the pole has a tangent plane exactly when those spokes are coplanar, and its normal is that plane’s. If the spokes are not coplanar, as for a cone’s apex, there is no single normal, and it is an error: “not an arbitrary pick among them”.

Check that against a sphere and a cone. A sphere’s next row near the pole lies in the tangent plane at the pole, because a circle’s control polygon leaves along the tangent: the spokes are coplanar. A cone’s next row is a circle further down, and the spokes form a cone. So the sphere’s pole has a normal, and the cone’s apex has none, by the same test, with no branch on the kind of surface.

What does a loop on such a face look like? Lecture 18 answers that. A collapsed row becomes a coedge that lies on a vertex, and the pcurve runs along the row.

What the hull does for surfaces

We close with the property that makes all this searchable. Each net cell’s hull encloses a piece of surface. Split in $u$ and $v$ and the net refines, and the distance from the net to the surface falls with the square of the cell’s size, as for curves. So a question like “could this surface contain this point?” comes down to a question about a finite array of numbers. That is surface_could_contain in core/geop-core-geometry/src/contains/surface.rs, with fat-line clipping in place of plain bisection, as Lecture 13 will explain. For the cheapest case there is NurbSurface::aabb.

Recap

Exercises

1. Write the net, weights and knots for a quarter of the unit cylinder of height 2, as a NURBS surface of degree (2,1).

2. For the quarter cylinder of exercise 1 at $u=0$, compute $S_u$, $S_v$, $n$, $E,F,G$ and $L,M,N$, and hence $K$ and $H$. Use $S_u(0)=2w_1(P_1-P_0)$, and find $S_{uu}$ from the circle: at parameter speed $\sigma$, a unit circle’s second derivative has normal component $\sigma^2$ towards the axis.

3. Show that for $F=0$, $\kappa_n$ along $u$ is $L/E$. Compare with curvature_radius in nurb_surface/curvature.rs, and say exactly when its result is a true principal curvature.

4. Why must a surface’s two knot vectors be separate, with a single net? A designer asks for a patch with an extra row of control points only in one part of it: what happens? (This is a real limit of tensor products. It is why trimmed faces and multiple patches exist.)

5. A revolved profile has one control point on the axis, and its other control points clear of it. Which row of the net collapses? Which partial vanishes at that row, and which does not? Compute $S_v$ at the pole for a sphere of radius 1, with $v$ the parameter along a quarter-arc profile from the pole to the equator. Does its direction depend on $u$?

6. Show that a tensor-product Bézier patch with an axis-aligned box has its normal bounded by the cross products of the difference nets. How would you bound the normal over a sub-box? Give a condition under which the patch cannot fold over onto itself.

Lecture 11: Evaluating relative to the patch

This lecture is about one line of arithmetic, and it is a model of how numerical debugging goes in this kernel. It starts with a real failure, descends to a specific division, and ends with a rule. We will compute it live.

A cylinder that overlapped itself

AGENTS.md, in “Evaluate relative to the patch”, records a failure in fixed-point arithmetic. A cylinder’s cap lying 50 along the axis came out, when evaluated over a parameter box, spanning $[38,94]$ in height. All the points of the cap are at height 50. Validation found the cylinder’s two caps overlapping. A quarter arc at $(1000,-500)$ gave a tangent 30 wide.

Today’s Scalar is one concrete f64 interval (core/geop-core-math/src/scalars/mod.rs), and the f64 format hides the rounding part of this. It does not hide the rest, as we are about to see: what went wrong is the formula, and it can be reproduced in doubles.

The formula

Recall from Lecture 8 that a point of a NURBS curve is $C(t)=A(t)/W(t)$, with $A$ and $W$ polynomial. Take the height $z$ of a curve lying in the plane $z=z_0$. Upstairs $A_z(t)=z_0W(t)$, so $C_z(t)=z_0$, exactly.

Now let $t$ be an interval, as it often is. Interval de Casteljau gives an enclosure of $A_z$ and of $W$, computed separately. Let $W(t)$ range over $[W_{\min},W_{\max}]$ on this box. Then $A_z$ encloses $[z_0W_{\min},z_0W_{\max}]$, and interval division gives

$$\frac{[z_0W_{\min},z_0W_{\max}]}{[W_{\min},W_{\max}]}=z_0\left[\frac{W_{\min}}{W_{\max}},\frac{W_{\max}}{W_{\min}}\right].$$

The numerator and denominator move together, exactly in step, but interval arithmetic has no way to know that. So the quotient has width about

$$\text{width}\approx 2\,\lvert z_0\rvert\,\frac{W_{\max}-W_{\min}}{W}.$$

It is proportional to the distance from the origin, and to how much the weight varies over the box (in practice more than the true $W_{\max}-W_{\min}$, since interval de Casteljau overestimates that range too). It has nothing to do with the geometry: a plane at height $z_0$ has no extent in $z$. Move the part away from the origin and the enclosure grows, though the part has not changed.

The same formula shows why this kernel’s arcs are the usual victims. Polynomial curves have $W=1$, and this effect is absent. On a quarter arc, $W$ runs from 1 at the span’s ends to about 0.85 in its middle.

Let’s put numbers to it. We computed an arc of radius 3 at $(1000,-500)$ at height 50, the configuration of the kernel’s test far_arc in core/geop-core-geometry/src/nurb_curve/evaluate.rs, with $t\in[0.33,0.34]$, by the demo below. The height’s enclosure comes out $[49.66,50.34]$, width $0.67$. A wide $t\in[0.1,0.9]$ gives $[22.2,112.4]$. In doubles! So that fixed-point number $[38,94]$ was not primarily about fixed point. The test a_cap_far_from_the_origin_stays_flat in nurb_surface/evaluate.rs records the original case.

The cure: a local origin

Why should the distance to the origin matter at all? There is nothing special about the origin. The patch has its own natural place: any control point of it.

So translate the patch to bring one of its control points to the origin, evaluate, and translate back. Pick $O=P_k$, as a sharp point. Upstairs, translation by $-O$ keeps the weights and changes the numerators: $A_i\mapsto A_i-w_iO$. Then

$$C(t)=O+\frac{\tilde A(t)}{W(t)},\qquad\tilde A=\sum_iN_i(t)\,w_i\,(P_i-O).$$

Now $\lvert\tilde A\rvert$ is at most about the patch’s size, however far the patch is from the origin. The same $A/W$ width formula now gives

$$\text{width}\approx 2\,\lvert\text{extent}\rvert\,\frac{W_{\max}-W_{\min}}{W},$$

proportional to the patch’s size, and zero along any direction the patch is flat in. Adding $O$ back costs a single rounding.

Which point is a free choice. And in the spirit of “a free choice still has to be a good one”, the kernel takes a control point at the nearer end of the span: nearer_end in core/geop-core-geometry/src/spline.rs. The doc comment on centered adds why that end is the best: at a pole, the homogeneous rows differ by their weights times the pole’s distance from the origin, and so do $A_u$ and $W_uS$, whose difference is the derivative along the row. “Relative to the pole, both are of the order of that distance, and so is their difference.”

One more point of detail, typical of this kernel. The origin is sharpened: at.sharpen(), since which origin to use is a free choice and a sharp one costs nothing. It is one division, “not a reciprocal and a product: a quotient the format holds, 2 / 1, is then the middle of its enclosure, and the point sits exactly on the origin, with nothing of its own left over for the derivatives to see.”

Now see it for yourself.

Things to notice:

1. With a box of width $0.01$, red rises in proportion to the distance; blue stays flat at about $0.065$, for distances up to $10^6$, while the true stretch is about $0.043$.
2. In the readout, the height of the relative form is narrower than $10^{-10}$ at any distance on the slider, where the direct form’s reaches $6.7\times10^2$ at $10^6$.
3. Change the box to “sharp”: red and blue come within a small factor of each other. Both are rounding only.

That last point is exactly the AGENTS.md remark: “in f64 the effect hides behind the format’s relative precision until the width is read”. A test that only checks that the enclosure contains the right value keeps passing, because a wide enclosure contains everything. That is why a test of geometry far from the origin must assert the width of what it gets:

#![allow(unused)]
fn main() {
assert!(p[2].width().to_f64() < 1e-6, "{p:?}");
}

is from an_arc_far_from_the_origin_is_as_wide_as_its_stretch in nurb_curve/evaluate.rs. And vertex.point[k].width().to_f64() < 1e-9 is in ops/geop-ops-surface/src/knit/tests.rs, among others.

That one rounding is not free

Should we then always evaluate relative to the patch? AGENTS.md’s next paragraph says no, with a beautifully tiny example:

A STEP edge lay 2e-31 above its plane face. Evaluated relative to a point 1.8e-15 up, its enclosure took in the plane, the importer measured no gap and did not widen the edge, and validation’s clipping, which reads the control points, found the face’s points off it.

Follow what happens. A number near $10^{-31}$ is held with a relative precision of $10^{-16}$, so to about $10^{-47}$. It is added to an origin $O=1.8\times10^{-15}$, whose spacing between doubles is about $4\times10^{-31}$. That single rounding is bigger than the value. The direct form never adds that origin back, so it has no such rounding.

So neither form dominates. Direct evaluation is tighter where the numbers are small to begin with; relative evaluation is tighter where the patch is far from the origin. Both are valid enclosures of one value. By the rule of Lecture 5: intersect them. In evaluate_into:

#![allow(unused)]
fn main() {
let absolute = de_boor(p, knots, local, piece, p);
let (moved, origin) = centered(local, near);
let relative = de_boor(p, knots, &moved, piece, p);
// The weights are not moved: one `W` for both.
let inv_w = inverse(absolute[D - 1], span)?;
for (c, o) in out.iter_mut().enumerate() {
    let near = origin[c].add(relative[c].mul(inv_w));
    let value = absolute[c].mul(inv_w).intersect(near);
    *o = if reached { o.union(value) } else { value };
}
}

Read that code as a summary of three lectures. The intersect combines two enclosures of the same point. The union covers the spans a wide parameter reaches, which are different points (Lecture 7). And the division by $W$ has one numerator for each form, with the same $W$, as translating leaves it unchanged.

NurbCurve::cartesian_derivatives in nurb_curve/tangent.rs does the same for derivatives: at(local), and at(&centered(local, near).0), intersected. Derivatives do not depend on where the curve lies, but their computation does.

Surfaces follow the pattern differently. NurbSurface::evaluate_at evaluates only relative to the corner nearest the point (corner_nearest), and does not form the absolute value at all.

Beyond the local origin: the mean value form

Moving the origin takes away the part of the width caused by distance. There is a second effect it does not fix. Look at the formula again: even relative to a control point, the width is $2\lvert\text{extent}\rvert\Delta W/W$, so the quotient widens with how far the point lies from the span’s control point times how much the weight changes. On a span of a rational surface, with a foot-point box of $(u,v)$, that can still be much too wide.

AGENTS.md’s “A bound is as local as what it bounds” gives the case. A shell’s copy of a traced edge, fitted through foot points evaluated over the boxes they project to, came out eight times as wide as the edge, and two offsets in a row went past the 1e-4 an entity may carry. The doc comment on NurbSurface::evaluate adds: “a cylinder evaluated at a foot point a micron wide came out ten microns wide along its axis”.

Enter a classic tool of interval analysis. By the mean value theorem, for $t$ in a box with middle $m$,

$$S(t)\in S(m)+S'(\text{box})\,(\text{box}-m).$$

This is the mean value form. Its width is $\lvert S'(\text{box})\rvert$ times the box’s width: the actual stretch of geometry, as long as the derivative’s enclosure is decent. It adds no factor of $\Delta W/W$. The price is computing the derivative over the box. A classical result of interval analysis is that the mean value form’s excess width shrinks as the square of the box’s width, whereas the direct evaluation’s shrinks only with its width. So for narrow boxes it wins, and for wide ones it can be worse, in a given coordinate: when $S'(\text{box})$ is wide, the product is large. In the demo, at a box of $0.0001$, green lies on the dashed line, at $0.01$ just above it, and at $0.1$ further off.

So NurbSurface::evaluate again intersects:

#![allow(unused)]
fn main() {
let direct = self.evaluate_at(u, v)?;
if u.is_sharp() && v.is_sharp() {
    return Ok(direct);
}
// Only ever tighter: where the partials over the box cannot be had,
// as at a pole or past the domain, the direct enclosure stands.
let (um, vm) = (u.midpoint(), v.midpoint());
let (Ok(middle), Ok((su, sv))) = (self.evaluate_at(um, vm), self.derivatives(u, v)) else {
    return Ok(direct);
};
let (du, dv) = (u.sub(um), v.sub(vm));
Ok(Vector3::from_array(std::array::from_fn(|c| {
    let mean = middle[c].add(su[c].mul(du)).add(sv[c].mul(dv));
    direct[c].intersect(mean)
})))
}

Note the structure: the improved form may only tighten, and any failure leaves the direct enclosure standing. That matches the refinement rule of Lecture 14: make it infallible, or callers start needing fallbacks of their own. Curves do not currently do this: evaluate_into has no mean value step.

Ratios of Bernstein coefficients

The same section of AGENTS.md has one more related trick, worth noting here. A ratio of two polynomials bounded by $\max(\text{numerator})/\min(\text{denominator})$ over a box suffers the same $\Delta W$ factor. For a polynomial $\sum B_ia_i$ over a positive polynomial $\sum B_id_i$, there is a sharper bound:

$$\frac{\sum_iB_ia_i}{\sum_iB_id_i}\le\max_i\frac{a_i}{d_i},$$

as the left side is a convex combination of the right-side ratios. This is free, and it took a vertex’s width from 6e-11 to 1.4e-12, according to AGENTS.md (ratio_bound, in the fitted offset). It is the same phenomenon: when correlated quantities are divided, take them together, not one at a time.

The general rule

We can now state it, in AGENTS.md’s words: “Wherever an interval formula subtracts or divides two large correlated quantities, move it to a local origin first.”

You met the phenomenon as $x-x$ in Lecture 4. Here it appears in its mature form. Ask of any interval formula: which quantities in it are in truth related, by how much, and is there a form of the computation where the relation is exploited?

• $A/W$ with $A\propto W$: subtract the local origin first, so that $A$ is small.
• $(1-\alpha)a+\alpha b$ with $a\approx b$: take $b-a$ first (Lecture 9).
• A sum of basis-function derivatives with $t$ in every term: difference the control points first (Lecture 7).
• A rational surface over a box: take its value at the middle, and add the derivative over the box times the half width (the mean value form).

None of those are approximations, and none use a tolerance. Each is a different exact expression for the same number, chosen so that outward rounding overestimates it less.

Why this matters so much

Every search in the following lectures is an interval computation over a box. If the enclosure of a surface over a box is $1000$ times wider than the piece of surface, then “could these overlap?” is true far too often, the search subdivides for nothing, and at the bottom it reports solutions where none are. A part’s distance from the world origin is not supposed to change anything, and now it doesn’t.

An exercise in humility, also from AGENTS.md: a plate placed at $(10000,-5000,2000)$ once took its extent as the greatest coordinate less the least over all three axes: 15080, for a plate 80 long. It is the same mistake in a different place: a question about the patch answered in terms of the world. Lecture 19 tells the rest.

Recap

Exercises

1. Take $P_0=(100,0)$, $P_1=(100,1)$, $P_2=(99,1)$ with weights $1,\frac{\sqrt2}2,1$, a unit quarter arc around $(99,0)$. For $t\in[0.4,0.6]$, compute by hand the enclosure of $W$, and the width of $x$ as $A_x/W$, ignoring rounding. Now do it relative to $P_0$. Compare with the actual $x$-extent of the arc over that box, using the formula $x-99=(a+b)/(a+b+e)$.

2. Write $C(t)=O+\tilde A/W$ for the point $O$ a sharp point near the patch but not a control point. Is the formula still valid? Which properties of the choice did the kernel’s nearer_end optimize, and what is lost with another $O$?

3. The demo at sharp $t$ shows both forms growing with distance. Explain that growth: where does it come from, and why is it unavoidable? At what distance is the absolute precision of an f64 coordinate $10^{-9}$?

4. Prove the mean value form: for $f$ continuously differentiable and $t$ in a box $B$ with middle $m$, $f(t)\in f(m)+f'(B)(B-m)$. What does the form give in the degenerate case $f(t)=z_0$ (constant height), and how does this compare with $A/W$?

5. Given the Bernstein coefficients of numerator $a=(2,1,3)$ and denominator $d=(1,2,1)$, bound $a(t)/d(t)$ on $[0,1]$ by (a) $\max a/\min d$, (b) $\max a_i/d_i$, and (c) the true maximum. Which rule does ratio_bound use, and why is it free?

6. Test design. Write the assertions a test should make for a cylinder of radius 2 placed with its axis through $(5000,-2000,800)$, after one boolean. Name at least three: what would you measure, and what do you compare it to?

Lecture 12: Subdivide to isolate

Part II gave us a representation of shape: a curve is a list of control points, a surface a grid of them. In Part III we ask that representation questions. Does this edge cross that face? Where? Does this curve pass through this point? Almost everything a boolean or a fillet does comes down to these.

They are hard because a curve has infinitely many points. To say “these two curves do not meet” we have to make a statement about all of them, and a program that runs for a finite time has to make it from finitely many numbers. This lecture shows how, using nothing beyond the convex hull property of Lecture 6 and the splitting of Lecture 9.

Why not just sample?

The obvious algorithm takes 100 points along each curve and looks for a pair that is close. It has two faults, and each one is fatal.

It cannot prove absence. Suppose the closest pair of samples is 0.01 apart. That says nothing about what lies between the samples. Two curves can cross and leave again between two neighbouring samples, and no number of samples rules that out.

It cannot say what “close” means. To decide that a pair is a crossing, it needs a threshold. Too small a threshold misses crossings that the sample spacing skipped. Too large a threshold reports curves that pass near each other. The threshold is an epsilon by another name, and Lecture 5 explained why those are not allowed.

The kernel’s rule is blunt: a fixed number of samples along a curve or over a surface proves nothing about what lies between them (AGENTS.md, “Do not decide from samples”). The rule is there because the code broke it. As AGENTS.md records, which way a loop runs was once decided from a polygon of eight points per pcurve, in five places with five sample counts, and a nearest point was found as the best of a few samples per span refined by Newton, so that an interference check could miss two faces touching. Both now come from the representation itself: the first from an integral (nurb_curve::loop_signed_area), the second from the search we study at the end of this lecture.

What we need is a test that talks about a whole piece at once.

The test: a box around a piece

Recall that a Bézier curve lies in the convex hull of its control points, and that the hull of a piece of the curve is found by splitting (the same holds for a NURBS curve with positive weights, Lecture 8). A convex hull is awkward to compute with, so take the coarser and cheaper thing, the axis-aligned box of the control points. Then:

If the boxes of two pieces are disjoint, the two pieces do not meet.

Read how much this says. It holds for all the infinitely many points of both pieces. It uses $O(n)$ numbers. And with interval coordinates it remains a proof: two boxes are “definitely disjoint” only if the comparison says so, and rounding is already enclosed.

The converse is false. Overlapping boxes do not mean the pieces meet, only that we have not yet ruled it out. So the test has two possible readings: no, or maybe. There is no yes, which is why the kernel’s containment functions are named curve_could_contain and surface_could_contain. The word “could” is part of the contract.

This is what splitting is for. A “maybe” is resolved by cutting the piece in two and asking each half again. The hull of a piece hugs the piece more closely the shorter it gets: for a piece of parameter width $h$, the distance from the curve to the chord between its ends is $O(h^2)$. So every “maybe” that is really a “no” turns into a “no” after a finite number of cuts, provided the two curves are a positive distance apart there.

Branch and bound

Put together, this is one of the oldest ideas in numerical search.

pending = [ the whole problem ]
while pending is not empty:
    take a piece from pending
    if the piece provably holds no solution:   drop it
    elif the piece is small enough to stop:    report it
    else:                                      cut it, add both halves to pending

The search’s central property is easy to prove: no solution is ever lost. A solution lies in the root piece. If it lies in a piece, that piece is not dropped (since dropping requires a proof that no solution lies there), and on cutting, the solution lies in one of the halves. So at the end every solution lies in a reported piece.

Note what is not promised. A reported piece need not hold a solution. It is a candidate. The kernel says so in as many words: in curve_curve_crossings a pair that converges is “a candidate, not an existence proof”.

For a pair of curves the “piece” is a pair of pieces, and cutting means halving one of the two. Which one? A free choice, and as AGENTS.md warns, a free choice still has to be a good one: cut the one whose box is larger, since parameter widths from differently scaled curves are not comparable.

Run the “three crossings” scene level by level. For the first few levels, nothing is rejected, because the boxes are large and overlap anyway. Then the pairs near no crossing start to fall away, and the alive set collapses onto three spots.

Counting the work

How many nodes does this cost? Three cases are worth doing carefully, since they explain much of what follows in Part III. All numbers in this section come from the demo’s algorithm, run at decreasing handoff sizes $h$ (the test book/scripts/demo_tests/subdivide-isolate.test.mjs checks the trends).

A transversal crossing. Near a crossing at a good angle, two pieces of size $h$ have overlapping boxes only if they lie within about $h$ of the crossing. So at every level a bounded number of pairs survive, regardless of how deep we are. Cost: a constant number of nodes per level, or $O(\log(1/h))$ overall.

Going 64 times finer took two times the nodes, and the leaves stayed constant. The width of what is reported fell in proportion to $h$.

A tangency. Now let the curves touch. Near the contact they separate as $\tfrac12\kappa x^2$, where $x$ is the distance along the common tangent and $\kappa$ the difference of the two curvatures. A box of a piece of size $h$ running obliquely across the axes has a thickness of about $h$, so two such boxes overlap until $\kappa x^2\approx h$, which is at $x\approx\sqrt{h/\kappa}$. Surviving leaves therefore cover a stretch of length $\sqrt{h}$, using pieces of length $h$, so there are $\sqrt{h}/h=h^{-1/2}$ of them.

Each quadrupling of the resolution doubles the leaves and only halves the width: the signature of $\sqrt{h}$. A tangency is hard not because anything breaks, but because the data only determines the contact to the square root of its precision. Lecture 15 takes this up. (The preset is tilted on purpose. If the common tangent runs along a coordinate axis, a box’s thickness across that axis is small near the contact, and leaves stay bounded. A search’s cost depends on how the problem sits in its coordinates. The kernel’s clip, too, measures across the other curve’s chord rather than along the axes, as Lecture 13 explains.)

An overlap. If the two curves share an arc, every pair of pieces along the arc that are near each other survives, at all depths. Then $O(1/h)$ leaves remain:

Leaves double per halving. Follow this to a kernel-sized handoff of $10^{-7}$ and the budget is gone. An overlap cannot be solved by search alone, and it cannot be recognised by counting either: read on.

Stopping: the handoff

When should a piece count as “small enough”? The naive answer is “when it is as narrow as I want the answer to be”. Subdivision gains one bit per cut, so a parameter accurate to machine precision would take about fifty levels. AGENTS.md puts the comparison plainly: about 50 levels against the roughly 7 that isolate a solution.

The better answer is that subdivision and accuracy are separate concerns.

• Subdivision is the global method. It finds every solution, separates them, and its proof of absence is the thing nothing else provides. Its convergence is linear and slow.
• Newton is the local method. It cannot find anything, but given an isolated solution it doubles the number of correct digits at every step (Lecture 14).

So min_subdivision_size is a handoff threshold, not an accuracy target: subdivide until the solutions are separated, then refine each one. In AGENTS.md’s words, the tolerance then stops influencing the final accuracy at all, which is a stronger form of the epsilon rule than making the search itself tolerance-free.

There is a second, subtler reason not to run subdivision to the bottom. At a very fine handoff, a transversal crossing surrounded by numerical width also leaves many leaves, so any signal one might read from the number of leaves disappears.

Width: where a search must stop anyway

A coordinate of the kernel is an interval. Its width is a statement about how well that coordinate is known. Cutting a piece halves its extent, the part of its box that comes from the curve actually going somewhere, but it does nothing about the width the control points carry. Once a piece’s extent is below its own width, further cuts multiply survivors without telling us anything new.

The kernel’s searches take this into account. In core/geop-core-geometry/src/fat_line.rs, stalled declares a box converged once every spatial extent is no longer definitely_greater than the larger of min_subdivision_size and carried_width, the width the data carries. Its documentation says why, and it is the right way to think about it: “resolving one object finer than another’s width answers nothing new — every piece within that width stays a candidate — and only multiplies them”.

Budgets, and why running out is an error

Every search has a second tunable: max_nodes, the number of pieces it may look at. What should a search return when it runs out?

The tempting answers are all wrong. “No solution” is wrong, because a solution may be in the pieces not yet visited. “These solutions” is wrong for the same reason: there may be more. An incomplete search has established neither, and a downstream decision built on either would be wrong in a way nothing can detect.

So the kernel gives the search a way to say so. core/geop-core-math/src/solvers/subdivide.rs is the whole of the driver, under a hundred lines:

pub enum Split<T> {
    Done,           // nothing more to search in it
    Restrict(T),    // the piece narrowed to where its solutions can be
    Halves(T, T),   // the piece cut in two, both to be searched
}

#[must_use]
pub enum Searched {
    Complete,
    Exhausted { pending: usize },
}

pub fn subdivide<T, E>(
    roots: impl IntoIterator<Item = T>,
    order: Order,
    max_nodes: usize,
    visit: impl FnMut(T, usize) -> Result<Split<T>, E>,
) -> Result<Searched, E>

This is branch and bound stripped to its skeleton. Which pieces hold nothing, when one is small enough and how to cut it are the caller’s visit; the driver promises only completeness, or says otherwise. Searched is #[must_use], and every caller with a budget it can run out of turns Exhausted into an error. In curve_curve_crossings, for example:

curve_curve_crossings: exhausted max_nodes=... with ... pairs pending; the result would be incomplete

This is the rule from AGENTS.md’s “Never use epsilons”: max_nodes and min_subdivision_size “are the only tunables allowed to affect how hard the search tries”. Change either, and the running time changes, and how wide an answer is; what a converged answer means stays the same. Their defaults live in core/geop-core-geometry/src/search.rs: MAX_NODES is 20 000, MIN_SUBDIVISION_SIZE is $10^{-7}$. The boolean’s RemeshParams (lib/geop-lib-booleans/src/remesh/mod.rs) gives its pairwise searches a looser curve_curve_min_subdivision_size of $10^{-4}$, as it documents, because a pairwise search empirically needs that to converge within its budget at all.

The usual order is breadth first (Order::BreadthFirst), level by level. That spends the budget evenly over the region, so one difficult corner does not eat the budget before the rest has been looked at.

The question of containment

The simplest search in the kernel has one curve and a point: could this curve pass through this point, and at which parameters? It is curve_could_contain in core/geop-core-geometry/src/contains/curve.rs, and surface_could_contain in contains/surface.rs is the same for a surface. Its signature is the whole contract:

pub fn curve_could_contain<const D: usize, const C: usize>(
    curve: &NurbCurve<D>,
    point: &Vector<C>,
    max_nodes: usize,
    min_subdivision_size: Scalar,
) -> GeopResult<Option<Scalar>>

• Ok(None): a proof that the curve misses the point.
• Ok(Some(t)): an enclosure of every parameter at which it could pass through. If the curve runs through the point twice, $t$ is the union of both places and is wide. It is a true statement, and the caller who needs them separately asks differently.
• Err(..): the search ran out. No statement about the point.

Its derivation is written up next to the code, in core/geop-core-geometry/src/contains/curve.md and contains/surface.md, and it is worth reading after the next lecture. The visit function does what the skeleton above says, with one addition: a piece is not only rejected or cut, it is also clipped, narrowed directly to where its solutions can lie (Split::Restrict). This is the subject of Lecture 13.

Here is the first test curve_could_contain applies to each piece: “Cheap prefilter against the cached bounding box.”

if !aabb_could_contain(&seg.aabb, point) {
    return Ok(Split::Done);
}

seg.aabb is the box of the control points the curve already carries (compute_aabb, in core/geop-core-geometry/src/aabb.rs). This is the box test of this lecture, in its plainest form, with a could in its name.

A worked example

Take $P_0=(0,0)$, $P_1=(1,2)$, $P_2=(3,2)$, $P_3=(4,0)$, the cubic of Lecture 6, and ask whether it could pass through $q=(2,1.7)$.

The whole curve’s box is $[0,4]\times[0,2]$, which contains $q$: maybe. Split at $\tfrac12$. The pieces, as computed in Lecture 6, have control points $(0,0),(\tfrac12,1),(\tfrac54,\tfrac32),(2,\tfrac32)$ and $(2,\tfrac32),(\tfrac{11}{4},\tfrac32),(\tfrac72,1),(4,0)$. Their boxes are $[0,2]\times[0,\tfrac32]$ and $[2,4]\times[0,\tfrac32]$. Since $1.7>1.5$, neither contains $q$. We have proved, with two splits’ worth of arithmetic, that the curve misses $q$. That is a statement about all the curve’s points. Notice that the original control polygon went up to $y=2$, so that one level of subdivision was what brought the bound below $q$.

The nearest point: branch and bound with a bound on the answer

The same skeleton answers a different question, and seeing it in that form is instructive. What is the least distance between two sets? It is what closest_point and the interference check need, and it was the second of the two AGENTS.md stories above. The answer now comes from least_distance in core/geop-core-geometry/src/nearest.rs. (A question that needs no least distance takes less: could_come_within in that file drops a pair definitely beyond reach before anything is asked of it, and the first points attained within it end the search. could_touch in core/geop-core-topology/src/distance.rs is that for a reach of zero. What the search cannot tell apart, at the handoff or at max_nodes, could be within it.)

Branch and bound over pairs of pieces (a span of a curve, a patch of a surface, a point) keeps two bounds.

• A lower bound for each pair, from the hulls: the two pieces are at least this far apart.
• An upper bound on the answer overall: a point of each piece that is actually attained gives a real distance, and the best so far is what the answer is at most.

A pair whose lower bound is definitely above the best attained distance holds nothing nearer, and it is dropped; any other pair is cut, its larger piece in two. Pairs are taken least lower bound first, from a priority queue, which is why this search does not use subdivide: the order is part of the method, and when the pair at the head of the queue is beyond the best distance, so is everything waiting. The result is itself an enclosure of the distance: from the least lower bound still standing to the attained distance.

Two details are worth keeping.

Which lower bound matters. Axis-aligned boxes give a lower bound that is off by the size of the piece, since an oblique piece has a fat box. Measured instead along the direction between the two pieces, the hulls give a bound that is off by the square of the piece’s size. “Near a minimum that keeps a handful of pairs alive at each level, rather than a number growing with every level”, as the module doc says. This is the same contrast we just saw in the tangent preset: what a search survives depends on how tightly its bound sits around the geometry.

The handoff again. The distance is flat at its least, so the search pins the point attaining it only to the square root of the distance’s precision. So here too subdivision isolates and Newton refines: NurbSurface::closest_point calls least_distance with NEAREST_HANDOFF ($10^{-3}$ of the size of what it searches), then Newton from each surviving pair. The search’s job is to prove that no basin was missed.

Recap

Exercises

1. In the demo, choose “three crossings”, set the handoff to 0.04, and step to the last level. Then drag one control point until two crossings merge and vanish. How do the leaves and clusters behave just before the crossings merge?

2. Prove termination. Assume two polynomial Bézier curves are at a distance $\delta>0$ from each other, and that the search halves the larger piece. Show that the search stops, and give an order of magnitude for its depth in terms of $\delta$ and the curves’ sizes.

3. A search reports a leaf whenever both pieces are below $h$, but it could equally test whether their boxes overlap by a little. Suppose a colleague proposes to drop a pair when the boxes overlap by less than $10^{-9}$, “since that can only be rounding”. Give a concrete way in which this loses a real crossing.

4. Write out the search for $p(t)=t^3-t+\tfrac14$ on $[0,1]$ by hand, using the Bernstein coefficients as the range bound, cutting at the midpoint, down to width $\tfrac18$. How many roots does it find? Which intervals survive? (On $[0,1]$ the coefficients are $\tfrac14,-\tfrac1{12},-\tfrac5{12},\tfrac14$.)

5. Why does least_distance not call subdivide? What would go wrong if its queue were first in, first out?

6. The “overlapping” preset is two pieces of one curve. Why can’t a search in principle distinguish it from two curves that are $10^{-30}$ apart? What can a kernel do instead?

Lecture 13: Fat-line clipping

The search of Lecture 12 has one weakness: it asks yes-or-no questions. A piece either survives or it does not, and a piece that survives gets cut exactly at its middle, wherever the solution is. Each cut gains one bit, no matter how much the control points know about where the solution lies.

They know a lot. Think of a curve crossing a nearly straight line. The control points say almost exactly where along the curve the crossing is. Throwing that away and halving is a waste. In 1990, Sederberg and Nishita showed how to read it off: the Bézier clipping method (Appendix A). This lecture derives it, and then the form Geop uses, which is the same idea in clothes better suited to intervals.

The fat line

Let $Q(s)$ be a Bézier curve with control points $Q_0,\dots,Q_m$. If the curve were straight, we would know its line. In general it is not, but the line $L$ through its end points $Q_0$ and $Q_m$ is a fair stand-in. Write $L$ in normal form, with unit normal $n$:

$$\mathrm{dist}(x) = n\cdot x + c,$$

the signed distance of a point $x$ from $L$. The control points are at distances $\mathrm{dist}(Q_j)$, and $\mathrm{dist}$ is affine, so the curve’s distances are a convex combination of the control points’ distances. Therefore the whole curve lies in the band

$$d_{\min}\le \mathrm{dist}(Q(s))\le d_{\max},\qquad d_{\min}=\min_j \mathrm{dist}(Q_j),\ d_{\max}=\max_j \mathrm{dist}(Q_j).$$

This band, a line with a thickness, is the fat line of $Q$. It is a better bound than a box. A box of an oblique curve is thick in both directions, whereas a fat line is aligned with the curve, and its thickness is of the order of the curve’s sag $\tfrac18 \kappa h^2$. It is second order in the piece’s length $h$, where a box’s thickness is first order. All of this lecture follows from that.

The distance function is a Bézier function

Now let $P(t)$ be the other curve, with control points $P_0,\dots,P_n$. Measure its distance from $L$:

$$d(t) = n\cdot P(t) + c = \sum_{i=0}^n B_i^n(t)\,d_i,\qquad d_i = n\cdot P_i + c.$$

We only used that the Bernstein polynomials sum to one. The distance along $P$ is a scalar polynomial in Bernstein form, whose coefficients are the control points’ distances. That is, computing it costs one dot product per control point.

We want the $t$ where $P(t)$ can be in the fat line of $Q$, that is where $d(t)\in[d_{\min},d_{\max}]$. Now comes the step that gives the method its power. Since $\sum_i \tfrac in B_i^n(t)=t$ (Bernstein polynomials have linear precision), the graph of $d$ is itself a Bézier curve in the $(t,d)$ plane:

$$\bigl(t,\ d(t)\bigr)=\sum_{i=0}^n B_i^n(t)\,\Bigl(\tfrac in,\ d_i\Bigr).$$

Its control points are $(\tfrac in, d_i)$. The convex hull property now applies to the graph. So: the graph lies in the convex hull of the points $(\tfrac in,d_i)$, and any $t$ where the graph is inside the band lies where the hull is inside the band. That is the clip: intersect the hull with the horizontal strip $d_{\min}\le d\le d_{\max}$, and read off the range of $t$. All $t$ outside that range are proved to hold no crossing.

A worked example

Suppose a cubic’s distances from a fat line are $d=(-1,\,-0.4,\,0.5,\,1)$, and the band is $[-0.1,0.1]$. The graph’s control points are $(0,-1)$, $(\tfrac13,-0.4)$, $(\tfrac23,0.5)$, $(1,1)$.

The upper hull runs from $(0,-1)$ to $(\tfrac23,0.5)$ with slope $2.25$, then on to $(1,1)$; the point $(\tfrac13,-0.4)$ lies below that first edge, whose height at $\tfrac13$ is $-0.25$. The lower hull runs from $(0,-1)$ to $(\tfrac13,-0.4)$, then from there to $(1,1)$ with slope $2.1$.

The leftmost point of the hull inside the band is where the upper hull reaches $d=-0.1$: $-1+2.25\,t=-0.1$, or $t=0.4$. The rightmost is where the lower hull reaches $d=0.1$: $-0.4+2.1(t-\tfrac13)=0.1$, or $t=\tfrac13+\tfrac{0.5}{2.1}\approx0.5714$. So the clip is

$$[0.4,\ 0.5714],$$

and 83% of the curve is proved not to meet the other one, in a single step. Bisection would have thrown away at most half, and that only if the other half’s box happened to be clear.

The iteration

One clip narrows $P$. Now restrict $P$ to $[t_0,t_1]$ by de Casteljau ($O(n^2)$, and exact as far as the arithmetic is), and swap roles. The restricted $P$ is shorter, so its fat line is much thinner, and $Q$ is clipped against that. Then $P$ again. Three things can happen to a pair.

• The intersection with the band is empty. The pair holds no crossing: rejected.
• The clip removes a good part of the piece. Restrict and go on. The usual rule, which is the kernel’s too, is that a clip is progress if it removes at least 20%.
• The clip barely removes anything. The piece holds several crossings, or a tangency, or the arithmetic has reached its limit. Halve the longer piece and go on with both halves. Clipping cannot separate two crossings on one piece, since the hull covers both of them, but bisection can.

So the search stays the branch and bound of Lecture 12. Clipping replaces only the “cut at the middle” with “cut where the solutions are”, wherever that is possible.

The readout under the demo compares the search with plain bisection on the same scene. For its default curves, with two crossings and a width of $10^{-9}$ in both parameters, the exact hull took 21 clips and one halving, the chord band 40 clips and 5 halvings. Plain bisection with boxes looked at 501 pairs of pieces.

Why quadratic

Follow a single transversal crossing. The curves in book/scripts/demo_tests/fat-line-clip.test.mjs cross once at a good angle, and here are the widths of the parameter ranges as the demo’s search narrows them, to a target width of $10^{-12}$:

Eight clips in all, against 479 pairs of boxes for plain bisection on the same two curves to the same width. The exponent roughly doubles per round trip: $10^{-2}$, $10^{-5}$, $10^{-12}$.

The argument is short. Let $Q$ have parameter width $h_Q$, so its fat line is $O(h_Q^2)$ thick. Let $P$ be nearly straight, of width $h_P$, so its distance function is nearly linear in $t$, with a slope $\sigma$ determined by the crossing angle, and departs from linear by $O(h_P^2)$. Then the zone where $d(t)$ can lie in the band has width

$$\frac{O(h_Q^2)+O(h_P^2)}{\sigma},$$

which is second order in the sizes. So one clip squares the widths, up to the constant, and it does so as long as $\sigma$ is not small: as long as the crossing is transversal. When $\sigma$ is small, the denominator gives that away. We meet this at tangencies in Lecture 15. At a tangency, clipping stalls and the search falls back to bisection, with all the $h^{-1/2}$ costs of Lecture 12.

That gives a useful comparison. Clipping is a solver with the convergence of Newton and the safety of subdivision: it cannot lose a root, and it does not need a start value. Its cost is rebuilding and restricting pieces, which is why the kernel still hands over to Newton at the end (Lecture 14).

From the published method to the kernel’s

Geop implements this differently in four ways. Each serves one of its rules.

1. No division: the rational case as a polynomial

A NURBS curve is a quotient, $C(t)=H(t)/W(t)$ in homogeneous form. Clipping is about polynomials. It would seem you need to divide, and dividing intervals is what makes them wide (Lecture 11). But weights are definitely positive (NurbCurve::try_new rejects anything else). So for a target coordinate $p_k$,

$$C_k(t)=p_k \iff g_k(t) := H_k(t)-p_k\,W(t)=0.$$

This is a polynomial spline again, with coefficients computed straight from the stored homogeneous coordinates:

$$d_{k,i}=P_{i,k}-p_k\,P_{i,w},$$

one multiplication and one subtraction per coefficient. That is the containment search in core/geop-core-geometry/src/contains/curve.rs, and the opening lines of its clip are nothing but this formula:

d.extend(
    seg.control_points
        .iter()
        .map(|q| q[k].sub(point[k].mul(q[D - 1]))),
);
let zeros = fat_line_zeros(xi, &d)?;

For two curves $A=H_A/W_A$ and $B=H_B/W_B$ the same trick, with cross multiplication, gives a spline in two independent parameters:

$$g_k(s,t)=H_{A,k}(s)\,W_B(t)-W_A(s)\,H_{B,k}(t),\qquad d_{k,ij}=P_{i,k}\,Q_{j,w}-P_{i,w}\,Q_{j,k}.$$

The derivation is in core/geop-core-geometry/src/intersection/curve_curve.md. Equating the two curves at one shared parameter would be a mistake here, since a crossing’s two parameters generally differ.

2. Greville abscissae: B-splines rather than Béziers

A segment in the kernel need not be a single Bézier span, and the abscissae $\tfrac in$ have a B-spline generalization. The Greville abscissae of degree $p$ on a knot vector $u$ are

$$\xi_i=\frac{u_{i+1}+\dots+u_{i+p}}{p}.$$

B-splines have linear precision too, $\sum_i\xi_iN_{i,p}(t)=t$, so the graph $(t,g(t))$ is a B-spline curve with control points $(\xi_i,d_i)$. On a Bézier segment over $[a,b]$ this reduces to $a+\tfrac in(b-a)$. See greville_abscissae in core/geop-core-geometry/src/fat_line.rs.

3. A chord band in place of the hull

Computing an exact hull with intervals is treacherous. As core/geop-core-geometry/src/contains/curve.md points out, a monotone-chain hull’s orientation tests are three-valued under intervals, and both keeping and dropping a doubtful vertex move the chain to the wrong place, which is unsound. An all-pairs form has no such decision, but costs $O(n^2)$. The kernel uses the cheaper one, called Clip B in that note, with all the control points’ uncertainty in one number.

Take the chord of the graph, from its first control point to its last, with slope $s$:

$$s=\frac{d_n-d_0}{\xi_n-\xi_0},\qquad \ell(t)=d_0+s\,(t-\xi_0).$$

Take the union of every control point’s signed vertical offset from it:

$$\Delta=\bigcup_i\bigl(d_i-\ell(\xi_i)\bigr).$$

All graph control points lie in the band $\ell(t)+\Delta$, so a zero of $g$ needs $\ell(t)\in-\Delta$, that is

$$\hat Z=\xi_0-\frac{d_0+\Delta}{s}.$$

This single interval expression is the clip. Evaluated in interval arithmetic it encloses every real choice of the inputs; losing the correlation between $d_0$, $s$ and $\Delta$ can only widen it. It is fat_line_zeros, and there is little more to it than this formula, with two safeguards. First, a test independent of the slope: on the domain, $\ell$ stays between $d_0$ and $d_n$, so if $(d_0\cup d_n)+\Delta$ misses zero there is no root. Second, if $s$ could be zero, the division fails and the axis gives no information, which is said honestly by returning the whole line (ENTIRE).

Note that $\Delta$ is the signed union, generally not symmetric, and contains zero. “Half a tolerance on each side” would be either too loose or, on the far side, unsound.

Let us apply this to the worked example. With the band included, the coefficients are the intervals $D_i=d_i-[-0.1,0.1]$:

$$D=\bigl([-1.1,-0.9],\ [-0.5,-0.3],\ [0.4,0.6],\ [0.9,1.1]\bigr).$$

The slope is $D_3-D_0=[1.8,2.2]$. The offsets $D_i-D_0-\xi_i\,s$ are $[-0.2,0.2]$, $[-\tfrac13,0.2]$, $[-\tfrac16,0.5]$, $[-0.4,0.4]$, so $\Delta=[-0.4,0.5]$. Then $D_0+\Delta=[-1.5,-0.4]$ and

$$\hat Z=-\frac{[-1.5,-0.4]}{[1.8,2.2]}=\Bigl[\frac{0.4}{2.2},\ \frac{1.5}{1.8}\Bigr]=[0.1818,\ 0.8333].$$

Compare the exact hull’s $[0.4,0.5714]$. The chord band is coarser, as it must be, since the band contains the hull. It is much cheaper, however, and it is still second order, so that with the same search only a few extra clips are needed: 8 either way on the one-crossing curves above, where the $s$ widths ran $1$, $0.23$, $9.9\times10^{-4}$, $2.6\times10^{-9}$.

4. The second parameter: union along one index, intersection across equations

The published method alternates: clip $P$ against $Q$’s fat line, then the other way. The kernel has one tensor of coefficients $d_{ij}$ in two parameters, and needs to clip each one. How?

Fix $t$. Then $g(s,t)=\sum_i\bigl(\sum_jd_{ij}M_j(t)\bigr)N_i(s)$, and the inner sum is a convex combination of the row $d_{i0},d_{i1},\dots$. So it lies in the row’s union

$$D^s_i=\bigcup_jd_{ij}.$$

For every $t$, the spline in $s$ has coefficients in $D^s_i$. A clip of the envelope $(\xi_i,D^s_i)$ therefore keeps every zero for every $t$. That is what clip_tensor does, in every direction.

Now look at what this is for two polynomial curves, measured in one direction $n$. The coefficient is $d_{ij}=n\cdot P_i-n\cdot Q_j$, so that

$$D^s_i=n\cdot P_i-\bigcup_j n\cdot Q_j.$$

That is the distance of $P_i$ minus a band: it is the fat line of $Q$ measured along $n$. If $n$ is the normal of $Q$’s chord, the band is exactly Sederberg and Nishita’s fat line. So the tensor formulation is the same method, and fat lines come out of it on their own. Which $n$? Any $n$ is valid, since $g_n=\sum_kn_kg_k$ vanishes wherever all $g_k$ do: a free choice, and a good one is across the other curve. clip in core/geop-core-geometry/src/intersection/curve_curve.rs takes the directions square to $B$’s chord, nearly independent of $t$ and so pinning $s$, then those square to $A$’s (chord, spanning in fat_line.rs). If the directions fail to span space, say for a closed segment without a chord or for two collinear curves, the axes are added.

Why not simply clip the coordinate axes? Because projections lose coupling. Two curves crossing on a diagonal, seen in the $x$ and $y$ axes alone, have boxes overlapping wherever their ranges do. A combination across the chord is nearly independent of the other parameter, so its envelope is thin, which is the whole point. curve_curve.md works through this example.

Two more rules deserve to be memorized, since they are correct for a reason that is not obvious.

Union over eliminated indices, intersection over equations. Each row of a tensor is an alternative possibility, so the rows are combined by union. Each equation is a constraint that must hold, simultaneously, so the clips are combined by intersect. Intersecting per-row clips would be wrong: as the doc of clip_tensor says, “a zero of the tensor need not be a zero of any row”. This is the contrast from Lecture 5: union when you must cover both possibilities, intersect when both constraints must hold at once.

intersect of disjoint enclosures is not empty. It returns an input. So disjointness is checked first, every time, with could_be_equal. Look at clip_tensor or clip for this pattern; it would be a nasty bug to forget it.

Restriction and the shared schedule

Curve containment, surface containment and the two intersection searches share one schedule for what to do with a clipped box, in fat_line.rs.

restriction decides whether to restrict. A direction is cut when its clip removed at least 20% of its domain (min_progress is $\tfrac45$) and one of the clip’s sharp outer bounds lies strictly inside the domain. Note what is and is not sharpened here, by the rules of Lecture 5. The restriction cuts at the clip’s outer bounds, Scalar::lower and Scalar::upper: any cut outside the enclosure loses nothing, and a cut through it would. A halving, by contrast, uses a sharpened midpoint, since splitting may happen anywhere.

stalled decides, for a box clipping could not shrink, whether it has converged or which direction to bisect, largest spatial extent first (parameter widths of differently scaled domains are not comparable). A direction with nothing left to give is skipped. Those skips have a history worth repeating. At a pole, where the solutions are a whole row and clipping closes in on it by a constant factor, it did so for thousands of steps, down to boxes of $10^{-280}$, which is why restriction stops once the spatial extents are within the width the data carries (carried_width).

And note how min_subdivision_size enters: it only bounds how far bisection goes. Clipping is never cut short by it, so a point merely within min_subdivision_size of a curve is still rejected if clipping can reject it. That makes a search’s answer more precise than its handoff, without any extra thought.

Recap

Exercises

1. Compute the clip for $d=(-2,\,-1,\,1,\,3)$ against the band $[-0.5,0.5]$, both by the exact hull and by the chord band. Check your answer in a calculator or by reasoning on the picture.

2. In the demo, try to make clipping stall. What kind of curves make it need many halvings? Make a pair that has two crossings close together. What do the first clips do?

3. Show that a Bernstein polynomial $\sum_i\tfrac inB_i^n(t)$ equals $t$, using the derivative formula of Lecture 6. Conclude the graph claim.

4. Why does fat_line_zeros test $(d_0\cup d_n)+\Delta$ before dividing, and why would a test of $d_0+\Delta$ alone be unsound? Construct an example.

5. Consider a pair of curves in 3-D that do not cross. Could the tensor clip fail to reject them, however far subdivided, if it only used the three coordinate axes? What about in-plane directions?

6. A colleague suggests a speedup: instead of the union over $j$ in clip_tensor, use the row at the middle of the other parameter, since “it is the most representative”. Show with a small tensor that this can reject a true solution.

Lecture 14: Newton to refine

Lecture 12 and Lecture 13 leave us with a small box known to hold every solution near it. For many purposes that box is too wide: a split parameter that wide makes a fat vertex and a fat sub-curve, and the width then grows (we saw why in Lecture 9). More subdivision would narrow the box by a bit per level. There is a better tool for the last stretch.

Newton in one variable

Near a point $x_k$, a smooth function looks like its tangent:

$$f(x)\approx f(x_k)+f'(x_k)(x-x_k).$$

Solve the right-hand side for zero, and call that the next guess:

$$x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)}.$$

For $f(x)=x^2-2$, starting from $x_0=1.5$:

The number of correct digits doubles every step. Here is why. Let $x^*$ be a root with $f'(x^*)\ne0$, and let $e_k=x_k-x^*$. Taylor’s theorem around $x_k$ gives

$$0=f(x^*)=f(x_k)-f'(x_k)e_k+\tfrac12f''(\zeta)e_k^2.$$

Divide by $f'(x_k)$, and rearrange using the definition of $x_{k+1}$:

$$e_{k+1}=\frac{f''(\zeta)}{2f'(x_k)}\,e_k^2.$$

So the error is squared, up to the constant $\lvert f''/2f'\rvert$, as long as $f'$ stays away from zero. Three caveats hide in this one line, and each of them has its own section below or in Lecture 15.

1. It is a statement about a neighbourhood of a root. Newton finds nothing: from far away it may go anywhere.
2. The constant blows up where $f'(x^*)$ is small. At a double root the method degrades to halving the error per step.
3. In a kernel the data has width, and so the answer does. The iteration must not pretend otherwise.

This is the division of labour from AGENTS.md’s “Subdivide to isolate, then Newton to refine”, which the kernel follows everywhere. Subdivision reaches machine accuracy in about 50 levels, of which about 7 isolate the solution. Newton then polishes it in a dozen iterations, and no more accuracy is asked of the subdivision than separating the solutions.

Newton on a system

For $F:\mathbb R^n\to\mathbb R^n$ the tangent is the Jacobian $J=\partial F/\partial x$, and each step solves a linear system:

$$J(x_k)\,\delta=-F(x_k),\qquad x_{k+1}=x_k+\delta.$$

The convergence proof is the same with norms. For a curve crossing a surface, $F(t,u,v)=C(t)-S(u,v)$ has three equations in three unknowns, with

$$J=\begin{bmatrix}C'(t)&-S_u(u,v)&-S_v(u,v)\end{bmatrix},\qquad \det J=C'\cdot(S_u\times S_v).$$

The determinant says the geometry out loud. It vanishes when the curve’s tangent lies in the surface’s tangent plane (a tangency), or where $S_u\times S_v=0$ (a pole). That is refine_crossing in core/geop-core-geometry/src/intersection/curve_surface.rs, which builds exactly this 3×3 system.

Two curves in 3-D give more equations than unknowns: three equations in $(s,t)$. The system is overdetermined, so refine_crossing in intersection/curve_curve.rs steps by least squares, solving the normal equations $(J^\top J)\delta=-J^\top F$. That is Gauss-Newton: it uses $J$ only.

Projection: the problem that teaches the most

Finding the point on a curve or surface nearest a given target is a kernel’s most common iterative computation. It is NurbSurface::project in core/geop-core-geometry/src/nurb_surface/project.rs, called on every sample of every fit_pcurve. So let us look carefully.

Take a curve $C(t)$ and a target $T$. The nearest point is a critical point of

$$\phi(t)=\tfrac12\lVert r(t)\rVert^2,\qquad r(t)=T-C(t),$$

so it satisfies $\phi'(t)=-r\cdot C'=0$: the residual is square to the curve. Now there are two reasonable steps.

Gauss-Newton linearizes the curve: $C(t+\delta)\approx C+C'\delta$, and minimizes $\lVert r-C'\delta\rVert^2$:

$$\delta_{\text{GN}}=\frac{r\cdot C'}{C'\cdot C'}.$$

This is the projection of the residual on the tangent. No second derivative is needed, which is its attraction.

Newton linearizes the condition $\phi'=0$: $\phi''=C'\cdot C'-r\cdot C''$, so

$$\delta_{\text{N}}=\frac{r\cdot C'}{C'\cdot C'-r\cdot C''}.$$

The difference is the term $r\cdot C''$. It is the residual’s length times the curvature, so it is small when the target is near the curve, and not small otherwise.

The circle, in closed form

Let $C(t)=R(\cos t,\sin t)$, and put the target at distance $\rho$ from the centre, at angle $\theta$. Let $e=\theta-t$ be the angle still to go (the foot point is $t=\theta$). Then $r\cdot C'=R\rho\sin e$, $C'\cdot C'=R^2$, and $r\cdot C''=-R\rho\cos e+R^2$. So

$$\delta_{\text{GN}}=\frac{\rho}{R}\sin e,\qquad \delta_{\text{N}}=\tan e.$$

The errors evolve as

$$e_{\text{GN}}'=e-\frac\rho R\sin e\approx\Bigl(1-\frac\rho R\Bigr)e,\qquad e_{\text{N}}'=e-\tan e\approx-\frac{e^3}{3}.$$

Write $\rho=R+d$, with $d$ the target’s distance outside the circle. Then Gauss-Newton multiplies the error by $-d/R$ per step.

That is the whole story of Gauss-Newton on a curved surface. Its convergence is linear, at the rate distance over radius of curvature. A target on the surface converges at once. A target half the radius away, with $-\tfrac12$ per step. A target exactly one radius off has a multiplier of $-1$: to first order it stays at its distance, swinging from side to side, and only the next term of the sine wears it down, too slowly for any iteration budget (from $e_0=0.5$ the error is still $0.24$ after 20 steps). Beyond that the foot point is repelling, and the iterates swing out until they settle in a cycle.

Newton, on the circle, is cubic: the second derivative here is so well behaved that even the second-order term cancels. On a general curve it is quadratic.

Let us check with numbers, with $R=1$, a target at $\rho=3$, so twice the radius off, and seed $e_0=0.2$. (The demo’s algo, and these claims, are tested in book/scripts/demo_tests/newton-foot.test.mjs.)

Gauss-Newton ends in a 2-cycle: $e\mapsto-e$ with $2e=3\sin e$, which gives $e\approx1.4958$. It will swing between the two sides for ever.

A bug this explains

AGENTS.md’s section “Newton needs the second derivatives where it is far from flat” has this from real life. NurbSurface::project used to step by Gauss-Newton, the first fundamental form alone. A pocket’s rim rounded with the radius its corner’s fill was rounded with put the ball’s centre at just such a distance, and meet, the iteration that moves a vertex onto offset surfaces, ran off along the fill’s end. As the code’s comment puts it: a target twice the radius off a cylinder “swung from side to side of its foot point for ever”.

The current project takes Newton’s step wherever its Hessian is definitely positive definite, and Gauss-Newton’s where not:

$$H=\begin{bmatrix}S_u\cdot S_u-r\cdot S_{uu}&S_u\cdot S_v-r\cdot S_{uv}\\S_u\cdot S_v-r\cdot S_{uv}&S_v\cdot S_v-r\cdot S_{vv}\end{bmatrix}.$$

Why fall back at all? Newton’s step targets any critical point. If $H$ is not positive definite, the iterate is not near a minimum of the distance, and Newton could step to a saddle or a maximum. Gauss-Newton’s matrix (the first fundamental form) is always positive semi-definite, so its step is a descent direction. In interval arithmetic “definitely positive definite” is a real test: h11.definitely_greater(Scalar::ZERO) and the determinant likewise.

Clamping, and poles

Two practical matters.

The domain. A surface’s domain is a rectangle. An iterate can step outside on its way in, and evaluating a B-spline outside its knots is not defined, so project clamps into domain_u and domain_v before each step. A foot point clamped at the edge means the nearest point is on the boundary (the true minimum may be beyond, on the surface carried on).

The pole. At a pole a row of control points collapses to a point, so $S_v=0$ there, and the first fundamental form is singular (see Lecture 10). A seed on the pole would then never move, whatever the target. As AGENTS.md’s “Degenerate parametrizations” notes: only the collapsed parameter stops moving the point. The other still moves it, so step along that one, down a meridian. Which meridian is a free choice, since at the pole the collapsed parameter says nothing, and a free choice should be a good one: of a few candidates, off_pole takes the one along which the step goes furthest towards the target.

The first seed has to come from a search

Newton converges to whichever foot point it is near. On a curved patch that may be a local minimum, or one clamped against a domain edge where the residual is merely square to the boundary. So NurbSurface::locate first proves a box with surface_could_contain, and only then refines with project. That is subdivide-then-Newton in one function. Similarly fit_pcurve seeds its first projection globally, unless its caller already knows the (u, v), and every later one from its predecessor.

And the same reasoning makes NurbSurface::closest_point call least_distance first, then refine each surviving patch (Lecture 12). The search proves no basin was missed. Newton’s job is the digits.

Honest refinement

So far this was classical numerics. Now the interval rules of Lecture 5 meet it, and all of this section is hard-won.

Suppose every number is an interval. Run $x_{k+1}=x_k-f(x_k)/f'(x_k)$ on intervals, and the iterates widen: $x_k$ appears twice, the two appearances are uncorrelated, and the width feeds forward. Left alone it compounds over a few iterations until evaluate fails: AGENTS.md measured this at 36 of 175 test scenes for the tracing corrector.

The way out is to notice that an iterate is a seed, and a seed is a free choice. Any point of the iterate serves equally well to start the next step from. So sharpen it: replace the interval with its midpoint, and take the next step from that one point. That removes the compounding at no cost.

But the last iterate is not a seed. It is the answer. Its width says how well the data pins down the root: the residual $F(x_k)$ inherits the data’s width, the step inherits the residual’s, and that is the honest statement of how precisely an uncertain target determines a foot point. Sharpen the last one too, and we return one sharp $(u,v)$ for a target that only ever determined a range of them. Then every later comparison against anything else derived from that uncertain point fails.

This is no theoretical nicety. NurbSurface::project used to sharpen every iteration including the final one. Fixing just this took validate_fast failures across the scene sweep from 21/175 to 1/175, as AGENTS.md records.

So keep two copies, one as the seed and the other as the answer. Here is the idea in full, from core/geop-core-math/src/solvers/iterate.rs:

pub fn iterate<X: Iterate, E>(
    seed: X,
    answer: X,
    steps: usize,
    repeat: Repeat,
    settle: impl Fn(&X) -> (X, X),
    step: impl FnMut(&X) -> Result<Step<X>, E>,
) -> Result<Iterated<X>, E>

Every Newton-type loop in the kernel runs through this: foot points, a ball rolled between faces, a vertex moved onto offset surfaces, a traced curve’s corrector. step takes a seed to the next iterate. settle turns that iterate into two values, the sharp seed of the next step and the honest answer it gave. The returned width is that of a single step from a sharp seed: bounded, not compounded, and not discarded.

A repeating seed is done

Newton with sharp seeds rarely hits an exact fixed point in floating point. It usually enters a cycle in the last bit, and it would run out its iterations going round. But the map from a sharp seed to the next is deterministic, so once a seed repeats exactly, all the rest is periodic and the end result known. iterate with Repeat::Cycle compares bit for bit (Iterate::same), and stops there, with the identical result, in less than half the iterations (AGENTS.md, “Count the work”). This is no convergence tolerance: the seed is bit for bit one seen before.

The refinement a search hands to

The generic version sits in core/geop-core-math/src/solvers/newton.rs:

pub fn refine<const N: usize>(
    start: [Scalar; N],
    domain: [(Scalar, Scalar); N],
    step: impl Fn([Scalar; N]) -> Option<[Scalar; N]>,
) -> [Scalar; N]

Read the design choices, each of which comes from a rule.

• It returns a value, not a Result. Refinement is infallible: a singular Jacobian at a tangential crossing or a pole, a step that comes back non-finite (over a singular Jacobian a step can be everything, and its middle is no number), or a refined box disjoint from the start, returns the incoming box unchanged, and an iterate leaving the domain is clamped back into it, since a foot point can step outside on its way in. “Refinement may only tighten, never fail, or callers start needing fallbacks of their own” (AGENTS.md).
• The result is start.intersect(x). The search proved the root lies in start, and Newton computed another enclosure of the same root. Both hold, so their intersection is valid and tighter than either: exactly the case Lecture 5 says intersection is for.
• Disjoint means reject. If Newton ended somewhere that does not overlap the proved box, it wandered off, and start is the answer.
• Clamp after sharpening. Clamping an interval that straddles a domain bound leaves it straddling, because its three-valued comparisons have no answer. Clamping a sharp value is exact.

ITERATIONS is 12. Quadratic convergence makes that plenty, and it cannot change which solutions are found.

The doc comment adds a warning: this is “plain, unverified Newton: it does not itself certify existence or uniqueness, so callers check that the refined box can still hold what they are looking for”. And they do: refine_crossing evaluates both shapes over the refined box and returns the incoming one unless they could coincide. AGENTS.md has a term for this habit: check an iterative answer against every condition it was asked to meet.

Refine the target first

A trap is worth describing separately. A refinement is only as good as its target. Refining a pcurve parameter against a $(u,v)$ that came straight from a search, itself min_subdivision_size wide, gives a wide residual, a wide step, and a result no narrower than it started. Silently, because the path that fails looks identical to success. The target has to be refined first (NurbSurface::project on the vertex point). The same applies to a vertex built by evaluating a wide crossing parameter: fix the parameter at the source, or everything downstream inherits its width.

Refine where the parameter is used

curve_surface_intersect deliberately does not refine everything it returns. Most callers only need to know where and how many, and refining changes results they already agree with. As AGENTS.md notes, doing it unconditionally broke offset_sphere_is_valid. So refining is done when the parameter becomes a split point, where width genuinely matters.

That leads to the kernel’s rule for cutting a curve at a vertex, which illustrates everything above in one function: split_parameter_at_point in core/geop-core-geometry/src/nurb_curve/refine.rs.

1. Refine the parameter with Newton on $C(t)=$ the vertex.
2. Refuse if it could be an end of the domain: the vertex is then an end of the curve, and there is nothing to cut.
3. Choose where to cut, a sharp middle of the refined enclosure, or a knot if one lies in it. That is legitimate: every parameter of that enclosure ends a piece on the curve and within the vertex, so it is a free choice.
4. Check the value you are about to use. Evaluate at the sharp parameter, and refuse unless the curve reaches the vertex there.

Step 4 comes from a bug that AGENTS.md describes under “Validate the value you are about to use, not a wider one”. Model::split_edge_at_vertex checked curve.evaluate(edge_t).could_be_equal(&vertex_point) on the wide edge_t, a fat box that genuinely contains the vertex, and then split at the sharpened edge_t. The check answered “is the vertex somewhere along this stretch?” and the split asked “cut exactly here”. A different question.

Interval Newton and Krawczyk: refining with a proof

All the above gives a narrow box, and a check that it could hold a root. It does not prove that a root exists. Sometimes we want a proof, and interval analysis has a neat way of getting one.

Take the mean value theorem. For a root $x^*$ in an interval $X$ and a point $m\in X$,

$$0=f(x^*)=f(m)+f'(\zeta)(x^*-m)\quad\text{for some }\zeta\in X.$$

Solving for the root, with $f'(\zeta)$ enclosed by $F'(X)$ over the interval:

$$x^*\in N(X)=m-\frac{f(m)}{F'(X)}.$$

This is the interval Newton operator, and it gives three outcomes.

• If $N(X)\cap X$ is empty, then $X$ holds no root.
• Otherwise every root in $X$ lies in $N(X)\cap X$: an honest contraction.
• If $N(X)\subset X$, then $X$ holds exactly one root: one exists by Brouwer’s fixed point theorem applied to the operator, and there is no second since $F'(X)$ does not contain 0.

Example: $f(x)=x^2-2$ on $X=[1.4,1.5]$. Take $m=1.45$, $f(m)=0.1025$, $F'(X)=[2.8,3]$. Then $f(m)/F'(X)=[0.03417,0.03661]$, and $N(X)=[1.41339,1.41583]\subset X$. We have just proved that $\sqrt2$ is in that narrow interval, using finitely many interval operations.

For systems, dividing by an interval matrix is awkward. Krawczyk’s form avoids it, with a fixed preconditioner $Y\approx J(m)^{-1}$:

$$K(X)=m-YF(m)+\bigl(I-YJ(X)\bigr)(X-m).$$

Again every root in $X$ lies in $K(X)$, and if $K(X)\subset X$ there is exactly one. $Y$ is a free choice: the proof is valid for any invertible $Y$, and only its tightness depends on how close it is to the true inverse. So it can be sharp. In the example, $Y=1/2.9$, $I-YJ(X)=[-0.0345,0.0345]$, and $K(X)=[1.41293,1.41638]\subset X$.

Recap

Exercises

1. Derive $\delta_{\text{N}}=\tan e$ for the circle, and show that $e-\tan e=-\tfrac{e^3}{3}+O(e^5)$. Then run the demo’s ellipse and show from its readout that the error there is squared rather than cubed.

2. What does Gauss-Newton do on a concave side, that is for a target inside the circle at $\rho<R$? For which $\rho$ is it fastest? What about $\rho=0$?

3. Why would a “convergence tolerance” be wrong in project? Suppose we stop once $\lVert\delta\rVert<10^{-12}$. Give one way this goes wrong for a small model and one for a large model.

4. Consider $f(x)=x^2$ with the double root at 0. Starting from $x_0=1$, what does Newton do? What does this say about refine_crossing at a tangency? Compare with the iteration counts in the readout of the tangency demo.

5. Do one Krawczyk step by hand for the system $x^2+y^2=1$, $x=y$, on $X=[0.6,0.8]^2$, with $m=(0.7,0.7)$, and $Y$ the inverse of $J(m)$. Does $K(X)\subset X$ prove a unique root?

6. newton::refine sharpens before clamping. Show why the other order can produce an iterate outside the domain.

Lecture 15: Tangency and overlap

The last three lectures assumed that two curves cross at a good angle. Real CAD models are full of the opposite. A fillet is tangent to the faces it blends, by definition. A hole drilled to the same radius as a boss has coincident walls. Two parts are mated with faces flush against each other. So a kernel’s hardest cases are not rare accidents. They are what users do on purpose.

This lecture is about the geometry of those cases, and what an honest kernel can and cannot say about them.

A crossing is an equation, and its condition number is the angle

Consider the simplest possible model. A curve $y=t^2-c$ crosses the line $y=0$, and we know its height only to within $\pm\varepsilon$: that is, a coordinate of the curve is an interval of half-width $\varepsilon$. Where can the crossing be?

$$t^2-c\in[-\varepsilon,\varepsilon]\iff t^2\in[c-\varepsilon,c+\varepsilon].$$

Far from tangency, when $c\gg\varepsilon$, there are two crossings. Each lies in an interval around $\pm\sqrt c$ of width

$$\sqrt{c+\varepsilon}-\sqrt{c-\varepsilon}\approx\frac{\varepsilon}{\sqrt c}.$$

The crossing angle has slope $2\sqrt c$, so this is the data’s whole uncertainty $2\varepsilon$ divided by the slope. A shallower crossing is less well determined, in proportion.

At tangency, with $c=0$, there is one contact and its enclosure is

$$t\in[-\sqrt\varepsilon,\sqrt\varepsilon],\qquad\text{width }2\sqrt\varepsilon.$$

A coordinate known to $10^{-12}$ pins the contact point to $2\times10^{-6}$. In double precision with $\varepsilon\approx10^{-16}$, half of the digits are gone before any algorithm starts. No method can recover them, because the data does not contain them. Every parameter in that range is equally a true contact point of some curve consistent with the data.

In between, when $c\approx\varepsilon$, there is a question that has no answer. With $c=10^{-13}$ and $\varepsilon=10^{-12}$: are there two crossings, or one, or none? All three are consistent with the data. An honest search returns one interval, wide, and it may hold anything.

Here are numbers from that search (checked against the exact set $t^2\in[c-\varepsilon,c+\varepsilon]$ in book/scripts/demo_tests/tangency-width.test.mjs).

Read each row. Each one is $\varepsilon/\sqrt c$ until $\varepsilon$ reaches $c$, where it merges into one interval of width about $2\sqrt{c+\varepsilon}$. That is the whole theory. What about geometry in more dimensions? A surface bending away from another along a curve is the same picture in a cross-section. Near any tangent contact, the two sides separate quadratically, and so the same square root appears.

What each tool does at a tangency

The search (Lecture 12). Leaves spread across a stretch of about $\sqrt h$, so there are $h^{-1/2}$ of them, all touching. They merge into one cluster, and the cluster’s width is $\sqrt h$, not $h$. See the tangent scene of Lecture 12’s demo. This is the same square root again, this time from the search’s own resolution rather than from the data.

The clip (Lecture 13). A clip divides by the slope of the distance function, and at a tangency that slope is about the size of the piece, not a constant. So the zone it keeps is $O(h^2)/O(h)=O(h)$, a constant fraction of the piece. Clipping stops being quadratic, it fails the 20% rule, and the search goes back to halving. In core/geop-core-geometry/src/intersection/curve_curve.rs, a comment on the converged branch of curve_curve_crossings describes the outcome:

// The pieces' whole domains, not the tighter clip: pieces the
// search could not separate — a tangency converges on a chain
// of adjacent ones — must merge into one unresolved cluster,
// and adjacent domains share an endpoint where clips need not touch.
solutions.insert((a.domain_as_scalar(), b.domain_as_scalar()));

The boxes are merged by union (a DisjointSet), never averaged. The answer is a wide box and the statement “somewhere in here”.

Newton (Lecture 14). For the double root of $t^2$, a step is $t\mapsto t-t^2/2t=t/2$. The error halves per step, down from squaring: in the demo’s readout, Newton from $t=0.7$ to a step of $10^{-14}$ takes 6 iterations for $c=0.25$, 11 for $c=10^{-4}$, and 46 at $c=0$. With intervals it fares worse: the Jacobian’s enclosure soon contains zero, the linear solve fails, and newton::refine returns the incoming box. This is correct. Refinement promised to tighten or leave alone, and here it leaves alone.

So at a tangency the tools all agree. They report a wide, honest box. What is then to be done with that box is up to the caller, and every caller has to be written with it in mind. Which is why a kernel should try to not make tangent contact where it can choose (“A tilt”, below), and try to make the topology answer where it cannot.

Overlap: why counting does not work

A tangency is a crossing that has gone soft. An overlap is a different thing: two curves running along the same arc, or a curve lying in a surface, so that the solution set is not points at all.

Lecture 12 showed that an overlap leaves $O(1/h)$ leaves. It is tempting to detect it that way: “if the search returns more than max_solutions answers, the curves must coincide”. That was indeed the code’s approach once, and it is wrong in both directions, which AGENTS.md lists under “Proximity is not membership”: the same shape of bug appears in find_coincident_pair reading a solution count as coincidence.

• A transversal crossing at a very fine handoff, among curves with width, also leaves many leaves. Lecture 12’s signal depends on stopping early.
• A short overlap, or a coarse handoff, leaves few leaves. So a partial overlap can be missed.
• Many isolated crossings, say a wavy curve along a line, look like an overlap by count, and are not.

And a search that hits its budget hasn’t answered at all.

The direct test

There is a much better way, based on a piece of algebra. A rational curve is an analytic function of its parameter. If two rational curve pieces coincide on any interval, they lie on the same algebraic curve throughout. So a stretch of overlap between two single-piece curves cannot just stop in the interior of both. It can only begin or end where one of the curves itself ends.

This means there are finitely many candidate places, and you know where they are:

• an end of curve $A$ that lies on curve $B$;
• an end of curve $B$ that lies on curve $A$.

These are point-containment questions, which we solved in Lecture 12. Sort the candidates along $A$. Between two neighbouring candidates, $A$ either lies on $B$ throughout, or nowhere in the interior. So one probe at the midpoint decides it.

That is find_overlaps in core/geop-core-geometry/src/intersection/coincidence.rs, and curve_curve_overlaps in curve_curve.rs is the candidate list. (The module states this reasoning for single-piece geometry, and the candidates are the ends of the curves as given. A curve of several spans is one analytic function only span by span.)

Let us do an example. Let $A$ be the straight curve from $(0,0)$ to $(1,1)$ with $s\in[0,1]$, and $B$ from $(0.5,0.5)$ to $(2,2)$.

1. Candidates. $A(0)$ is not on $B$. $A(1)=(1,1)$ is on $B$. $B$’s start $(0.5,0.5)$ is on $A$ at $s=0.5$. $B$’s end $(2,2)$ is not on $A$. So the candidates are $s=0.5$ and $s=1$.
2. Probe $s=0.75$: the point $(0.75,0.75)$ lies on $B$. So $[0.5,1]$ is an overlap.
3. The crossing search then runs on the rest of $A$, here $[0,0.5]$, where it is entitled to assume no overlap.

For a curve against a surface, curve_surface_overlaps adds the surface’s boundary. A stretch of a curve lying in a patch can end at the curve’s own ends, or where it crosses one of the patch’s four boundary curves, and those are found by the curve-curve search again.

The result is the variant of Intersections we saw: Found (no overlap, isolated crossings only), or Coincident, which lists the overlaps’ ends, then crossings on the rest, and then points spread evenly over the overlaps.

The module’s documentation is honest about the limit: “The midpoint probe is the one place this is a test rather than a proof: a curve touching the other object exactly at two candidates and the midpoint between them, while leaving it in between, is read as an overlap.” Note how little that leaves. It is one specific algebraic coincidence, stated in the code, not a tolerance.

Overlap needs the same direction

The test of the midpoint, as stated, asks whether a point lies on the other curve. Enclosures can answer yes to that for two curves that plainly cross. A helix met a cylinder’s rim at about three degrees, in the example square_thread, and validation called the two edges coincident along an arc. Each carried the width its construction left it, each end was found on the other curve at a parameter of its own, and those two ends lay no further apart than a shared vertex’s width. Probed between them, on so short a stretch, one curve “lies on” the other at any angle, as far as enclosures tell.

The geometry has more to say: two curves lying on each other run the same way. So curve_curve_overlaps also asks that the tangents could be parallel there, same_way, every 2-by-2 minor of the two tangents’ components having to enclose zero. A caller who gave a reach above zero is exempt, since a curve within reach of another counts as one with it whichever way it runs. The test edges_meeting_at_a_shallow_angle_do_not_overlap in curve_curve.rs states what must hold, with two straight edges 1e-10 wide meeting at three degrees, and two lying on each other for comparison. curve_surface_overlaps has no such test.

Note the shape of this fix. It adds no tolerance: it checks another condition an overlap has to meet, in the same three-valued arithmetic.

Asking the right question about membership

After a boolean starts to split faces, one more trap appears. Splitting a face leaves both halves on one surface, so “this edge lies on this face’s surface” no longer means it belongs to this face. An edge already imprinted into one half still lies on the surface of the half it does not belong to, and an earlier find_coincident_pair re-imprinted it, splitting again, forever. The fix is to ask the question about the face’s trim, whether the edge lies inside it (face_contains, Lecture 19), and to first rule out the exact topological case, edge_is_boundary_of_face. AGENTS.md tells this story in “Make the invariant a type, not a convention”.

This is the general form of the rule we keep meeting: proximity is not membership. When you select something because it is near or on, ask what property you are about to assume it has, and test that.

Tangent branches

The subtlest configuration in this lecture arises when a revolve and an extrude are built from one sketch. The revolved surface at its seam and the extruded wall both contain the profile and the direction perpendicular to the sketch. So they touch tangentially all along the shared profile edge. Their intersection contains that edge, and it may also contain other curves that leave the edge transversally. Where such a branch leaves, an intersection curve ends, and so a vertex is needed. But the edge lies on the other surface throughout, so no piercing search can single out that point.

core/geop-core-geometry does not solve this; lib/geop-lib-booleans/src/remesh/remesh_tangent_branches.rs does, and its module doc has the derivation. Let $s$ run along the edge $E$ and $h$ across it, in the common tangent plane. One surface stands off the other by a height $f(s,h)$. Tangency along $E$ means that $f$ and $\partial f/\partial h$ both vanish at $h=0$. So $f$ has a double factor:

$$f(s,h)=h^2\,g(s,h).$$

The surfaces meet where $f=0$, that is on $\{h=0\}\cup\{g=0\}$. The first set is $E$ itself. The second is the other branches, and they meet $E$ where

$$g(s,0)=\tfrac12\frac{\partial^2f}{\partial h^2}(s,0)=0.$$

The second derivative of the height difference across $E$ is the difference of the two surfaces’ normal curvatures across $E$. Along $E$’s own direction the normal curvatures agree, as both surfaces contain $E$, and they share $E$‘s normal. So the points are those where the two surfaces’ mean curvatures agree.

A geometric puzzle has turned into a scalar function along an edge, and a search finds its sign changes (sign_changes in core/geop-core-math/src/solvers/bisect.rs): subdivide along the edge to isolate them, then bisect. This is a fine example of what to do with a degeneracy. Do not widen a tolerance until the search stops failing; find the exact equation that holds on the degenerate set, factor out what you already know, and search the remainder.

A related search is in core/geop-core-geometry/src/intersection/surface_surface.rs. A boolean finds the curves in which two faces meet by tracing from where edges pierce faces (Lecture 16), and a closed loop lying inside both faces is pierced by no edge. So turning_points looks for points where the curve turns back along a chosen direction $d$, solving $S_a=S_b$ together with $(n_a\times n_b)\cdot d=0$. Where the normals could be parallel, that equation is vacuous, so those pairs are set aside as tangent, and refine_node tells apart isolated nodes (where branches cross) from contacts along a curve. Its doc is exact about what remains: “a loop shrinking into a tangency is the one case left to the caller to notice, as a face the boolean cannot classify”.

Do not make the degenerate case yourself

So far, tangency arrived from outside. Often we cause it ourselves. AGENTS.md has a whole sequence of stories on this, and they generalize to one question: did something our code chose put the entities there?

Name the entities before you name the cause

A revolve cut “up to next” failed with a vertex too wide to validate, and was explained to the user as a near-tangency. The user looked at the scene: there was none. Naming the vertex showed that it lay where the drilled hole’s rim crossed the revolve’s end face, and the end face was there because the operation’s own fallback had stopped the tool flat at the first point of contact, an existing corner. The bug was a logic error in our own code, not numerics.

“Near-tangent crossing” sounds plausible and cannot be checked. Before you state a cause, resolve the failing vertex or edge to its entities by name (names_mentioned, used by the regression tests’ assert_builds_valid), and test the explanation against them. Lecture 42 has more of this.

Never place a cut on a corner that already exists

The same story generalizes. When an operation chooses where a new face goes, landing it exactly on an existing vertex or edge asks the boolean to cross an edge at its own end point. That configuration is degenerate, and no enclosure stays tight through it. So choose positions the geometry actually defines. For up to next, that means going up to the next face from where the tool enters the target (or all the way round, or through all), rather than stopping flat at a sampled contact point.

A tool’s own curves must not lie in the solid’s faces

A blend tool meets the solid where the boolean finds it. Any curve of the tool lying in a face of the solid, crossing that face’s boundary, asks the boolean to resolve an overlap rather than a crossing, and it exhausts its search. AGENTS.md lists four instances, all from geop-ops-fillet. Here is the most instructive one. A straight chamfer built as a tool once used its chord’s ends on the faces: exact lines lying in them. Extending the chord past the faces, as the swept chamfer always did, made it cross them instead. In other words, the tool was made to cross rather than touch, at the cost of nothing, because the extra part of the tool is outside the solid and is dropped.

That is a general trick. A boolean is robust on transversal crossings and fragile on contacts, and a tool is yours to shape. Where a patch ends is a free choice. And the last of the four instances adds a warning about how to make it: reaching a third of a patch’s size past a section only moved its side onto a stepped block’s chamfer contact, because “a fraction of one design dimension lands on another”. The cap now reaches past the extent of every edge of the model (FlatPatch::past), where nothing can lie along it.

An approximate tangency dips through the face

Here is the case with the most interesting cure. A surface only approximately tangent to a face, say a rolling-ball blend skinned through exact stations, is tangent at each and interpolated in between. It then crosses the face as often as not between the stations, by a hair. A boolean then has to cut the face along each such crossing: near-duplicate vertices next to the contact, “nothing to split” errors. As AGENTS.md puts it: tangency to within the approximation’s error is no tangency to the boolean.

The cure, in ops/geop-ops-fillet/src/rolling.rs, is to leave each face at a millionth of a radian towards the ball (TILT there, equal to DEVIATION, $10^{-6}$), and to check between stations that it really leaves on that side. A tangency is thereby turned into a very shallow but definite crossing, along the contact curve alone. The boolean can do that, as Lecture 13’s clip has a slope to divide by, if a small one.

And an exact tangency needs no tilt, and must not get one: AGENTS.md records that, tilted, a revolved handle’s rim showed two silhouettes per face where the crease ran, and a draft no longer saw the fillet as tangent. Where the untilted skin meets the ball between stations to rounding (span_is_exact), it is kept untilted. Lecture 29 develops all this.

Recap

Exercises

1. Derive the width $\varepsilon/\sqrt c$ and the merging point. For $\varepsilon=10^{-12}$, at what $c$ do the two enclosures of the demo merge? Check with the demo.

2. The surfaces $z=\tfrac12\kappa_1x^2$ and $z=\tfrac12\kappa_2x^2$ touch along the $y$ axis. Given data uncertainty $\varepsilon$ in $z$, how far from the axis could the two be said to meet? What does a containment search on their pieces see, and what would make refine_crossing on a curve in one of them give up?

3. Show by example that the midpoint probe can be fooled: construct two curves that coincide at the two candidates and the midpoint, but are not the same between. Why is this not a practical problem for two lines? For a line and a conic?

4. Go through $f(s,h)=h^2g(s,h)$ for two cylinders of radii $a$ and $b$, axes parallel to $y$, touching along a line, one with radius changing linearly along $s=y$, so that $b=b_0+\gamma y$, with $b_0=a-\gamma y_0$. Where does a transversal branch leave the common line, approximately? (You can ignore that the second is then only a cone.)

5. A colleague fixes a “nothing to split” error by moving a fillet’s boundary 1e-9 off the face. Explain why that differs from a tilt, and give a configuration where it fails.

6. Look at the tangent scene in Lecture 12’s demo and the demo above. Both have square roots. What is the same and what differs?

Lecture 16: Tracing and pcurves

An edge crosses a face at a point. A face crosses a face along a curve. And we have no formula for that curve: two bicubic patches meet in a curve of very high algebraic degree, and it has no NURBS form of any low degree. Yet the boolean’s result needs an edge there, with a NURBS curve for it, and a $(u,v)$ curve on each face.

So the curve has to be traced: followed numerically, step by step, from one end, then fitted. Tracing is where Lectures 12 to 15 come together in the kernel’s most delicate loop. As the developer notes in DEVELOPERS.md record, an earlier generation of the kernel got stuck here: “Booleans looked promising but tracing was too difficult, too many edge cases”. This lecture goes through those cases one at a time, because each has a reason.

The equations

Let $S_a(u_a,v_a)$ and $S_b(u_b,v_b)$ be two surfaces. A point on both satisfies

$$S_a(u_a,v_a)-S_b(u_b,v_b)=0,$$

three equations in four unknowns. Where the two surfaces are not tangent, the solution set is locally a curve (the implicit function theorem: the $3\times4$ Jacobian has full rank when the two tangent planes differ).

In which direction does the curve go? It lies in both tangent planes, so its tangent is square to both normals:

$$T=\frac{n_a\times n_b}{\lVert n_a\times n_b\rVert},\qquad n_a=\partial_uS_a\times\partial_vS_a,\ \ n_b=\partial_uS_b\times\partial_vS_b.$$

Where the normals are parallel, $T$ does not exist. This is the tangent contact of Lecture 15: either the surfaces touch there, or two branches of the curve cross. A trace cannot go through such a point on its own account.

The curve’s ends are the first thing to know. Every intersection curve inside two trimmed faces ends where it leaves one of them. It leaves a face through that face’s boundary, so at a point where an edge of one solid pierces a face of the other. The boolean’s earlier stages have already found those piercing points and put a vertex at each (Lecture 21). So every open curve is known at both ends before tracing begins: find_tracing_start_points lists them. Closed loops have no ends; they are found by the different search turning_points of the last lecture.

Predictor and corrector

From a point $p_k$ on the curve, with tangent $T_k$ and a step length $\ell$, the predictor is the obvious guess:

$$\hat p=p_k+\ell\,T_k.$$

That point lies on neither surface, by about $\tfrac12\kappa\ell^2$ for a curve of curvature $\kappa$. The corrector brings it back. There are two ways to do so.

Alternating projection. Project $\hat p$ onto $S_a$, and onto $S_b$. The trouble, as the doc of predictor_corrector_step in core/geop-core-geometry/src/intersection/branch.rs explains: separate projections each enforce their own surface while ignoring the other, so they “land a full prediction-error apart and stay there no matter how good the predictor is”. Alternating between them converges slowly (linearly, at a rate set by the angle between the surfaces).

One square system. Solve for all four parameters at once. The three equations above are short of one, and the missing one is what pins the free direction. Add a fourth: the point must lie on the plane through the prediction, normal to the direction of march.

$$F(u_a,v_a,u_b,v_b)=\begin{bmatrix}S_a(u_a,v_a)-S_b(u_b,v_b)\\T_k\cdot\bigl(S_a(u_a,v_a)-\hat p\bigr)\end{bmatrix}=0.$$

The Jacobian is square:

$$J=\begin{bmatrix}\partial_uS_a&\partial_vS_a&-\partial_uS_b&-\partial_vS_b\\T_k\cdot\partial_uS_a&T_k\cdot\partial_vS_a&0&0\end{bmatrix},$$

where each $\partial S$ in the top row is a 3-vector, giving three rows. This is an ordinary Newton iteration on a square system (Lecture 14), and it lands on the curve at the plane the predictor stepped to. The corrector “may move freely onto the curve but not drift along it”, as the code says.

Note the role of the fourth row in the stopping test. The corrector is handed $(u,v)$ that already lie on both surfaces. If only the first three residuals were checked, it would report success without moving, “leaving the march pinned at its first step forever”. The fourth residual says that the point has actually advanced.

On this scene with cylinder offset $a=0.35$, radius $0.5$, and a sphere of radius 1, the demo’s algorithm gives (checked in book/scripts/demo_tests/tracing.test.mjs):

So longer strides cost corrector iterations and then fail outright. The curve’s length here is about 3.3, and it bends with a radius of the order of the cylinder’s, 0.5: a stride of that size leaves the region where Newton converges.

Worse than failing is succeeding at the wrong place. With offset $a=0.495$ and stride 0.2 the demo’s march does not fail. Both surfaces are smooth and the corrector converges at every step. But the loops pass within 0.2 of each other near $(1,0,0)$, since $a+r=0.995$ leaves them a gap of $2\sqrt{1-0.995^2}\approx0.2$. The march hops from the upper loop to the lower one, and then goes round that loop until it runs out of its 400 steps, never coming back to its start. A corrector that converges has only found a point of the intersection. It says nothing about whether it is on the same branch.

Sharp and honest, in the same loop

The corrector is the place where the sharpening rules of Lecture 5 had to be worked out, and it contains three instances.

The prediction is sharp, once. A prediction is a seed and a bias, free to be any point along the step, and it fixes the plane the corrector solves against. Likewise the plane’s normal: any plane crossing the curve pins the free direction equally well. Left wide, the normal is the one wide row of an otherwise sharp Jacobian. Where the two surfaces meet at a shallow angle, the normalized cross product of nearly parallel normals is wide by percents, enough to put zero inside a pivot of a nearly singular system. So both are sharpened, right at the top of the function, and “Nothing downstream of this line sharpens”.

Two copies of the iterate. iterate keeps a sharp seed for the next step, and the same value with its width as the answer. The sharp seed is necessary: left wide, 20 Newton steps compound width until evaluate fails outright, measured at 36 of 175 scenes. The honest answer is equally necessary, and here is the case that taught it.

predictor_corrector_step once returned sharp $(u,v)$. A traced curve that ran exactly along a face’s trim boundary therefore came back with a sharp $(u,v)$ sitting $6\times10^{-17}$ off that boundary, and face_contains classified it Inside, since could_be_equal had no width to work with. The direction check that exists to reject such a curve waved it through. The curve was spliced in as a duplicate of the boundary, carving a zero-area sliver, which surfaced much later as a face a boolean could not classify, since it had no interior point. As AGENTS.md notes, there was no geometric error in this at all. What sharpening destroyed was the statement of uncertainty every three-valued comparison depends on.

The start of each step, too. The $(u,v)$ handed in from the previous step are that step’s honest answer, and if used as the seed unsharpened, they compound through the march. The test cube_minus_box_with_wide_corners, on a tool whose corners are only known to $10^{-12}$, reached a whole patch’s width within a dozen steps that way.

Union at the end. When the corrector is done, it has two enclosures of one point, from the two surfaces. They are combined by union, the rule of Lecture 5: an average would claim more precision than either had, and hide how far apart they were.

The solve

The linear system is solved by damped normal equations:

$$(J^\top J+\lambda I)\,\delta=-J^\top F,\qquad\lambda=10^{-12}\ (\texttt{CORRECTOR\_DAMPING}).$$

The damping only keeps the step finite where $J$ loses rank (nearly parallel surfaces, a pole), at the cost of a shorter step that the next iteration continues. It cannot change what a converged answer means, since what decides convergence is the interval test on the residuals. And $J^\top J+\lambda I$ is sharpened before it is solved. That is legitimate, by the logic of Krawczyk’s preconditioner $Y$ from Lecture 14: the matrix only steers the step, and the residual $F$ alone decides where the fixed point is. The step’s honest width is the residual’s, carried through a sharp matrix. Left as intervals, the elimination would amplify the surfaces’ width by the system’s condition, which $J^\top J$ squares.

A corrector that fails to bring the two surfaces onto a common point within CORRECTOR_ITERATIONS (20) reports an error, never a best guess.

Which way, how far

The first step

From a start vertex, a curve leaves along $+T$ or $-T$. Exactly one of those (for a typical piercing point) runs into the shared interior of the two trimmed faces. To decide, trace_one_side tries a step each way, and asks face_contains on each face whether the result is Inside. Neither direction being viable is a legitimate outcome: the start points are found against untrimmed surfaces, so two cubes meeting at a corner produce a start point from which no curve really leaves.

Returning quietly carries a risk, which the source code names. A trace that should have produced a curve and didn’t looks identical to this. So it is guarded from outside: box_cylinder_blind_hole_remesh_succeeds asserts that tracing produced edges shared between both solids’ faces, which a dead trace cannot fake.

The first step has a trap of its own, which AGENTS.md calls “the invariant has one more condition”: a trace must stop at every vertex on its curve. If a known vertex lies nearer than the trial step, the full step lands past it. Where the half step toward it was rejected, because the curve runs along a boundary up to that vertex, the full step was accepted, and the trace jumped the vertex and spliced a curve partially duplicating an existing edge. So the trial now goes half way to any near vertex first, and the full step is not allowed to decide a direction past a vertex known to lie on the curve.

The stride

adaptive_step_size in lib/geop-lib-booleans/src/remesh/remesh_edges_x_faces.rs takes the smaller of two lengths.

• A sixteenth of the smaller surface (STEPS_ACROSS_PATCH). The curve lies on both surfaces, so that is the scale of features it can pass, and the stride grows and shrinks with the part. A fixed 0.1 took 400 steps along a 40 mm cut, and was also more than a cut 0.05 wide.
• The arc that takes a sixty-fourth of a full turn (STEPS_PER_REVOLUTION) around the tighter surface. A step has to be short compared to how fast the curve turns, or the straight-line predictor leaves the surface and the corrector drags it somewhere else entirely.

A surface that runs straight the way the curve does, curvature_radius_along returning None, imposes no limit. There is a good story behind this exception.

Sized only by the tightest bend, the trace misbehaves at the apex of a cone. There the tightest bend is the circle round the axis, which shrinks with the distance to the apex. A curve running along a generator into it, a straight line, took a step a fraction of its distance to the apex each time and never got there: 359 steps, spaced geometrically, for a straight line. The interpolation through those points then widened exponentially along the chain, to 111 in the control points, and the splice failed on a “corner” it could not decide.

Going further, sizing the step by the bend along the curve for every surface, was tried and is wrong: where the surfaces bend a curve a little here it may bend sharply a stride on (the blades of turbine_blade, whose pattern’s trace left its patch). So the rule is narrow: no limit from a surface that does not bend the curve at all, the tightest bend from any that does. AGENTS.md’s lesson: look at the number of steps a trace took. A line that took hundreds went wrong before it was fitted.

When to stop

A trace ends when it reaches a known vertex. This takes more care than it sounds.

Proximity is not membership. candidate_within chose the nearest vertex within one stride. But a traced curve lies on both of the surfaces it traces between, so its terminating vertex must lie on both too, and an oversized or overhanging solid puts plenty of unrelated vertices within a stride of the curve’s end. Adopting one of those ended the edge at a point on neither face. fit_pcurve then had nothing to project onto and silently returned the nearest foot point, about $8\times10^{-3}$ away, which finally surfaced as 32 validate_fast errors about pcurve endpoints not matching their 3-D points. The fix asks for what the geometry asserts: surface_could_contain on both surfaces.

AGENTS.md adds a diagnostic that paid for itself. Reporting whether the point lies on the surface at all separates “the pcurve names the wrong $(u,v)$” from “the edge was spliced onto a face it does not belong to”. The error text is identical, the fixes are opposite. Measuring it: 32 failures, every one with exactly one endpoint on-surface and one off.

Ahead of us, on this stride. A candidate must also lie between the current point and the plane the corrector lands on. One definitely behind is not it: on narrow_groove, a trace leaving the cut’s inner corner found, one stride later, the vertex 0.05 behind its start, nearer than anything ahead, and spliced a spur out of the wall to it. One definitely past the plane is the next stride’s to reach.

Before and after each stride. A vertex is looked for within the stride ball first. It matters where the curve ends at the edge of a patch the stride would have to leave. Then, once more within the chord actually taken: where the curve bends away from $T$, it meets the plane farther than $\ell$ away, and a vertex can lie ahead of the plane but outside the ball. Checking only the ball stepped over such a vertex on chained_differences_block_with_two_slots_and_a_sphere.

Never at a guess. “An intersection curve of two faces has to end at a vertex.” If a trace exhausts max_trace_steps, the result is an error naming the last point and how many steps were taken. If a stride fails on the way, the error’s context says which vertex on both surfaces lies nearest within four strides, so the next failure reports itself: AGENTS.md’s rule of adding context to the real path instead of debugging it from outside.

The same curve, twice

Every curve has two ends, both start points, so it would be traced twice. The second trace would re-derive the curve and splice a duplicate edge between the same two vertices, carving a zero-area sliver along the first.

The tempting fix samples the finished curve, and requires some sample strictly inside the face. AGENTS.md’s “Ask a topological question topologically” says this works and is wrong: “a sampling heuristic standing in for a question with an exact answer”. The right question is is this curve already represented here?, answered from the topology: an edge joining the same two vertices. In the current code (trace_one_side), the check first requires the same two end vertices, then asks one more thing, whether that edge’s curve could contain the trace’s middle point (curve_could_contain), which separates an arc from its complement sharing both ends. Ties go to the lowest id, not to whatever a hash map yields first.

What makes such a cheap test sufficient is a pipeline invariant: the edge × edge stage has already put a vertex at every coincidence, so no curve can partially overlap an edge. A curve either fully duplicates one or does not touch it. Where an invariant like that holds, reach for the exact test it enables.

Fitting the curve

The march leaves a list of points with $(u,v)$ on both surfaces. An edge needs a curve. Interpolating a cubic through the points (NurbCurve::interpolate, Piegl and Tiller’s method) passes through each exactly, and drifts from the true curve between them. For the interval kernel this drift is no detail: “otherwise a point on the true curve tests as not on its interpolant, and every containment or intersection question asked of it is answered for the wrong curve”.

So the kernel measures it, and encloses it. NurbCurve::interpolate_enclosing takes, besides the points, for each interval some true points of the curve between them, located by the same corrector (true_point_fractions says at which fractions, because a single midpoint under-measures a hump whose peak can lie elsewhere). The control points are then widened until the curve holds these. The result has the honest width of the fit, and its error is a property of the data structure, not a guess about it.

fit_branch in core/geop-core-geometry/src/intersection/branch.rs wraps this with NurbCurve::fit_enclosing: every fit has at least MIN_FITTED_LEGS legs, and while the width is more than FIT_DRIFT ($10^{-5}$) and halving the legs still halves the width, the legs are split finer. Once halving no longer halves it, the width is the points’ own, and the coarser fit is kept. A gear’s involute flank, whose curvature radius grows tenfold from root to tip, left a cubic drifting a hundred times the usual, as the code comment records; this is what handles it.

Since cubic interpolation is smooth to second order, its drift shrinks as $h^4$ on smooth data. A good test of a trace is thus to read the width: trace_one_side records its widest marched point and the widest true point between them, which tells where a wide curve got wide.

Pcurves

An edge of a B-rep has one 3-D curve, and, for each face it bounds, a curve in that face’s $(u,v)$ domain: its pcurve (Lecture 17). The face’s trim is defined by those. Answering “is this $(u,v)$ inside the face?” takes a 2-D question about pcurves, not a 3-D one about edges.

A new traced edge has both faces’ $(u,v)$ at each marched point, and a 3-D curve. Those marched values are not what the pcurve is made from. trace_one_side splices the edge into each face with splice_edge_into_face, and that calls NurbSurface::fit_pcurve (core/geop-core-geometry/src/nurb_surface/fit_pcurve.rs), which builds a pcurve from the edge’s curve alone, on any surface: sample the curve, project each sample onto the surface with Newton, seeding each projection from the previous one’s result so the walk stays continuous, and interpolate. Enough detail matters here to go through the pieces.

The first seed is found globally. Newton from anywhere settles on whichever local foot point is nearest. So the start is either pin_start, the face’s authoritative $(u,v)$, or locate. A curve that does not start on the surface is an error: there is no trace of it to fit.

The ends are pinned. Where a curve’s end is at a place the face already has a coedge for, that coedge’s own pcurve end is what the new one has to match, exactly: a boundary loop requires could_be_equal continuity from one pcurve’s end to the next one’s start. An independently projected $(u,v)$ lands a hair away. And the clamped interpolant passes exactly through its end data, so the ends are free to pin. (splice_edge_into_face goes further, pinning each end of a new loop to the arriving end’s $(u,v)$ intersected with the vertex’s projection: both enclose the same point. Using the end alone let width spread around the loop until a face split came out degenerate.)

Use the structure you know. The kernel’s arcs and helices are rational quadratics, only $C^1$ where their spans meet. fit_pcurve once fitted one $C^2$ cubic through 48 samples of the whole curve. Across those joints, its drift shrank only as $h^2$, and a helix of pitch 1 came within 15% of the $10^{-4}$ an entity may carry (MAX_WIDTH, Lecture 20). The joints are known, they are the knots of the source curve. So it now fits each smooth piece on its own and joins them.

Then adapt. A piece is sampled more densely while its pcurve is wider than a target, between MIN_PIECE_SAMPLES and MAX_PIECE_SAMPLES. A pcurve still drifting at the maximum has a trace that is not smooth where the knots say nothing, say at a surface’s seam: “it is returned as it is, wide and honest”. As AGENTS.md concludes: use the structure you know before you sample blindly, and let the sample count follow the width it produces.

The answer keeps its width. The seeds of project are sharp, the projected $(u,v)$ is not, because these values enter the fit and are later compared against the edge’s 3-D points. Narrowing them would claim a precision the projection did not have.

Recap

Exercises

1. In the demo, find (roughly) the longest stride that closes at offset $a=0.35$, and at $a=0.45$. Relate them to the curve’s tightest bend. Then read off the kernel’s stride for a sphere of radius 1: what would STEPS_PER_REVOLUTION give?

2. Compute the prediction’s distance from the sphere for stride $\ell$ on a great circle of a unit sphere, starting on it. How does the corrector’s first residual scale with $\ell$? Why does this suggest the corrector iteration counts in the table?

3. Why does predictor_corrector_step plant its plane at $\hat p$, and not at $p_k$? Why is a plane better than a sphere of radius $\ell$ around $p_k$?

4. A trace ran from vertex $v$ to a vertex $w$, and a second trace from $w$ found a curve to $v$. Under what circumstance are they two different branches, and how does trace_one_side tell?

5. The kernel’s pcurves are fitted to the edge’s 3-D curve, even though the march has $(u,v)$ for both surfaces at each point. Argue for this choice. (Hint: what is compared with what, later?)

6. Let a cubic interpolate samples of a curve at spacing $h$. Its drift on a smooth stretch is $O(h^4)$, and across a jump in curvature $O(h^2)$. Estimate how many times fewer samples are needed if you know where the jump is, for a target drift $10^{-5}$ and a curve of unit size.

Lecture 17: The boundary representation

Part III gave us curves and surfaces we can intersect with proof. But a curve is not a part. A part has an inside, and it has faces you can pick, edges you can round and a volume you can weigh. Something has to say how pieces of surface fit together into a closed skin, and that something is the subject of Part IV. It is much less mathematical than what came before and every bit as easy to get wrong.

What should a solid be, in a computer?

A solid is a set of points in space. A computer cannot store a set of infinitely many points, so it stores a description of it. There are four serious candidates.

Voxels divide space into small cubes and mark each as in or out. Booleans are trivial, one logical operation per voxel, and that is the end of the good news. A face is never flat, a dimension is never exact, and memory grows with the cube of the resolution.

Triangle meshes store the skin as flat triangles. Every graphics system uses them, and Geop’s rasterizer (core/geop-core-rasterize) produces them for drawing. But a mesh of a cylinder is not a cylinder. It has no radius you can read and no tangent you can trust, and a fillet between two faces needs to know what those faces are.

Constructive solid geometry (CSG) stores a tree: this box, minus that cylinder, union this sphere. The description is exact and compact. But the faces and edges the result has are nowhere in it, so you cannot pick “that edge” and round it without first computing what the tree means.

Boundary representation (B-rep) stores the skin itself: exact pieces of surface, each trimmed by exact curves, and a record of which piece meets which. A point is inside if it is enclosed by the skin. This is what every mechanical CAD system of consequence uses, because it answers the questions a designer asks: there is an entity for the top face and an entity for the edge around the hole, each with exact geometry.

The price is that a B-rep has two kinds of content that must agree. Geometry says where things are: surfaces, curves, points. Topology says what is connected to what: this edge bounds these two faces, this vertex ends these three edges. Most of what makes a boolean hard is keeping the two in agreement while both are being changed. This lecture introduces the data structure. Lecture 18 changes it. Lecture 20 checks the agreement.

Geop is also a little CSG-like at a higher level: a part is a program of steps (Lecture 2), and each step’s output is a B-rep, not a tree.

The entities, from the bottom up

Take a cube. It has 8 corners, 12 straight sides and 6 square faces. The B-rep has entities for those, and needs a few more besides.

A vertex is a point in space. In Geop, Vertex (core/geop-core-topology/src/vertex.rs) holds exactly one field, point: Vector3, and the three coordinates are interval Scalars, as everywhere. So a vertex is an enclosure of a point (see Lecture 4), and two vertices are the same point only if they are the same entity.

An edge is an arc of a curve between two vertices. Edge (core/geop-core-topology/src/edge.rs) holds curve: Curve3, start_vertex and end_vertex. Note a simplification that many kernels do not make: there is no window [t0, t1] into a longer curve. The curve is trimmed so that its own domain is exactly the edge, start_vertex at the start of the domain and end_vertex at the end. Splitting an edge splits its curve.

A face is a piece of a surface. Face (core/geop-core-topology/src/face.rs) holds a surface: NurbSurface3D, one outer boundary, a list of holes, and the shell it belongs to. It means: everything of the surface inside the outer boundary and outside every hole.

How is a boundary given? You might think: by a list of edges. That is where the naive design fails, and understanding why is the key to the whole structure.

Why an edge is not enough: the coedge

Consider the edge between the top and front faces of the cube. It belongs to two faces. Walking the boundary of the top face counter-clockwise (seen from outside), you pass along it in one direction. Walking the boundary of the front face, also counter-clockwise, you pass along it in the other direction. Moreover, each of the two faces has its own coordinates: the same 3-D curve appears as one curve in the (u, v) plane of the top face and as a different curve in the (u, v) plane of the front face.

So there are facts that belong neither to the edge, of which there is one, nor to the face, which has many sides. They belong to one side of the edge, on one face. This is what a coedge is. (Other texts call it a half-edge, or a fin, or a use of an edge.) Coedge (core/geop-core-topology/src/coedge.rs) has:

The pcurve deserves a moment. We met pcurves in Lecture 16. A curve in parameter space is the only way to say exactly which part of a surface a face is. Its existence means every edge of a solid is represented three times: once in 3-D, and once in each of its two faces’ parameter planes. They must agree, and they must agree to within their widths. That is the most-checked fact in the kernel.

Loops

A face’s boundary is one or more closed cycles of coedges. Such a cycle is a loop. In Geop a loop is not an entity at all. A face’s boundary is a BoundaryType (core/geop-core-topology/src/boundary.rs), and BoundaryType::Loop(id) names one coedge of the cycle, the anchor. Following next from the anchor walks the loop until it returns. Model::iterate_loop_coedges does this. There is nothing else to keep consistent.

(BoundaryType has a second variant, Vertex(id), for a boundary that is a bare point, a face with no edges yet or another body touching a face at a point. Then there is no pcurve, only the vertex. A face always has an outer boundary, so “no boundary” is not a state the type can even express. We return to this.)

Now use the demo. It is a real half-edge structure, built from a list of faces the way a BodySpec describes a body. Each arrow is a coedge, drawn shifted to its left, into the material of the face it bounds.

Notice what you never saw: a pointer from a coedge to its opposite. Geop does not store one. The opposite coedge is found through the edge, by Model::coedges_of_edge, because a stored pointer would be one more thing that could disagree with the rest.

Why a doubly linked list rather than an array in each face? Because the operations of the next lecture are all of the form “cut this ring here and join it to that ring there”. On linked lists these are four pointer assignments. On arrays, they are reorderings and index fix-ups everywhere.

Shell, solid, sheet, wire, body

A shell is a connected set of faces. Shell (core/geop-core-topology/src/shell.rs) holds faces and solid: Option<SolidId>.

A solid is a region of space bounded by closed shells. Solid (core/geop-core-topology/src/solid.rs) holds shells. It usually has one. It has more if it consists of separate lumps, or has a void inside. What is interesting is how lumps and voids are told apart: they are not. The doc comment says it directly: “a point is inside the solid when it is inside an odd number of them, so a shell is a lump of its own or a void in one by where it lies, not by its place in shells”. A void’s faces point into it, out of the material, and the volume it encloses comes out negative.

A sheet is a shell that belongs to no solid: solid is None. Its faces stand on their own, bounding nothing, and its border edges are used by one face only. A sheet is what a surface operation builds, and what a split cuts a solid with.

A wire (core/geop-core-topology/src/wire.rs) holds edges and vertices that bound no face: curves in space. A sweep’s path can be one.

Body (core/geop-core-topology/src/body.rs) is a plain enum of Solid(SolidId) and Sheet(ShellId): what operations act on. Remesh and imprint (Lecture 21) take bodies because all they ask of one is its faces, edges and vertices.

The model is the arena that holds it all:

Solid ──▶ Shell(s) ──▶ Face(s) [+Surface] ──▶ Coedge(s) [+Curve2] ──▶ Edge [+Curve3] ──▶ Vertex
Wire ──▶ Edge(s) [+Curve3], Vertex(es)

Orientation has no flag

Which side of a face is outside? In many kernels, a face carries a boolean, “the normal is reversed”. Geop has no such flag. The material side follows from the surface alone: the normal is $S_u\times S_v$, and the convention is that it points out of the solid. A face that has to be turned around has its surface turned around.

The conventions, all of which are checked (core/geop-core-topology/src/validation/face_orientation.rs, and the third by validation/manifold_check.rs):

1. The normal $S_u\times S_v$ points out of the material.
2. In the (u, v) plane, an outer loop runs counter-clockwise and a hole runs clockwise. Equivalently, the material lies to the left of every coedge, for outer loops and holes alike.
3. Two coedges of an edge have opposite senses.

Are these independent? Check with the top of a cube. If $(u, v)$ are $(x, y)$, the normal is $+z$, up and out. An outer loop counter-clockwise in $(u, v)$ is counter-clockwise seen from above, with the top of the square to the left. Fine. Now the neighbouring face must traverse the shared edge the other way, for it too to have its material on its left when seen from its own outside: fold the two faces flat along the edge to see it. So (3) follows from (1) and (2) on a closed orientable surface.

Model::reverse_face (core/geop-core-topology/src/edit/reverse_face.rs) shows what turning a face around amounts to. It mirrors the surface in u (NurbSurface::reverse_u), applies the same mirror to every pcurve, and reverses each loop: swap next and prev, run each pcurve the other way, flip each sense. The edges, the vertices and the coedges on the other side of the edges are not touched, which is what lets a boolean flip one operand’s faces while they stay glued to the other’s.

Make the invariant a type

Face once had a flat Vec of boundaries with “index 0 is the outer loop” as an unwritten rule. It now has outer: BoundaryType and holes: Vec<BoundaryType>, and lookups return a BoundaryIndex, either Outer or Hole(i). AGENTS.md calls this “not cosmetic”, and this is why.

An edge joining two points of a face’s boundary means different things according to which boundaries the points lie on. Two different holes merge. A hole bridged to the outer loop stops being a boundary. Two points of the outer loop cut the face in two. With a flat list, these looked like “same index or different index”: two cases, handled as two cases, which silently produced faces carrying several disjoint patches. As a type, match (loop_a, loop_b) over BoundaryIndex is four cases, and the compiler insists on all of them. We work through them in the next lecture.

The general lesson goes beyond this one type. Any time a comment says “by convention, the first element is special” or “this is only set if”, ask whether the type could make the wrong state unrepresentable. Each such state a type rules out is a state nothing else ever has to handle.

Counting: Euler and Poincaré

There is a global check on all these counts. For a polyhedron with $V$ vertices, $E$ edges and $F$ faces, topologically a sphere,

$$V-E+F=2.$$

For a cube: $8-12+6=2$. The reason is worth having in your fingers. Remove one face and flatten the rest onto the plane. Now repeatedly remove things: an edge between two faces (removes one edge and one face, leaves $V-E+F$ unchanged), or a dangling edge with its end vertex (removes one vertex and one edge, unchanged). Eventually one vertex remains, for which $V-E+F=1$. Add back the face we removed first, for 2.

Two generalizations cover everything a B-rep can be. A handle, like the hole through a ring, lowers the count by 2 per handle. A face with a hole loop is not simply connected, so the argument above needs one extra edge to bridge the hole to the rest of the boundary, which adds one to $E$ without adding a face. So for $S$ shells, total genus $G$, and $L$ hole loops over all faces,

$$V-E+F-L=2(S-G).$$

Choose the block with a through hole in the demo. It reads $16-24+10-2=0$, so one shell of genus 1, as a square ring has.

A worked example from the real kernel

The example box_with_drill_hole is a 2 by 2 by 1 box with a blind hole of radius 0.4 and depth 0.5 in its top. Running it (geop 0.12.0, list_elements over geop mcp, capped) gives 16 vertices, 24 edges and 11 faces:

The circle is four quarter arcs because a sketched circle is four rational arcs, and a sweep turns each into a wall of its own (see Lecture 8 on why a quarter arc is a good piece, and Lecture 26 on the sweep). Now $V-E+F=16-24+11=3$. The top face of the box has the rim as a hole loop, so $L=1$ and $V-E+F-L=2$: one shell, genus 0. A blind hole is no handle. Correct.

Compare cross_drilled_shaft, a stepped shaft with a hole straight through it from side to side (the same run): 26 vertices, 50 edges, 24 faces, so $V-E+F=0$. A through hole makes genus 1, and $2(1-1)=0$, so $L=0$: no face there has a hole loop. That is true because the hole’s rim runs across the shaft’s quarter seams, so each bite is part of the outer loop of the quarter it is in. The formula told you something about how the kernel cut the faces. We return to this example in Lecture 21.

Saying a body whole

How do you make a B-rep? You could insert entities one at a time, keeping a partly built body valid after each step. The operations of the next lecture do that, because they are editing what exists. But a sweep comes in many shapes (one span or several, closed around an axis or not, with poles, holes or caps), and each would need its own sequence of steps. It is far simpler said whole: these vertices, these edges between them, these faces, each bounded by these coedges.

That is BodySpec (core/geop-core-topology/src/build.rs). It is plain data, entities referring to each other by index in lists: vertices, edges (EdgeSpec: a curve, a start, an end), faces (FaceSpec: a surface, an outer loop, holes), shells, and a flag solid. A loop is a list of CoedgeSpec, each saying what it runs along (CoedgeOn::Edge(index, sense) or CoedgeOn::Vertex(index)) and giving its pcurve.

Model::build_body takes a spec, runs BodySpec::check, and only then inserts anything. The check is entirely combinatorial. Every edge runs between vertices that exist. Every coedge ends at the vertex where the next one starts, so that every loop closes. Every face is in exactly one shell. Every vertex is used. And the key test: an edge is used an even number of times, as often one way as the other, with exactly one use allowed only for a sheet; used twice in the same sense, or by one face in a solid, it is an error that names the edge.

This is a worthwhile design to imitate. A construction builds a description, the description is checked before the model is touched, and a failed check leaves nothing half done. Model::body_spec goes the other way, describing faces of a model as a copy, and says where each entity came from (SpecSources). This pair is what geop_lib_sweep’s builders and the operations that copy or rebuild faces are written with. A description with no faces is a wire.

Exercises

1. A triangular prism has how many vertices, edges, faces and coedges? Verify $V-E+F=2$, then choose the prism in the demo to check, and walk one of its triangular loops.

2. Geop’s sphere is eight faces meeting at two poles (lib/geop-lib-sweep/src/shapes/sphere.rs). Count $V$, $E$, $F$, the real coedges and the degenerate ones. Check the Euler formula. Why does a degenerate coedge leave the count alone?

3. A disc with two round holes is extruded. Each circle is a single closed edge from a vertex to itself, ignoring how Geop would cut it into quarters. Count $V,E,F,L$ and find $G$.

4. Show that a coedge’s pcurve, its edge’s curve and the face’s surface being all exact does not by itself prevent $S(p(t))\ne C(t)$. What would have to hold, and what is the kernel’s stance on “hold” given that everything is an enclosure?

5. Rules (1) and (2) of Orientation has no flag imply (3). State this as a precise claim about two faces sharing a straight edge in the plane $z=0$, one to each side, with both surfaces parametrized so that the normal is $+z$.

6. Why does Model not store an opposite pointer on Coedge, and what is the cost of that choice? How would you decide whether to add one?

Lecture 18: Poles, seams and splicing

Lecture 17 described a B-rep at rest. Now we change one. Everything a boolean does to topology comes down to two edits: split an edge at a vertex, and splice an edge into a face. Before that, we have to deal with the two places where a surface’s coordinates stop being a faithful picture of the surface.

The trouble with coordinates

A pcurve lives in the (u, v) plane. Face containment (Lecture 19) counts crossings of a ray with the pcurves in that plane, and so relies on a fact that sounds too obvious to state: a face’s loops are closed polygons in (u, v). Two kinds of surface break the picture.

A closed surface has a seam. Take a whole cylinder as one patch, $u\in[0,1)$ round the axis. The lines $u=0$ and $u=1$ are one line in space. A loop around the wall is closed in space and is no closed curve in (u, v). It runs from $(0,v)$ to $(1,v)$. The classical remedy is a seam edge, one edge that a single face uses twice, once on each side of the parameter domain. Then the outer loop goes along the bottom, up the seam, back along the top, and down the seam. It works, and it infects everything: an edge whose two coedges are on one face, a point on the seam with two (u, v), a surface evaluated past its domain’s end by anyone who walks across. And a NURBS patch cannot close on itself smoothly and stay a single simple tensor-product patch unless it is periodic.

A pole collapses a row. Take a sphere with $v$ running pole to pole. The whole row $v=0$, every $u$, is one point in space. Here the surface is perfectly smooth and only the coordinates are singular. As AGENTS.md puts it, the surface is fine, the coordinates are not.

Geop’s way with seams: there are none

Every surface in Geop is a plain patch over a rectangle, and nothing closes on itself. A full turn of a revolve is cut into four spans, each a face of its own, so a cylinder is four quarters; lib/geop-lib-sweep/src/revolve.rs does this, and the wall of the drilled box in Lecture 17 was four faces for this reason. A “seam” in Geop is an ordinary edge between two ordinary faces.

STEP files, from systems which do use seams, have periodic faces with one edge used twice. The importer takes them out, cutting such a face along meridians into sectors, each on a patch of its own (ops/geop-ops-step/src/import/body.rs; see Lecture 35). So the rest of the kernel does not need a seam path at all. It is the same bargain as Everything is a NURBS: one more piece of geometry, and one whole class of branch gone.

A pole needs a coedge with no edge

A pole cannot be cut away: a cone has an apex and a sphere has two poles, whichever way you cut it. So look at one octant of Geop’s sphere (lib/geop-lib-sweep/src/shapes/sphere.rs). It is a face over the unit square whose row at the pole has collapsed. Three of its sides are real arcs. The fourth is a single point in space.

Now consider the loop. In space it goes round three arcs: down a meridian, along the equator, back up. It is closed. But in (u, v), one meridian ends at $(0,1)$, say, and the other starts at $(1,1)$. Both corners are the pole. If the loop jumped from one to the other, a ray in (u, v) could pass through the gap uncounted, and the parity of the crossings would be wrong for points near it. It is exactly the kind of bug no single-coedge check can see, and core/geop-core-topology/src/validation/pcurve_loop_continuity.rs exists to catch it.

So the loop needs a piece of boundary that runs along the row in (u, v) while staying at one point in space. That is CoedgeGeometry::Vertex(id): “a degenerate loop segment that sits at a single, already-existing vertex the whole way”. It has a pcurve, along the row. It has no edge, and so no opposite: there is nothing to share between two faces. It adds no vertex and no edge, so it leaves the Euler count alone. In a BodySpec, it is CoedgeOn::Vertex.

Picture the octant’s outer loop, with its pole at $v=1$:

(0,1) ◀── Vertex coedge along the row ── (1,1)     all of this row is the pole
  │                                        ▲
  ▼  meridian 1                  meridian 2│
(0,0) ───────── equator arc ─────────▶ (1,0)

Four coedges, closed in (u, v); three edges, closed in space.

One vertex at many (u, v)

A pole then behaves unlike any other vertex: a loop passes it at several places in parameter space. Which corner an arriving edge belongs to is told by its (u, v), not by the direction it leaves the vertex in. This caused a real bug: pinning an arriving edge to the first coedge that happened to end at the vertex bends the edge to a corner at the row’s end, wherever it actually arrives.

The fix (AGENTS.md, A pole is one vertex at many (u, v)) is in splice_edge_into_face. It first fits the new edge’s pcurve leaving its end at a pole unpinned. Then split_pole_row looks at where the pcurve ended. If that is strictly inside the row’s coedge, it cuts the row in two there, at a sharp $u$ (where to cut is a free choice), and the new edge is spliced between the two halves. So the row gets a corner exactly where the edge arrives.

Degenerate parametrizations have other consequences, for projection and for Newton, which we saw in Lecture 14: a seed on a pole never leaves it unless you step along a meridian. Together these give the general lesson. Whenever an algorithm reads a coordinate near a pole, ask whether what it wants is a fact about the coordinates or about the surface.

Edit one: splitting an edge at a vertex

Suppose an edge passes through a point where a vertex is wanted. Model::split_edge_at_vertex (core/geop-core-topology/src/edit/split_edge_at_vertex.rs) makes two edges joined at that vertex. The existing edge keeps the first piece, and a new edge holds the second, “in the same direction”. Every coedge along it is split to match, and the two pieces are linked into the loop in place of the one.

It sounds simple. The code’s comments say why it is not. Three parameters are involved: $t$ on the 3-D curve, and one on each pcurve.

Where to cut is a free choice, and a good one is a sharp one. The search that found the crossing returned an enclosure of the parameter. The cut goes at the sharp middle of its Newton refinement, NurbCurve::split_parameter_at_point (core/geop-core-geometry/src/nurb_curve/refine.rs), which then checks that the curve reaches the vertex at that parameter. AGENTS.md’s section Width along a split edge gives the numbers behind this. Cutting instead at the whole refined enclosure carried the parameter’s width through Boehm’s formula ($\alpha=(t-a)/(b-a)$) into the pieces’ control points and knots. The next cut’s parameter, found on the wider piece, came out wider. The width grew threefold per cut: a quarter arc cut forty times went from $3\times10^{-15}$ to $0.08$. Cut at a sharp point, each cut adds one subdivision’s rounding.

Validate the value you are about to use. The check is done at the sharpened parameter, not the wide one. A wide box genuinely contains the vertex, but “is the vertex somewhere along this stretch?” and “does cutting exactly here give that vertex?” are different questions.

Cut the pcurves at that same point. The point where the pieces will meet, $C(t)$, is taken first. Each coedge’s face then locates that point on its surface (NurbSurface::locate), not the vertex: “Cut each at a point of the vertex of its own and the pieces of the edge and of its pcurves end apart, by as much as the vertex is wide”. In the case that exposed this, 1.4e-13 apart, where a cubic crossed a plane at its inflection.

Refine the target, then the parameter. The (u, v) has to be pinned down first. A pcurve parameter refined against a target only as sharp as the containment search’s tolerance gives a wide residual, a wide step, and a result no narrower than it started, silently. That is why it is surface.locate here: located, not merely isolated.

Everything is checked before anything is changed, so a pcurve that cannot be split does not leave a half-updated model.

Do not reconcile a mismatch by moving geometry. AGENTS.md records pinning the split curve’s endpoint onto the vertex: it looked like it worked for the 3-D curve, and applied to the pcurve it pinned a control point to a box so wide the pcurve lost its tangent, taking a sweep of test scenes from 0/175 failing to 12/175. Fix the parameter, not the geometry.

Edit two: splicing an edge into a face

Model::splice_edge_into_face (core/geop-core-topology/src/edit/splice_edge_into_face.rs) takes an edge that already exists, usually because two faces are meant to share it, and makes it part of one face’s boundary. It adds two new coedges on that face, fwd and rev, one each way along the edge. (So each face gets a pair. A finished manifold edge has one coedge on each of two faces; what a splice imprints into a face is a cut with a side for each piece of it. When the face is split, the two end up on two different faces. In a boolean, faces on the other side of the cut are dropped.)

What happens depends only on which of the edge’s end vertices are already on the face’s boundary and, if both are, on which boundaries.

The last four all begin with the same rewiring. Let a be a coedge arriving at the edge’s start, b a coedge arriving at its end. Then four assignments:

let a_next = a.next;
let b_next = b.next;
link(a, fwd);      // arrive at the start, run the edge forward
link(fwd, b_next); // leave the end along what used to follow b
link(b, rev);      // arrive at the end, run the edge back
link(rev, a_next); // leave the start along what used to follow a

These are the actual lines, abbreviated. Think through what they do to a ring. If a and b are on one ring $R$, write $R=a\,X\,b\,Y$, with $X$ the coedges from a_next to b, and $Y$ those from b_next round to a. After the rewiring, $a\to\text{fwd}\to Y\to a$ is one ring, and $X\to b\to\text{rev}\to X$ is another. One ring became two. If a and b are on two different rings, the same four assignments give one ring: a walk goes a, fwd, all round the second ring, b, rev, then round the first back to a. Two rings became one. All of the case analysis follows from this simple fact: the same four pointer moves either divide a ring or join two.

Try every case. In the demo, click two vertices of the boundary, or click in the middle of the face for a free point (it snaps to a half-unit grid). It runs the same case analysis and the same rewiring on a small real half-edge structure, and reports what it found. The total area of the faces is shown too. It never changes: a splice rewires boundaries, and moves no material.

The six cases

Neither end on the face. The new edge floats in the interior of the face. fwd and rev link to each other to make a ring of two coedges, and that ring is pushed onto holes. It encloses no area. It happens when an edge lying in the face is imprinted before anything else has reached it: the next edge imprinted beside it attaches to this ring, as we will see in the next lectures.

One end on the face. splice_spur inserts fwd and rev after the arriving coedge a: out and straight back. The ring gets longer by two and its area does not change. face.holes and the set of faces are unaffected. The code insists that the free end “has to lie strictly inside” the face, with a test by face_contains, and returns an error rather than a loop that runs out of the face and back. That check exists because of the test scene narrow_groove, whose failure surfaced far downstream in classification.

Both on two different holes. Rings join. Set the boundary at hole $i$ to the new ring and remove hole $j$.

A hole and the outer loop. Rings join. The hole “stops being a separate boundary and its coedges become part of the one ring that now bounds the face”. The face has a slit from its outer boundary to what was the hole. It is a valid face, with a loop that runs along that slit twice.

Both on the outer loop. One ring becomes two, both bounding material. This is the case the boolean is built on. split_face keeps one ring in the existing face and makes a new face from the other, on the same surface, in the same shell (new_face_from_ring). Then each of the old face’s holes must go with whichever piece it now lies in: take a point on the hole, and ask face_contains. “Inside” moves it to the new face, anything else stays. A hole with no usable point stays put “rather than being dropped: leaving it on the wrong face is a recoverable error, losing it silently fills in a hole that should exist”.

Before doing that, it computes the signed area of both rings. “A half of zero area means the new edge runs along the stretch of boundary between its own endpoints.” It returns the error DEGENERATE_SPLIT rather than building a sliver, because a zero-area face is structurally perfect and every other check accepts it. (We see the bug that led here in Lecture 21.)

Both on one hole. One ring becomes two, but which is which? Let the hole be the region in the middle, and let the edge run from one point of its boundary through material to another. Together with one arc of the hole, the edge encloses a patch of material, now cut off. Together with the other arc, it goes around everything: the hole, plus that patch. The first ring bounds material, and is a new face. The second still has material outside it, and stays a hole. To tell them apart: “both rings pass through the same vertices and contain the same edge, so no purely topological test can tell them apart”. Only the winding can: the ring of positive signed area bounds material. If the area cannot be decided, the code reports it, rather than guessing.

(AGENTS.md’s short summary says two points of one hole “divide it into two holes”, and says only the outer split creates a face; and the doc table above splice_edge_into_face still says “divided into two holes”. The code is what is described here, and its own comment agrees: dividing a hole is a case that creates a face.)

Which corner?

A vertex can appear more than once on a boundary. The root of a spur is passed twice, once on each side of the spur; so is a point where a face touches itself. Then there are several arriving coedges, each with its own corner, and the new edge belongs in the one it leaves into. choose_corner decides it from the rule that the material fills “the angle swept counter-clockwise from the way the boundary leaves to the way it arrived from”. The test scene three_inscribed_bores_remesh_imprints_every_arc shows the failure. An arc spliced as a spur from a point of a bore’s rim, and the rim arc spliced at that same point a moment later, landed the spur in the wrong piece of the wall. Since the corner choice can move the pin, the new pcurve is then refitted.

Pins

The new pcurve has to start exactly where the arriving coedge ends in (u, v), for the loop to stay closed there. pin takes the arriving end, intersected with the vertex’s projection onto the surface. Using the end alone let interpolate_enclosing’s width spread around the loop until a split came out degenerate. Both values enclose the one (u, v) the vertex is at, so their intersection is honest and tighter. AGENTS.md uses this as its example of combine two enclosures of the same value: union when you must cover both possibilities, intersect when both constraints must hold, average never (Lecture 5).

Splitting a face invalidates “on this surface”

There is a subtle consequence of the face-split case. Both halves share one surface, so every edge already imprinted into the original is still “on the surface” of the half it does not belong to. Before faces had proper topology, an imprinting loop asked curve_surface_intersect whether an edge was on a face’s surface, found the edge on the other half too, imprinted it there, split that half, and so on, forever.

Faces with real outer loops and holes can answer the right question: is this edge inside this face’s trim? That is face_contains, in find_coincident_pair of lib/geop-lib-booleans/src/remesh/remesh_edges_x_faces.rs. AGENTS.md’s Proximity is not membership generalizes this: being near, or on an untrimmed surface, is not belonging to the face. The surface of a face runs on past its trim.

Exercises

1. Show that every case of the splice table preserves $V-E+F-L$.

2. In the demo, merge the two holes by connecting vertices (2,1) and (4,1). Then connect (2,3) and (4,3). What case is the second, what are the rings, and what is the area of the new face?

3. Write out, in terms of $R=a\,X\,b\,Y$, what happens when a and b are the same coedge: an edge from a vertex back to itself.

4. new_face_from_ring decides which of the old face’s holes move to the new face by one point of each hole. Why is a single point enough?

5. Why does split_pole_row take a sharp $u$ at which to cut the row, when the end of the pcurve is an enclosure? Is this the “free choice” of AGENTS.md’s rule on sharpening?