search

JuliaGraphics / Luxor.jl

Simple drawings using vector graphics; Cairo "for tourists!"
http://juliagraphics.github.io/Luxor.jl/
Other
588 stars 71 forks source link

Adding bspline function in polygons.jl (evaluates npoints B-spline curve based on control points) #305

Closed j-adel closed 5 months ago

j-adel commented 6 months ago

B-Spline function with clamped/unclamped options using De Boor's Algorithm. It's implemented in a pretty efficient way. However, I'm sure there's plenty of ways to speed it up using vectorization/parallelization. I think this was the objective of the polyfit() experimental function? I didn't want to replace it so I leave that decision to the maintainer.

But yeah I've been working with B-Splines lately and thought it would be useful to have for this library! I tried to make it similar to the coding style and inputs/outputs of the library functions however any changes in that regard might be a good idea. I think the only commonly used feature of splines is making them closed (i.e. a loop) but I haven't gotten around to that yet. The other kind of spline that would be crucial for the library is a Catmull-Rom spline. With several spline types, I was thinking maybe it can be a unified function curve() or spline() which takes as input the spline type, sorta like how Processing does it?

Let me know if you think any changes are necessary for the PR. This is my first significant PR actually so I'm guessing I'm missing a lot of stuff lol.

cormullion commented 6 months ago

This is great. Thanks very much!

Other things to do:

cormullion commented 6 months ago

As for

The other kind of spline that would be crucial for the library is a Catmull-Rom spline.

perhaps just use https://github.com/JeffreySarnoff/CatmullRom.jl ?

hyrodium commented 6 months ago

Thank you for the great PR!

It's implemented in a pretty efficient way. However, I'm sure there's plenty of ways to speed it up using vectorization/parallelization.

Yeah, the performance can be improved. Here are comparisons with my package, BasicBSpline.jl:

julia> using Luxor, BasicBSpline, BenchmarkTools

julia> ps = [Point(randn(), randn()) for _ in 1:8]  # control points
8-element Vector{Point}:
 Point(-0.7994182143543495, 0.9390793075737545)
 Point(-0.22753695196095303, -0.7882913054541987)
 Point(1.2845720190835004, 0.5623233646865434)
 Point(-1.8355529584037287, 0.8438730777812197)
 Point(0.5510575401065174, 1.402746732606732)
 Point(0.45539845640833243, -1.265031994706883)
 Point(-0.1908845986894419, -0.10172057455987063)
 Point(-1.386438336193102, -0.29019694915633065)

julia> P = BSplineSpace{3}(KnotVector([0,0,0,0,1,2,3,4,5,5,5,5]))
BSplineSpace{3, Int64, KnotVector{Int64}}(KnotVector([0, 0, 0, 0, 1, 2, 3, 4, 5, 5, 5, 5]))

julia> C = BSplineManifold(ps, P)
BSplineManifold{1, (3,), Point, Int64, Tuple{BSplineSpace{3, Int64, KnotVector{Int64}}}}(Point[Point(-0.7994182143543495, 0.9390793075737545), Point(-0.22753695196095303, -0.7882913054541987), Point(1.2845720190835004, 0.5623233646865434), Point(-1.8355529584037287, 0.8438730777812197), Point(0.5510575401065174, 1.402746732606732), Point(0.45539845640833243, -1.265031994706883), Point(-0.1908845986894419, -0.10172057455987063), Point(-1.386438336193102, -0.29019694915633065)], (BSplineSpace{3, Int64, KnotVector{Int64}}(KnotVector([0, 0, 0, 0, 1, 2, 3, 4, 5, 5, 5, 5])),))

