Merge pull request #305 from j-adel/patch-2

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

docs/src/howto/polygons.md

@@ -714,6 +714,28 @@ nothing # hide
 
 ![offset poly](../assets/figures/polyfit.png)
 
+There is also the [`polybspline`](@ref) function, which constructs a B-spline curve based on the control points and desired degree with the option of being unclamped.
+
+```@example
+using Luxor, Random # hide
+Drawing(600, 250, "../assets/figures/polyfit.png") # hide
+origin() # hide
+background("white") # hide
+Random.seed!(42) # hide
+
+pts = [Point(x, rand(-100:100)) for x in -280:30:280] # hide
+setopacity(0.7) # hide
+sethue("red") # hide
+prettypoly(pts, () -> circle(O, 5, action = :fill)) # hide
+sethue("darkmagenta") # hide
+poly(polybspline(pts, 200; degree = 5), action = :stroke)
+
+finish() # hide
+nothing # hide
+```
+
+![B-Spline](../assets/figures/polybspline.png)
+
 ## Converting paths to polygons
 
 You can convert the current path to an array of polygons, using [`pathtopoly`](@ref).

src/Luxor.jl

@@ -115,7 +115,7 @@ export Drawing,
     rotatepoint, perpendicular, crossproduct,
     dotproduct, determinant3, distance, prettypoly, polysmooth, polysplit,
     poly, simplify, polycentroid, polysortbyangle, polyhull,
-    polysortbydistance, offsetpoly, polyfit, currentpoint,
+    polysortbydistance, offsetpoly, polyfit, currentpoint, polybspline,
     hascurrentpoint, getworldposition, anglethreepoints, polyperimeter, polydistances, polyportion,
     polyremainder, nearestindex, polyarea, polysample,
     insertvertices!, polymove!, polyscale!, polyrotate!, polyreflect!, @polar, polar, strokepreserve, fillpreserve,

src/polygons.jl

@@ -775,6 +775,81 @@ function polyfit(plist::Array{Point,1}, npoints = 30)
     return resultpoly
 end
 
+"""
+    polybspline(controlpoints::Array{Point,1}, npoints; degree=3, clamped=true)
+
+Generate a B-spline curve from a given set of control points.
+
+# Arguments
+- `controlpoints::Array{Point,1}`: An array of control points that define the B-spline.
+- `npoints=100`: The number of points to generate on the B-spline curve.
+- `degree=3`: The degree of the B-spline. Default is 3.
+- `clamped=true`: A boolean to indicate if the B-spline is clamped. Default is true.
+
+# Returns
+- An array of points on the B-spline curve.
+"""
+function polybspline(controlpoints::Array{Point,1}, npoints=30;degree=3, clamped=true)
+    nCP::Int64 = length(controlpoints)
+    nCP == 0 && error("Error: controlpoints array cannot be empty.")
+    npoints <= 0 && error("Error: npoints must be greater than zero.")
+    degree <= 0 && error("Error: degree must be greater than zero.")
+    degree >= nCP && error("Error: degree cannot be greater than the number of control points.")
+    points = Array{Point,1}(undef, npoints)
+    T = Array{Float64,1}(undef, nCP+degree+1)
+    if clamped
+        T[1:degree] .= 0.
+        for i=degree+1:nCP+1
+            T[i] = (i - degree - 1) / (nCP - degree)
+        end
+        T[nCP+2:end] .= 1.
+    else
+        for i=1:nCP+degree+1
+            T[i] = (i - 1) / (nCP + degree)
+        end
+    end
+
+    """De Boor's algorithm for B-spline evaluation.
+        from https://en.wikipedia.org/wiki/De_Boor%27s_algorithm
+        Arguments
+        ---------
+        k: Index of knot interval that contains x.
+        x: Position.
+        t: Array of knot positions, needs to be padded as described above.
+        c: Array of control points.
+        p: Degree of B-spline.
+    """
+    function deBoor(k::Int64,x::Float64,t::Array{Float64,1},c::Array{Point,1},p::Int64)::Point
+        d = Array{Point,1}(undef, p+1)
+        for j=1:p+1
+            d[j] = c[j+k-p]
+        end
+        @inbounds for r=1:p
+            for j=p+1:-1:r+1
+                alpha=(x-t[j+k-p])/(t[j+1+k-r] - t[j+k-p])
+                d[j]=(1-alpha)*d[j-1]+alpha*d[j]
+            end
+        end
+        return d[p+1]
+    end
+
+    @inbounds for i=0:npoints-1
+        t=i/(npoints-1)
+        if !clamped
+            t =t*(T[nCP+1]-T[degree+1])+T[degree+1]
+        end
+        k=1 #index of knot interval
+        while k < nCP
+            if t<T[k+1]
+                break
+            end
+            k += 1
+        end
+        points[i+1]=deBoor(k-1,t,T,controlpoints,degree)
+    end
+    return points
+end
+
 """
     pathtopoly()
 

test/polybspline-test.jl

@@ -0,0 +1,33 @@
+using Luxor
+
+using Test
+
+using Random
+Random.seed!(42)
+
+function test_polybspline()
+    controlpoints = [Point(rand(), rand()) for _ in 1:10]
+
+    resultClamped = polybspline(controlpoints, 10;clamped = true)
+    resultUnClamped = polybspline(controlpoints, 10;clamped = false)
+    # Check if the result is an array of the correct length
+    @test length(resultClamped) == 10
+    @test length(resultUnClamped) == 10
+
+    # Check if the result is an array of Points
+    @test all(x -> x isa Point, resultClamped)
+    @test all(x -> x isa Point, resultUnClamped)
+
+    # Check that all the values in the result points are non NaN
+    @test all(x -> !isnan(x.x) && !isnan(x.y), resultClamped)
+    @test all(x -> !isnan(x.x) && !isnan(x.y), resultUnClamped)
+
+    # Check that all the result points lie within the bounding box of the control points
+    min_x, max_x = extrema(p.x for p in controlpoints)
+    min_y, max_y = extrema(p.y for p in controlpoints)
+
+    @test all(min_x <= p.x <= max_x && min_y <= p.y <= max_y for p in resultClamped)
+    @test all(min_x <= p.x <= max_x && min_y <= p.y <= max_y for p in resultUnClamped)
+end
+
+test_polybspline()

test/runtests.jl

@@ -40,6 +40,7 @@ function run_all_tests()
         include("point-intersection.jl")
         include("polycross-test.jl")
         include("polyfit-test.jl")
+        include("polybspline-test.jl")
         include("polysample.jl")
         include("polygon-centroid-sort-test.jl")
         include("polygon-test.jl")