refactor: improve BSplines with simpler parameterization

README.md

@@ -55,12 +55,11 @@ corresponding to `(u,t)` pairs.
   - `QuadraticSpline(u,t)` - A quadratic spline interpolation.
   - `CubicSpline(u,t)` - A cubic spline interpolation.
   - `AkimaInterpolation(u, t)` - Akima spline interpolation provides a smoothing effect and is computationally efficient.
-  - `BSplineInterpolation(u,t,d,pVec,knotVec)` - An interpolation B-spline. This is a B-spline which hits each of the data points. The argument choices are:
+  - `BSplineInterpolation(u,t,d,knotVec)` - An interpolation B-spline. This is a B-spline which hits each of the data points. The argument choices are:
 
       + `d` - degree of B-spline
-      + `pVec` - Symbol to Parameters Vector, `pVec = :Uniform` for uniform spaced parameters and `pVec = :ArcLen` for parameters generated by chord length method.
       + `knotVec` - Symbol to Knot Vector, `knotVec = :Uniform` for uniform knot vector, `knotVec = :Average` for average spaced knot vector.
-  - `BSplineApprox(u,t,d,h,pVec,knotVec)` - A regression B-spline which smooths the fitting curve. The argument choices are the same as the `BSplineInterpolation`, with the additional parameter `h<length(t)` which is the number of control points to use, with smaller `h` indicating more smoothing.
+  - `BSplineApprox(u,t,d,h,knotVec)` - A regression B-spline which smooths the fitting curve. The argument choices are the same as the `BSplineInterpolation`, with the additional parameter `h<length(t)` which is the number of control points to use, with smaller `h` indicating more smoothing.
   - `CubicHermiteSpline(du, u, t)` - A third order Hermite interpolation, which matches the values and first (`du`) order derivatives in the data points exactly.
   - `PCHIPInterpolation(u, t)` - a type of `CubicHermiteSpline` where the derivative values `du` are derived from the input data in such a way that the interpolation never overshoots the data.
   - `QuinticHermiteSpline(ddu, du, u, t)` - A fifth order Hermite interpolation, which matches the values and first (`du`) and second (`ddu`) order derivatives in the data points exactly.

docs/src/clarification.md

@@ -9,7 +9,7 @@ using DataInterpolations, Plots
 using DataInterpolations: derivative
 u = [2.0, 1.0, 5.0, 4.0, 5.0, 4.0]
 t = [0.0, 2.0, 3.5, 4.0, 5.0, 6.5]
-bspline = BSplineInterpolation(u, t, 0, :Uniform, :Uniform; extrapolation_left=ExtrapolationType.Extension,  extrapolation_right=ExtrapolationType.Extension)
+bspline = BSplineInterpolation(u, t, 0, :Uniform; extrapolation_left=ExtrapolationType.Extension,  extrapolation_right=ExtrapolationType.Extension)
 plot(bspline)
 ```
 
@@ -19,7 +19,7 @@ Thus, the plot for B-Spline interpolation does not appear the same as the plot f
 # Derivative behavior of quadratic B-Spline 
 
 ```@example interpclarity