julia> C.(range(domain(P), length=21))
21-element Vector{Point}:
 Point(-0.7994182143543495, 0.9390793075737545)
 Point(-0.35288449596438315, 0.051968983954802477)
 Point(0.0612559346484744, -0.18664398382581504)
 Point(0.34154483103380284, -0.029266470466635504)
 Point(0.38652394674118207, 0.2715946493338039)
 Point(0.1400172011610212, 0.5107818139356997)
 Point(-0.27302282295295666, 0.6725347139341815)
 Point(-0.6823613770061979, 0.7884423529831135)
 Point(-0.9177637124041496, 0.890093734736359)
 Point(-0.8581722393237627, 0.99468809056206)
 Point(-0.5792380030280005, 1.0618655626854696)
 Point(-0.205789207551331, 1.0368765210461206)
 Point(0.1373459430717789, 0.864971335583544)
 Point(0.35033275819609594, 0.5204195808298326)
 Point(0.4333106007333227, 0.09356764968732044)
 Point(0.41141234698439616, -0.29621886034909783)
 Point(0.309770873250253, -0.5295743517845274)
 Point(0.13797423332130898, -0.5305523328668529)
 Point(-0.1565688090541049, -0.3968827348150751)
 Point(-0.6419943126381785, -0.26971459459097424)
 Point(-1.386438336193102, -0.29019694915633065)

julia> Luxor.bspline(ps, 20)  # The same results as previous one (except for floating point errors)
21-element Vector{Any}:
 Point(-0.7994182143543495, 0.9390793075737545)
 Point(-0.35288449596438315, 0.0519689839548025)
 Point(0.061255934648474314, -0.18664398382581504)
 Point(0.3415448310338027, -0.029266470466635615)
 Point(0.3865239467411821, 0.27159464933380395)
 Point(0.14001720116102095, 0.5107818139356994)
 Point(-0.27302282295295643, 0.6725347139341815)
 Point(-0.6823613770061976, 0.7884423529831134)
 Point(-0.9177637124041498, 0.8900937347363591)
 Point(-0.8581722393237624, 0.9946880905620599)
 Point(-0.5792380030280007, 1.0618655626854698)
 Point(-0.2057892075513309, 1.0368765210461206)
 Point(0.13734594307177872, 0.8649713355835442)
 Point(0.3503327581960958, 0.5204195808298324)
 Point(0.4333106007333225, 0.09356764968732088)
 Point(0.41141234698439616, -0.29621886034909783)
 Point(0.30977087325025293, -0.5295743517845275)
 Point(0.1379742333213091, -0.5305523328668529)
 Point(-0.15656880905410495, -0.39688273481507497)
 Point(-0.6419943126381779, -0.26971459459097435)
 Point(-1.386438336193102, -0.29019694915633065)

julia> @benchmark C.(range(domain(P), length=21))
BenchmarkTools.Trial: 10000 samples with 23 evaluations.
 Range (min … max):  950.913 ns … 219.537 μs  ┊ GC (min … max): 0.00% … 98.23%
 Time  (median):       1.001 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):     1.156 μs ±   2.916 μs  ┊ GC (mean ± σ):  3.50% ±  1.39%

  ▂▆██▆▄▂    ▁▂▂▁▁▁▃▄▄▃▂▁▁    ▁ ▁▁▁▁▂▂▂▂▂▂▁ ▁▁                  ▂
  ████████████████████████████████████████████████▇▇▇▆▆▅▆▆▆▅▄▅▅ █
  951 ns        Histogram: log(frequency) by time       1.79 μs <

 Memory estimate: 576 bytes, allocs estimate: 4.

julia> @benchmark Luxor.bspline(ps, 20)  # I have commented out `println("t: ",t," k: ",k)` in this PR.
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  16.361 μs …  3.592 ms  ┊ GC (min … max): 0.00% … 97.29%
 Time  (median):     18.024 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   22.123 μs ± 66.233 μs  ┊ GC (mean ± σ):  6.41% ±  2.16%

  ▂▇██▆▅▄▃▄▄▄▃▂▁▁     ▁▁▁▁▂▂▂▂▃▄▃▃▂▂▂▂▁▁                      ▂
  █████████████████████████████████████████▇▇▆▆▆▆▅▆▆▆▆▅▅▅▄▄▄▅ █
  16.4 μs      Histogram: log(frequency) by time      39.5 μs <

 Memory estimate: 24.31 KiB, allocs estimate: 908.

