Bézier polynomial curve regression

Sharks. I never saw that coming.

June 09, 2026 · John Peach

Introduction

How often have you said to yourself, “Gosh, I sure wish I knew an easy way to approximate a $C^1$ curve with an $n$th-degree Bézier polynomial!” Several times a day would be my guess. Well, there is a way! In fact, just a little bit of linear algebra, orthogonal distance regression, and the Akaike criterion is about all you need. The linear algebra gives us a starting guess, orthogonal distance regression refines the control points by letting the curve slide its closest points onto the data, and the Akaike criterion tells us how many control points are actually worth keeping.

What is this $C^1$ stuff? It’s a way to describe how smooth a function or curve is. Smoothness depends on continuity and differentiability. It’s often possible to look at the plot of a function and see that it’s at least $C^1$. If the function is continuous and doesn’t have any sharp bends anywhere, you’ve probably got a $C^1$ or higher function. The functions we’ll study will be like this, that is, without corners and are continuous everywhere. The $C^n$ background section near the end of the article gives the precise definitions, and we’ll confirm that the fitted Bézier result inherits this smoothness.

Why Bézier polynomials? The advantage is that using only a few control points, the entire curve can be defined. Instead of using hundreds of points along the curve, you can generate the same curve from a handful of control points, and the Akaike criterion will tell how many of these control points are needed.

Bézier curves are derived from Bernstein polynomials, but were not applied until Paul de Casteljau derived a numerically stable method now called de Casteljau’s algorithm which he used for computer aided design at the French automaker Citroën. Pierre Bézier independently discovered the polynomials and used them to design car bodies at Renault. We used Bézier curves in Yacht Design with Mathematics. My interest in the inverse problem that is, recovering control points from curve data, started with a 2002 paper by Borges and Pastva.

Bézier curves are now applied in illustrator programs like Adobe Illustrator and Inkscape, in font design, in many Computer-Aided Design (CAD) software programs, and for robotics path planning. By manipulating the locations of the control points, you can quickly modify the shape of the resulting curve. The more difficult inverse problem, which we’ll solve here, is to reconstruct the locations of the control points from a given curve.

Program Outline

The Julia file BezierCurveFit.jl fits Bézier polynomials to a set of points $Q = (q_x,q_y)$ describing the curve to be fit. The polynomials are parameterized by $t \in [0,1]$, so that $q_x(i) = f_x(t_i)$ for some polynomial $f_x$ at $t_i$. Similarly, $q_y(i) = f_y(t_i)$. Since we don’t know the values for $t$, our first task is to estimate them from the measured points in $Q$ but, as we’ll see, if the $t_i$ aren’t known, we won’t get a good fit. The real trick is to stop treating $t$ as known. Instead of guessing $t$ and then solving for the control points, we let an optimizer adjust both the control points and the $t$ values at the same time. That joint optimization is called orthogonal distance regression (ODR) Orthogonal distance regression minimizes the perpendicular (shortest) distance from each data point to the fitted curve, rather than the vertical distance minimized by ordinary least squares. When data points can have errors in both $x$ and $y$ — as digitized curve points do — perpendicular distance is the geometrically correct thing to minimize., and it’s what turns a mediocre fit into a good one.

The curve points in $Q$ can be found by multiplying a matrix of coefficients, $A$, and the control points $P$, $Q = AP$. The inverse also works, so $P = A^{-1}Q$, (using some special tricks described below). To create $A$, we need to know how many control points there are, and we need to generate the Vandermonde matrix from $t$.

If we knew $t$ exactly, we could perfectly generate $A$, and we’d know where the control points $P$ were. Errors in our estimate of $t$ cause fairly large errors when we calculate the estimated control points, $P_{est} = A^{-1}Q$, and consequently errors in locations of the estimated curve points, $Q_{est} = AP_{est}$. The points in $Q_{est}$ are interpolated to the nearest approximation of points in $Q$.

Using some optimization techniques, we can iterate on the estimates of the locations of the control points until the error between $Q_{est}$ and $Q$ is as small as is computationally possible. Critically, we let the parameter values $t$ move during this iteration too. Each $t_i$ slides until $B(t_i)$ arrives on the point of the curve closest to the measured $Q_i$. Minimizing that perpendicular (geometric) distance, rather than the distance at a fixed, and possibly wrong $t$, is the defining feature of orthogonal distance regression.

In most cases, the curve points $Q$ weren’t generated using Bézier polynomials, so the approximating polynomials in this process depend on the number and locations of the control points $P$. Starting with a small number of control points, we can use this procedure to get $Q_{est}$, then increase the number of control points to generate a new $Q_{est}$, repeating until the fit is good enough. “Good enough” is exactly the judgment the Akaike criterion The Akaike Information Criterion (AIC) was introduced by statistician Hirotugu Akaike in 1974. It balances model fit against complexity by adding a penalty proportional to the number of free parameters — the mathematical formalization of Occam’s Razor. Lower AIC is better; what matters is differences in AIC between models, not the absolute value. automates. More control points always reduce the residual error, so we need a penalty that punishes extra parameters and stops us from overfitting.

