1 /*
2 * Copyright (c) 2002, 2017 Jens Keiner, Stefan Kunis, Daniel Potts
3 *
4 * This program is free software; you can redistribute it and/or modify it under
5 * the terms of the GNU General Public License as published by the Free Software
6 * Foundation; either version 2 of the License, or (at your option) any later
7 * version.
8 *
9 * This program is distributed in the hope that it will be useful, but WITHOUT
10 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
12 * details.
13 *
14 * You should have received a copy of the GNU General Public License along with
15 * this program; if not, write to the Free Software Foundation, Inc., 51
16 * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17 */
18
19 #include "infft.h"
20
bspline_help(const INT k,const R x,R * scratch,const INT j,const INT ug,const INT og,const INT r)21 static inline void bspline_help(const INT k, const R x, R *scratch, const INT j,
22 const INT ug, const INT og, const INT r)
23 {
24 INT i; /* Row index of the de Boor scheme */
25 INT idx; /* Index in scratch */
26 R a; /* Alpha of de Boor scheme */
27
28 /* computation of one column */
29 for (i = og + r - k + 1, idx = og; idx >= ug; i--, idx--)
30 {
31 a = (x - (R)i) / ((R)(k - j));
32 scratch[idx] = (K(1.0) - a) * scratch[idx - 1] + a * scratch[idx];
33 }
34 }
35
36 /* Evaluate the cardinal B-Spline B_{n-1} supported on [0,n]. */
Y(bsplines)37 R Y(bsplines)(const INT k, const R _x)
38 {
39 const R kk = (R)k;
40 R result_value;
41 INT r;
42 INT g1, g2; /* boundaries */
43 INT j, idx, ug, og; /* indices */
44 R a; /* Alpha of de Boor scheme*/
45 R x = _x;
46 R scratch[k];
47
48 result_value = K(0.0);
49
50 if (K(0.0) < x && x < kk)
51 {
52 /* Exploit symmetry around k/2, maybe. */
53 if ( (kk - x) < x)
54 {
55 x = kk - x;
56 }
57
58 r = (INT)LRINT(CEIL(x) - K(1.0));
59
60 /* Do not use the explicit formula x^k / k! for first interval! De Boor's
61 * algorithm is more accurate. See https://github.com/NFFT/nfft/issues/16.
62 */
63
64 for (idx = 0; idx < k; idx++)
65 scratch[idx] = K(0.0);
66
67 scratch[k-r-1] = K(1.0);
68
69 /* Bounds of the algorithm. */
70 g1 = r;
71 g2 = k - 1 - r;
72 ug = g2;
73
74 /* g1 <= g2 */
75
76 for (j = 1, og = g2 + 1; j <= g1; j++, og++)
77 {
78 a = (x + (R)(k - r - og - 1)) / ((R)(k - j));
79 scratch[og] = (K(1.0) - a) * scratch[og-1];
80 bspline_help(k, x, scratch, j, ug + 1, og - 1, r);
81 a = (x + (R)(k - r - ug - 1)) / ((R)(k - j));
82 scratch[ug] = a * scratch[ug];
83 }
84
85 for (og-- ; j <= g2; j++)
86 {
87 bspline_help(k, x, scratch, j, ug + 1, og, r);
88 a = (x + (R)(k - r - ug - 1)) / ((R)(k - j));
89 scratch[ug] = a * scratch[ug];
90 }
91
92 for(; j < k; j++)
93 {
94 ug++;
95 bspline_help(k, x, scratch, j, ug, og, r);
96 }
97
98 result_value = scratch[k-1];
99 }
100
101 return(result_value);
102 }
103