I'm not sure where the performance difference comes from, but one reason would be type-instability as commented in https://github.com/JuliaGraphics/Luxor.jl/pull/305#discussion_r1544298801. (Note that BasicBSpline.jl does not use vectorization/parallelization explicitly.)

As a maintainer of BasicBSpline.jl, I'm happy to help updating the code with BasicBSpline.jl. But there are more B-spline packages in Julia ecosystem, so another package might be better choice. One possible concern with BasicBSpline.jl would be increased load-time.

julia> @time using Luxor
  0.567440 seconds (546.01 k allocations: 39.472 MiB, 11.75% gc time, 3.79% compilation time)

julia> @time using BasicBSpline
  0.347235 seconds (334.66 k allocations: 17.003 MiB)

In any cases, avoiding reinventing the wheel is advisable for better interoperability with other packages and enhanced performance.

hyrodium commented 6 months ago

A B-spline curve can be split into Bezier curves, so I think adding a method to generate a BezierPath instance would be useful. Here is an example with BasicBSpline.jl.

julia> ps  # Defined randomly as in the previous comment
8-element Vector{Point}:
 Point(-0.7994182143543495, 0.9390793075737545)
 Point(-0.22753695196095303, -0.7882913054541987)
 Point(1.2845720190835004, 0.5623233646865434)
 Point(-1.8355529584037287, 0.8438730777812197)
 Point(0.5510575401065174, 1.402746732606732)
 Point(0.45539845640833243, -1.265031994706883)
 Point(-0.1908845986894419, -0.10172057455987063)
 Point(-1.386438336193102, -0.29019694915633065)

julia> P = BSplineSpace{3}(KnotVector([0,0,0,0,1,2,3,4,5,5,5,5]))
BSplineSpace{3, Int64, KnotVector{Int64}}(KnotVector([0, 0, 0, 0, 1, 2, 3, 4, 5, 5, 5, 5]))

julia> C = BSplineManifold(ps, P)  # Same curve definition as in the previous comment
BSplineManifold{1, (3,), Point, Int64, Tuple{BSplineSpace{3, Int64, KnotVector{Int64}}}}(Point[Point(-0.7994182143543495, 0.9390793075737545), Point(-0.22753695196095303, -0.7882913054541987), Point(1.2845720190835004, 0.5623233646865434), Point(-1.8355529584037287, 0.8438730777812197), Point(0.5510575401065174, 1.402746732606732), Point(0.45539845640833243, -1.265031994706883), Point(-0.1908845986894419, -0.10172057455987063), Point(-1.386438336193102, -0.29019694915633065)], (BSplineSpace{3, Int64, KnotVector{Int64}}(KnotVector([0, 0, 0, 0, 1, 2, 3, 4, 5, 5, 5, 5])),))

julia> P2 = BSplineSpace{3}(unique(knotvector(P)) * 4)
BSplineSpace{3, Int64, KnotVector{Int64}}(KnotVector([0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5]))

julia> C2 = refinement(C, P2)
BSplineManifold{1, (3,), Point, Int64, Tuple{BSplineSpace{3, Int64, KnotVector{Int64}}}}(Point[Point(-0.7994182143543495, 0.9390793075737545), Point(-0.22753695196095303, -0.7882913054541987), Point(0.5285175335612737, -0.11298397038382763), Point(0.3865239467411823, 0.27159464933380395), Point(0.3865239467411823, 0.27159464933380395), Point(0.24453035992109085, 0.6561732690514355), Point(-0.7955112992413191, 0.7500231734163277), Point(-0.9177637124041496, 0.8900937347363591), Point(-0.9177637124041496, 0.8900937347363591), Point(-1.0400161255669802, 1.0301642960563906), Point(-0.24447929273023133, 1.2164555143315614), Point(0.1373459430717789, 0.8649713355835442), Point(0.1373459430717789, 0.8649713355835442), Point(0.5191711788737892, 0.513487156835527), Point(0.4872848176410608, -0.37577241893567803), Point(0.309770873250253, -0.5295743517845274), Point(0.309770873250253, -0.5295743517845274), Point(0.13225692885944523, -0.6833762846333767), Point(-0.1908845986894419, -0.10172057455987063), Point(-1.386438336193102, -0.29019694915633065)], (BSplineSpace{3, Int64, KnotVector{Int64}}(KnotVector([0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5])),))

