RasmahRasmah

NURBS: Curves, Surfaces, and Curve Fitting

using Rasmah, LinearAlgebra

Overview

A NURBS (non-uniform rational B-spline) is the industry-standard exact curve and surface representation behind STEP and every parametric CAD kernel. Rasmah’s NurbsCurve / NurbsSurface implement the Level-3a core — de Boor evaluation, degree elevation, knot insertion, and circle/line conversion — plus the Level-3b least-squares curve fitting used to turn a traced surface–surface intersection back into an exact curve.

A degree-p NURBS curve is C(t) = Σ N_{i,p}(t) w_i P_i / Σ N_{i,p}(t) w_i, where N_{i,p} are the Cox–de Boor basis functions over a knot vector, P_i are control points, and w_i are weights. Rational weights are what let a NURBS represent a circle exactly — a plain polynomial B-spline can only approximate one.

Building a curve by hand

NurbsCurve(degree, knots, control_points, weights) stores everything directly, and curve_point(c, t) evaluates it. A straight line is a degree-1 curve with the endpoint weights 1:

line = NurbsCurve(1, [0.0, 0.0, 1.0, 1.0], [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]], [1.0, 1.0])

println("line(0.0) = ", curve_point(line, 0.0))
println("line(0.5) = ", curve_point(line, 0.5))
println("line(1.0) = ", curve_point(line, 1.0))
line(0.0) = [0.0, 0.0, 0.0]
line(0.5) = [0.5, 0.0, 0.0]
line(1.0) = [1.0, 0.0, 0.0]

A rational quadratic (degree 2) with the middle weight cos(π/4) ≈ 0.707 sweeps out a quarter circle exactly:

w = cos(π / 4)
arc = NurbsCurve(
    2,
    [0.0, 0.0, 0.0, 1.0, 1.0, 1.0],
    [[1.0, 0.0, 0.0], [1.0, 1.0, 0.0], [0.0, 1.0, 0.0]],
    [1.0, w, 1.0],
)

println("arc(0.0) = ", round.(curve_point(arc, 0.0); digits = 4))
println("arc(0.5) = ", round.(curve_point(arc, 0.5); digits = 4))
println("arc(1.0) = ", round.(curve_point(arc, 1.0); digits = 4))
arc(0.0) = [1.0, 0.0, 0.0]
arc(0.5) = [0.7071, 0.7071, 0.0]
arc(1.0) = [0.0, 1.0, 0.0]

curve_point(arc, 0.5) is [√2/2, √2/2, 0] — exactly on the unit circle, which a degree-2 polynomial Bézier could never hit.

Converting curves to NURBS

to_nurbs converts any supported analytic curve to its exact NURBS form:

ln = to_nurbs(LineCurve([0.0, 0.0, 0.0], [1.0, 2.0, 0.0]))
println("line → degree ", ln.degree, ", ", length(ln.control_points), " control points")

quarter = to_nurbs(
    CircleArc([0.0, 0.0, 0.0], [0.0, 0.0, 1.0], 1.0, 0.0, π / 2, [1.0, 0.0, 0.0]),
)
println("quarter arc → degree ", quarter.degree, ", ", length(quarter.control_points), " control points")
println("  weights: ", round.(quarter.weights; digits = 4))
line → degree 1, 2 control points
quarter arc → degree 2, 3 control points
  weights: [1.0, 0.7071, 1.0]

A full circle becomes four 90° rational arcs stitched together — 9 control points, uniform interior knots:

circ = to_nurbs(CircleCurve(1.0, 0.0))
println("circle → degree ", circ.degree, ", ", length(circ.control_points), " control points")
println("  t = 0.0: ", round.(curve_point(circ, 0.0); digits = 4))
println("  t = 0.5: ", round.(curve_point(circ, 0.5); digits = 4))
println("  t = 1.0: ", round.(curve_point(circ, 1.0); digits = 4))   # closes back on itself
circle → degree 2, 9 control points
  t = 0.0: [1.0, 0.0, 0.0]
  t = 0.5: [-1.0, 0.0, 0.0]
  t = 1.0: [1.0, 0.0, 0.0]

Knot insertion

knot_insert(c, u) inserts a knot at parameter u (Boehm’s algorithm). The shape is unchanged but the control net gains a point — the mechanism used to split or trim a curve at an arbitrary parameter:

c2 = knot_insert(circ, 0.25)
println("after one insertion:  ", length(c2.control_points), " control points")

c3 = knot_insert(circ, 0.25; times = 2)
println("after two insertions: ", length(c3.control_points), " control points")

println("shape preserved: ", curve_point(circ, 0.5) ≈ curve_point(c2, 0.5) ≈ curve_point(c3, 0.5))
after one insertion:  10 control points
after two insertions: 11 control points
shape preserved: true

Degree elevation

degree_elevate(c) raises the degree (Prautzsch’s algorithm, The NURBS Book A5.9), again without changing the shape:

c4 = degree_elevate(circ; times = 2)
println("degree ", circ.degree, " → ", c4.degree, ", ", length(circ.control_points), " → ",
    length(c4.control_points), " control points")
println("shape preserved: ", curve_point(circ, 0.125) ≈ curve_point(c4, 0.125))
degree 2 → 4, 9 → 17 control points
shape preserved: true

NURBS surfaces