The outline for the curve fitting function is:

• Read in measured curve points $Q$ (an $n \times 2$ array from WebPlotDigitizer, described below).
• Estimate the parameter vector $t$ from the distances between successive points in $Q$ (the chord length). chord_length_param(Q) $\rightarrow$ t.
• Construct the Bernstein basis matrix $A$. bernstein_matrix(t, m) $\rightarrow$ A (m = degree, so $A$ has m+1 columns).
• Warm-start the control points by endpoint-constrained least squares. init_control_points(Q, t, m) $\rightarrow$ P.
• Jointly optimize P and t by orthogonal distance regression. fit_bezier_odr(Q, m) $\rightarrow$ P_fit, t_fit.
• Select the best degree $m$ by residual-elbow detection (with AIC fallback). select_degree_robust(Q) $\rightarrow$ m_best, P_best.
• Optionally enforce a monotone parameterization. monotone_reproject(Q, P_fit) $\rightarrow$ t_monotone.

The first half of this article builds up the linear algebra behind bernstein_matrix and init_control_points; the second half explains the orthogonal distance regression in fit_bezier_odr, the degree selection in select_degree_robust, and the reprojection in monotone_reproject.

Bézier Polynomials

In the article Yacht Design with Mathematics we used Bézier polynomials to approximate the curves of a boat hull. The Bézier polynomial is defined as

The $P_i$ array represents the locations of $n+1$ control points and can have any $(x,y)$ coordinates you like. Notice that the index starts at $i=0$ which is why there are $n+1$ rather than $n$ control points. (In the code sections later, we follow the Julia convention and call the degree $m$, so the curve has $m+1$ control points — the math here uses $n$ to keep the standard textbook notation.) The parameter $t$ goes between $0$ and $1$, and the symbol

is the binomial coefficient, and $k! = k \cdot (k-1) \cdot (k-2) \cdots 1$ is $k$ factorial. The binomial coefficient tells how many ways you can choose $i$ things from a set of $n$, and it also gives the coefficients in the expansion of $(x+y)^n$. You can check the expansion coefficients by typing

expand((x+y)^4)

in Wolfram Alpha

and checking that the same coefficients appear in the vector,

n = 4
[binomial(n,k) for k = 0:n]'
[1, 4, 6, 4, 1]

when you run it in Julia. You could have gotten each coefficient by typing binomial(4,0), binomial(4,1), …, binomial(4,4), but the notation above is a shortcut to generate them all at once, called list comprehension.

Another way to get the coefficients is to use Pascal’s Triangle. Each row represents the $n$th set of expansion coefficients. The zeroth row is the single $1$ at the top, so the expansion for $(x+y)^4$ is the fifth row down, $1,4,6,4,1$. If you want to explore more about Pascal’s Triangle, read Zaphod Beeblebrox’s Brain and the Fifty-ninth row of Pascal’s Triangle by Andrew Granville.

The location of the points $P_i$ controls the shape of the curve. Here’s a curve created with five control points that looks a little like a shark fin drawn with GeoGebra. Notice that the control points are not in order. Place five points anywhere on the canvas, and then type

curve[x(A)*(1−t)^4+4*x(B)*(1−t)^3*t+6*x(C)*(1−t)^2*t^2+4*x(D)*(1−t)*t^3+x(E)*t^4,y(A)*(1−t)^4+4*y(B)*(1−t)^3*t+6*y(C)*(1−t)^2*t^2+4*y(D)*(1−t)*t^3+y(E)*t^4,t,0,1]

in the Input: dialog box.

You can find any point on the curve by varying $t$ between $0$ and $1$, and using the formula for $B(t)$ above. You can do this in a spreadsheet:

In the first column of the spreadsheet I entered the point labels and in the next column is the index, $i$, running from $0$ to $4$. Below the main table, I entered two labeled cells, $n$ and $t$. The value for $n$ should be one less than the number of points since indexing begins at zero. You can put any value for $t$ so long as it’s in the range $[0,1]$. The third column is the binomial expansion coefficient for $n$ and $i$. Excel uses =combin(n,a2).

The third and fourth columns, “D” and “E” contain the powers of $t^i$ and $(1-t)^{(n-i)}$. In the next two columns, put in the $x,y$ coordinates for each of the five points. Right click on point A and click on “Object Properties” to open a window on the right side of GeoGebra. Under the “Basic” tab you should see the name of the point and the value. Copy the values for the $x$ and $y$ coordinates into the columns labelled P_x and P_y in the spreadsheet. If you’ve put point A at the origin, GeoGebra will have Definition: Intersect(xAxis, yAxis), but just enter zeros.

