Unexpected squiggles in the spline

[Kevin-Mattheus-Moerman]

@jipolanco

I tried to make a horribly irregularly sampled curve to test the even sampling approach discussed over at https://github.com/jipolanco/BSplineKit.jl/issues/85. Using the below code where I created some sorted unique random points from 0 up to 2*pi I unfortunately see unexpected squiggles in the spline. Am I doing something wrong?

using GeometryBasics
using GLMakie
using LinearAlgebra
using BSplineKit
using Random

Random.seed!(16)

# Define raw curve
np = 50
t = sort(unique(rand(np)*2*pi))  # t = range(0,2.0*pi,np)
V = [GeometryBasics.Point{3, Float64}(tt,3*sin(tt),0.0) for tt ∈ t]

T = range(0.0,1.0,length(V)) # Parameterise with unit total length
S = interpolate(T, V, BSplineOrder(4),Natural()) # Create interpolator

TT = range(0.0,1.0,200) 
Vi = S.(TT)

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

ax1 = Axis(fig[1, 1])
hp1 = scatter!(ax1, V,markersize=15,color=:red)
hp2 = lines!(ax1, V,linewidth=3,color=:red)

ax2 = Axis(fig[1, 2])
hp1 = scatter!(ax2, Vi,markersize=20,color=:black)
hp2 = lines!(ax2, Vi,linewidth=2,color=:red)

fig

[jipolanco]

It does look horrible! I think the problem is that, unlike in #85, here you're not parametrising the curve by the (approximate) arc length, which can lead to this kind of artefacts. I mean, this line:

T = range(0.0,1.0,length(V)) # Parameterise with unit total length

creates a uniform parametrisation, but your input points are strongly non-uniformly spaced.

If instead, following your initial example in #85, I do the following:

T = pushfirst!(cumsum(norm.(diff(V))),0.0) # Along curve distance from start
T ./= last(T)  # normalise to [0, 1] range

then I get something which looks nice:

I hope this answers the question, I'm not sure I understood correctly what you're trying to do.

[Kevin-Mattheus-Moerman]

Right! Of course, my mistake. After I include that suggestion it all works. I was only trying to pass it something tricky and very non-uniform.

Updated solution here:

function evenly_sample(V,n)    
    m = length(V)
    T = pushfirst!(cumsum(norm.(diff(V))),0.0) # Along curve distance from start
    T ./= last(T)  # normalise to [0, 1] range
    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

Thanks for your help!