Optimizing Spline:SolveLength()

[ecurtiss]

Currently, little effort has been put into making the numerical integration performed in Spline:SolveLength() faster and more accurate. The first iteration simply sampled points and summed the distances between adjacent points.

local sum = 0
local prev = self:SolvePosition(0)
for i = 1, N do
    local pos = self:SolvePosition(i / N)
    sum += (pos - prev).Magnitude
    prev = pos
end

return sum

This changed in v0.4.0 when we switched to doing five iterations of Simpson's 3/8ths rule.

return Integrate.Simp38Comp(function(x)
    return self:SolveVelocity(x).Magnitude
end, a or 0, b or 1, 5)

Today I compared Simpson's rule with Gaussian quadrature. I did this by creating N random splines, and for each spline, I performed M iterations of Simpson's rule with 1 <= i <= 64 intervals and M iterations of Gaussian quadrature with 2 <= i <= 63 basis functions. Then I averaged the results for each spline, and finally averaged the results over all splines. Try it for yourself: SolveLengthComparison.zip

In the following plots, the x-axis is the amount of intervals/basis functions used. Note that the bottom plot uses a log scale.

We see from this that Gaussian quadrature is always faster and more accurate than Simpson's rule. Furthermore, Gaussian quadrature performs no better after about 20 functions. In order to compromise between speed and accuracy, I think 10 basis functions should be used.

[ecurtiss]

Here is a benchmark using a 5-interval composite Simpson's 3/8th rule vs a 10-point Gauss rule.

local CatRom = require(game.ReplicatedStorage.CatRom)
local CatRom2 = require(game.ReplicatedStorage.CatRom2)

local N = 1e3

return {
    ParameterGenerator = function()
        local points = table.create(4)
        for i = 1, 4 do
            points[i] = Vector3.new(math.random(), math.random(), math.random()) * 10
        end
        return points, math.random(), math.random()
    end,

    Functions = {
        Gauss = function(Profiler, points, a, b)
            local spline = CatRom.new(points).splines[2]
            for i = 1, N do
                spline:SolveLength(a, b)
            end
        end,

        Simp = function(Profiler, points, a, b)
            local spline = CatRom2.new(points).splines[2]
            for i = 1, N do
                spline:SolveLength(a, b)
            end
        end,
    }
}

[ecurtiss]

Resolved by v0.5.0.