Any number raised to the $0^{th}$ power is $1$, but Excel incorrectly assumes the answer is undefined. If this happens with your spreadsheet, replace the cells with $1$‘s. In the example, this would be in cells D2 and E6. Or, you could use the open-source LibreOffice suite which correctly calculates zero powers.

Using the formula,

we can calculate Q_x and Q_y for each point, which are the products of columns C,D,E, and F for Q_x and C,D,E, and G for Q_y. Then, sum the columns to get the coordinates for $Q$. In GeoGebra, add a new point and rename it $Q$. Right click on “Q” and enter the calculated coordinates in the Value box. It should fall somewhere on the shark fin plot. If you change the value of $t$ in the spreadsheet, you’ll get a new position for $Q$ somewhere on the plot.

After choosing the locations for control points $A,B,C,D$, and $E$ you can generate the shark fin curve by varying the value of $t$ between $0$ and $1$.

Bézier Linear Algebra

Let’s look at $B(t)$ again,

Suppose we have $5$ control points, $\{ P_0, P_1, P_2, P_3,P_4 \}$, so $n = 4$ since the index begins at zero. If we write out the terms of $B(t)$ we get

Each row of the matrix $R$ is binomial expansion of $n-k, \; k = 0, \dots, n-1$. Both $\circ$ and $\odot$ are used here for element-wise (Hadamard) multiplication of two matrices with the same dimensions — the distinction in what follows is purely structural: $\circ$ multiplies two full matrices entry-by-entry, while $\odot$ broadcasts a column vector across all columns of a matrix (the Julia dot-broadcast operation). The sign matrix $S$ with alternating $\pm1$‘s changes every other sign of the values in the matrix $R$. The binomial expansion $\begin{bmatrix} 1 & 4 & 6 & 4 & 1 \end{bmatrix}^T$ multiplies each column of $R \circ S$ by the appropriate coefficient, and is the first row of $R$, denoted by $R[1,:]^T$. Multiplying the coefficient matrices together gives

Since we know how many control points there are, we can quickly create the coefficient matrix $C$. Multiply the control points array $P$ by the coefficient matrix to get

The last step is to create the vector of $t$ values between $0$ and $1$ using the Julia range function:

t = range(0, 1, num = nPts)

Setting nPts = 100 gives 100 points along the curve. For each value of $t$ we need $\begin{bmatrix} t^4 & t^3 & t^2 & t & 1 \end{bmatrix}$, which is a Vandermonde matrix, but flipped left and right. For any input vector $v$, with exponents up to $n$, this code returns the Vandermonde matrix, $T$, which has size $nPts \times (n+1)$:

VM = Matrix{Float64}(undef, length(v), n + 1);
for k in 0:n
    VM[:, k + 1] = v .^ (n - k)
end

The Bézier curve points can be found by multiplying $B = T \times C \times P$, which will be $nPts \times 2$. The first column contains the $x$-coordinates and the second the $y$-coordinates. To check that the math is correct, I plotted the Bézier curve of the shark fin in Julia and then overlaid it onto the GeoGebra curve. The red dashed line is the curve generated by GeoGebra, and the solid blue is the Julia curve.

But, suppose you had the plot already, and wanted to know where the control points should be placed to reproduce the curve? We’d like to be able to reverse the process to get the control points.

Curve point coordinates

To get the coordinates of the curve points in GeoGebra first turn off extraneous points. On the left side of the window, under Algebra is a list of the control points with small dots to the left of the label. Click on each one to hide the point in the graph. Save the graph as an image using File $\rightarrow$ Export $\rightarrow$ Graphics View as Picture (png, eps).

Now, go to the WebPlotDigitizer site and click on Launch Now! Under File select Load Image(s) and Browse to find your image of the curve. The default type is 2D (X,Y) Plot which is right for our plot. Click Align Axes at the bottom and Proceed on the next pop-up window. Click on a known $x$-coordinate somewhere on the left side and then a second known coordinate to the right. Repeat this for two known $y$-coordinates. You can use the origin for both if you like. If you need to adjust points, click on the point and use the arrow keys to move it. When the points are in the correct positions, click on Complete! The next dialog asks for the coordinates of these points.

Under Automatic Extraction, click on the blue square to the right of Foreground Color, then choose Color Picker and click on a point of the curve, then Done. Click Run and you should see points that WebPlotDigitizer has selected along the curve. If any points don’t seem to be on the curve, click Adjust Point (S) and use the arrow keys to move them back onto the curve.

On the left side under Dataset click View Data, Sort by: Nearest Neighbor and then Download .CSV. The Nearest Neighbor sorting option puts the points in order along the length of the curve. Name the file something like SharkFinCurve.csv and save it. Open the file with a text editor such as Notepad++. Bézier polynomials go through the first and last control points exactly, so if WebPlotDigitizer didn’t pick up these points, add them in to the list of curve points. Do this by entering $t=0$ and $t=1$ in the spreadsheet and copying the summed values for Q_x and Q_y for each.