julia> BezierPath([BezierPathSegment(controlpoints(C2)[1+4i:4+4i]) for i in 0:4])
5-element BezierPath:
 [Point(-0.7994182143543495, 0.9390793075737545), Point(-0.22753695196095303, -0.7882913054541987), Point(0.5285175335612737, -0.11298397038382763), Point(0.3865239467411823, 0.27159464933380395)]
 [Point(0.3865239467411823, 0.27159464933380395), Point(0.24453035992109085, 0.6561732690514355), Point(-0.7955112992413191, 0.7500231734163277), Point(-0.9177637124041496, 0.8900937347363591)]
 [Point(-0.9177637124041496, 0.8900937347363591), Point(-1.0400161255669802, 1.0301642960563906), Point(-0.24447929273023133, 1.2164555143315614), Point(0.1373459430717789, 0.8649713355835442)]
 [Point(0.1373459430717789, 0.8649713355835442), Point(0.5191711788737892, 0.513487156835527), Point(0.4872848176410608, -0.37577241893567803), Point(0.309770873250253, -0.5295743517845274)]
 [Point(0.309770873250253, -0.5295743517845274), Point(0.13225692885944523, -0.6833762846333767), Point(-0.1908845986894419, -0.10172057455987063), Point(-1.386438336193102, -0.29019694915633065)]

The curves C and C2 can be visualized with Plots.jl. They are equivalent to each other as a curve, but the control points are different. C2 has more control points, and the positions are equal to the control points of Bezier curves.

using StaticArrays, Plots
# We need to convert the type of control points from `Luxor.Point` to `StaticArrays.SVector{2, Float64}` to plot the curves.
plot(BSplineManifold(SVector{2}.(Tuple.(controlpoints(C))), bsplinespaces(C)), xlims=(-2,2), ylims=(-2,2), label="original B-spline curve")
plot(BSplineManifold(SVector{2}.(Tuple.(controlpoints(C2))), bsplinespaces(C2)), xlims=(-2,2), ylims=(-2,2), label="B-spline curve with Bezier control points")

C C2

An example script to plot the BezierPath instance with Luxor.jl:

# Invert around the y-axis because Luxor and Plots have different coordinate systems.
# Magnify the points by `*100` to update the size of the curve.
ps2 = [Point(p.x, -p.y)*100 for p in controlpoints(C2)]
bezpath = BezierPath([BezierPathSegment(ps2[1+4i:4+4i]) for i in 0:4])
@svg begin
    sethue(1.0, 0.0, 0.0)
    drawbezierpath(bezpath, action = :stroke, close=true)
    sethue(0.5, 0.5, 0.5)
    circle.(ps2, 3, action = :fill)
    poly(ps2, action = :stroke, close=false)
end

luxor-drawing-220339_049

cormullion commented 6 months ago

@hyrodium Thanks for your great analysis!

>julia> @time using Luxor
  0.567440 seconds (546.01 k allocations: 39.472 MiB, 11.75% gc time, 3.79% compilation time)

julia> @time using BasicBSpline
  0.347235 seconds (334.66 k allocations: 17.003 MiB)

– this is mostly StaticArrays.jl?

