hyrodium / basicbspline.jl Goto Github PK

View Code? Open in Web Editor NEW

91.0 91.0 2.0 374.51 MB

Basic (mathematical) operations for B-spline functions and related things with julia

Home Page: https://hyrodium.github.io/BasicBSpline.jl

License: MIT License

Julia 100.00%

b-spline fitting julia nurbs refinement

basicbspline.jl's People

Contributors

Stargazers

Watchers

Forkers

tangwang-ustc romanwsgit

basicbspline.jl's Issues

autogeneration of index.md

Currently, README.md and docs/src/index.md are duplicated.
The difference is just paths to images.

`changebasis` sometimes returns incorrect nonzero element

julia> using BasicBSpline

julia> p = 3
3

julia> P1 = BSplineSpace{p}(KnotVector(1:8))
BSplineSpace{3, Float64}(KnotVector([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]))

julia> P2 = BSplineSpace{p}(KnotVector(2,3,3,4,5,6,7,8))
BSplineSpace{3, Float64}(KnotVector([2.0, 3.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]))

julia> changebasis(P1,P2)
4×4 Matrix{Float64}:
  0.666667      2.04041e-16  -1.41585e-16  3.30839e-16
  0.333333      1.0          -1.18601e-15  5.25506e-15
 -1.73591e-15   9.12098e-16   1.0          4.36688e-15
  9.68114e-16  -3.80046e-16   3.74332e-16  1.0

add DOI

https://guides.github.com/activities/citable-code/

use Meshes.jl instead of GeometryBasics.jl

Use `Threads` for faster fitting

Add supports for multi-dimension

Currently, some code treats only one or two dimension.

Make `BSplineManifold` function-like object

Comparison with other B-spline Julia packages

I'm not much familiar with the usage of the following packages, but it would be nice to have a detailed comparison in the documentation.

Reduce inner constructor

https://docs.julialang.org/en/v1/manual/constructors/#man-inner-constructor-methods

It is good practice to provide as few inner constructor methods as possible: only those taking all arguments explicitly and enforcing essential error checking and transformation. Additional convenience constructor methods, supplying default values or auxiliary transformations, should be provided as outer constructors that call the inner constructors to do the heavy lifting. This separation is typically quite natural.

`changebasis` fails for non-float `BSplineSpace`

julia> P = BSplineSpace{3}(KnotVector{Int}(1:12))
BSplineSpace{3, Int64}(KnotVector([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]))

julia> changebasis(P,P)
ERROR: InexactError: Int64(NaN)
Stacktrace:
  [1] Int64
    @ ./float.jl:812 [inlined]
  [2] convert
    @ ./number.jl:7 [inlined]
  [3] fill!(A::SubArray{Int64, 2, Matrix{Int64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Vector{Int64}}, false}, x::Float64)
    @ Base ./multidimensional.jl:1062
  [4] copyto!
    @ ./broadcast.jl:921 [inlined]
  [5] materialize!
    @ ./broadcast.jl:871 [inlined]
  [6] materialize!
    @ ./broadcast.jl:868 [inlined]
  [7] _changebasis_R(P::BSplineSpace{0, Int64}, P′::BSplineSpace{0, Int64})
    @ BasicBSpline ~/.julia/dev/BasicBSpline/src/_ChangeBasis.jl:22
  [8] _changebasis_R(P::BSplineSpace{1, Int64}, P′::BSplineSpace{1, Int64}) (repeats 3 times)
    @ BasicBSpline ~/.julia/dev/BasicBSpline/src/_ChangeBasis.jl:26
  [9] changebasis(P::BSplineSpace{3, Int64}, P′::BSplineSpace{3, Int64})
    @ BasicBSpline ~/.julia/dev/BasicBSpline/src/_ChangeBasis.jl:162
 [10] top-level scope
    @ REPL[14]:1

`changebasis(::UniformBSplineSpace, ::UniformBSplineSpace)` can be faster

add method `_array2arrayofvector`

`expandspace` sometimes doesn't produce expanded space

julia> using BasicBSpline

