@@ -5,14 +5,14 @@ Cubic{BC<:Flag}(::BC) = Cubic{BC}()
5
5
Assuming uniform knots with spacing 1, the `i`th piece of cubic spline
6
6
implemented here is defined as follows.
7
7
8
- y_i(x) = ϕm p(x-i) + ϕ q(x-i) + ϕp q(1- (x-i)) + ϕpp p(1 - (x-i))
8
+ y_i(x) = cm p(x-i) + c q(x-i) + cp q(1- (x-i)) + cpp p(1 - (x-i))
9
9
10
10
where
11
11
12
- p(☆ ) = 1/6 * (1-☆ )^3
13
- q(☆ ) = 2/3 - ☆ ^2 + 1/2 ☆ ^3
12
+ p(δx ) = 1/6 * (1-δx )^3
13
+ q(δx ) = 2/3 - δx ^2 + 1/2 δx ^3
14
14
15
- and the values `ϕX ` for `X ⋹ {m, _, p, pp}` are the pre-filtered coefficients.
15
+ and the values `cX ` for `X ∈ {m, _, p, pp}` are the pre-filtered coefficients.
16
16
17
17
For future reference, this expands out to the following polynomial:
18
18
@@ -24,6 +24,11 @@ When we derive boundary conditions we will use derivatives `y_0'(x)` and
24
24
"""
25
25
Cubic
26
26
27
+ """
28
+ `define_indices_d` for a cubic b-spline calculates `ix_d = floor(x_d)` and
29
+ `fx_d = x_d - ix_d` (corresponding to `i` `and `δx` in the docstring for
30
+ `Cubic`), as well as auxiliary quantities `ixm_d`, `ixp_d` and `ixpp_d`
31
+ """
27
32
function define_indices_d {BC} (:: Type{BSpline{Cubic{BC}}} , d, pad)
28
33
symix, symixm, symixp = symbol (" ix_" ,d), symbol (" ixm_" ,d), symbol (" ixp_" ,d)
29
34
symixpp, symx, symfx = symbol (" ixpp_" ,d), symbol (" x_" ,d), symbol (" fx_" ,d)
@@ -39,6 +44,14 @@ function define_indices_d{BC}(::Type{BSpline{Cubic{BC}}}, d, pad)
39
44
end
40
45
end
41
46
47
+ """
48
+ `define_indices_d` for a cubic, periodic b-spline calculates `ix_d = floor(x_d)`
49
+ and `fx_d = x_d - ix_d` (corresponding to `i` and `δx` in the docstring entry
50
+ for `Cubic`), as well as auxiliary quantities `ixm_d`, `ixp_d` and `ixpp_d`.
51
+
52
+ If any `ixX_d` for `x ∈ {m, p, pp}` (note: not `c_d`) should fall outside of
53
+ the data interval, they wrap around.
54
+ """
42
55
function define_indices_d (:: Type{BSpline{Cubic{Periodic}}} , d, pad)
43
56
symix, symixm, symixp = symbol (" ix_" ,d), symbol (" ixm_" ,d), symbol (" ixp_" ,d)
44
57
symixpp, symx, symfx = symbol (" ixpp_" ,d), symbol (" x_" ,d), symbol (" fx_" ,d)
@@ -55,15 +68,16 @@ padding{BC<:Flag}(::Type{BSpline{Cubic{BC}}}) = Val{1}()
55
68
padding (:: Type{BSpline{Cubic{Periodic}}} ) = Val {0} ()
56
69
57
70
"""
58
- In this function we assume that `fx_d = x-ix_d` and we produce `cX_d` for
59
- `X ⋹ {m, _, p, pp}` such that
71
+ In `coefficients` for a cubic b-spline we assume that `fx_d = x-ix_d`
72
+ and we define `cX_d` for `X ⋹ {m, _, p, pp}` such that
60
73
61
74
cm_d = p(fx_d)
62
75
c_d = q(fx_d)
63
76
cp_d = q(1-fx_d)
64
77
cpp_d = p(1-fx_d)
65
78
66
- where `p` and `q` are defined in the docstring entry for `Cubic`.
79
+ where `p` and `q` are defined in the docstring entry for `Cubic`, and
80
+ `fx_d` in the docstring entry for `define_indices_d`.
67
81
"""
68
82
function coefficients {C<:Cubic} (:: Type{BSpline{C}} , N, d)
69
83
symm, sym = symbol (" cm_" ,d), symbol (" c_" ,d)
@@ -82,15 +96,16 @@ function coefficients{C<:Cubic}(::Type{BSpline{C}}, N, d)
82
96
end
83
97
84
98
"""
85
- In this function we assume that `fx_d = x-ix_d` and we produce `cX_d` for
86
- `X ⋹ {m, _, p, pp}` such that
99
+ In `gradient_coefficients` for a cubic b-spline we assume that `fx_d = x-ix_d`
100
+ and we define `cX_d` for `X ⋹ {m, _, p, pp}` such that
87
101
88
102
cm_d = p'(fx_d)
89
103
c_d = q'(fx_d)
90
104
cp_d = q'(1-fx_d)
91
105
cpp_d = p'(1-fx_d)
92
106
93
- where `p` and `q` are defined in the docstring for `Cubic`.
107
+ where `p` and `q` are defined in the docstring entry for `Cubic`, and
108
+ `fx_d` in the docstring entry for `define_indices_d`.
94
109
"""
95
110
function gradient_coefficients {C<:Cubic} (:: Type{BSpline{C}} , d)
96
111
symm, sym, symp, sympp = symbol (" cm_" ,d), symbol (" c_" ,d), symbol (" cp_" ,d), symbol (" cpp_" ,d)
147
162
# Prefiltering #
148
163
# ------------ #
149
164
165
+ """
166
+ `Cubic`: continuity in function value, first and second derivatives yields
167
+
168
+ 2/3 1/6
169
+ 1/6 2/3 1/6
170
+ 1/6 2/3 1/6
171
+ ⋱ ⋱ ⋱
172
+ """
150
173
function inner_system_diags {T,C<:Cubic} (:: Type{T} , n:: Int , :: Type{C} )
151
174
du = fill (convert (T, SimpleRatio (1 , 6 )), n- 1 )
152
175
d = fill (convert (T, SimpleRatio (2 , 3 )), n)
@@ -155,8 +178,8 @@ function inner_system_diags{T,C<:Cubic}(::Type{T}, n::Int, ::Type{C})
155
178
end
156
179
157
180
"""
158
- `Flat` `OnGrid` amounts to setting `y_0 '(x) = 0` at `x = 0`. Applying this
159
- condition gives:
181
+ `Cubic{ Flat} ` `OnGrid` amounts to setting `y_1 '(x) = 0` at `x = 1`.
182
+ Applying this condition yields
160
183
161
184
-cm + cp = 0
162
185
"""
@@ -173,10 +196,18 @@ function prefiltering_system{T,TC}(::Type{T}, ::Type{TC}, n::Int,
173
196
end
174
197
175
198
"""
176
- `Flat` `OnCell` amounts to setting `y_0'(x) = 0` at `x = -1/2`. Applying this
177
- condition gives:
199
+ `Cubic{Flat}`, `OnCell` amounts to setting `y_1'(x) = 0` at `x = 1/2`.
200
+ Applying this condition yields
201
+
202
+ -9/8 cm + 11/8 c - 3/8 cp + 1/8 cpp = 0
203
+
204
+ or, equivalently,
178
205
179
206
-9 cm + 11 c -3 cp + 1 cpp = 0
207
+
208
+ (Note that we use `y_1'(x)` although it is strictly not valid in this domain; if we
209
+ were to use `y_0'(x)` we would have to introduce new coefficients, so that would not
210
+ close the system. Instead, we extend the outermost polynomial for an extra half-cell.)
180
211
"""
181
212
function prefiltering_system {T,TC} (:: Type{T} , :: Type{TC} , n:: Int ,
182
213
:: Type{Cubic{Flat}} , :: Type{OnCell} )
@@ -198,10 +229,14 @@ function prefiltering_system{T,TC}(::Type{T}, ::Type{TC}, n::Int,
198
229
end
199
230
200
231
"""
201
- `Line` `OnCell` amounts to setting `y_0 ''(x) = 0` at `x = - 1/2`. Applying this
202
- condition gives:
232
+ `Cubic{ Line} ` `OnCell` amounts to setting `y_1 ''(x) = 0` at `x = 1/2`.
233
+ Applying this condition yields
203
234
204
235
3 cm -7 c + 5 cp -1 cpp = 0
236
+
237
+ (Note that we use `y_1'(x)` although it is strictly not valid in this domain; if we
238
+ were to use `y_0'(x)` we would have to introduce new coefficients, so that would not
239
+ close the system. Instead, we extend the outermost polynomial for an extra half-cell.)
205
240
"""
206
241
function prefiltering_system {T,TC} (:: Type{T} , :: Type{TC} , n:: Int ,
207
242
:: Type{Cubic{Line}} , :: Type{OnCell} )
@@ -223,13 +258,10 @@ function prefiltering_system{T,TC}(::Type{T}, ::Type{TC}, n::Int,
223
258
end
224
259
225
260
"""
226
- `Line` `OnGrid` amounts to setting `y_0 ''(x) = 0` at `x = 0 `. Applying this
261
+ `Cubic{ Line} ` `OnGrid` amounts to setting `y_1 ''(x) = 0` at `x = 1 `. Applying this
227
262
condition gives:
228
263
229
- 1 cm -2 c + 1 cp = 0 = 0
230
-
231
- This is the same system as `Quadratic{Line}`, `OnGrid` so we reuse the
232
- implementation
264
+ 1 cm -2 c + 1 cp = 0
233
265
"""
234
266
function prefiltering_system {T,TC} (:: Type{T} , :: Type{TC} , n:: Int ,
235
267
:: Type{Cubic{Line}} , :: Type{OnGrid} )
@@ -247,6 +279,14 @@ function prefiltering_system{T,TC}(::Type{T}, ::Type{TC}, n::Int,
247
279
Woodbury (lufact! (Tridiagonal (dl, d, du), Val{false }), specs... ), zeros (TC, n)
248
280
end
249
281
282
+ """
283
+ `Cubic{Periodic}` `OnGrid` closes the system by looking at the coefficients themselves
284
+ as periodic, yielding
285
+
286
+ c0 = c(N+1)
287
+
288
+ where `N` is the number of data points.
289
+ """
250
290
function prefiltering_system {T,TC,GT<:GridType} (:: Type{T} , :: Type{TC} , n:: Int ,
251
291
:: Type{Cubic{Periodic}} , :: Type{GT} )
252
292
dl, d, du = inner_system_diags (T,n,Cubic{Periodic})
@@ -260,14 +300,11 @@ function prefiltering_system{T,TC,GT<:GridType}(::Type{T}, ::Type{TC}, n::Int,
260
300
end
261
301
262
302
"""
263
- The free boundary condition for either `OnGrid` or `OnCell` makes sure the
264
- interpoland has a continuous third derivative at the second-to-outermost cell
265
- boundary: `y_0'''(1) = y_1'''(1)` and `y_{n-1}'''(n) = y_n'''(n)`. Applying this
266
- condition gives:
303
+ `Cubic{Free}` `OnGrid` and `Cubic{Free}` `OnCell` amount to requiring an extra
304
+ continuous derivative at the second-to-last cell boundary; this means
305
+ `y_1'''(2) = y_2'''(2)`, yielding
267
306
268
307
1 cm -3 c + 3 cp -1 cpp = 0
269
-
270
- This is the same system as `Quadratic{Free}` so we reuse the implementation
271
308
"""
272
309
function prefiltering_system {T,TC,GT<:GridType} (:: Type{T} , :: Type{TC} , n:: Int ,
273
310
:: Type{Cubic{Free}} , :: Type{GT} )
0 commit comments