NurbsSurface is the tensor-product analogue: two knot vectors, a grid of control points, and a grid of weights, evaluated by surface_point(s, u, v):

knots = [0.0, 0.0, 1.0, 1.0]
cps = [[u, v, u * v] for u in [0.0, 1.0], v in [0.0, 1.0]]
surf = NurbsSurface(1, 1, knots, knots, cps, ones(2, 2))

println("surface(0, 0)     = ", surface_point(surf, 0.0, 0.0))
println("surface(0.5, 0.5) = ", surface_point(surf, 0.5, 0.5))
println("surface(1, 1)     = ", surface_point(surf, 1.0, 1.0))
surface(0, 0)     = [0.0, 0.0, 0.0]
surface(0.5, 0.5) = [0.5, 0.5, 0.25]
surface(1, 1)     = [1.0, 1.0, 1.0]

This bilinear patch maps the unit square to the ruled surface z = u·v.

Least-squares curve fitting

fit_nurbs(points, degree) fits a non-rational NURBS to a polyline by global least squares — chord-length parameters, an averaging knot vector, and endpoint interpolation. It is the Level-3b step that turns a sampled curve (a traced intersection, a digitised profile) back into a compact exact curve:

helix(t) = [cos(2π * t), sin(2π * t), t]
samples = [helix(t) for t in range(0, 1; length = 31)]

fit = fit_nurbs(samples, 3; nctrl = 8)
println("fit: degree ", fit.degree, ", ", length(fit.control_points), " control points")

err = maximum(norm(curve_point(fit, t) - helix(t)) for t in range(0, 1; length = 201))
println("max error over the helix: ", round(err; digits = 4))

# the first and last control points are pinned to the data endpoints
println("endpoints interpolated: ", curve_point(fit, 0.0) ≈ samples[1] && curve_point(fit, 1.0) ≈ samples[end])
fit: degree 3, 8 control points
max error over the helix: 0.0063
endpoints interpolated: true

Raising nctrl trades smoothness for fidelity — nctrl = 12 drops the max error to ~1e-4, and nctrl = length(samples) reproduces the samples to machine precision.

A real use: reverse-engineering an airfoil

Fitting is the bridge from measured to exact. A digitized or scanned profile arrives as a cloud of points; fit_nurbs compresses it into a curve that any downstream feature can consume. The classic profile is the NACA 0012 airfoil — a 12%-thick symmetric section defined by the thickness formula y_t = 5t\,(0.2969\sqrt{x} - 0.1260x - 0.3516x^2 + 0.2843x^3 - 0.1015x^4). Sample it as if it had been digitized, then fit:

t = 0.12; c = 0.200                      # 12% thickness, 200 mm chord
yt(x) = 5 * t * c * (0.2969 * sqrt(x) - 0.1260 * x - 0.3516 * x^2 + 0.2843 * x^3 - 0.1015 * x^4)

xs = collect(range(0.0, 1.0; length = 61))
top = [[x * c, yt(x), 0.0] for x in xs]
bot = [[x * c, -yt(x), 0.0] for x in reverse(xs[1:end-1])]
samples = vcat(top, bot)                # 121 "digitized" points

fit = fit_nurbs(samples, 3; nctrl = 20)
println("fit: degree ", fit.degree, ", ", length(fit.control_points), " control points")
fit: degree 3, 20 control points

The engineering check on a fitted profile: does it reproduce the thickness the section is named for? At quarter-ish chord, x/c = 0.3, the NACA formula gives a thickness ratio of about 0.500 — and so does the fitted curve:

function fitted_y_at(fit, xc)
    best = Inf; by = 0.0
    for t0 in range(0, 1; length = 2001)
        p = curve_point(fit, t0)
        if abs(p[1] - xc) < best
            best = abs(p[1] - xc); by = p[2]
        end
    end
    by
end
println("thickness at x/c = 0.3: fit ", round(abs(fitted_y_at(fit, 0.3 * c)) / (c * t); digits = 4),
        "  vs NACA ", round(yt(0.3) / (c * t); digits = 4))
thickness at x/c = 0.3: fit 0.5004  vs NACA 0.5001

The sampled outline then becomes geometry: extruded into a wing section whose volume should be the closed-form NACA area (≈ 0.6853 · t · c²) times the span:

wing = extrude(polygon_2d([[p[1], p[2]] for p in samples]), 0.200)   # 200 mm span
println("volume (extruded) = ", round(ustrip(uconvert(u"mm^3", volume(wing))); digits = 0), " mm³")
println("analytic ≈ 0.6853 t c² × span = ", round(0.6853 * t * c^2 * 0.200 * 1e9; digits = 0), " mm³")
volume (extruded) = 656213.0 mm³
analytic ≈ 0.6853 t c² × span = 657888.0 mm³
vtk(surface_mesh(wing; n = 64))

The polygon’s shoelace area under-samples the smooth profile slightly — the two volumes agree to within a quarter of a percent, which is the difference between a point cloud and the curve that fits it.

Differentiation flows through

Evaluation is type-generic, so AD flows through the de Boor evaluation — the derivative of a NURBS curve is smooth within each knot span:

c = to_nurbs(CircleCurve(1.0, 0.0))
dx = derivative(t -> curve_point(c, t)[1], 0.25)
println("∂(curve.x)/∂t at t = 0.25: ", round(dx; digits = 3))
∂(curve.x)/∂t at t = 0.25: -5.657

Next steps