Finding the Control Points

The equation for the Bézier polynomial is the sum of coefficients times the control point coordinates. We could write it as

A more compact way to write it would be as the multiplication of two vectors (multiply corresponding entries, sum the whole thing),

If we use $t = 0.2$ the result is the point $Q = [0.826948, 1.072299]$. We could do this multiplication for many different values of $t$ along the curve and get the matrix multiplication,

In the example, there are five control points $P_i$ and WebPlotDigitizer generated $k = 84$ steps along the curve. The matrix on the left is $84$ rows by $5$ columns, and the vector of $Q$‘s is $84$ points long.

Compactifying a bit more, we can write this system of equations as $A P = Q$ where $A$ is the matrix of coefficients on the left, $P$ is the vector of control points, and $Q$ is the vector of measured points along the curve. We need to find a vector $P$ that makes the difference between the estimated points and the measured points as small as possible. That is, we want to minimize the error $\lVert AP-Q \rVert$, but it’s slightly more convenient to minimize the square, $(AP-Q)^2$ which gives another solution.

If we expand the square we get

recalling that $(AB)^T = B^T A^T$ for compatible matrices. To minimize the error we can take the derivative, and set the result to zero,

The solution for $P$ is simply,

In Julia, the numerically stable way to evaluate this least-squares solution is the backslash operator, which avoids forming and inverting $A^TA$ explicitly:

P = bernstein_matrix(t, m) \ Q

so you just need to get $A$ and then the calculation is easy. There’s a minor flaw with this method, though. The $Q$ points along the curve seem to be uniformly spaced, but that doesn’t mean the values of $t_i$ are. We don’t know what $t_i$ should be for each $Q_i$, but we know that $t_0 = 0$ and $t_n = 1$, so we could calculate the cumulative distance along the curve and normalize by the total length to get approximate values for each $t_i.$ Keep this flaw in mind — it’s the whole reason we’ll need orthogonal distance regression. The least-squares solution here isn’t our final answer; it’s a warm start that gets the optimizer into the right neighborhood.

Building the A matrix

The “Bézier Linear Algebra” section above expanded the basis into monomials — a coefficient matrix $C$ times a Vandermonde matrix in $t$ — which is great for seeing where the binomial coefficients come from. In code, though, it’s both shorter and numerically safer to evaluate the Bernstein basis directly from its definition,

where $m$ is the polynomial degree (the same thing the math above called $n$), so there are $m+1$ control points. This is exactly bernstein_matrix in BezierCurveFit.jl:

function bernstein_matrix(t::AbstractVector{Float64}, m::Int)::Matrix{Float64}
    n = length(t)
    B = zeros(Float64, n, m + 1)
    for j = 0:m
        c = binomial(m, j)
        for i = 1:n
            B[i, j + 1] = c * (t[i]^j) * ((1.0 - t[i])^(m - j))
        end
    end
    return B
end

From here on I’ll switch to the code’s convention: m is the degree, and the curve has m + 1 control points. The returned matrix is $n \times (m+1)$ — one row per measured point, one column per control point — and evaluating the curve is then just the matrix product bernstein_matrix(t, m) * P, which is what the helper bezier_eval does.

We still don’t know the true values for $t$, though. The simplest guess is to space them evenly between $0$ and $1$, but a better guess measures the distances between the points in $Q$, the array returned by WebPlotDigitizer:

function chord_length_param(Q::Matrix{Float64})::Vector{Float64}
    n = size(Q, 1)
    d = [norm(Q[i+1, :] - Q[i, :]) for i in 1:(n-1)]
    s = cumsum([0.0; d])
    return s ./ s[end]
end

This is the chord-length parameterization The chord-length parameterization assigns $t$ values proportional to the cumulative straight-line (chord) distances between successive data points, as a proxy for arc length. It distributes parameter values roughly uniformly along the curve, avoiding the clustering that would result from spacing $t$ evenly when the data points themselves are unevenly spaced.. Each d[i] is the Euclidean distance between successive points in $Q$,

The cumulative sum s is the running arc length with a $0$ prepended, and dividing by the total length s[end] normalizes everything into $[0,1]$ with $t_1 = 0$ and $t_n = 1$.

The result is:

There’s very little difference between even spacing and chord length for densely sampled curves. Using the basis matrix at these estimated $t$ values, we can take a first stab at the control points with an ordinary least-squares solution, and reconstruct the curve from them:

A     = bernstein_matrix(t, m)   # n × (m+1)
P_est = A \ Q                    # least-squares solve, (m+1) × 2
Q_est = A * P_est

Julia’s backslash operator \ solves the overdetermined system $A P = Q$ in the least-squares sense. Plotting the estimated curve points gives a result that is not very shark-fin-ish:

The problem is that even small errors in the estimates of $t$ for each of the measured points in $Q$ can result in large errors in the estimated control points $P_{est}$.

But, notice that the blue curve generated from the estimated control points $P$ overlays the fit generated by GeoGebra, which uses a different $t$ vector from the one our least-squares solve used. The implication is that we should be able to match the original shark-fin curve if we can somehow maneuver the control points into place and, just as importantly, let the $t$ values move too.

A better warm start with endpoint-constrained least squares

Before handing the problem to an optimizer, we can do better than the raw least-squares solution. A Bézier curve passes exactly through its first and last control points, and we already know where the curve starts and ends — those are just the first and last measured points, $Q_1$ and $Q_n$. So we can fix $P_0 = Q_1$ and $P_m = Q_n$ and solve the least-squares problem for only the interior control points. The function init_control_points does this:

function init_control_points(Q::Matrix{Float64}, t::Vector{Float64}, m::Int)::Matrix{Float64}
    if m == 1
        return [Q[1, :]'; Q[end, :]']
    end

    A       = bernstein_matrix(t, m)          # (n, m+1)
    A_fixed = A[:, [1, m+1]]                  # endpoint basis columns
    A_free  = A[:, 2:m]                       # interior basis columns (n, m-1)
    rhs     = Q - A_fixed * Q[[1,end], :]     # subtract fixed contributions
    P_free  = A_free \ rhs                    # (m-1, 2), overdetermined LS

    return vcat(Q[[1],:], P_free, Q[[end],:])
end

The trick is to move the known endpoint contributions to the right-hand side. We split the basis matrix into its two endpoint columns A_fixed and its interior columns A_free. The endpoints contribute A_fixed * [Q₁; Qₙ] to every fitted point, so we subtract that off and solve A_free \ rhs for the interior points alone. Stacking the fixed endpoints back around the solved interior points gives an $(m+1)\times 2$ control-point matrix that already honors the endpoints and sits close to the right answer. This is the warm start, and now we’ll refine it.

Orthogonal Distance Regression

Ordinary least squares minimizes

with the $t_i$ held fixed at whatever we guessed. The trouble, as the table above showed, is that a wrong $t_i$ drags the control points to the wrong place. Orthogonal distance regression (ODR) fixes this by minimizing over the control points and the parameter values at the same time,

subject to the endpoint constraints $P_0 = Q_1$, $P_m = Q_n$, $t_1 = 0$, and $t_n = 1$. Because each interior $t_i$ is now free to move, at the optimum $B(t_i)$ lands on the point of the curve closest to $Q_i$. The residual we minimize is therefore the perpendicular, or geometric, distance from each data point to the curve which is exactly the closest-point problem we’ll work out by hand in the next section. ODR does it for us, for every point, on every iteration.

The Julia package Odrpack.jl wraps the Fortran ODRPACK95 solver, whose engine is a trust-region Levenberg–Marquardt minimizer designed for exactly this problem (Zwolak, Boggs & Watson, Algorithm 869, ACM TOMS 33(4), 2007). It minimizes the general weighted form

Here $\beta$ holds the model parameters (our flattened control points), $x$ is the explanatory variable (our parameter vector $t$), $y$ is the response (our data $Q$), and $\delta_i$ is the adjustment the solver makes to each $x_i$. The weight $w_{x,i}$ penalizes moving $t$; driving it toward zero makes $t$ effectively free and recovers pure geometric distance. Setting it to exactly zero is numerically degenerate in ODRPACK95, so the code uses a tiny value, $w_x = 10^{-10}$ making the parameters effectively free, while keeping the solver well-conditioned.

ODRPACK calls a model function with the signature f!(x, beta, y) that writes the curve points into y in place. Ours just reshapes beta into control points and evaluates the curve:

function bezier_fcn!(x::AbstractVector{Float64},
                     beta::Vector{Float64},
                     y::Matrix{Float64})
    y .= bezier_eval(x, beta_to_P(beta))
    return nothing
end

The control points live in beta as an interleaved vector $[P_{0x}, P_{0y}, P_{1x}, P_{1y}, \dots]$, so we convert back and forth with a reshape:

beta_to_P(beta) = Matrix(reshape(beta, 2, :)')   # interleaved vector → (p × 2)
P_to_beta(P)    = vec(P')                          # (p × 2) → interleaved vector

Putting the warm start, the constraints, and the model together gives fit_bezier_odr:

function fit_bezier_odr(Q::Matrix{Float64}, m::Int;
                        tol::Float64=1e-10,
                        reproject::Bool=false)
    n = size(Q, 1)



    t_init = chord_length_param(Q)
    P_init = init_control_points(Q, t_init, m)
    beta_init = P_to_beta(P_init)



    fix_b = falses(length(beta_init))
    fix_b[1:2] .= true       # fix P₀ = Q[1,:]
    fix_b[end-1:end] .= true # fix Pₘ = Q[end,:]



    fix_x_arr = falses(n)
    fix_x_arr[1] = true      # fix t₁ = 0
    fix_x_arr[end] = true    # fix tₙ = 1



    result = odr_fit(
        bezier_fcn!,
        t_init,
        Q,
        beta_init;
        fix_beta=fix_b,
        fix_x=fix_x_arr,
        jac_beta! = bezier_jac_beta!,
        jac_x! = bezier_jac_x!,
        weight_x = fill(1e-10, n)
    )



    P_fit = beta_to_P(result.beta)
    t_fit = vec(result.xplusd)



    if reproject
        t_fit = monotone_reproject(Q, P_fit; tol=tol)
    end



    return P_fit, t_fit, result
end

The fix_beta mask freezes the first and last control points; fix_x freezes the first and last $t$ values at $0$ and $1$. Everything else, the interior control points and every interior $t_i,$ is free to move. The optimized parameter values come back as result.xplusd (the initial x plus the fitted adjustments $\delta$).

Analytic Jacobians

We could let ODRPACK estimate the derivatives by finite differences, but the Bernstein basis has analytic derivatives, so supplying them is faster and more accurate. The curve is linear in the control points, so the Jacobian with respect to $\beta$ is just the basis itself. Each coordinate of the output depends only on the matching coordinate of each control point, which is why the entries land on alternating slots:

function bezier_jac_beta!(x, beta, jb)
    m = (length(beta) ÷ 2) - 1
    n = length(x)
    B = bernstein_matrix(x, m)
    fill!(jb, 0.0)
    for j in 1:(m + 1)
        for i in 1:n
            jb[i, 2*j - 1, 1] = B[i, j]
            jb[i, 2*j, 2]     = B[i, j]
        end
    end
    return nothing
end

The Jacobian with respect to $t$ needs the derivative of the curve, which uses the derivative of the Bernstein basis,

with the convention $B_{m-1,-1} \equiv 0$ and $B_{m-1,m} \equiv 0$. That recurrence is bernstein_deriv, and bezier_jac_x! multiplies it by the control points to get $dB/dt$ at each point. With both Jacobians supplied, the solver recovers control points that put the curve right on top of the data, because the optimizer now has the freedom to slide the parameter values into place.

Choosing the Degree

We’ve been assuming we know how many control points to use, but in practice we won’t. Adding control points (raising the degree $m$) can only lower the residual error, and a high enough degree will fit the curve through every measured point even if it’s noisy. This is called overfitting, and looking at the residual error alone will never tell you that the curve has been over fit. We need a score that rewards a good fit but penalizes extra parameters.

Degree selection with AIC

The Akaike Information Criterion (AIC) gives this score for good fits,

where $k$ is the number of free parameters and $\hat{L}$ is the maximized likelihood. For a model with Gaussian residuals and unknown variance, the log-likelihood term reduces to $n\ln(\text{SSE}/n)$ plus a constant (after profiling out the variance), where SSE is the sum of squared residuals. The constant is the same for every degree, so it drops out of the comparison, leaving the form select_degree_aic uses:

A degree-$m$ curve has $m+1$ control points, but the two endpoints are fixed, leaving $m-1$ free interior control points. Each has an $x$ and a $y$ coordinate, so there are $2(m-1)$ free parameters, and the penalty is $2k = 4(m-1)$, which is the final term above. We pick the degree with the lowest AIC:

function select_degree_aic(Q::Matrix{Float64};
                           m_min::Int=2, m_max::Int=8,
                           tol::Float64=1e-10)
    n = size(Q, 1)
    m_range = m_min : min(m_max, n-2)
    aic = Dict{Int,Float64}()
    best_aic = Inf
    m_best = m_min
    P_best = zeros(0, 2)
    t_best = zeros(0)



    for m in m_range
        P_fit, t_fit, _ = fit_bezier_odr(Q, m; tol=tol)
        SSE = sum((bezier_eval(t_fit, P_fit) .- Q).^2)
        aic[m] = n * log(max(SSE, eps()) / n) + 4*(m-1)
        if aic[m] < best_aic
            best_aic = aic[m]
            m_best = m
            P_best = P_fit
            t_best = t_fit
        end
    end
    return m_best, P_best, t_best, aic
end

The loop fits every degree from $2$ up to $\min(8, n-2)$ by ODR, computes the SSE at the fitted parameters, scores it, and keeps the best. The max(SSE, eps()) guard keeps the logarithm finite if a fit is essentially perfect.

Note the specific form of the penalty: a degree-$m$ curve has $m+1$ control points, but the two endpoints are fixed, leaving $m-1$ free interior control points. Each has an $x$ and a $y$ coordinate, so there are $k = 2(m-1)$ free parameters, and substituting into the standard $2k$ penalty gives $4(m-1)$, the final term in the formula above.

There are three further caveats worth keeping in mind:

• The latent $t$ values aren’t counted as parameters. ODR estimates $n-2$ interior $t_i$ alongside the control points, which could be thought of as additional degrees of freedom. The usual defense is that they are nuisance parameters bound one-to-one to the data points rather than free model structure, so they don’t inflate model complexity the way control points do.
• No term for the variance parameter. Strictly, estimating $\sigma^2$ adds one parameter, so some authors write the penalty as $2(k+1)$. That shifts every AIC by the same constant, so it doesn’t change which degree wins.
• Small samples. When $n$ is not large relative to $k = 2(m-1)$, the bias-corrected $\text{AIC}_c = \text{AIC} + \frac{2k(k+1)}{n-k-1}$ is safer and penalizes high degrees harder. For the shark fin with dozens of curve points and a handful of control points, the plain AIC and $\text{AIC}_c$ will usually agree, but it’s worth checking on sparse data.

The function select_degree_aic works well when the data is noisy, because noise creates a floor on SSE that the penalty can eventually overcome, and the AIC minimum selects a sensible degree. For the shark-fin curve with $n=30$ points and Gaussian noise $\sigma=0.05$, the AIC minimum is $m^*=4$, which is exactly the true degree.

When AIC misleads: the SSE floor problem

On clean, digitized data the story is different. The shark-fin CSV (SharkFinCurve.csv, $n=84$ points, no measurement noise) produces the following table when select_degree_aic sweeps $m = 2, \ldots, 10$:

Once the degree reaches 4, SSE drops to numerical zero and $R^2 = 1.000$. The log-likelihood term saturates, and AIC keeps improving only because the penalty $4(m-1)$ grows more slowly than the marginal gain from floating-point noise. The same failure occurs with $\text{AIC}_c$ and BIC. All three criteria prefer the maximum degree in the search range, wherever that happens to be set.

The real decision boundary is visible in the max residual column: a 58× drop between $m=3$ (0.142) and $m=4$ (0.0025), after which further improvement is less than 30% across six additional degrees. AIC can’t see this because it operates on SSE, which is already near zero. The geometric residual is what carries the signal.

Automatic degree selection: select_degree_robust

select_degree_robust in DegreeSelection.jl automates the elbow argument above. Its algorithm has three steps.

Step 1 — Sweep. Call fit_bezier_odr for every $m$ in $[m_{\min}, \min(m_{\max}, n-2)]$, recording SSE, AIC, and the maximum geometric residual $r_{\max}(m) = \max_i \lVert B(t_i; P) - Q_i \rVert$ at each degree.

Step 2 — Elbow detection. For each consecutive pair of degrees, compute the improvement ratio

The elbow is the first $m$ where $\rho(m) > \tau$ (default $\tau = 10$). If $r_{\max}(m) < \varepsilon_\text{machine}$ the ratio is treated as $\infty$, so a numerically exact fit always triggers the elbow. If no ratio exceeds $\tau,$ which happens on noisy data where the residual declines smoothly, the function falls back to the AIC minimum, matching the behavior of select_degree_aic.

Step 3 — Plateau selection. Among all degrees $m \geq m_\text{elbow}$ where

return the smallest such $m$.

function select_degree_robust(Q::Matrix{Float64};
                              m_min::Int    = 2,
                              m_max::Int    = 12,
                              τ::Float64   = 10.0,
                              ε::Float64   = 2.0,
                              tol::Float64 = 1e-10)
    sweep   = fit_sweep(Q, m_min, m_max, tol)
    degrees = sort(collect(keys(sweep)))



    m_elbow = _detect_elbow(sweep, τ)



    if m_elbow === nothing
        m_elbow = degrees[argmin([sweep[m].AIC for m in degrees])]
        println("[select_degree_robust] No elbow found; " *
                "falling back to AIC minimum at m=$m_elbow")
    end



    m_best = _plateau_select(sweep, m_elbow, ε)



    println("[select_degree_robust] Elbow at m=$m_elbow. " *
            "Selected m=$m_best (smallest in plateau). " *
            "AIC=$(round(sweep[m_best].AIC; digits=2))")



    # diagnostics Dict omitted for brevity — see DegreeSelection.jl
    return m_best, sweep[m_best].P, sweep[m_best].t
end

On the shark-fin curve, the function prints:

[select_degree_robust] Elbow at m=4 (ratio=57.8). Selected m=4 (smallest in plateau). AIC=-1245.85

The ratio 57.8 is the measured $\rho(4)$: the max residual fell by nearly 58× when the degree increased from 3 to 4.

The sharp drop between $m=3$ and $m=4$ is the elbow that select_degree_robust detects. AIC continues to improve above $m=4$ because SSE has already reached numerical zero, so the penalty no longer dominates. The dashed line marks the selected degree $m^*=4$.

The Akaike Information Criterion isn’t the only model selection tool. The Bayesian Information Criterion

replaces the fixed penalty $2k = 4(m-1)$ with $k\ln(n) = 2(m-1)\ln(n)$, which grows with sample size and consistently favors lower degrees than AIC. A good overview of both is in Choosing the Best Model: A Friendly Guide to AIC and BIC. On the clean shark-fin data, BIC and $\text{AIC}_c$ both suffer the same SSE-floor failure as AIC, and all three prefer the maximum degree in the sweep range.

All three information criteria continue improving past the elbow at $m=4$ because SSE is at numerical zero for $m \geq 4$. The residual panel (log scale) shows the true signal: a 58× drop at $m=4$, with less than 30% further improvement across six additional degrees. The vertical gray line marks the elbow.

To summarize the selection strategy:

• Noisy data ($\sigma > 0$, SSE declines smoothly): use select_degree_aic. The AIC penalty will dominate once the degree exceeds the signal in the data, and the minimum is meaningful.
• Clean or near-clean data (digitized curves, CAD traces, noise-free simulations): use select_degree_robust. AIC will overfit; the geometric elbow will not.
• When in doubt: call select_degree_robust with default parameters. On noisy data it falls back to the AIC minimum automatically, so it is strictly more general than select_degree_aic.

Monotone reparameterization

The ODR routine is free to move each $t_i$ wherever it lowers the distance, and nothing in the objective forces the $t_i$ to stay in order. Usually they do, but on a curve that doubles back on itself, two data points can be assigned $t$ values that are out of sequence. If you need a strictly increasing parameterization for resampling, arc-length work, or just tidiness monotone_reproject re-solves for $t$ point by point, holding each one to the right of its predecessor:

function monotone_reproject(Q::Matrix{Float64}, P::Matrix{Float64};
                             tol::Float64=1e-10)::Vector{Float64}
    n = size(Q, 1)
    t_proj    = zeros(n)
    t_proj[n] = 1.0
    for i in 2:n-1
        t_lo = t_proj[i-1]
        dist2(t) = sum((bezier_eval([t], P)[1, :] .- Q[i, :]) .^ 2)
        t_proj[i] = _golden_section(dist2, t_lo, 1.0; tol=tol)
    end
    return t_proj
end

For each interior point it minimizes the squared distance to the curve over the interval $[t_{i-1}, 1]$ with a golden-section search, which guarantees the new $t_i$ can’t fall below $t_{i-1}$. The control points are left untouched and only the parameterization changes so the SSE may increase slightly relative to the unconstrained ODR fit. Notice that this per-point closest-point search is exactly the geometric problem of the next section, solved numerically one point at a time. ODR was solving the same problem for all points jointly.

Note: The two sections that follow (“Finding the closest point by hand” and “Cubic and Quadratic Spline Interpolation”) are geometric background. None of it is needed to run the fit — ODR handles the closest-point problem automatically. These sections are here for readers who want to understand what the optimizer is minimizing, or who need to compute closest points independently of the fitting pipeline.

Finding the closest point by hand

Rather than trying to find the optimal $t$ vector, suppose we wanted to find the best values for the control point locations using the $t$ we have. The locations of the points in $Q_{est}$ are determined by the positions of the control points $P$, so the goal is to move the control points around until the difference between $Q_{est_i}$ and $Q_i^*$ is minimized. A couple of ways to minimize these distances are to minimize the sum of the errors, or to minimize the worst case, or maximum error. This is the manual version of the orthogonal distance that ODR minimizes automatically.

Perhaps the easiest way to find the closest point on the curve $Q$ to an estimated point $Q_{est_i}$ would be to draw straight lines between successive points in $Q$ and then get the minimum distance between $Q_{est_i}$ and the nearest line segment. This will be perpendicular to the line segment and passing through $Q_{est_i}$.

This works well so long as $Q_{est}$ can be reached by a perpendicular segment (Case 1). But as shown in Case 2, if $Q_{est}$ is near node $Q_i$, and between the endpoint perpendiculars, the method fails. A sharper turn at $Q_i$ leaves more room for an undefined condition, but in most cases, we’ll have lots of points along the length of $Q$, and the turns will be fairly shallow. If $Q_{est}$ falls in the undefined region, the easiest solution would be to simply measure the distance between $Q_{est}$ and $Q_i$.

To find $Q^*$ in Case 1, lets parameterize the segment from $Q_{i-1}$ to $Q_i$ by $\tau$ where $0 \leq \tau \leq 1$, so the segment becomes

If $\tau = 0$, then you’re at $Q_{i-1}$ and if $\tau = 1$ you’ve arrived at $Q_i$. The point $Q^*$ will be at some value of $\tau$ where the vector $Q_{est} - Q^*$ is perpendicular to the vector $Q_i - Q_{i-1}$, or equivalently, the dot product is zero,