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?