-bspline = BSplineInterpolation(u, t, 2, :Uniform, :Uniform; extrapolation_left=ExtrapolationType.Extension,  extrapolation_right=ExtrapolationType.Extension)
+bspline = BSplineInterpolation(u, t, 2, :Uniform; extrapolation_left=ExtrapolationType.Extension,  extrapolation_right=ExtrapolationType.Extension)
 plot(t->derivative(bspline, t))
 ```
 

docs/src/index.md

@@ -25,12 +25,11 @@ corresponding to `(u,t)` pairs.
   - `QuadraticSpline(u,t)` - A quadratic spline interpolation.
   - `CubicSpline(u,t)` - A cubic spline interpolation.
   - `AkimaInterpolation(u, t)` - Akima spline interpolation provides a smoothing effect and is computationally efficient.
-  - `BSplineInterpolation(u,t,d,pVec,knotVec)` - An interpolation B-spline. This is a B-spline that hits each of the data points. The argument choices are:
+  - `BSplineInterpolation(u,t,d,knotVec)` - An interpolation B-spline. This is a B-spline that hits each of the data points. The argument choices are:
 
       + `d` - degree of B-spline
-      + `pVec` - Symbol to Parameters Vector, `pVec = :Uniform` for uniformly spaced parameters, and `pVec = :ArcLen` for parameters generated by the chord length method.
       + `knotVec` - Symbol to Knot Vector, `knotVec = :Uniform` for uniform knot vector, `knotVec = :Average` for average spaced knot vector.
-  - `BSplineApprox(u,t,d,h,pVec,knotVec)` - A regression B-spline which smooths the fitting curve. The argument choices are the same as the `BSplineInterpolation`, with the additional parameter `h<length(t)` which is the number of control points to use, with smaller `h` indicating more smoothing.
+  - `BSplineApprox(u,t,d,h,knotVec)` - A regression B-spline which smooths the fitting curve. The argument choices are the same as the `BSplineInterpolation`, with the additional parameter `h<length(t)` which is the number of control points to use, with smaller `h` indicating more smoothing.
   - `CubicHermiteSpline(du, u, t)` - A third order Hermite interpolation, which matches the values and first (`du`) order derivatives in the data points exactly.
   - `PCHIPInterpolation(u, t)` - a type of `CubicHermiteSpline` where the derivative values `du` are derived from the input data in such a way that the interpolation never overshoots the data.
   - `QuinticHermiteSpline(ddu, du, u, t)` - a fifth order Hermite interpolation, which matches the values and first (`du`) and second (`ddu`) order derivatives in the data points exactly.

docs/src/methods.md

@@ -124,10 +124,10 @@ that every data point is taken into account for each point of the curve.
 The interpolating B-spline is the version which hits each of the points. This
 method is described in more detail [here](https://pages.mtu.edu/%7Eshene/COURSES/cs3621/NOTES/INT-APP/CURVE-INT-global.html).
 Let's plot a cubic B-spline (3rd order). Since the data points are not close to
-uniformly spaced, we will use the `:ArcLen` and `:Average` choices:
+uniformly spaced, we will use the `:Average` knot vector:
 
 ```@example tutorial
-A = BSplineInterpolation(u, t, 3, :ArcLen, :Average)
+A = BSplineInterpolation(u, t, 3, :Average)
 plot(A)
 ```
 
@@ -138,7 +138,7 @@ is a least square approximation. This has a natural effect of smoothing the
 data. For example, if we use 4 control points, we get the result:
 
 ```@example tutorial
