Arc length for a spline

General Usage
, ,
21

You can, but splatting a large array is problematic too.

I’ve been meaning to add an API for this so that you at least don’t need to call the undocumented QuadGK.Segment constructor, e.g. a new method for alloc_segbuf that lets you pass an array of points. There’s already a documented API quadgk(f, knots; kws...) (no splatting), see below.

22

I think you already did. Looks like quadgk(f, knots; kws...) (no splatting) should work too: QuadGK.jl/src/api.jl at master · JuliaMath/QuadGK.jl · GitHub

You called the helper to_segbuf rather than adding a method to alloc_segbuf, but since users can pass the array to quadgk directly that doesn’t seem important.

1 Like
23

Doh, I forgot! It’s even documented.

2 Likes
24

I guess even more improvements might come once that PR gets cleared and merged (thank you for that!), but here are the corrected benchmarks:

julia> @btime method1($spl, extrema(t)...)
  5.940 ms (60 allocations: 9.16 MiB)
70.44577023545651

julia> @btime method2($spl, extrema(t)...)
  1.424 ms (51046 allocations: 1.51 MiB)
70.44577026541424

julia> @btime method2a($spl)
  134.311 μs (4983 allocations: 151.69 KiB)
70.44577024751408

julia> @btime method3($spl, $(gauss(15))...)
  83.455 μs (3354 allocations: 101.50 KiB)
70.44577024750694

julia> @btime method4($spl, order=4, rtol=1e-5)
  69.240 μs (2467 allocations: 75.44 KiB)
70.44577023146906

These have been compared so the relative accuracy is ~1e-10 (and I tried to fix method4 following your discussion, hope I didn’t mess it up):

function method1(spl, t1, t2)
    tl = LinRange(t1, t2, 100_000)
    xy = spl(tl)
    Δxy = xy[:,2:end] - xy[:,1:end-1]
    sum(norm.(eachcol(Δxy)))
end
function method2(spl, t1, t2)
    res, _ = quadgk(t -> norm(derivative(spl, t)), t1, t2, rtol = 1e-8)
    res
end
function method2a(spl)
    knots = Dierckx.get_knots(spl)
    s = 0.0
    for (k1, k2) in zip(knots[1:end-1], knots[2:end])
        res, _ = quadgk(t -> norm(derivative(spl, t)), k1, k2)
        s += res
    end
    return s
end
function method3(spl, x_gauss, w_gauss)
    knots = Dierckx.get_knots(spl)
    integral = sum(eachindex(knots)[2:end]) do i
        t1, t2 = knots[i-1], knots[i]
        scale = (t2 - t1) / 2
        sum(zip(x_gauss, w_gauss)) do ((xg, wg))
            t = (xg + 1) * scale + t1 # map from (-1,1) to (t1,t2)
            norm(derivative(spl, t)) * wg
        end * scale
    end
end
function method4(spl; kws...)
    knots = Dierckx.get_knots(spl)
    s, _ = quadgk(t -> norm(derivative(spl, t)), knots; kws...)
    return s
end



row = df[5, :]
spl = row.spl
t = row.t
s, _ = quadgk(t -> norm(derivative(spl, t)), extrema(t)..., rtol = 1e-9)
rel(x) = (x - s)/s

rel(method1(spl, extrema(t)...))
rel(method2(spl, extrema(t)...))
rel(method2a(spl))
rel(method3(spl, gauss(15)...))
rel(method4(spl, order=4, rtol=1e-5))
25

For future reference, if you need a higher-level interface for integrals over geometries, including BezierCurve and ParametrizedCurve, consider the MeshIntegrals.jl package developed and maintained by @mike.ingold , @JoshuaLampert and @kylebeggs

It relies on QuadGK.jl, FastGaussQuadrature.jl and HCubature.jl as backends, and can facilitate the lives of developers who need to perform integration over curves, surfaces or volumes.

4 Likes
26

I think you need more parens @btime method3($spl, $(gauss(15))...), otherwise you are only interpolating the function name.

27

Oops, you’re right, I edited it now.

28

Just to address your excellent point: I need a smoothing spline, hence the use of Dierckx.jl (and not Interpolations.jl for instance).

In other news, I’ve updated the winning implementation, “method4”, to allow for bounded arc-lengths:

function arclength(spl, t1, t2; kws...)
    knots = get_knots(spl)
    filter!(t -> t1 < t < t2, knots)
    pushfirst!(knots, t1)
    push!(knots, t2)
    s, _ = quadgk(t -> norm(derivative(spl, t)), knots; kws...)
    return s
end