Skip to content

Commit 76c545b

Browse files
ChrisRackauckasChrisRackauckas
authored
Merge pull request #349 from ashutosh-b-b/bb/bspline_array
BSplineInterpolation: Add support for higher order arrays
2 parents c6c0090 + 26b109f commit 76c545b

File tree

5 files changed

+357
-5
lines changed

5 files changed

+357
-5
lines changed

‎src/derivatives.jl‎

Lines changed: 47 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -164,6 +164,30 @@ function _derivative(A::BSplineInterpolation{<:AbstractVector{<:Number}}, t::Num
164164
ducum * A.d * scale
165165
end
166166

167+
function _derivative(
168+
A::BSplineInterpolation{<:AbstractArray{<:Number, N}}, t::Number, iguess) where {N}
169+
# change t into param [0 1]
170+
ax_u = axes(A.u)[1:(end - 1)]
171+
t < A.t[1] && return zeros(size(A.u)[1:(end - 1)]...)
172+
t > A.t[end] && return zeros(size(A.u)[1:(end - 1)]...)
173+
idx = get_idx(A, t, iguess)
174+
n = length(A.t)
175+
scale = (A.p[idx + 1] - A.p[idx]) / (A.t[idx + 1] - A.t[idx])
176+
t_ = A.p[idx] + (t - A.t[idx]) * scale
177+
sc = t isa ForwardDiff.Dual ? zeros(eltype(t), n) : A.sc
178+
spline_coefficients!(sc, A.d - 1, A.k, t_)
179+
ducum = zeros(size(A.u)[1:(end - 1)]...)
180+
if t == A.t[1]
181+
ducum = (A.c[ax_u..., 2] - A.c[ax_u..., 1]) / (A.k[A.d + 2])
182+
else
183+
for i in 1:(n - 1)
184+
ducum = ducum +
185+
sc[i + 1] * (A.c[ax_u..., i + 1] - A.c[ax_u..., i]) /
186+
(A.k[i + A.d + 1] - A.k[i + 1])
187+
end
188+
end
189+
ducum * A.d * scale
190+
end
167191
# BSpline Curve Approx
168192
function _derivative(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, iguess)
169193
# change t into param [0 1]
@@ -185,6 +209,29 @@ function _derivative(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, ig
185209
ducum * A.d * scale
186210
end
187211

212+
function _derivative(
213+
A::BSplineApprox{<:AbstractArray{<:Number, N}}, t::Number, iguess) where {N}
214+
# change t into param [0 1]
215+
ax_u = axes(A.u)[1:(end - 1)]
216+
t < A.t[1] && return zeros(size(A.u)[1:(end - 1)]...)
217+
t > A.t[end] && return zeros(size(A.u)[1:(end - 1)]...)
218+
idx = get_idx(A, t, iguess)
219+
scale = (A.p[idx + 1] - A.p[idx]) / (A.t[idx + 1] - A.t[idx])
220+
t_ = A.p[idx] + (t - A.t[idx]) * scale
221+
sc = t isa ForwardDiff.Dual ? zeros(eltype(t), A.h) : A.sc
222+
spline_coefficients!(sc, A.d - 1, A.k, t_)
223+
ducum = zeros(size(A.u)[1:(end - 1)]...)
224+
if t == A.t[1]
225+
ducum = (A.c[ax_u..., 2] - A.c[ax_u..., 1]) / (A.k[A.d + 2])
226+
else
227+
for i in 1:(A.h - 1)
228+
ducum = ducum +
229+
sc[i + 1] * (A.c[ax_u..., i + 1] - A.c[ax_u..., i]) /
230+
(A.k[i + A.d + 1] - A.k[i + 1])
231+
end
232+
end
233+
ducum * A.d * scale
234+
end
188235
# Cubic Hermite Spline
189236
function _derivative(
190237
A::CubicHermiteSpline{<:AbstractVector{<:Number}}, t::Number, iguess)

‎src/interpolation_caches.jl‎

Lines changed: 170 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -597,7 +597,7 @@ struct BSplineInterpolation{uType, tType, pType, kType, cType, scType, T, N} <:
597597
end
598598

599599
function BSplineInterpolation(
600-
u, t, d, pVecType, knotVecType; extrapolate = false, assume_linear_t = 1e-2)
600+
u::AbstractVector, t, d, pVecType, knotVecType; extrapolate = false, assume_linear_t = 1e-2)
601601
u, t = munge_data(u, t)
602602
n = length(t)
603603
n < d + 1 && error("BSplineInterpolation needs at least d + 1, i.e. $(d+1) points.")
@@ -665,6 +665,79 @@ function BSplineInterpolation(
665665
u, t, d, p, k, c, sc, pVecType, knotVecType, extrapolate, assume_linear_t)
666666
end
667667

668+
function BSplineInterpolation(
669+
u::AbstractArray{T, N}, t, d, pVecType, knotVecType; extrapolate = false,
670+
assume_linear_t = 1e-2) where {T, N}
671+
u, t = munge_data(u, t)
672+
n = length(t)
673+
n < d + 1 && error("BSplineInterpolation needs at least d + 1, i.e. $(d+1) points.")
674+
s = zero(eltype(u))
675+
p = zero(t)
676+
k = zeros(eltype(t), n + d + 1)
677+
l = zeros(eltype(u), n - 1)
678+
p[1] = zero(eltype(t))
679+
p[end] = one(eltype(t))
680+
681+
ax_u = axes(u)[1:(end - 1)]
682+
683+
for i in 2:n
684+
s += √((t[i] - t[i - 1])^2 + sum((u[ax_u..., i] - u[ax_u..., i - 1]) .^ 2))
685+
l[i - 1] = s
686+
end
687+
if pVecType == :Uniform
688+
for i in 2:(n - 1)
689+
p[i] = p[1] + (i - 1) * (p[end] - p[1]) / (n - 1)
690+
end
691+
elseif pVecType == :ArcLen
692+
for i in 2:(n - 1)
693+
p[i] = p[1] + l[i - 1] / s * (p[end] - p[1])
694+
end
695+
end
696+
697+
lidx = 1
698+
ridx = length(k)
699+
while lidx <= (d + 1) && ridx >= (length(k) - d)
700+
k[lidx] = p[1]
701+
k[ridx] = p[end]
702+
lidx += 1
703+
ridx -= 1
704+
end
705+
706+
ps = zeros(eltype(t), n - 2)
707+
s = zero(eltype(t))
708+
for i in 2:(n - 1)
709+
s += p[i]
710+
ps[i - 1] = s
711+
end
712+
713+
if knotVecType == :Uniform
714+
# uniformly spaced knot vector
715+
# this method is not recommended because, if it is used with the chord length method for global interpolation,
716+
# the system of linear equations would be singular.
717+
for i in (d + 2):n
718+
k[i] = k[1] + (i - d - 1) // (n - d) * (k[end] - k[1])
719+
end
720+
elseif knotVecType == :Average
721+
# average spaced knot vector
722+
idx = 1
723+
if d + 2 <= n
724+
k[d + 2] = 1 // d * ps[d]
725+
end
726+
for i in (d + 3):n
727+
k[i] = 1 // d * (ps[idx + d] - ps[idx])
728+
idx += 1
729+
end
730+
end
731+
# control points
732+
sc = zeros(eltype(t), n, n)
733+
spline_coefficients!(sc, d, k, p)
734+
c = (sc \ reshape(u, prod(size(u)[1:(end - 1)]), :)')'
735+
c = reshape(c, size(u)...)
736+
sc = zeros(eltype(t), n)
737+
BSplineInterpolation(
738+
u, t, d, p, k, c, sc, pVecType, knotVecType, extrapolate, assume_linear_t)
739+
end
740+
668741
"""
669742
BSplineApprox(u, t, d, h, pVecType, knotVecType; extrapolate = false)
670743
@@ -738,7 +811,7 @@ struct BSplineApprox{uType, tType, pType, kType, cType, scType, T, N} <:
738811
end
739812

740813
function BSplineApprox(
741-
u, t, d, h, pVecType, knotVecType; extrapolate = false, assume_linear_t = 1e-2)
814+
u::AbstractVector, t, d, h, pVecType, knotVecType; extrapolate = false, assume_linear_t = 1e-2)
742815
u, t = munge_data(u, t)
743816
n = length(t)
744817
h < d + 1 && error("BSplineApprox needs at least d + 1, i.e. $(d+1) control points.")
@@ -827,6 +900,101 @@ function BSplineApprox(
827900
u, t, d, h, p, k, c, sc, pVecType, knotVecType, extrapolate, assume_linear_t)
828901
end
829902

903+
function BSplineApprox(
904+
u::AbstractArray{T, N}, t, d, h, pVecType, knotVecType; extrapolate = false,
905+
assume_linear_t = 1e-2) where {T, N}
906+
u, t = munge_data(u, t)
907+
n = length(t)
908+
h < d + 1 && error("BSplineApprox needs at least d + 1, i.e. $(d+1) control points.")
909+
s = zero(eltype(u))
910+
p = zero(t)
911+
k = zeros(eltype(t), h + d + 1)
912+
l = zeros(eltype(u), n - 1)
913+
p[1] = zero(eltype(t))
914+
p[end] = one(eltype(t))
915+
916+
ax_u = axes(u)[1:(end - 1)]
917+
918+
for i in 2:n
919+
s += √((t[i] - t[i - 1])^2 + sum((u[ax_u..., i] - u[ax_u..., i - 1]) .^ 2))
920+
l[i - 1] = s
921+
end
922+
if pVecType == :Uniform
923+
for i in 2:(n - 1)
924+
p[i] = p[1] + (i - 1) * (p[end] - p[1]) / (n - 1)
925+
end
926+
elseif pVecType == :ArcLen
927+
for i in 2:(n - 1)
928+
p[i] = p[1] + l[i - 1] / s * (p[end] - p[1])
929+
end
930+
end
931+
932+
lidx = 1
933+
ridx = length(k)
934+
while lidx <= (d + 1) && ridx >= (length(k) - d)
935+
k[lidx] = p[1]
936+
k[ridx] = p[end]
937+
lidx += 1
938+
ridx -= 1
939+
end
940+
941+
ps = zeros(eltype(t), n - 2)
942+
s = zero(eltype(t))
943+
for i in 2:(n - 1)
944+
s += p[i]
945+
ps[i - 1] = s
946+
end
947+
948+
if knotVecType == :Uniform
949+
# uniformly spaced knot vector
950+
# this method is not recommended because, if it is used with the chord length method for global interpolation,
951+
# the system of linear equations would be singular.
952+
for i in (d + 2):h
953+
k[i] = k[1] + (i - d - 1) // (h - d) * (k[end] - k[1])
954+
end
955+
elseif knotVecType == :Average
956+
# NOTE: verify that average method can be applied when size of k is less than size of p
957+
# average spaced knot vector
958+
idx = 1
959+
if d + 2 <= h
960+
k[d + 2] = 1 // d * ps[d]
961+
end
962+
for i in (d + 3):h
963+
k[i] = 1 // d * (ps[idx + d] - ps[idx])
964+
idx += 1
965+
end
966+
end
967+
# control points
968+
c = zeros(eltype(u), size(u)[1:(end - 1)]..., h)
969+
c[ax_u..., 1] = u[ax_u..., 1]
970+
c[ax_u..., end] = u[ax_u..., end]
971+
q = zeros(eltype(u), size(u)[1:(end - 1)]..., n)
972+
sc = zeros(eltype(t), n, h)
973+
for i in 1:n
974+
spline_coefficients!(view(sc, i, :), d, k, p[i])
975+
end
976+
for k in 2:(n - 1)
977+
q[ax_u..., k] = u[ax_u..., k] - sc[k, 1] * u[ax_u..., 1] -
978+
sc[k, h] * u[ax_u..., end]
979+
end
980+
Q = Array{eltype(u), N}(undef, size(u)[1:(end - 1)]..., h - 2)
981+
for i in 2:(h - 1)
982+
s = zeros(eltype(sc), size(u)[1:(end - 1)]...)
983+
for k in 2:(n - 1)
984+
s = s + sc[k, i] * q[ax_u..., k]
985+
end
986+
Q[ax_u..., i - 1] = s
987+
end
988+
sc = sc[2:(end - 1), 2:(h - 1)]
989+
M = transpose(sc) * sc
990+
Q = reshape(Q, prod(size(u)[1:(end - 1)]), :)
991+
P = (M \ Q')'
992+
P = reshape(P, size(u)[1:(end - 1)]..., :)
993+
c[ax_u..., 2:(end - 1)] = P
994+
sc = zeros(eltype(t), h)
995+
BSplineApprox(
996+
u, t, d, h, p, k, c, sc, pVecType, knotVecType, extrapolate, assume_linear_t)
997+
end
830998
"""
831999
CubicHermiteSpline(du, u, t; extrapolate = false, cache_parameters = false)
8321000

‎src/interpolation_methods.jl‎

Lines changed: 36 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -197,6 +197,25 @@ function _interpolate(A::BSplineInterpolation{<:AbstractVector{<:Number}},
197197
ucum
198198
end
199199

200+
function _interpolate(A::BSplineInterpolation{<:AbstractArray{T, N}},
201+
t::Number,
202+
iguess) where {T <: Number, N}
203+
ax_u = axes(A.u)[1:(end - 1)]
204+
t < A.t[1] && return A.u[ax_u..., 1]
205+
t > A.t[end] && return A.u[ax_u..., end]
206+
# change t into param [0 1]
207+
idx = get_idx(A, t, iguess)
208+
t = A.p[idx] + (t - A.t[idx]) / (A.t[idx + 1] - A.t[idx]) * (A.p[idx + 1] - A.p[idx])
209+
n = length(A.t)
210+
sc = t isa ForwardDiff.Dual ? zeros(eltype(t), n) : A.sc
211+
nonzero_coefficient_idxs = spline_coefficients!(sc, A.d, A.k, t)
212+
ucum = zeros(eltype(A.u), size(A.u)[1:(end - 1)]...)
213+
for i in nonzero_coefficient_idxs
214+
ucum = ucum + (sc[i] * A.c[ax_u..., i])
215+
end
216+
ucum
217+
end
218+
200219
# BSpline Curve Approx
201220
function _interpolate(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, iguess)
202221
t < A.t[1] && return A.u[1]
@@ -213,6 +232,23 @@ function _interpolate(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, i
213232
ucum
214233
end
215234

235+
function _interpolate(
236+
A::BSplineApprox{<:AbstractArray{T, N}}, t::Number, iguess) where {T <: Number, N}
237+
ax_u = axes(A.u)[1:(end - 1)]
238+
t < A.t[1] && return A.u[ax_u..., 1]
239+
t > A.t[end] && return A.u[ax_u..., end]
240+
# change t into param [0 1]
241+
idx = get_idx(A, t, iguess)
242+
t = A.p[idx] + (t - A.t[idx]) / (A.t[idx + 1] - A.t[idx]) * (A.p[idx + 1] - A.p[idx])
243+
sc = t isa ForwardDiff.Dual ? zeros(eltype(t), A.h) : A.sc
244+
nonzero_coefficient_idxs = spline_coefficients!(sc, A.d, A.k, t)
245+
ucum = zeros(eltype(A.u), size(A.u)[1:(end - 1)]...)
246+
for i in nonzero_coefficient_idxs
247+
ucum = ucum + (sc[i] * A.c[ax_u..., i])
248+
end
249+
ucum
250+
end
251+
216252
# Cubic Hermite Spline
217253
function _interpolate(
218254
A::CubicHermiteSpline{<:AbstractVector{<:Number}}, t::Number, iguess)

‎test/derivative_tests.jl‎

Lines changed: 38 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -186,6 +186,44 @@ end
186186
:Uniform,
187187
:Uniform],
188188
name = "BSpline Approx (Uniform, Uniform)")
189+
190+
f3d(t) = [sin(t) cos(t);
191+
0.0 cos(2t)]
192+
193+
t3d = 0.1:0.1:1.0 |> collect
194+
u3d = cat(f3d.(t3d)...; dims = 3)
195+
test_derivatives(BSplineInterpolation;
196+
args = [u3d, t3d,
197+
2,
198+
:Uniform,
199+
:Uniform],
200+
name = "BSpline Interpolation (Uniform, Uniform): AbstractArray"
201+
)
202+
203+
test_derivatives(BSplineInterpolation;
204+
args = [u3d, t3d,
205+
2,
206+
:ArcLen,
207+
:Average],
208+
name = "BSpline Interpolation (Arclen, Average): AbstractArray"
209+
)
210+
211+
test_derivatives(BSplineApprox;
212+
args = [u3d, t3d,
213+
3,
214+
4,
215+
:Uniform,
216+
:Uniform],
217+
name = "BSpline Approx (Uniform, Uniform): AbstractArray")
218+
219+
test_derivatives(BSplineApprox;
220+
args = [u3d, t3d,
221+
3,
222+
4,
223+
:ArcLen,
224+
:Average],
225+
name = "BSpline Approx (Arclen, Average): AbstractArray"
226+
)
189227
end
190228

191229
@testset "Cubic Hermite Spline" begin

0 commit comments

Comments
 (0)