Error with smoothing spline for 3D points

[Kevin-Mattheus-Moerman]

@jipolanco I tried testing the new smoothing splines using the below. I essentially create the points for a noisy circle and tried smoothing it. The below code works if V is instead a scalar quantity, but if I let V be a 3D point I get the error. Did you perhaps only implement the smoothing for 1D interpolation? Any help would be appreciated. Thanks.

using Comodo
using Comodo.GeometryBasics
using Comodo.GLMakie
using Comodo.BSplineKit

n = 25
t = collect(range(0.0,(2*pi)-(2*pi)/n,n))
r = 1.0 .+ 0.5.*rand(n)

V = Vector{Point{3,Float64}}(undef,n)
for (i,tt) in enumerate(t)
    V[i] = Point{3,Float64}(r[i]*cos(tt),r[i]*sin(tt),0.0)
end
 
λ = 0.1  # smoothing parameter
# S_fit_natural = BSplineKit.fit(t, V, λ)  # for comparison, compute a natural spline (no implied periodicity)
S_fit_periodic = BSplineKit.fit(t, V, λ, BSplineKit.Periodic(2*pi))

m = 200
t_fit = range(0.0,(2*pi)-(2*pi)/m,m)
V_fit = S_fit_periodic.(t_fit)

# Visualization
fig = Figure(size = (800,800))
ax = Axis3(fig[1, 1],aspect = :data)

hp1 = scatter!(ax, V,markersize=15,color=:red)
hp2 = lines!(ax,V_fit,linewidth=5,color=:black)

Legend(fig[1, 2],[hp1,hp2],["Raw","Smoothing spline"])

fig

ERROR: MethodError: no method matching +(::Float64, ::Point{3, Float64})
For element-wise addition, use broadcasting with dot syntax: scalar .+ array
The function `+` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  +(::Any, ::Any, ::Any, ::Any...)
   @ Base operators.jl:596
  +(::ChainRulesCore.ZeroTangent, ::Any)
   @ ChainRulesCore ~/.julia/packages/ChainRulesCore/U6wNx/src/tangent_arithmetic.jl:99
  +(::Any, ::ChainRulesCore.NotImplemented)
   @ ChainRulesCore ~/.julia/packages/ChainRulesCore/U6wNx/src/tangent_arithmetic.jl:25
  ...

Stacktrace:
  [1] _At_or_Ac_mul_B!(tfun::typeof(transpose), C::Vector{…}, A::SparseArrays.SparseMatrixCSC{…}, B::Vector{…}, α::Bool, β::Bool)
    @ SparseArrays ~/.julia/juliaup/julia-1.11.3+0.x64.linux.gnu/share/julia/stdlib/v1.11/SparseArrays/src/linalg.jl:111
  [2] spdensemul!
    @ ~/.julia/juliaup/julia-1.11.3+0.x64.linux.gnu/share/julia/stdlib/v1.11/SparseArrays/src/linalg.jl:61 [inlined]
  [3] generic_matvecmul!
    @ ~/.julia/juliaup/julia-1.11.3+0.x64.linux.gnu/share/julia/stdlib/v1.11/SparseArrays/src/linalg.jl:54 [inlined]
  [4] _mul!
    @ ~/.julia/juliaup/julia-1.11.3+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:73 [inlined]
  [5] mul!
    @ ~/.julia/juliaup/julia-1.11.3+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:70 [inlined]
  [6] mul!
    @ ~/.julia/juliaup/julia-1.11.3+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:253 [inlined]
  [7] fit(xs::Vector{…}, ys::Vector{…}, λ::Float64, order_in::BSplineOrder{…}, bc::Periodic{…}; weights::Nothing)
    @ BSplineKit.SplineInterpolations ~/.julia/packages/BSplineKit/Sxz1p/src/SplineInterpolations/smoothing.jl:178
  [8] fit
    @ ~/.julia/packages/BSplineKit/Sxz1p/src/SplineInterpolations/smoothing.jl:116 [inlined]
  [9] fit(x::Vector{Float64}, y::Vector{Point{3, Float64}}, λ::Float64, bc::Periodic{Float64})
    @ BSplineKit.SplineInterpolations ~/.julia/packages/BSplineKit/Sxz1p/src/SplineInterpolations/smoothing.jl:186
 [10] top-level scope
    @ ~/DATA/Julia/PROJECTS/Comodo_dev/dev/Comodo/examples/demo_smoothing_spline.jl:20
Some type information was truncated. Use `show(err)` to see complete types.

[Kevin-Mattheus-Moerman]

It works if I do it per coordinate like this:

n = 250
t = collect(range(0.0,(2*pi)-(2*pi)/n,n))
r = 1.0 .+ 0.5.*rand(n)

V = Vector{Point{3,Float64}}(undef,n)
for (i,tt) in enumerate(t)
    V[i] = Point{3,Float64}(r[i]*cos(tt),r[i]*sin(tt),0.0)
end

λ = 1  # smoothing parameter

m = 200
t_fit = range(0.0,(2*pi)-(2*pi)/m,m)
X = [v[1] for v in V]
Y = [v[2] for v in V]
S_fit_periodic_X = BSplineKit.fit(t, X, λ, BSplineKit.Periodic(2*pi))
S_fit_periodic_Y = BSplineKit.fit(t, Y, λ, BSplineKit.Periodic(2*pi))

V_fit = Vector{Point{3,Float64}}(undef,m)
for (i,tt) in enumerate(t_fit)
    V_fit[i] = Point{3,Float64}(S_fit_periodic_X.(tt),S_fit_periodic_Y.(tt),0.0)
end

I also have a separate question on the smoothing for these periodic functions. I tried to set $\lambda$ rather high e.g. $\lambda=20$. The hope was that it would approximate a circle, however it seems to shrink a lot as $\lambda$ is increased.

[jipolanco]

Did you perhaps only implement the smoothing for 1D interpolation?

Oops, thanks for reporting this. That's right, I didn't think about parametric splines. This is easy to fix, I'll tag a new version soon.

[jipolanco]

I also have a separate question on the smoothing for these periodic functions. I tried to set λ rather high e.g. λ = 20 . The hope was that it would approximate a circle, however it seems to shrink a lot as λ is increased.

That's an interesting question! I hadn't thought about it either. It seems like, for large $\lambda$, the optimal solution tends to the constant solution $x(t) = x_0$, $y(t) = y_0$. This kind of makes sense as it indeed minimises the curve derivatives... It can also be understood if one separately considers the smoothing of the $(t_i, x_i)$ or the $(t_i, y_i)$ data. In both cases, the result should tend to a straight line, which moreover needs to have zero slope to be periodic.

This is what I obtain with $\lambda = 20$:

I guess the issue is that separately considering each Cartesian coordinate is not the best way of dealing with parametric curves. Or at least, one should not minimise the derivatives of each component but something else. So the question is what's the relevant objective function for smoothing parametric curves.

[Kevin-Mattheus-Moerman]

@jipolanco Could the objective function be defined using curvature? E.g. the derivative of curvature.