julia> k = KnotVector(0,0,0,1,1,1,2)
KnotVector([0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 2.0])

julia> P = BSplineSpace{2}(k)
BSplineSpace{2, Float64}(KnotVector([0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 2.0]))

julia> P ⊑ expandspace(P, p₊=1)
false

equality `==` for `BSplineSpace` is not defined

julia> P1 = BSplineSpace{2}(KnotVector(1,2))
BSplineSpace{2, Float64}(KnotVector([1.0, 2.0]))

julia> P2 = BSplineSpace{2}(KnotVector(1,2))
BSplineSpace{2, Float64}(KnotVector([1.0, 2.0]))

julia> P1 == P2
false

Not just a issue: Visualization of B-Spline curves in Jupyter Notebook

Hi,
Two days ago I came across your package, and I'm very excited because it is more mathematically oriented than Dierkx or Interpolations. It is very useful for illustrating B-spline manifold properties in a course of Geometric Modeling or Approximation Theory.
I worked on a few examples of curves and displayed them instantly in Jupyter Notebook, instead of saving them as png files. They can be saved as raster or vector images with PlotlyJS.jl, too.
Now I'm working on B-spline surfaces.

`⊊(::BSplineSpace, ::BSplineSpace)` is not supported

julia> using BasicBSpline

julia> P = BSplineSpace{1}(KnotVector(1:8))
BSplineSpace{1, Float64}(KnotVector([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]))

julia> P
BSplineSpace{1, Float64}(KnotVector([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]))

julia> P ⊆ P
true

julia> P ⊊ P
ERROR: MethodError: no method matching length(::BSplineSpace{1, Float64})
Closest candidates are:
  length(::Union{Base.KeySet, Base.ValueIterator}) at ~/julia/julia-1.7.1/share/julia/base/abstractdict.jl:58
  length(::Union{ArrayInterface.BidiagonalIndex, ArrayInterface.TridiagonalIndex, ArrayInterface.BandedBlockBandedMatrixIndex, ArrayInterface.BandedMatrixIndex, ArrayInterface.BlockBandedMatrixIndex}) at ~/.julia/packages/ArrayInterface/TCj9U/src/array_index.jl:209
  length(::Union{LinearAlgebra.Adjoint{T, <:Union{StaticArrays.StaticVector{<:Any, T}, StaticArrays.StaticMatrix{<:Any, <:Any, T}}}, LinearAlgebra.Diagonal{T, <:StaticArrays.StaticVector{<:Any, T}}, LinearAlgebra.Hermitian{T, <:StaticArrays.StaticMatrix{<:Any, <:Any, T}}, LinearAlgebra.LowerTriangular{T, <:StaticArrays.StaticMatrix{<:Any, <:Any, T}}, LinearAlgebra.Symmetric{T, <:StaticArrays.StaticMatrix{<:Any, <:Any, T}}, LinearAlgebra.Transpose{T, <:Union{StaticArrays.StaticVector{<:Any, T}, StaticArrays.StaticMatrix{<:Any, <:Any, T}}}, LinearAlgebra.UnitLowerTriangular{T, <:StaticArrays.StaticMatrix{<:Any, <:Any, T}}, LinearAlgebra.UnitUpperTriangular{T, <:StaticArrays.StaticMatrix{<:Any, <:Any, T}}, LinearAlgebra.UpperTriangular{T, <:StaticArrays.StaticMatrix{<:Any, <:Any, T}}, StaticArrays.StaticVector{<:Any, T}, StaticArrays.StaticMatrix{<:Any, <:Any, T}, StaticArrays.StaticArray{<:Tuple, T}} where T) at ~/.julia/dev/StaticArrays/src/abstractarray.jl:1

Switch to GitHub Actions from Travis

register v0.1.3

@JuliaRegistrator register

Add `issqsubset` for `BSplineDerivativeSpace`

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

use GeometryBasics.jl for control points

bug in refinement algorithm