julia-1.10> @time_imports using BasicBSpline
               ┌ 5.5 ms SuiteSparse_jll.__init__() 
     26.5 ms  SuiteSparse_jll 76.05% compilation time
               ┌ 5.5 ms SparseArrays.CHOLMOD.__init__() 93.91% compilation time
    154.0 ms  SparseArrays 3.36% compilation time
     12.3 ms  IntervalSets
      0.5 ms  IntervalSets → IntervalSetsRandomExt
     10.3 ms  Preferences
      0.6 ms  PrecompileTools
      1.0 ms  StaticArraysCore
    225.0 ms  StaticArrays
    356.7 ms  BasicBSpline

In any cases, avoiding reinventing the wheel is advisable for better interoperability with other packages and enhanced performance.

Yes. Due to Luxor.jl's great age (10 years old this year!) it was never very good at interoperability with other packages - it predates most of the cool and useful ones. Even StaticArrays.jl was a risky and costly dependency back when it was first developed. Recently I've started trying use external packages more - for example, replacing some of the polygon operations with PolygonAlgorithms.jl. And with the new package dependencies/extensions, we can use more external packages without incurring the cost of loading. Perhaps BasicBSplines would work as an extension?

A B-spline curve can be split into Bezier curves, so I think adding a method to generate a BezierPath instance would be useful.

Yes. I've forgotten how they work, now, but I think more Bezier handling is always good.

j-adel commented 6 months ago

perhaps just use https://github.com/JeffreySarnoff/CatmullRom.jl ?

Maybe... I can try it out and see if it's easily compatible with Luxor. I don't think implementing the basic algorithm will take too long though. Depends on how extensive we want the feature set for the splines and whether you're ok with adding dependencies.

hyrodium commented 6 months ago

– this is mostly StaticArrays.jl?

Yes! There is a plan to replace the dependency with StaticArraysCore.jl, but I haven't finished yet.

Perhaps BasicBSplines would work as an extension?

I think adding a constructor method BezierPath(::BasicBSpline.BSplineManifold{1, (3,), Point}) in an extension is enough and we don't have to add a new function like polybspline here.

julia> C2  # The same definition as in my previous comment
BSplineManifold{1, (3,), Point, Int64, Tuple{BSplineSpace{3, Int64, KnotVector{Int64}}}}(Point[Point(-0.7994182143543495, 0.9390793075737545), Point(-0.22753695196095303, -0.7882913054541987), Point(0.5285175335612737, -0.11298397038382763), Point(0.3865239467411823, 0.27159464933380395), Point(0.3865239467411823, 0.27159464933380395), Point(0.24453035992109085, 0.6561732690514355), Point(-0.7955112992413191, 0.7500231734163277), Point(-0.9177637124041496, 0.8900937347363591), Point(-0.9177637124041496, 0.8900937347363591), Point(-1.0400161255669802, 1.0301642960563906), Point(-0.24447929273023133, 1.2164555143315614), Point(0.1373459430717789, 0.8649713355835442), Point(0.1373459430717789, 0.8649713355835442), Point(0.5191711788737892, 0.513487156835527), Point(0.4872848176410608, -0.37577241893567803), Point(0.309770873250253, -0.5295743517845274), Point(0.309770873250253, -0.5295743517845274), Point(0.13225692885944523, -0.6833762846333767), Point(-0.1908845986894419, -0.10172057455987063), Point(-1.386438336193102, -0.29019694915633065)], (BSplineSpace{3, Int64, KnotVector{Int64}}(KnotVector([0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5])),))