-A = BSplineApprox(u, t, 3, 4, :ArcLen, :Average)
+A = BSplineApprox(u, t, 3, 4, :Average)
 plot(A)
 ```
 

src/derivatives.jl

@@ -249,94 +249,102 @@ function _derivative(A::CubicSpline{<:AbstractVector}, t::Number, iguess)
 end
 
 function _derivative(A::BSplineInterpolation{<:AbstractVector{<:Number}}, t::Number, iguess)
-    # change t into param [0 1]
-    t < A.t[1] && return zero(A.u[1])
-    t > A.t[end] && return zero(A.u[end])
-    idx = get_idx(A, t, iguess)
+    if A.d == 0
+        return isempty(searchsorted(A.t, t)) ? zero(A.u[1]) : typed_nan(A.u)
+    end
     n = length(A.t)
-    scale = (A.p[idx + 1] - A.p[idx]) / (A.t[idx + 1] - A.t[idx])
-    t_ = A.p[idx] + (t - A.t[idx]) * scale
     sc = t isa ForwardDiff.Dual ? zeros(eltype(t), n) : A.sc
-    spline_coefficients!(sc, A.d - 1, A.k, t_)
+    spline_coefficients!(sc, A.d - 1, A.k, t)
     ducum = zero(eltype(A.u))
     if t == A.t[1]
-        ducum = (A.c[2] - A.c[1]) / (A.k[A.d + 2])
+        denom = A.k[A.d + 2] - A.k[2]
+        ducum = denom != 0 ? (A.c[2] - A.c[1]) / denom : zero(eltype(A.u))
     else
         for i in 1:(n - 1)
-            ducum += sc[i + 1] * (A.c[i + 1] - A.c[i]) / (A.k[i + A.d + 1] - A.k[i + 1])
+            denom = A.k[i + A.d + 1] - A.k[i + 1]
+            if denom != 0
+                ducum += sc[i + 1] * (A.c[i + 1] - A.c[i]) / denom
+            end
         end
     end
-    return ducum * A.d * scale
+    return ducum * A.d
 end
 
 function _derivative(
         A::BSplineInterpolation{<:AbstractArray{<:Number}}, t::Number, iguess
     )
-    # change t into param [0 1]
+    if A.d == 0
+        return isempty(searchsorted(A.t, t)) ? zero(A.u[:, 1]) :
+            typed_nan(A.u) .* A.u[:, 1]
+    end
     ax_u = axes(A.u)[1:(end - 1)]
-    t < A.t[1] && return zeros(size(A.u)[1:(end - 1)]...)
-    t > A.t[end] && return zeros(size(A.u)[1:(end - 1)]...)
-    idx = get_idx(A, t, iguess)
     n = length(A.t)
-    scale = (A.p[idx + 1] - A.p[idx]) / (A.t[idx + 1] - A.t[idx])
-    t_ = A.p[idx] + (t - A.t[idx]) * scale
     sc = t isa ForwardDiff.Dual ? zeros(eltype(t), n) : A.sc
-    spline_coefficients!(sc, A.d - 1, A.k, t_)
+    spline_coefficients!(sc, A.d - 1, A.k, t)
     ducum = zeros(size(A.u)[1:(end - 1)]...)
     if t == A.t[1]
-        ducum = (A.c[ax_u..., 2] - A.c[ax_u..., 1]) / (A.k[A.d + 2])
+        denom = A.k[A.d + 2] - A.k[2]
+        if denom != 0
+            ducum = (A.c[ax_u..., 2] - A.c[ax_u..., 1]) / denom
+        end
     else
         for i in 1:(n - 1)
-            ducum = ducum +
-                sc[i + 1] * (A.c[ax_u..., i + 1] - A.c[ax_u..., i]) /
-                (A.k[i + A.d + 1] - A.k[i + 1])
+            denom = A.k[i + A.d + 1] - A.k[i + 1]
+            if denom != 0
+                ducum = ducum +
+                    sc[i + 1] * (A.c[ax_u..., i + 1] - A.c[ax_u..., i]) / denom
+            end
         end
     end
-    return ducum * A.d * scale
+    return ducum * A.d
 end
 # BSpline Curve Approx
 function _derivative(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, iguess)
-    # change t into param [0 1]
-    t < A.t[1] && return zero(A.u[1])
-    t > A.t[end] && return zero(A.u[end])
-    idx = get_idx(A, t, iguess)
-    scale = (A.p[idx + 1] - A.p[idx]) / (A.t[idx + 1] - A.t[idx])
-    t_ = A.p[idx] + (t - A.t[idx]) * scale
+    if A.d == 0
+        return isempty(searchsorted(A.t, t)) ? zero(A.u[1]) : typed_nan(A.u)
+    end
     sc = t isa ForwardDiff.Dual ? zeros(eltype(t), A.h) : A.sc
-    spline_coefficients!(sc, A.d - 1, A.k, t_)
+    spline_coefficients!(sc, A.d - 1, A.k, t)
     ducum = zero(eltype(A.u))
     if t == A.t[1]
-        ducum = (A.c[2] - A.c[1]) / (A.k[A.d + 2])
+        denom = A.k[A.d + 2] - A.k[2]
+        ducum = denom != 0 ? (A.c[2] - A.c[1]) / denom : zero(eltype(A.u))
     else
         for i in 1:(A.h - 1)
-            ducum += sc[i + 1] * (A.c[i + 1] - A.c[i]) / (A.k[i + A.d + 1] - A.k[i + 1])
+            denom = A.k[i + A.d + 1] - A.k[i + 1]
+            if denom != 0
+                ducum += sc[i + 1] * (A.c[i + 1] - A.c[i]) / denom
+            end
         end
     end
-    return ducum * A.d * scale
+    return ducum * A.d
 end
 
 function _derivative(
         A::BSplineApprox{<:AbstractArray{<:Number}}, t::Number, iguess
     )
-    # change t into param [0 1]
+    if A.d == 0
+        return isempty(searchsorted(A.t, t)) ? zero(A.u[:, 1]) :
+            typed_nan(A.u) .* A.u[:, 1]
+    end
     ax_u = axes(A.u)[1:(end - 1)]
-    t < A.t[1] && return zeros(size(A.u)[1:(end - 1)]...)
-    t > A.t[end] && return zeros(size(A.u)[1:(end - 1)]...)
-    idx = get_idx(A, t, iguess)
-    scale = (A.p[idx + 1] - A.p[idx]) / (A.t[idx + 1] - A.t[idx])
-    t_ = A.p[idx] + (t - A.t[idx]) * scale
     sc = t isa ForwardDiff.Dual ? zeros(eltype(t), A.h) : A.sc
-    spline_coefficients!(sc, A.d - 1, A.k, t_)
+    spline_coefficients!(sc, A.d - 1, A.k, t)
     ducum = zeros(size(A.u)[1:(end - 1)]...)
     if t == A.t[1]
-        ducum = (A.c[ax_u..., 2] - A.c[ax_u..., 1]) / (A.k[A.d + 2])
+        denom = A.k[A.d + 2] - A.k[2]
+        if denom != 0
+            ducum = (A.c[ax_u..., 2] - A.c[ax_u..., 1]) / denom
+        end
     else
         for i in 1:(A.h - 1)
-            ducum += sc[i + 1] * (A.c[ax_u..., i + 1] - A.c[ax_u..., i]) /
-                (A.k[i + A.d + 1] - A.k[i + 1])
+            denom = A.k[i + A.d + 1] - A.k[i + 1]
+            if denom != 0
+                ducum += sc[i + 1] * (A.c[ax_u..., i + 1] - A.c[ax_u..., i]) / denom
+            end
         end
     end
-    return ducum * A.d * scale
+    return ducum * A.d
 end
 # Cubic Hermite Spline
 function _derivative(

src/integrals.jl

@@ -287,11 +287,86 @@ end
 function _integral(A::LagrangeInterpolation, idx::Number, t1::Number, t2::Number)
     throw(IntegralNotFoundError())
 end
-function _integral(A::BSplineInterpolation, idx::Number, t1::Number, t2::Number)
-    throw(IntegralNotFoundError())
+# Evaluate the antiderivative of a B-spline at a point.
+# The antiderivative of a degree-d B-spline is a degree-(d+1) B-spline with
+# extended knot vector and coefficients derived from the original.
+function _bspline_antiderivative_val(c, d, k, t_eval)
+    nc = length(c)
+    dp1 = d + 1
+    # Antiderivative coefficients: C[1] = 0, C[i+1] = C[i] + c[i]*(k[i+d+1]-k[i])/(d+1)
+    T = promote_type(eltype(c), eltype(k), typeof(t_eval))
+    C = zeros(T, nc + 1)
+    for i in 1:nc
+        C[i + 1] = C[i] + c[i] * (k[i + dp1] - k[i]) / dp1
+    end
+    # Extended knot vector: prepend k[1], append k[end]
+    nk = length(k)
+    k_ext = zeros(eltype(k), nk + 2)
+    k_ext[1] = k[1]
+    for i in 1:nk
+        k_ext[i + 1] = k[i]
+    end
+    k_ext[end] = k[end]
+    # Evaluate degree-(d+1) B-spline with coefficients C on k_ext
+    sc = zeros(T, nc + 1)
+    nonzero = spline_coefficients!(sc, dp1, k_ext, t_eval)
+    result = zero(T)
+    for i in nonzero
+        result += sc[i] * C[i]
+    end
+    return result
 end
-function _integral(A::BSplineApprox, idx::Number, t1::Number, t2::Number)
-    throw(IntegralNotFoundError())
+
+function _integral(
+        A::BSplineInterpolation{<:AbstractVector{<:Number}}, idx::Number, t1::Number, t2::Number
+    )
+    return _bspline_antiderivative_val(A.c, A.d, A.k, t2) -
+        _bspline_antiderivative_val(A.c, A.d, A.k, t1)
+end
+function _integral(A::BSplineApprox{<:AbstractVector{<:Number}}, idx::Number, t1::Number, t2::Number)
+    return _bspline_antiderivative_val(A.c, A.d, A.k, t2) -
+        _bspline_antiderivative_val(A.c, A.d, A.k, t1)
+end
+
+# Override integral to bypass the hasfield(:I) check in the generic method.
+# The antiderivative is computed on the fly, so no cached I field is needed.
+const _BSplineTypes = Union{
+    BSplineInterpolation{<:AbstractVector{<:Number}},
+    BSplineApprox{<:AbstractVector{<:Number}},
+}
+
+function integral(A::_BSplineTypes, t::Number)
+    return integral(A, first(A.t), t)
+end
+
+function integral(A::_BSplineTypes, t1::Number, t2::Number)
+    t1 == t2 && return zero(eltype(A.u))
+    t1 > t2 && return -integral(A, t2, t1)
+
+    total = zero(eltype(A.u))
+    lo = t1
+    hi = t2
+
+    if lo < first(A.t)
+        if hi <= first(A.t)
+            return _extrapolate_integral_left(A, lo) - _extrapolate_integral_left(A, hi)
+        end
+        total += _extrapolate_integral_left(A, lo)
+        lo = first(A.t)
+    end
+
+    if hi > last(A.t)
+        if lo >= last(A.t)
+            return _extrapolate_integral_right(A, hi) - _extrapolate_integral_right(A, lo)
+        end
+        total += _extrapolate_integral_right(A, hi)
+        hi = last(A.t)
+    end
+
+    total += _bspline_antiderivative_val(A.c, A.d, A.k, hi) -
+        _bspline_antiderivative_val(A.c, A.d, A.k, lo)
+
+    return total
 end
 
 # Cubic Hermite Spline