## 2-dim B-spline manifold
p = 2 # degree of polynomial
k = Knots(1:8) # knot vector
P = FastBSplineSpace(p, k) # B-spline space
rand_a = [rand(2) for i in 1:dim(P), j in 1:dim(P)]
a = [[2 * i - 6.5, 2 * j - 6.5] for i in 1:dim(P), j in 1:dim(P)] + rand_a # random generated control points
M = BSplineSurface([P, P], a) # Define B-spline manifold
save_png("docs/src/img/2dim.png", M, unitlength = 50)

## Refinement
p₊ = [1,1]
M′ = refinement(M, k₊ = k₊, p₊ = p₊)
save_png("docs/src/img/2dim_refinement.png", M′, unitlength = 50)

Add dependency on `RecipesBase`

Fitting for arbitrary knot vector

Add more recipes for plotting

In #183, RecipesBase was added in [deps] table, but the currently supported recipe is only B-spline curve in 2-dimensional space.

We need more recipes for the following types:

BSplineManiofold{1} in 3-dim space
BSplineManiofold{2}
RationalBSplineManiofold
BSplineSpace

register v0.1.4

@JuliaRegistrator register

Add support for NURBS (`RationalBSplineManifold`)

Like this:

struct RationalBSplineManifold{Dim,Deg,T,S<:Tuple,Dim₊₁} <: AbstractRationalBSplineManifold{Dim,Deg}
    bsplinespaces::S
    controlpoints::Array{T,Dim₊₁}
    weights::Array{T,Dim}
end

struct CustomRationalBSplineManifold{Dim,Deg,C,S<:Tuple} <: AbstractRationalBSplineManifold{Dim,Deg}
    bsplinespaces::S
    controlpoints::Array{C,Dim}
    weights::Array{T,Dim}
end

Support `intersection(::BSplineSpace, ::BSplineSpace)`?

In BasicBSpline, BSplineSpace (AbstractBSplineSpace) structs are regarded as a linear space, so we can compute their intersection and +.

I'm not sure we have practical benefits for the features.

Error on `changebasis`

This should be fixed:

julia> using BasicBSpline

julia> p1 = 3
3

julia> k1 = KnotVector(1:8)
KnotVector([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0])

julia> P1 = BSplineSpace{p1}(k1)
BSplineSpace{3, Float64}(KnotVector([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]))

julia> domain(P1)
4.0..5.0

julia> P2 = expandspace(P1, p₊=1, k₊=KnotVector(4.95))
BSplineSpace{4, Float64}(KnotVector([1.0, 2.0, 3.0, 4.0, 4.0, 4.95, 5.0, 5.0, 6.0, 7.0, 8.0]))

julia> P1 ⊑ P2
true

julia> changebasis(P1,P2)
ERROR: LinearAlgebra.SingularException(6)
Stacktrace:
  [1] checknonsingular
    @ ~/julia/julia-1.7.1/share/julia/stdlib/v1.7/LinearAlgebra/src/factorization.jl:19 [inlined]
  [2] checknonsingular
    @ ~/julia/julia-1.7.1/share/julia/stdlib/v1.7/LinearAlgebra/src/factorization.jl:21 [inlined]
  [3] #lu!#146
    @ ~/julia/julia-1.7.1/share/julia/stdlib/v1.7/LinearAlgebra/src/lu.jl:82 [inlined]
  [4] lu(A::LinearAlgebra.Adjoint{Float64, Matrix{Float64}}, pivot::LinearAlgebra.RowMaximum; check::Bool)
    @ LinearAlgebra ~/julia/julia-1.7.1/share/julia/stdlib/v1.7/LinearAlgebra/src/lu.jl:279
  [5] lu (repeats 2 times)
    @ ~/julia/julia-1.7.1/share/julia/stdlib/v1.7/LinearAlgebra/src/lu.jl:278 [inlined]
  [6] \(A::LinearAlgebra.Adjoint{Float64, Matrix{Float64}}, B::LinearAlgebra.Adjoint{Float64, Matrix{Float64}})
    @ LinearAlgebra ~/julia/julia-1.7.1/share/julia/stdlib/v1.7/LinearAlgebra/src/generic.jl:1142
  [7] /(A::Matrix{Float64}, B::Matrix{Float64})
    @ LinearAlgebra ~/julia/julia-1.7.1/share/julia/stdlib/v1.7/LinearAlgebra/src/generic.jl:1150
  [8] _changebasis_sim(P1::BSplineSpace{4, Float64}, P2::BSplineSpace{4, Float64})
    @ BasicBSpline ~/.julia/dev/BasicBSpline/src/_ChangeBasis.jl:124
  [9] _changebasis_I(P::BSplineSpace{3, Float64}, P′::BSplineSpace{4, Float64})
    @ BasicBSpline ~/.julia/dev/BasicBSpline/src/_ChangeBasis.jl:148
 [10] changebasis(P::BSplineSpace{3, Float64}, P′::BSplineSpace{4, Float64})
    @ BasicBSpline ~/.julia/dev/BasicBSpline/src/_ChangeBasis.jl:156
 [11] top-level scope
    @ REPL[8]:1