julia> path = BezierPath([BezierPathSegment(controlpoints(C2)[1+4i:4+4i]) for i in 0:4])  # This instance will be obtained via the method `BezierPath(::BasicBSpline.BSplineManifold{1, (3,), Point})`.
5-element BezierPath:
 [Point(-0.7994182143543495, 0.9390793075737545), Point(-0.22753695196095303, -0.7882913054541987), Point(0.5285175335612737, -0.11298397038382763), Point(0.3865239467411823, 0.27159464933380395)]
 [Point(0.3865239467411823, 0.27159464933380395), Point(0.24453035992109085, 0.6561732690514355), Point(-0.7955112992413191, 0.7500231734163277), Point(-0.9177637124041496, 0.8900937347363591)]
 [Point(-0.9177637124041496, 0.8900937347363591), Point(-1.0400161255669802, 1.0301642960563906), Point(-0.24447929273023133, 1.2164555143315614), Point(0.1373459430717789, 0.8649713355835442)]
 [Point(0.1373459430717789, 0.8649713355835442), Point(0.5191711788737892, 0.513487156835527), Point(0.4872848176410608, -0.37577241893567803), Point(0.309770873250253, -0.5295743517845274)]
 [Point(0.309770873250253, -0.5295743517845274), Point(0.13225692885944523, -0.6833762846333767), Point(-0.1908845986894419, -0.10172057455987063), Point(-1.386438336193102, -0.29019694915633065)]

julia> bezierpathtopoly(path; steps=3)  # A bezier path can be converted into a vector of points. I'm not sure how to include the last end point. (See `[1:20]` in the next evaluation)
20-element Vector{Point}:
 Point(-0.7994182143543495, 0.9390793075737545)
 Point(-0.3528844959643831, 0.05196898395480254)
 Point(0.06125593464847445, -0.18664398382581515)
 Point(0.341544831033803, -0.02926647046663583)
 Point(0.3865239467411823, 0.27159464933380395)
 Point(0.14001720116102112, 0.5107818139356995)
 Point(-0.27302282295295655, 0.6725347139341815)
 Point(-0.6823613770061978, 0.7884423529831135)
 Point(-0.9177637124041496, 0.8900937347363591)
 Point(-0.8581722393237626, 0.99468809056206)
 Point(-0.5792380030280007, 1.0618655626854698)
 Point(-0.205789207551331, 1.0368765210461206)
 Point(0.1373459430717789, 0.8649713355835442)
 Point(0.35033275819609594, 0.5204195808298328)
 Point(0.43331060073332284, 0.09356764968732056)
 Point(0.4114123469843964, -0.2962188603490976)
 Point(0.309770873250253, -0.5295743517845274)
 Point(0.13797423332130898, -0.5305523328668528)
 Point(-0.15656880905410486, -0.396882734815075)
 Point(-0.6419943126381785, -0.2697145945909743)

julia> C.(range(domain(bsplinespaces(C)[1]), length=21))[1:20]  # The same results as previous one (except for floating point errors)
20-element Vector{Point}:
 Point(-0.7994182143543495, 0.9390793075737545)
 Point(-0.35288449596438315, 0.051968983954802477)
 Point(0.0612559346484744, -0.18664398382581504)
 Point(0.34154483103380284, -0.029266470466635504)
 Point(0.38652394674118207, 0.2715946493338039)
 Point(0.1400172011610212, 0.5107818139356997)
 Point(-0.27302282295295666, 0.6725347139341815)
 Point(-0.6823613770061979, 0.7884423529831135)
 Point(-0.9177637124041496, 0.890093734736359)
 Point(-0.8581722393237627, 0.99468809056206)
 Point(-0.5792380030280005, 1.0618655626854696)
 Point(-0.205789207551331, 1.0368765210461206)
 Point(0.1373459430717789, 0.864971335583544)
 Point(0.35033275819609594, 0.5204195808298326)
 Point(0.4333106007333227, 0.09356764968732044)
 Point(0.41141234698439616, -0.29621886034909783)
 Point(0.309770873250253, -0.5295743517845274)
 Point(0.13797423332130898, -0.5305523328668529)
 Point(-0.1565688090541049, -0.3968827348150751)
 Point(-0.6419943126381785, -0.26971459459097424)

The changes can be found in this commit. I don't think I can edit this PR but you can merge the commit to this PR. I can create a new one which I think is not preferred?

I'm not sure I understand what you mean. You have already added the commit 2095bbb to this PR, right?

j-adel commented 6 months ago

