Repository: hyrodium/BasicBSpline.jl
Pull request #63 (2020-12-24 13:31:18)

Feature/tangents
solve #58 partly


Patch:
diff --git a/Project.toml b/Project.toml
index 38cc1d06..3e836d01 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,16 +4,14 @@ authors = ["hyrodium <hyrodium@gmail.com>"]
 version = "0.1.1"
 
 [deps]
-EllipsisNotation = "da5c29d0-fa7d-589e-88eb-ea29b0a81949"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
 [compat]
-julia = "1"
-IntervalSets = "0.5"
 FastGaussQuadrature = "0.4"
-EllipsisNotation = "0.4, 1"
+IntervalSets = "0.5"
+julia = "1"
 
 [extras]
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
diff --git a/src/BasicBSpline.jl b/src/BasicBSpline.jl
index d1cd8b77..0daedd7f 100644
--- a/src/BasicBSpline.jl
+++ b/src/BasicBSpline.jl
@@ -1,7 +1,6 @@
 module BasicBSpline
 
 using LinearAlgebra
-using EllipsisNotation
 using IntervalSets
 using FastGaussQuadrature
 
@@ -30,6 +29,7 @@ include("_FastBSplineManifold.jl")
 include("_Refinement.jl")
 include("_Integral.jl")
 include("_Fitting.jl")
+include("_Geometry.jl")
 
 if VERSION < v"1.2.0"
     <=(x) = Fix2(<=, x)
diff --git a/src/_BSplineManifold.jl b/src/_BSplineManifold.jl
index 22ddd5df..3b5019d7 100644
--- a/src/_BSplineManifold.jl
+++ b/src/_BSplineManifold.jl
@@ -84,10 +84,15 @@ Calculate the mapping of B-spline manifold for given parameter.
 """
 function (M::BSplineManifold)(t::AbstractVector{<:Real})
     Ps = M.bsplinespaces
-    𝒂 = M.controlpoints
+    a = M.controlpoints
     d = length(Ps)
-    d̂ = size(𝒂)[end]
-    return [sum(bsplinebasis(Ps, t) .* 𝒂[.., i]) for i in 1:d̂]
+    d̂ = size(a)[end]
+    N = prod(dim.(Ps))
+
+    B = bsplinebasis(Ps, t)
+    B_flat = reshape(B,N)
+    a_flat = reshape(a,N,d̂)
+    return [sum(B_flat .* a_flat[:,i]) for i in 1:d̂]
 end
 
 @doc raw"""
diff --git a/src/_FastBSplineManifold.jl b/src/_FastBSplineManifold.jl
index 1ab27df7..81ef1cd7 100644
--- a/src/_FastBSplineManifold.jl
+++ b/src/_FastBSplineManifold.jl
@@ -33,7 +33,7 @@ struct BSplineCurve{p1} <: AbstractBSplineManifold
     end
 end
 function BSplineCurve(P1,::AbstractBSplineSpace, a::AbstractArray{<:AbstractVector{<:Real},1})
-    a′ = arrayofvector2array(a)
+    a′ = _arrayofvector2array(a)
     return BSplineCurve(P1, a′)
 end
 function BSplineCurve(Ps::AbstractVector{<:AbstractBSplineSpace}, a::AbstractArray{<:Real})
@@ -65,7 +65,7 @@ struct BSplineSurface{p1,p2} <: AbstractBSplineManifold
     end
 end
 function BSplineSurface(P1,::AbstractBSplineSpace, P2,::AbstractBSplineSpace, a::AbstractArray{<:AbstractVector{<:Real},2})
-    a′ = arrayofvector2array(a)
+    a′ = _arrayofvector2array(a)
     return BSplineSurface(P1, P2, a′)
 end
 function BSplineSurface(Ps::AbstractVector{<:AbstractBSplineSpace}, a::AbstractArray{<:Real})
