Re-implement arc length reparametrization

[ecurtiss]

The ability to get uniformly spaced points along a spline was removed in v0.4.0 because it was too slow to be realistic to use. Today I investigated two faster methods.

The idea is that instead of solving for the position at t in [0, 1], we instead solve for the position such that the percent arc length is t. The formula for arc length is

and we want to solve for x for a given L. In our case, this integral cannot be evaluated analytically, but if we suppose that it has an antiderivative F, then by the Fundamental Theorem of Calculus,

At t = 0, the arc length should be 0, so F(0) = 0. This gives

or alternatively

which is the form of a root finding problem!

Furthermore, we can also find the derivative of the function

This allows us to perform Newton's method.

Algorithm comparison

I compared 3 algorithms: naive (sum distances between points until you reach the desired length), bisection with Gaussian quadrature, and Newton's method with Gaussian quadrature. Try it for yourself: CatRomReparameterizeBenchmark.zip. Newton is the clear winner.

Empirically, Newton tends to take 2-5 iterations whereas bisection can take 15-20.

[ecurtiss]

I found a paper that claims that Newton's method should be combined with bisection in order to avoid cases where the x-value falls outside of the desired interval.

Using a 5-point integration rule, Newton tends to diverge more around 0.5 (especially when the curve is symmetric) but less in [0.99, 1]. Using a 10-point rule, neither diverge and their speeds are the same. However, Newton fails horribly when the tension is cranked up (in particular where curvature is high), so we will use the hybrid method instead.

The methods for reference:

function Spline:ReparameterizeNewton(s: number)
    if s == 0 or s == 1 then
        return s
    end

    -- Newton's method
    local integrand = function(x)
        return self:SolveVelocity(x).Magnitude
    end

    local t = s
    for i = 1, MAX_NEWTON_ITERATIONS do
        local f =  GaussLegendre.Ten(integrand, 0, t) / self.length - s
        if math.abs(f) < EPSILON then
            return t
        end
        local g = self:SolveVelocity(t).Magnitude / self.length
        t -= f / g
    end

    warn("Failed to reparameterize; falling back to input")
    return s
end

function Spline:ReparameterizeHybrid(s: number)
    if s == 0 or s == 1 then
        return s
    end

    -- Hybrid of Newton's method and bisection
    -- https://www.geometrictools.com/Documentation/MovingAlongCurveSpecifiedSpeed.pdf
    local t = s
    local lower = 0
    local upper = 1
    local iterate = 0

    local v = function(x)
        return self:SolveVelocity(x).Magnitude
    end

    while iterate < MAX_NEWTON_ITERATIONS do
        iterate += 1

        local f =  GaussLegendre.Ten(v, 0, t) / self.length - s

        if math.abs(f) < EPSILON then
            return t
        end

        local g = self:SolveVelocity(t).Magnitude / self.length
        local candidate = t - f / g

        if f > 0 then
            -- Solution is below the current t
            upper = t
            t = candidate <= lower and (upper + lower) / 2 or candidate
        else
            -- Solution is above the current t
            lower = t
            t = candidate >= upper and (upper + lower) / 2 or candidate
        end
    end

    -- Failed to reparameterize; return input as "good enough"
    warn("Failed to reparameterize; falling back to input")
    return s
end

[ecurtiss]

Another idea is to approximate the parameter -> arc length parameter function. I tried approximating the function using a basis of Legendre polynomials. Sometimes it's spot-on, sometimes it's quite bad.

local Polynomial = require(script.Parent.Polynomial)
local GaussLegendre = require(script.Parent.GaussLegendre)

-- Precompute 20-point Gaussian quadrature evaluations
local function Precompute(f)
    return {
        [-0.07652652113349734] = f(-0.07652652113349734),
        [0.07652652113349734]  = f(0.07652652113349734),
        [-0.22778585114164507] = f(-0.22778585114164507),
        [0.22778585114164507]  = f(0.22778585114164507),
        [-0.37370608871541955] = f(-0.37370608871541955),
        [0.37370608871541955]  = f(0.37370608871541955),
        [-0.5108670019508271]  = f(-0.5108670019508271),
        [0.5108670019508271]   = f(0.5108670019508271),
        [-0.636053680726515]   = f(-0.636053680726515),
        [0.636053680726515]    = f(0.636053680726515),
        [-0.746331906460150]   = f(-0.746331906460150),
        [0.746331906460150]    = f(0.746331906460150),
        [-0.839116971822218]   = f(-0.839116971822218),
        [0.839116971822218]    = f(0.839116971822218),
        [-0.912234428251326]   = f(-0.912234428251326),
        [0.912234428251326]    = f(0.912234428251326),
        [-0.963971927277913]   = f(-0.963971927277913),
        [0.963971927277913]    = f(0.963971927277913),
        [-0.9931285991850949]  = f(-0.9931285991850949),
        [0.9931285991850949]   = f(0.9931285991850949),
    }
end

local LegendrePolynomials = {
    Polynomial.new({1}),
    Polynomial.new({0, 1}),
    Polynomial.new({-1, 0, 3}) / 2,
    Polynomial.new({0, -3, 0, 5}) / 2,
    Polynomial.new({3, 0, -30, 0, 35}) / 8,
    Polynomial.new({0, 15, 0, -70, 0, 63}) / 8,
    Polynomial.new({-5, 0, 105, 0, -315, 0, 231}) / 16,
    Polynomial.new({0, -35, 0, 315, 0, -693, 0, 429}) / 16,
    Polynomial.new({35, 0, -1260, 0, 6930, 0, -12012, 0, 6435}) / 128,
    Polynomial.new({0, 315, 0, -4620, 0, 18018, 0, -25740, 0, 12155}) / 128,
    Polynomial.new({-63, 0, 3465, 0, -30030, 0, 90090, 0, -109395, 0, 46189}) / 256
}

-- Normalize
for i, poly in ipairs(LegendrePolynomials) do
    LegendrePolynomials[i] = poly / math.sqrt(2 / (2 * i - 1))
end

local LegendrePrecomputes = {}
for i, poly in ipairs(LegendrePolynomials) do
    LegendrePrecomputes[i] = Precompute(poly)
end

local function ApproximateSpline(spline)
    local splinePrecompute = Precompute(function(x)
        return spline:ReparameterizeHybrid(x/2 + 0.5)
    end)
    local approximation = Polynomial.new()
    for i, poly in ipairs(LegendrePolynomials) do
        local polyPrecompute = LegendrePrecomputes[i]
        local inner = GaussLegendre.Twenty(function(x)
            return polyPrecompute[x] * splinePrecompute[x]
        end, -1, 1)
        approximation += poly * inner
    end
    return approximation
end

return ApproximateSpline

[ecurtiss]

I tried another approach, this time piecewise linearly interpolating the parameterization function. This gives passable results at significantly faster speeds. The accuracy can be tuned by varying the number of intervals of the interpolating function. In the future, the approximation can be improved by using more intervals in areas of higher curvature and less intervals in areas of lower curvature. This would likely require a binary search to find the correct interval, as they would not be evenly spaced. Alternatively, a quadratic approximation with uniform intervals might work.

[ecurtiss]

Resolved by v0.5.0.