@hyrodium First of all thank you for the feedback! I didn't see your BasicBSpline.jl package when I wrote this code, seems really good! I just wanted to write one myself and then thought I should submit it to Luxor. I am not the most experienced with package load times and pre-compilation so I'm unsure what is the best option for Luxor. I could spend some time to make polybspline run more efficiently to somewhat match that of BasicBSpline.

I'm not sure I understand what you mean. You have already added the commit 2095bbb to this PR, right?

Ah good, yes that's the one. As I said this my first ever PR I don't have that much experience with them😅

j-adel commented 6 months ago

Has a decision been made for going with this PR or instead using BasicBSpline.jl? I'm willing to spend some time to get the function up to speed with the package and doing the remainder documentation/package complementary tasks stated in this comment. I have to say I was quite excited to make my first ever major code contribution and to one of my favorite packages, non other. However, the quality and speed of the package should take priority so if going with the BasicBSpline is significantly better, then of course that's what you should go with.

cormullion commented 6 months ago

Yes, Bsplines are cool. Happy to have this function added, once the PR is ready to merge. Also happy to have more package extensions linking to specialized packages for other curves - although someone needs to write them first... 😂

I think it's good to have different ways of doing things here; the end result is to put graphics on drawings, so various approaches with different strengths and weaknesses with regard to performance or flexibility (assuming they're mathematically correct) is fine by me.

[Contributing to open source occasionally leads to disappointment as well as excitement - I've spent quite a few hours of my life on PRs that went nowhere. You get used to it! 😀]

j-adel commented 5 months ago

Ok the PR is now updated with significantly improved performance of the algorithm (mainly through eliminating a lot of type instabilities). Also added package requirements listed in https://github.com/JuliaGraphics/Luxor.jl/pull/305#issuecomment-2027005194 except for ChangeLog. I've added a suite of tests that I think if passed would mean the spline is very likely to be still working. The documentation I left quite light and similar to polyfit however I did not clarify the difference between them (perhaps polyfit should be deprecated as polybspline advertises the same purpose?) As for the Changelog, perhaps @cormullion could add it as I'm not exactly sure what to write there.

Following the test made by @hyrodium here 

julia> @benchmark polybspline(ps,20)
BenchmarkTools.Trial: 10000 samples with 10 evaluations.
 Range (min … max):  1.300 μs … 61.330 μs  ┊ GC (min … max): 0.00% … 93.01%
 Time  (median):     1.400 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   1.565 μs ±  1.412 μs  ┊ GC (mean ± σ):  2.03% ±  2.28%

   ▄█     
  ▂██▆▃▂▂▃▃▃▂▂▂▁▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▁
  1.3 μs         Histogram: frequency by time        3.45 μs <

 Memory estimate: 3.05 KiB, allocs estimate: 22.

 julia> @benchmark C.(range(domain(P), length=21))
BenchmarkTools.Trial: 10000 samples with 25 evaluations.
 Range (min … max):  928.000 ns … 90.380 μs  ┊ GC (min … max): 0.00% … 97.73%
 Time  (median):     968.000 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):     1.042 μs ±  1.541 μs  ┊ GC (mean ± σ):  2.47% ±  1.69%

  ▃█▇▄▂▂▂  ▁▄▃▁                                                ▂
  ████████▇████▇▇▆▆▆▅▃▄▁▄▄▃▃▁▁▃▆▆▇▆▆▆▅▄▆▅▆▅▅▅▆▄▅▄▅▅▅▅▃▃▁▃▄▄▅▅▅ █
  928 ns        Histogram: log(frequency) by time      2.12 μs <

 Memory estimate: 576 bytes, allocs estimate: 4.

I couldn't get the performance to match that of BasicBSpline which is currently ~150% the speed of polybspline. I've tried all the known suspects for code speed in Julia but I don't I can improve upon it without StaticArrays.jl or the new 1.11 Memory type.

cormullion commented 5 months ago

@j-adel Thanks for your contribution!