Interpolation of closed parametric curves

[Kevin-Mattheus-Moerman]

@jipolanco is there a way to use matched start and end slopes? This is useful for interpolating curves that should be deemed closed.

For instance, in this example the left might be a circle sample with only 5 points. It would be nice to be able to interpolate it more densely but such that the fact that the curve is closed is taken into account.

Of course a naive try would be to add a matched end=start point, however the trained eye would spot that the slopes do not match.

So that leads me back to my original question, can a boundary condition be added to match the start and end slopes? Or is there a better way to say that the curve is closed (I suppose it can be viewed as periodic?).

Thanks.

[jipolanco]

Hi @Kevin-Mattheus-Moerman, as you suggest, this is something which should be allowed by the Periodic boundary conditions (docs here).

However, I just realised that interpolation of parametric curves (using SVectors for the input points) currently doesn't work when using periodic conditions. I'll soon push a fix for this.

After the fix is released, one will be able to do the following:

using BSplineKit
using StaticArrays
using CairoMakie

θs = range(0, 2π; length = 6)[1:5]  # input angles in [0, 2π)
points = [SVector(cos(θ), sin(θ)) for θ ∈ θs]
S = interpolate(θs, copy(points), BSplineOrder(4), Periodic(2π))  # only cubic splines work for now

θs_interp = range(0, 2π; length = 100)
S_interp = S.(θs_interp)

fig = Figure()
ax = Axis(fig[1, 1]; aspect = DataAspect())
scatterlines!(ax, points; color = :red)
scatterlines!(ax, S_interp; color = :black)
fig

[jipolanco]

This should be fixed in the upcoming 0.17.5 version (for cubic splines only). Let me know if it works for you.

I'm reopening this issue since orders other than cubic (BSplineOrder(4)) still don't work.

[Kevin-Mattheus-Moerman]

@jipolanco thanks. Let me know when the release is ready. I'll test it out. The context is evenly_sample in COMODO.

Can the natural option be combined with periodic?

This is the untested prototype which I'll start testing once this functionality is available, let me know if you have any tips:

function evenly_sample(V::Vector{Point{ND,TV}}, n::Int64; rtol = 1e-8, niter = 1,spline_order=4, close_loop=false) where ND where TV<:Real
    m = length(V)
    LL = curve_length(V; close_loop=close_loop) # Initialise as along curve (multi-linear) distance
    LL ./= last(LL) # Normalise    
    if close_loop == true
        S = BSplineKit.interpolate(LL[1:end-1], V, BSplineOrder(spline_order), BSplineKit.Periodic(1.0)) # Create interpolator
    else
        S = BSplineKit.interpolate(LL, V, BSplineOrder(spline_order), BSplineKit.Natural()) # Create interpolator
    end
    for _ in 1:niter
        dS = Derivative() * S  # spline derivative
        L = similar(LL) # Initialise spline length vector
        L[1] = 0
        for i in 2:m
            # Compute length of segment [i-1, i]
            segment_length, _ = quadgk(LL[i-1], LL[i]; rtol) 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
        if close_loop == true
            S = BSplineKit.interpolate(L[1:end-1], V, BSplineOrder(spline_order), BSplineKit.Periodic(1.0)) # Create interpolator
        else
            S = BSplineKit.interpolate(L, V, BSplineOrder(spline_order), BSplineKit.Natural()) # Create interpolator
        end
    end
    # Even range for curve distance 
    if close_loop == true
        l = range(0.0, 1.0 - 1.0/n, n) 
    else
        l = range(0.0, 1.0, n) 
    end
    return S.(l), S # Evaluate interpolator at even distance increments
end

[jipolanco]

It's available in the latest release, you should be able to try it out now. Let me know if you have any issues.

Can the natural option be combined with periodic?

No, and that wouldn't make much sense either. Mathematically they are different and incompatible conditions. You really want to use Periodic if you want the curve and its derivatives to be continuous at the "endpoint".

[Kevin-Mattheus-Moerman]

@jipolanco thanks, testing it now. I've noticed the need to copy(points) above in your example. That wasn't needed for the Natural constraint. Is that intentional? Would be nice to avoid the need to copy (I see that without it will alter the points).

[Kevin-Mattheus-Moerman]

@jipolanco Cool it works, see implementation of evenly_sample (and related demos) in Comodo. If you had time to check to see if I did that right that would be great.

Let me know once other orders work, the linear one (order 2 I believe) is particularly relevant.

Thanks for fixing this.

[jipolanco]

I've noticed the need to copy(points) above in your example. That wasn't needed for the Natural constraint. Is that intentional? Would be nice to avoid the need to copy (I see that without it will alter the points).

The current implementation of periodic cubic spline interpolations reuse the input data to avoid allocations (this is mentioned in the docstring). But I agree that this is not very convenient, and it probably makes sense to change this...

Cool it works, see implementation of evenly_sample (and related demos) in Comodo. If you had time to check to see if I did that right that would be great.

I quickly looked at it and it looks good to me!

Let me know once other orders work, the linear one (order 2 I believe) is particularly relevant.

Order 2 in particular should be actually quite easy to fix, I'll see what I can do.

[jipolanco]

Linear splines should be fixed in the upcoming v0.17.6.

[Kevin-Mattheus-Moerman]

@jipolanco thanks. I've just updated the whole thing again. Rather than normalising the parameterisation e.g. so that L goes from 0-1, I let is be un-normalised so it goes from 0 - curve length. This way it was easier to implement a new feature I added, i.e. must_points. These are points the spline must visit.

For the images below in the middle the curve is not closed and there are no must points. For the image on the right I ask for a closed curve and also have the green points as must points.