Question: How to obtain the spline curve length for a given domain?

[Kevin-Mattheus-Moerman]

@jipolanco thanks for developing this package. I've only started testing it but it looks great.

TLDR: Is it possible to easily obtain the length of a spline curves for a given range? So ideally I'd set up an interpolator a bit like S = interpolate..., and then ask what is the length of S between x1 and x2. Do you have something implemented for that? I suppose we need something like:

$$L_{a-b} = \int_a^b |f'(t)| dt $$

Below is a longer explanation of what I am trying to use this for.

I am trying to use interpolation to evenly sample a non-evenly sampled ND input curve. Below I have an example of a first attempt. Here I parameterise the raw curve using the "along curve distance". I then evenly sample the distance from 0 to the final end distance (length) of the curve and feed that to one of your spline based interpolation methods. The below picture shows this works reasonably well. However, with this setup the distances are computed for a linear fit to the curve, i.e. straight line segments connecting the points. It would be better to fit a spline, then get the length of the spline and resample it evenly in this length/distance space.

Any suggestions? Thanks!

using GeometryBasics
using GLMakie
using LinearAlgebra
using BSplineKit

# Define raw curve
t = range(0,2.0*pi,20)
V = [GeometryBasics.Point{3, Float64}(tt,3.0*sin(tt),0.0) for tt ∈ t]

# Evenly sample curve
n = 50; # Number of points for spline 

L = pushfirst!(cumsum(norm.(diff(V))),0.0) # Along curve distance from start
l = range(0.0,maximum(L),n) #Even along curve distance 

S = interpolate(L, V, BSplineOrder(4),Natural()) # Create interpolator

Vi = S.(l) # Evaluate interpolator at even distance increments

# Visualization
fig = Figure(size = (1200,800))

ax1 = Axis3(fig[1, 1],aspect = :data)
hp1 = scatter!(ax1, V,markersize=25,color=:red)
hp2 = lines!(ax1, V,linewidth=3,color=:red)

ax2 = Axis3(fig[1, 2],aspect = :data)
hp1 = scatter!(ax2, Vi,markersize=20,color=:black)
hp2 = lines!(ax2, Vi,linewidth=3,color=:black)
hp2 = lines!(ax2, V,linewidth=2,color=:red)

fig

[jipolanco]

Hi, I completely understand what you're trying to do as I have been dealing with similar issues recently (parametrisation of 3D curves, and how to achieve something close to arc length parametrisation). Thanks for providing a nice minimal example!

In BSplineKit is no built-in way of doing this. To get the total length of the spline $S(t)$ the easiest would be to use numerical integration to estimate the length:

$$L = \int_a^b |S'(t)| \, \mathrm{d}t$$

where $S'(t)$ is the derivative of the spline with respect to its parameter $t$ (approximately the arc length here). In your example, $a$ and $b$ would be respectively L[begin] and L[end], which you can also get via a = first(knots(S)) and b = last(knots(S)). One can use something like QuadGK.jl to perform the numerical integration:

using QuadGK

Sp = Derivative() * S  # spline derivative

a, b = first(knots(S)), last(knots(S))
L_total, err = quadgk(a, b; rtol = 1e-8) do t
    norm(Sp(t))
end

If I run this on your example, I get L_total = 13.973590503382189 which, as expected, is slightly larger than your initial estimation of the total length using the straight segment approximation (L[end] = 13.925687625112573).

If what you want is to reparametrise the curve based on the actual length of each segment, you could also build a new L vector obtained from integrating the length segment-by-segment:

L_new = fill!(similar(L), 0)
for i ∈ eachindex(L_new)[2:end]
    # Compute length of segment [i, i + 1].
    ta = L[i - 1]
    tb = L[i]
    segment_length, err = quadgk(ta, tb; rtol = 1e-8) do t
        norm(Sp(t))  # integrate |S'(t)| in segment [i, i + 1]
    end
    L_new[i] = L_new[i - 1] + segment_length
end

L_new - L  # note: L_new[i] > L[i] for all i since the actual lengths are larger than approximating by straight segments

# We can now interpolate a new spline using the new lengths
S = interpolate(L_new, V, BSplineOrder(4), Natural())

If you really care about precision, you could do this iteratively multiple times until L_new converges to some constant.

[Kevin-Mattheus-Moerman]

@jipolanco wow! This is fantastic. Yeah I did think of an iterative approach too (for the linear distance metric, would be interesting to see if both always properly converge or if for instance the linear one could jump back and forth and not quite converge all the time). I'll test what you proposed now. I'll shortly publish a basic package featuring this as well. I can copy the above and work it in, or would you like to officially submit the above as a pull request to my project? That way you are down as a contributor so it is clear who gets the credit for the above (and perhaps co-author if more follows e.g. on a paper on the package)? Let me know what you prefer.

Thanks a lot for your help, and especially for the speed at which you responded! :partying_face:

[jipolanco]

That sounds interesting! I'd be happy to contribute to your package (and to a possible paper), so let me know when the package is ready to accept contributions :)

[Kevin-Mattheus-Moerman]

Below is the full prototype code implementing your solution.

function evenly_sample(V,n)    
    m = length(V)
    T = curve_length(V) # Initialise as along curve (multi-linear) distance
    T ./= last(T) # Normalise
    S = interpolate(T, V, BSplineOrder(4),Natural()) # Create interpolator
    dS = Derivative() * S  # spline derivative
    L = zeros(Float64,m) # Initialise spline length vector
    for i ∈ 2:m
        # Compute length of segment [i-1, i]
        segment_length, _ = quadgk(T[i-1], T[i]; rtol = 1e-8) do t
            norm(dS(t))  # integrate |S'(t)| in segment [i, i + 1]
        end        
        L[i] = L[i - 1] + segment_length
    end    
    L ./= last(L) # Normalise to 0-1 range
    l = range(0.0,1.0,n) #Even allong curve distance 
    S = interpolate(L, V, BSplineOrder(4),Natural()) # Create interpolator
    return S.(l),S # Evaluate interpolator at even distance increments
end

Note sure if iterating this will be a solution to get a pure homogeneous point spacing. I think this gives us a homogeneous points spacing along the spline. Hence if the spline has significant and varying curvature between the points, e.g. when a relatively large point spacing is chosen, then the distance may not be constant for the points obtained, but still correct from "the spline point of view".

Like in the below zoomed in region, if the spacing causes us to "skip" a bit with curvature, then the spacing of the result (black) deviates from the distance between the points according to the full spline (red), hence it appears to be wrong but... I think this is correct and just something to be aware of.

Having a truly even spacing in the black dotted result requires a different approach I think.

[Kevin-Mattheus-Moerman]

Question was answered so can be closed.

[Kevin-Mattheus-Moerman]

@jipolanco My project now lives here: https://github.com/COMODO-research/Comodo.jl. The current prototype of what we spoke about is here: https://github.com/COMODO-research/Comodo.jl/blob/main/wip/wip_evenly_sample_curve.jl. If you want to add an improved version to the functions.jl file feel free.

Note, I am relatively new to Julia, so you may see a lot of "spinach in my teeth" over at that brand new project :) (it started as a Julia version of GIBBON).

[jipolanco]

Your new project looks seriously impressive (and GIBBON as well)!

I'll be happy to contribute with the spline part. I'll submit a PR soon with something equivalent to your evenly_sample function.