@@ -98,7 +98,7 @@ struct BSplineSolid{p1,p2,p3} <: AbstractBSplineManifold
     end
 end
 function BSplineSolid(P1,::AbstractBSplineSpace, P2,::AbstractBSplineSpace, P3,::AbstractBSplineSpace, a::AbstractArray{<:AbstractVector{<:Real},3})
-    a′ = arrayofvector2array(a)
+    a′ = _arrayofvector2array(a)
     return BSplineSolid(P1, P2, P3, a′)
 end
 function BSplineSolid(Ps::AbstractVector{<:AbstractBSplineSpace}, a::AbstractArray{<:Real})
@@ -115,7 +115,7 @@ function BSplineSolid{q1,q2,q3}(Ps::AbstractVector{<:AbstractBSplineSpace}, a::A
     return BSplineSolid(Ps, a)
 end
 
-function arrayofvector2array(a::AbstractArray{<:AbstractVector{R}, d})::Array{R,d+1} where d where R<:Real
+function _arrayofvector2array(a::AbstractArray{<:AbstractVector{R}, d})::Array{R,d+1} where d where R<:Real
     d̂ = length(a[1])
     s = size(a)
     N = prod(s)
@@ -124,19 +124,9 @@ function arrayofvector2array(a::AbstractArray{<:AbstractVector{R}, d})::Array{R,
     return a′
 end
 
-function array2arrayofvector(a::AbstractArray{R,d})::Array{Vector{R}, d-1} where d where R<:Real
-    s = size(a)
-    d̂ = s[end]
-    N = s[1:end-1]
-    a_flat = reshape(a,prod(N),d̂)
-    a_vec = [a_flat[i,:] for i in 1:prod(N)]
-    a′ = reshape(a_vec,N)
-    return a′
-end
-
 for fname in (:BSplineCurve, :BSplineSurface, :BSplineSolid, :FastBSplineManifold, :BSplineManifold)
     @eval function $fname(Ps::AbstractVector{<:AbstractBSplineSpace}, a::AbstractArray{<:AbstractVector{<:Real}})
-        a′ = arrayofvector2array(a)
+        a′ = _arrayofvector2array(a)
         return $fname(Ps, a′)
     end
 end
diff --git a/src/_Geometry.jl b/src/_Geometry.jl
new file mode 100644
index 00000000..99222c7c
--- /dev/null
+++ b/src/_Geometry.jl
@@ -0,0 +1,82 @@
+# Geometry
+
+"""
+tangentvectors on B-spline manifold
+NOTICE: these methods are experimental
+"""
+tangentvectors
+
+function tangentvectors(M::BSplineCurve{p1}, t1::Real) where {p1}
+    P1, = bsplinespaces(M)
+    a = controlpoints(M)
+    n1 = dim(P1)
+    k1 = P1.knots
+    j1 = _knotindex(P1, t1)
+    b1′ = _bsplinebasis′(P1,t1,j1)
+    d̂ = size(a)[end]
+
+    S1 = Array{Float64}(undef, d̂)
+    f1 = j1-p1
+    l1 = j1
+    for xx in 1:d̂
+        S1[xx] = b1′[1]*a[f1,xx]
+        for i1 in f1+1:l1
+            S1[xx] += b1′[i1-f1+1]*a[i1,xx]
+        end
+    end
+
+    return S1
+end
+
+
+function tangentvectors(M::BSplineSurface{p1,p2}, t1::Real,t2::Real) where {p1} where {p2}
+    P1, P2 = bsplinespaces(M)
+    a = controlpoints(M)
+    n1, n2 = dim(P1), dim(P2)
+    k1, k2 = P1.knots, P2.knots
+    j1, j2 = _knotindex(P1, t1), _knotindex(P2, t2)
+    b1 = _bsplinebasis(P1,t1,j1)
+    b2 = _bsplinebasis(P2,t2,j2)
+    b1′ = _bsplinebasis′(P1,t1,j1)
+    b2′ = _bsplinebasis′(P2,t2,j2)
+    d̂ = size(a)[end]
+
+    S1 = Array{Float64}(undef, d̂)
+    S2 = Array{Float64}(undef, d̂)
+    f1, f2 = j1-p1, j2-p2
+    l1, l2 = j1, j2
+
+    for xx in 1:d̂
+        S2[xx] = b2[1]*a[f1,f2,xx]
+        for i2 in f2+1:l2
+            S2[xx] += b2[i2-f2+1]*a[f1,i2,xx]
+        end
+        S1[xx] = b1′[1]*S2[xx]
+        for i1 in f1+1:l1
+            S2[xx] = b2[1]*a[i1,f2,xx]
+            for i2 in f2+1:l2
+                S2[xx] += b2[i2-f2+1]*a[i1,i2,xx]
+            end
+            S1[xx] += b1′[i1-f1+1]*S2[xx]
+        end
+    end
+    X1 = copy(S1)
+
+    for xx in 1:d̂
+        S2[xx] = b2′[1]*a[f1,f2,xx]
+        for i2 in f2+1:l2
+            S2[xx] += b2′[i2-f2+1]*a[f1,i2,xx]
+        end
+        S1[xx] = b1[1]*S2[xx]
+        for i1 in f1+1:l1
+            S2[xx] = b2′[1]*a[i1,f2,xx]
+            for i2 in f2+1:l2
+                S2[xx] += b2′[i2-f2+1]*a[i1,i2,xx]
+            end
+            S1[xx] += b1[i1-f1+1]*S2[xx]
+        end
+    end
+    X2 = copy(S1)
+
+    return (X1,X2)
+end


Test patch:
diff --git a/test/runtests.jl b/test/runtests.jl
index fd07a776..d0b1eb8c 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -358,25 +358,6 @@ using Random
         @test norm(a_fit - a_ref) < ε
     end
 
-    @testset "array <-> array of vector" begin
-        Random.seed!(42)
-
-        A = rand(3,4)
-        A′ = BasicBSpline.array2arrayofvector(A)
-        A′′ = BasicBSpline.arrayofvector2array(A′)
-        @test A == A′′
-
-        B = rand(3,4,5)
-        B′ = BasicBSpline.array2arrayofvector(B)
-        B′′ = BasicBSpline.arrayofvector2array(B′)
-        @test B == B′′
-
-        C = rand(3,4,5,6)
-        C′ = BasicBSpline.array2arrayofvector(C)
-        C′′ = BasicBSpline.arrayofvector2array(C′)
-        @test C == C′′
-    end
-
     @testset "Fitting-surface_R" begin
         Random.seed!(42)
 
@@ -402,7 +383,7 @@ using Random
         a_ref = M′.controlpoints
 
         a_tmp = fittingcontrolpoints(M, [P1′, P2′])
-        a_fit = BasicBSpline.arrayofvector2array(a_tmp)
+        a_fit = BasicBSpline._arrayofvector2array(a_tmp)
 
         @test norm(a_fit - a_ref) < ε
     end
@@ -432,7 +413,7 @@ using Random
         a_ref = M′.controlpoints
 
         a_tmp = fittingcontrolpoints(M, [P1′, P2′])
-        a_fit = BasicBSpline.arrayofvector2array(a_tmp)
+        a_fit = BasicBSpline._arrayofvector2array(a_tmp)
 
         @test norm(a_fit - a_ref) < ε
     end
@@ -469,7 +450,7 @@ using Random
         a_ref = M′.controlpoints
 
         a_tmp = fittingcontrolpoints(M, [P1′, P2′, P3′])
-        a_fit = BasicBSpline.arrayofvector2array(a_tmp)
+        a_fit = BasicBSpline._arrayofvector2array(a_tmp)
 
         @test norm(a_fit - a_ref) < ε
     end
@@ -506,7 +487,7 @@ using Random
         a_ref = M′.controlpoints
 
         a_tmp = fittingcontrolpoints(M, [P1′, P2′, P3′])
-        a_fit = BasicBSpline.arrayofvector2array(a_tmp)
+        a_fit = BasicBSpline._arrayofvector2array(a_tmp)
 
         @test norm(a_fit - a_ref) < 1.0e-6
     end