`changebasis` for `DerivativeBSplineSpace`

Faster mapping

https://en.wikipedia.org/wiki/De_Boor%27s_algorithm#Optimizations

Add `UniformKnotVector` <: `AbstractKnotVector`

Stop promotion `KnotVector([1,2,3])` to `KnotVector([1.0,2.0,3.0])`

julia> KnotVector([1,2,3])
KnotVector([1.0, 2.0, 3.0])

struct BSplineDerivativeSpace
    bsplinespace::BSplineSpace  # original B-spline space
    order::Int  # r-th order derivation
end

`innerproduct_R(::UniformBSplineSpace)` can be calculated with Euler's triangle

Here's a sample code.

using BasicBSpline

k = KnotVector(1:12)
P0 = BSplineSpace{0}(k)
P1 = BSplineSpace{1}(k)
P2 = BSplineSpace{2}(k)
P3 = BSplineSpace{3}(k)
P4 = BSplineSpace{4}(k)

a01, = BasicBSpline.innerproduct_R(P0)
a11,a12 = BasicBSpline.innerproduct_R(P1)
a21,a22,a23 = BasicBSpline.innerproduct_R(P2)
a31,a32,a33,a34 = BasicBSpline.innerproduct_R(P3)
a41,a42,a43,a44,a45 = BasicBSpline.innerproduct_R(P4)

b0 = [rationalize(a01,tol=1e-12)]
b1 = [rationalize(a11,tol=1e-12),rationalize(a12,tol=1e-12)]
b2 = [rationalize(a21,tol=1e-12),rationalize(a22,tol=1e-12),rationalize(a23,tol=1e-12)]
b3 = [rationalize(a31,tol=1e-12),rationalize(a32,tol=1e-12),rationalize(a33,tol=1e-12),rationalize(a34,tol=1e-12)]
b4 = [rationalize(a41,tol=1e-12),rationalize(a42,tol=1e-12),rationalize(a43,tol=1e-12),rationalize(a44,tol=1e-12),rationalize(a45,tol=1e-12)]

The vectors b0 to b4 seem to be able to be calculated with Euler's triangle.
This will make innter_product_R faster and more precise.

add method for tangent vector

Fix type promotion in `bsplinebasis`

julia> k = KnotVector{Int}(1:12)
KnotVector([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12])

julia> P = BSplineSpace{3}(k)
BSplineSpace{3, Int64}(KnotVector([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]))

julia> bsplinebasis(P,1,1)
ERROR: InexactError: Int64(-0.5)
Stacktrace:
 [1] Int64
   @ ./float.jl:812 [inlined]
 [2] _d
   @ ~/.julia/dev/BasicBSpline/src/_BSplineBasis.jl:3 [inlined]
 [3] macro expansion
   @ ~/.julia/dev/BasicBSpline/src/_BSplineBasis.jl:110 [inlined]
 [4] bsplinebasis(P::BSplineSpace{3, Int64}, i::Int64, t::Int64)
   @ BasicBSpline ~/.julia/dev/BasicBSpline/src/_BSplineBasis.jl:110
 [5] top-level scope
   @ REPL[11]:1