From 70e21d8b19ede01f9f62070229b13a27f480b825 Mon Sep 17 00:00:00 2001 From: Kevin-Mattheus-Moerman Date: Tue, 9 Jul 2024 16:41:44 +0100 Subject: [PATCH] Updates to evenly sample functions and demos to allow for must points --- examples/demo_curve_length.jl | 26 ++++++++ examples/demo_evenly_sample.jl | 24 +++++-- examples/demo_evenly_space.jl | 31 +++++++-- examples/demo_filletcurve_extrude.jl | 4 +- examples/demo_loftlinear.jl | 2 +- examples/demo_regiontrimesh.jl | 8 +-- examples/demo_revolvecurve.jl | 6 +- examples/demo_sweeploft.jl | 26 ++++---- src/functions.jl | 99 +++++++++++++++++++++------- test/runtests.jl | 65 ++++++++++++------ 10 files changed, 213 insertions(+), 78 deletions(-) create mode 100644 examples/demo_curve_length.jl diff --git a/examples/demo_curve_length.jl b/examples/demo_curve_length.jl new file mode 100644 index 0000000..893fb41 --- /dev/null +++ b/examples/demo_curve_length.jl @@ -0,0 +1,26 @@ +using Comodo +using GeometryBasics +using GLMakie +using BSplineKit +using QuadGK: quadgk # For numerical integration +using LinearAlgebra + +testCase = 2 + +# Define input curve +if testCase == 1 + V = Point{3,Float64}[[0,0,0],[1.0,0,0],[1,1,0],[0,1,0]] + close_loop = false + c = 3 +elseif testCase == 2 + n = 250 + r = 2*pi + V = [Point{3,Float64}(r*cos(t),r*sin(t),0.0) for t in range(0,2*pi,n)] + close_loop = true + c = 2*pi*r +end + +L = curve_length(V; close_loop = close_loop) +d = last(L) + +println("Max. length: " * string(d) * ", theoretical/true: " * string(c)) diff --git a/examples/demo_evenly_sample.jl b/examples/demo_evenly_sample.jl index 6f47fb7..165b25f 100644 --- a/examples/demo_evenly_sample.jl +++ b/examples/demo_evenly_sample.jl @@ -4,24 +4,36 @@ using GLMakie # This demo shows how to evenly sample a curve using spline interpolations. -testCase = 2 +testCase = 3 # Define input curve if testCase == 1 n = 5 t = range(0,2π-2π/n,n) V = [GeometryBasics.Point{3, Float64}(3.0*cos(tt),3.0*sin(tt),0.0) for tt ∈ t] + + # Evenly sample curve + n = 50; # Number of points for spline + Vi1 = evenly_sample(V, n; spline_order=4, close_loop = false) # Returns points and spline interpolation object + Vi2 = evenly_sample(V, n; close_loop = true) # Returns points and spline interpolation object elseif testCase == 2 n = 4*7+4 rFun(t) = sin(4*t)+2 V = circlepoints(rFun,n) -end -# Evenly sample curve -n = 50; # Number of points for spline -Vi1, S = evenly_sample(V, n; niter = 10, close_loop = false) # Returns points and spline interpolation object + # Evenly sample curve + n = 75; # Number of points for spline + Vi1 = evenly_sample(V, n; spline_order=4, close_loop = false, niter = 10) # Returns points and spline interpolation object + Vi2 = evenly_sample(V, n; spline_order=4, close_loop = true, niter = 10) # Returns points and spline interpolation object +elseif testCase == 3 + V = batman(25; symmetric = true, dir=:acw) + + # Evenly sample curve + n = 50; # Number of points for spline + Vi1 = evenly_sample(V, n; niter = 10, spline_order=4, close_loop = false) # Returns points and spline interpolation object + Vi2 = evenly_sample(V, n; niter = 10, spline_order=4, close_loop = true) # Returns points and spline interpolation object +end -Vi2, S = evenly_sample(V, n; niter = 10, close_loop = true) # Returns points and spline interpolation object # Visualization diff --git a/examples/demo_evenly_space.jl b/examples/demo_evenly_space.jl index 613bc78..8425ce0 100644 --- a/examples/demo_evenly_space.jl +++ b/examples/demo_evenly_space.jl @@ -4,31 +4,45 @@ using GLMakie # This demo shows how to evenly space a curve using spline interpolations. -testCase = 2 +testCase = 1 # Define input curve if testCase == 1 n = 5 t = range(0,2π-2π/n,n) V = [GeometryBasics.Point{3, Float64}(3.0*cos(tt),3.0*sin(tt),0.0) for tt ∈ t] + + must_points = [1,2,3,4,5] + + # Evenly sample curve + Vi1 = evenly_space(V, 0.65; close_loop = false, spline_order = 4) # Returns points and spline interpolation object + Vi2 = evenly_space(V, 0.65; close_loop = true, spline_order = 4,must_points = must_points) # Returns points and spline interpolation object elseif testCase == 2 n = 4*7+4 rFun(t) = sin(4*t)+2 V = circlepoints(rFun,n) -end + must_points = [3,7,11,15,19,23,27,31] + + # Evenly sample curve + Vi1 = evenly_space(V, 0.5; close_loop = false, niter=10, spline_order = 4) # Returns points and spline interpolation object + Vi2 = evenly_space(V, 0.3; close_loop = true, niter=10, spline_order = 4,must_points = must_points) # Returns points and spline interpolation object +elseif testCase == 3 + V = batman(25; symmetric = true, dir=:acw) -# Evenly sample curve -n = 25; # Number of points for spline -Vi1 = evenly_space(V, 0.3; niter = 10, close_loop = false) # Returns points and spline interpolation object + must_points = [7,10,13,15,18,21] -Vi2 = evenly_space(V, 0.3; niter = 10, close_loop = true) # Returns points and spline interpolation object + # Evenly sample curve + Vi1 = evenly_space(V, 0.15; close_loop = false, niter=10, spline_order = 4) # Returns points and spline interpolation object + Vi2 = evenly_space(V, 0.1; close_loop = true, niter=10, spline_order = 4, must_points = must_points) # Returns points and spline interpolation object +end # Visualization - markersize1 = 10 markersize2 = 15 +markersize3 = 20 + linewidth = 2 fig = Figure(size = (1200,500)) @@ -45,6 +59,9 @@ hp4 = lines!(ax2, Vi1,linewidth=linewidth,color=:black) ax3 = Axis3(fig[1, 3],aspect = :data,azimuth=-pi/2,elevation=pi/2, title="evenly spaced, closed") hp1 = scatter!(ax3, V,markersize=markersize2,color=:red) +if !isnothing(must_points) + hp5 = scatter!(ax3, V[must_points],markersize=markersize3,color=:green) +end hp2 = lines!(ax3, V,linewidth=linewidth,color=:red) hp3 = scatter!(ax3, Vi2,markersize=markersize1,color=:black) hp4 = lines!(ax3, Vi2,linewidth=linewidth,color=:black) diff --git a/examples/demo_filletcurve_extrude.jl b/examples/demo_filletcurve_extrude.jl index 280782e..8cee3fd 100644 --- a/examples/demo_filletcurve_extrude.jl +++ b/examples/demo_filletcurve_extrude.jl @@ -31,7 +31,7 @@ rMax = 0.5 #collect(range(0.5,5,length(V))) n = 20 h =2.0 VC = filletcurve(V; rMax=rMax, constrain_method = :max, n=n, close_loop = close_loop, eps_level = 1e-6) -VC,_ = evenly_sample(VC,50; spline_order = 2) +VC = evenly_sample(VC,50; spline_order = 2) Fe,Ve = extrudecurve(VC; extent=h, direction=:both, close_loop=false,face_type=:quad) @@ -59,7 +59,7 @@ hSlider1 = Slider(fig[2, 1], range = stepRange1, startvalue = 0,linewidth=30) on(hSlider1.value) do stepIndex1 VC = filletcurve(V; rMax=stepIndex1, constrain_method = :max, n=n, close_loop = close_loop, eps_level = 1e-6) - VC,_ = evenly_sample(VC,53; spline_order = 2) + VC = evenly_sample(VC,53; spline_order = 2) indPlot = collect(1:length(VC)) if close_loop == true push!(indPlot,1) diff --git a/examples/demo_loftlinear.jl b/examples/demo_loftlinear.jl index 64b2011..7aacbeb 100644 --- a/examples/demo_loftlinear.jl +++ b/examples/demo_loftlinear.jl @@ -12,7 +12,7 @@ height = 2.5 rFun(t) = r + 0.5.* sin(4.0*t) V2 = circlepoints(rFun,nc;dir=:acw) V2 = [GeometryBasics.Point{3, Float64}(v[1],v[2],height) for v ∈ V2] -V2,_ = evenly_sample(V2,nc) +V2 = evenly_sample(V2,nc) ## Loft from curve 1 to curve 2 num_steps = 11 # Uneven works best for "tri" face type close_loop = true diff --git a/examples/demo_regiontrimesh.jl b/examples/demo_regiontrimesh.jl index 6b4997a..a33cbff 100644 --- a/examples/demo_regiontrimesh.jl +++ b/examples/demo_regiontrimesh.jl @@ -31,7 +31,7 @@ elseif testCase == 3 r = 1.0 rFun(t) = r + 0.5.*sin(6*t) V = circlepoints(rFun,n; dir=:acw) - V,_ = evenly_sample(V, n) + V = evenly_sample(V, n) pointSpacing = pointspacingmean(V) VT = (V,) @@ -41,17 +41,17 @@ elseif testCase == 4 rFun1(t) = 12.0 + 3.0.*sin(5*t) n1 = 120 V1 = circlepoints(rFun1,n1) - V1,_ = evenly_sample(V1,n1) + V1 = evenly_sample(V1,n1) rFun2(t) = 10.0 + 2.0.*sin(5*t) n2 = 100 V2 = circlepoints(rFun2,n2) - V2,_ = evenly_sample(V2,n2) + V2 = evenly_sample(V2,n2) rFun3(t) = 2.0 + 0.5 *sin(5*t) n3 = 50 Vp = circlepoints(rFun3,n3) - Vp,_ = evenly_sample(Vp,n3) + Vp = evenly_sample(Vp,n3) V3 = [Point{3,Float64}(v[1],v[2]+6,v[3]) for v in Vp] V4 = [Point{3,Float64}(v[1]-5,v[2]+1,v[3]) for v in Vp] diff --git a/examples/demo_revolvecurve.jl b/examples/demo_revolvecurve.jl index 3e82b69..2bbfd33 100644 --- a/examples/demo_revolvecurve.jl +++ b/examples/demo_revolvecurve.jl @@ -13,7 +13,7 @@ if testCase == 1 nc = 15 t = range(0,2*π,nc) Vc = [Point3{Float64}(2.0+tt,0.0,sin(tt)) for tt ∈ t] - Vc,_ = evenly_sample(Vc, nc) + Vc = evenly_sample(Vc, nc) n = Vec{3, Float64}(0.0,0.0,1.0) num_steps = 31 close_loop = false @@ -22,7 +22,7 @@ elseif testCase == 2 nc = 15 t = range(0,2*π,nc) Vc = [Point3{Float64}(2.0+tt,0.0,sin(tt)) for tt ∈ t] - Vc,_ = evenly_sample(Vc, nc) + Vc = evenly_sample(Vc, nc) n = normalizevector(Vec{3, Float64}(0.0,-1.0,1.0)) num_steps = 31 close_loop = false @@ -31,7 +31,7 @@ elseif testCase == 3 nc = 16 t = range(2.0*π-(2.0*π/nc),0,nc) Vc = [Point3{Float64}(2.0+cos(tt),0.0,sin(tt)) for tt ∈ t] - Vc,_ = evenly_sample(Vc, nc) + Vc = evenly_sample(Vc, nc) n = Vec{3, Float64}(0.0,0.0,1.0) num_steps = 25 close_loop = true diff --git a/examples/demo_sweeploft.jl b/examples/demo_sweeploft.jl index 62cbec3..5570f08 100644 --- a/examples/demo_sweeploft.jl +++ b/examples/demo_sweeploft.jl @@ -17,19 +17,19 @@ if testCase == 1 P[4 ] = GeometryBasics.Point{3, Float64}( 1.0, 1.0, 1.0) Vc = nbezier(P,nc) # Get Bezier fit points Vc = [vc.*10 for vc in Vc] - Vc,Sc = evenly_sample(Vc, nc) + Vc = evenly_sample(Vc, nc) # Define section curves np = 50 # Number of section points cFun(t) = 2.0 + 0.5.*sin(3*t) V1 = circlepoints(cFun,np; dir=:acw) - V1,_ = evenly_sample(V1, np) + V1 = evenly_sample(V1, np) Q = RotXYZ(0.0,0.5*π,0.0) # Define a rotation tensor using Euler angles V1 = [Q*v for v ∈ V1] # Rotate the coordinates cFun(t) = 2.0 + 0.5*sin(3*t) V2 = circlepoints(cFun,np; dir=:acw) - V2,_ = evenly_sample(V2, np) + V2 = evenly_sample(V2, np) V2 = [v2 .+ Vc[end] for v2 ∈ V2] elseif testCase == 2 # Define guide curve @@ -41,19 +41,19 @@ elseif testCase == 2 P[4 ] = GeometryBasics.Point{3, Float64}( 3.0, 0.0, 0.0) Vc = nbezier(P,nc) # Get Bezier fit points Vc = [vc.*10 for vc in Vc] - Vc,Sc = evenly_sample(Vc, nc) + Vc = evenly_sample(Vc, nc) # Define section curves np = 35 # Number of section points cFun(t) = 2.0 + 0.5.*sin(3*t) V1 = circlepoints(cFun,np; dir=:acw) - V1,_ = evenly_sample(V1, np) + V1 = evenly_sample(V1, np) Q = RotXYZ(0.0,0.5*π,0.0) # Define a rotation tensor using Euler angles V1 = [(Q*v) .+ Vc[1] for v ∈ V1] # Rotate the coordinates cFun(t) = 2.0 + 0.5*sin(3*t) V2 = circlepoints(cFun,np; dir=:acw) - V2,_ = evenly_sample(V2, np) + V2 = evenly_sample(V2, np) Q = RotXYZ(0.0,0.5*π,0.0) # Define a rotation tensor using Euler angles V2 = [(Q*v) .+ Vc[end] for v ∈ V2] elseif testCase == 3 @@ -68,7 +68,7 @@ elseif testCase == 3 # Section 1 cFun(t) = 1.5 + 0.5.*sin(3*t) V1 = circlepoints(cFun,np; dir=:cw) - V1,_ = evenly_sample(V1, np) + V1 = evenly_sample(V1, np) Q = RotXYZ(0.5*π,0.0,0.0) # Define a rotation tensor using Euler angles V1 = [Q*v for v ∈ V1] # Rotate the coordinates @@ -90,7 +90,7 @@ elseif testCase == 3 # Section 2 cFun(t) = 4 + 1.5*sin(5*t) V2 = circlepoints(cFun,np; dir=:cw) - V2,_ = evenly_sample(V2, np) + V2 = evenly_sample(V2, np) Q1 = RotXYZ(0.5*π,0.0,0.0) # Define a rotation tensor using Euler angles Q2 = RotXYZ(0.0,-0.25*π,0.0) # Define a rotation tensor using Euler angles Q = Q2*Q1 @@ -119,7 +119,7 @@ elseif testCase == 4 rc = rb/3 r = [rb + rc*sin(ff*tt) for tt in t] Vc = [GeometryBasics.Point{3, Float64}(r[i]*cos(t[i]),r[i]*sin(t[i]),t[i]+rc*cos(ff*t[i])) for i in eachindex(t)] - Vc,_ = evenly_sample(Vc, nc) + Vc = evenly_sample(Vc, nc) # Define section curves np = 150 # Number of section points @@ -127,7 +127,7 @@ elseif testCase == 4 # Section 1 cFun(t) = rc/3 + rc/5 .* sin(3*t) V1 = circlepoints(cFun,np; dir=:cw) - V1,_ = evenly_sample(V1, np) + V1 = evenly_sample(V1, np) Q = RotXYZ(0.5*π,0.0,0.0) # Define a rotation tensor using Euler angles V1 = [Q*v for v ∈ V1] # Rotate the coordinates @@ -149,7 +149,7 @@ elseif testCase == 4 # Section 2 cFun(t) = rc/4 + rc/6*sin(5*t) V2 = circlepoints(cFun,np; dir=:cw) - V2,_ = evenly_sample(V2, np) + V2 = evenly_sample(V2, np) Q1 = RotXYZ(0.5*π,0.0,0.0) # Define a rotation tensor using Euler angles Q2 = RotXYZ(0.0,-0.25*π,0.0) # Define a rotation tensor using Euler angles Q = Q2*Q1 @@ -181,7 +181,7 @@ elseif testCase ==5 # Section 1 cFun(t) = 3 + 0.25.*sin(3*t) V1 = circlepoints(cFun,np; dir=:cw) - V1,_ = evenly_sample(V1, np) + V1 = evenly_sample(V1, np) V1_ori = deepcopy(V1) # Ensure section is orthogonal to guide curve @@ -203,7 +203,7 @@ elseif testCase ==5 # Section 2 cFun(t) = 5 + 0.5*sin(5*t) V2 = circlepoints(cFun,np; dir=:cw) - V2,_ = evenly_sample(V2, np) + V2 = evenly_sample(V2, np) V2_ori = deepcopy(V2) # Q1 = RotXYZ(0.5*π,0.0,0.0) # Define a rotation tensor using Euler angles diff --git a/src/functions.jl b/src/functions.jl index dec1e58..b547d53 100644 --- a/src/functions.jl +++ b/src/functions.jl @@ -3482,41 +3482,51 @@ curvature exists for the B-spline between two adjacent data points then the spacing between points in the output may be non-uniform (despite the allong B-spline distance being uniform). """ -function evenly_sample(V::Vector{Point{ND,TV}}, n::Int; rtol = 1e-8, niter = 1,spline_order=4, close_loop=false) where ND where TV<:Real +function evenly_sample(V::Vector{Point{ND,TV}}, n::Int; rtol=1e-8, niter=1, spline_order=4, close_loop=false) where ND where TV<:Real + + S,_,D = make_geospline(V; rtol=rtol, niter=niter, spline_order=spline_order, close_loop=close_loop) + + # Even range for curve distance + if close_loop == true + l_end = D - D/n + else + l_end = D + end + l = range(0.0, l_end, n) + + return S.(l) # Evaluate interpolator at even distance increments +end + +function make_geospline(V::Vector{Point{ND,TV}}; rtol = 1e-8, niter = 10, spline_order=4, close_loop=false) where ND where TV<:Real LL = curve_length(V) # Initialise as along curve (multi-linear) distance if close_loop == true - LL ./= (last(LL)+ norm(V[1]-V[end])) # Normalise - bc = BSplineKit.Periodic(1.0) # Use periodic bc for closed curves + D = last(LL)+ norm(V[1]-V[end]) + bc = BSplineKit.Periodic(D) # Use periodic bc for closed curves else - LL ./= last(LL) # Normalise + D = last(LL) # Normalise bc = BSplineKit.Natural() # Otherwise use natural end S = BSplineKit.interpolate(LL, deepcopy(V), BSplineOrder(spline_order), bc) # Create interpolator + L = zeros(eltype(LL),length(LL)) # Initialise spline length vector @inbounds for _ in 1:niter - dS = Derivative() * S # spline derivative - L = zeros(eltype(LL),length(LL)) # Initialise spline length vector + dS = BSplineKit.Derivative() * S # spline derivative @inbounds for i in 2:lastindex(LL) # Compute length of segment [i-1, i] L[i] = L[i - 1] + integrate_segment_(dS,LL[i-1], LL[i],rtol) end - - # Normalise + if close_loop == true - L ./= (last(L) + integrate_segment_(dS,LL[end], 1.0,rtol)) + D = last(L) + integrate_segment_(dS,LL[end], D,rtol) + bc = BSplineKit.Periodic(D) else - L ./= last(L) - end + D = last(L) + end S = BSplineKit.interpolate(L, deepcopy(V), BSplineOrder(spline_order), bc) # Create interpolator + LL = L 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 + return S,L,D end function integrate_segment_(dS,l1,l2,rtol) @@ -3526,12 +3536,55 @@ function integrate_segment_(dS,l1,l2,rtol) return segment_length end -function evenly_space(V::Vector{Point{ND,TV}}, pointSpacing=nothing; rtol = 1e-8, niter = 1, close_loop=false) where ND where TV<:Real +function evenly_space(V::Vector{Point{ND,TV}}, pointSpacing=nothing; rtol = 1e-8, niter = 1, spline_order=4, close_loop=false, must_points=nothing) where ND where TV<:Real if isnothing(pointSpacing) pointSpacing = pointspacingmean(V) + end + + if !isnothing(must_points) + sort!(must_points) # sort the indices + if first(must_points)==1 # Check if first is 1 + popfirst!(must_points) # remove start, is already a must point + end + end + + if isnothing(must_points) || isempty(must_points) + n = 1 + ceil(Int,last(curve_length(V; close_loop=close_loop))/pointSpacing) + Vn = evenly_sample(V,n; rtol=rtol, niter=niter, spline_order=spline_order, close_loop=close_loop) + else + # Check must point set + m = length(V) + if close_loop == false && last(must_points)!=m # Check if last needs to be added + push!(must_points,m) + else + push!(must_points,1) # Add first to close over curve + end + + # Construct length parameterised spline + S,L,D = make_geospline(V; rtol=rtol, niter=niter, spline_order=spline_order, close_loop=close_loop) + + Vn = Vector{eltype(V)}() + l1 = 0.0 # Effectively makes first point a must point + for i in eachindex(must_points) + j = must_points[i] + if j==1 # i.e. last step when close_loop == true + l2 = D + else + l2 = L[j] + end + n = 1 + ceil(Int,(l2-l1)/pointSpacing) + l = range(l1,l2,n) + Vn_now = S.(l) + + if j==1 # i.e. last step when close_loop == true + append!(Vn,Vn_now[1:end-1]) + else + append!(Vn,Vn_now) + end + + l1 = deepcopy(l2) + end end - n = ceil(Int,maximum(curve_length(V; close_loop=close_loop))/pointSpacing) - Vn,_ = evenly_sample(V,n; rtol=rtol, niter=niter, close_loop=close_loop) return Vn end @@ -3782,14 +3835,14 @@ function batman(n::Int; symmetric = false, dir=:acw) else m = ceil(Int,n/2)+1 end - V,_ = evenly_sample([Point{3, Float64}(x[i]/22,y[i]/22,0.0) for i in eachindex(x)],m) + V = evenly_sample([Point{3, Float64}(x[i]/22,y[i]/22,0.0) for i in eachindex(x)],m) V2 = [Point{3, Float64}(-V[i][1],V[i][2],0.0) for i in length(V)-1:-1:2] append!(V,V2) else V = [Point{3, Float64}(x[i]/22,y[i]/22,0.0) for i in eachindex(x)] V2 = [Point{3, Float64}(-V[i][1],V[i][2],0.0) for i in length(V)-1:-1:2] append!(V,V2) - V,_ = evenly_sample(V,n) + V = evenly_sample(V,n) end if dir==:cw reverse!(V) diff --git a/test/runtests.jl b/test/runtests.jl index 48ff4b2..49f64dc 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4065,7 +4065,7 @@ end Vt = [GeometryBasics.Point{3,Float64}(x, 2.0*x, 0.0) for x in Xt] # Resample using evenly_sample - Vi, S = evenly_sample(V, n; niter=5) + Vi = evenly_sample(V, n; niter=5) @test sum(norm.(Vi-Vt)) < tol_level # Correct and even spacing @test typeof(V) == typeof(Vi) # Did not manipulate input type @@ -4084,7 +4084,7 @@ end Vt = [GeometryBasics.Point{3, Float64}(r*cos(t),r*sin(t),0) for t in range(0.0,2.0*π,n)] # Resample using evenly_sample - Vi, S = evenly_sample(V, n; niter=5) + Vi = evenly_sample(V, n; niter=5) @test sum(norm.(Vi-Vt)) < tol_level # Correct and even spacing @test typeof(V) == typeof(Vi) # Did not manipulate input type @@ -4096,7 +4096,7 @@ end nc = 6 V = circlepoints(r,nc) n = 25 - Vi, S = evenly_sample(V, n; niter=5,close_loop = true) + Vi = evenly_sample(V, n; niter=5,close_loop = true) @test typeof(V) == typeof(Vi) # Did not manipulate input type @test length(Vi) == n # Correct length end @@ -4108,24 +4108,51 @@ end r = 3.25 nc = 1000 V = circlepoints(r,nc) - Vi, S = evenly_sample(V, nc; niter=5,close_loop = true) - dS = Derivative() * S # spline derivative - L = Comodo.integrate_segment_(dS,0,1,1e-8) + S,_,D = Comodo.make_geospline(V; rtol=1e-8, niter=10, spline_order=4, close_loop=true) + dS = BSplineKit.Derivative() * S # spline derivative + L = Comodo.integrate_segment_(dS,0,D,1e-8) @test isapprox(L,2*pi*r,atol=eps_level) end -@testset "evenly_space" verbose = true begin - eps_level = 1e-3 - np = 10 - t = range(0.0,2.0*π,np) # Parameterisation metric - V = [GeometryBasics.Point{3, Float64}(cos(t[i]),sin(t[i]),0.0) for i in eachindex(t)] +@testset "Comodo.make_geospline" begin + eps_level = 1e-6 + r = 3.25 + nc = 10 + V = circlepoints(r,nc) + S,L,D = Comodo.make_geospline(V; rtol=1e-8, niter=10, spline_order=4, close_loop=true) - Vn = evenly_space(V) - @test isapprox(pointspacingmean(Vn),pointspacingmean(V),atol=eps_level) + @test isa(S,SplineInterpolation) + @test isapprox(V,S.(L),atol=eps_level) - pointSpacing = pointspacingmean(V) - Vn = evenly_space(V,pointSpacing) + S,L,D = Comodo.make_geospline(V; rtol=1e-8, niter=10, spline_order=4, close_loop=false) + + @test isa(S,SplineInterpolation) + @test isapprox(V,S.(L),atol=eps_level) +end + +@testset "evenly_space" verbose = true begin + eps_level = 1e-1 + eps_level2 = 1e-6 + + pointSpacing=0.05 + + n = 5 + r = 3.0 + t = range(0,2π-2π/n,n) + V = [GeometryBasics.Point{3, Float64}(r*cos(tt),r*sin(tt),0.0) for tt ∈ t] + + Vn = evenly_space(V, pointSpacing; close_loop = false, spline_order = 4) + @test isapprox(pointspacingmean(Vn),pointSpacing,atol=eps_level) + + must_points = [1] # This tests and empty must_points vector + Vn = evenly_space(V, pointSpacing; close_loop = true, spline_order = 4, must_points = must_points) # Returns points and spline interpolation object + @test isapprox(pointspacingmean(Vn),pointSpacing,atol=eps_level) + + must_points = [1,2,3] + Vn = evenly_space(V, pointSpacing; close_loop = true, spline_order = 4, must_points = must_points) # Returns points and spline interpolation object @test isapprox(pointspacingmean(Vn),pointSpacing,atol=eps_level) + d = mindist(V,Vn) + @test isapprox(sum(d[must_points]),0.0,atol=eps_level) r = 1000 nc = 1000 @@ -4133,7 +4160,7 @@ end pointSpacing = 1 Vn = evenly_space(V,pointSpacing; close_loop = true) @test isapprox(pointspacingmean(Vn),pointSpacing, atol=eps_level) - @test length(Vn) == ceil(2*pi*r/pointSpacing) + @test length(Vn) == ceil(2*pi*r/pointSpacing)+1 end @@ -4175,17 +4202,17 @@ end P[4 ] = GeometryBasics.Point{3, Float64}( 1.0, 1.0, 1.0) Vc = nbezier(P,nc) # Get Bezier fit points Vc = [vc.*10 for vc in Vc] - Vc,Sc = evenly_sample(Vc, nc) + Vc = evenly_sample(Vc, nc) # Define section curves np = 20 # Number of section points V1 = circlepoints(2.0,np; dir=:acw) - V1,_ = evenly_sample(V1, np) + V1 = evenly_sample(V1, np) Q = RotXYZ(0.0,0.5*π,0.0) # Define a rotation tensor using Euler angles V1 = [Q*v for v in V1] # Rotate the coordinates V2 = circlepoints(2.0,np; dir=:acw) - V2,_ = evenly_sample(V2, np) + V2 = evenly_sample(V2, np) V2 = [v2 .+ Vc[end] for v2 in V2] @testset "reversed" begin