NFFT  3.3.0
bspline.c
1 /*
2  * Copyright (c) 2002, 2015 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 /* $Id: util.c 3483 2010-04-23 19:02:34Z keiner $ */
20 
21 #include "infft.h"
22 
23 static inline void bspline_help(const INT k, const R x, R *scratch, const INT j,
24  const INT ug, const INT og, const INT r)
25 {
26  INT i; /* Row index of the de Boor scheme */
27  INT idx; /* Index in scratch */
28  R a; /* Alpha of de Boor scheme */
29 
30  /* computation of one column */
31  for (i = og + r - k + 1, idx = og; idx >= ug; i--, idx--)
32  {
33  a = (x - (R)i) / ((R)(k - j));
34  scratch[idx] = (K(1.0) - a) * scratch[idx - 1] + a * scratch[idx];
35  }
36 }
37 
38 /* Evaluate the cardinal B-Spline B_{n-1} supported on [0,n]. */
39 R Y(bspline)(const INT k, const R _x, R *scratch)
40 {
41  const R kk = (R)k;
42  R result_value;
43  INT r;
44  INT g1, g2; /* boundaries */
45  INT j, idx, ug, og; /* indices */
46  R a; /* Alpha of de Boor scheme*/
47  R x = _x;
48 
49  result_value = K(0.0);
50 
51  if (K(0.0) < x && x < kk)
52  {
53  /* Exploit symmetry around k/2, maybe. */
54  if ( (kk - x) < x)
55  {
56  x = kk - x;
57  }
58 
59  r = (INT)LRINT(CEIL(x) - K(1.0));
60 
61  /* Explicit case for first interval. */
62  if (r == 0)
63  {
64  result_value = K(1.0);
65  for (j = 0; j < k - 1; j++)
66  result_value *= x/((R)(j+1));
67  return result_value;
68  }
69 
70 
71  for (idx = 0; idx < k; idx++)
72  scratch[idx] = K(0.0);
73 
74  scratch[k-r-1] = K(1.0);
75 
76  /* Bounds of the algorithm. */
77  g1 = r;
78  g2 = k - 1 - r;
79  ug = g2;
80 
81  /* g1 <= g2 */
82 
83  for (j = 1, og = g2 + 1; j <= g1; j++, og++)
84  {
85  a = (x + (R)(k - r - og - 1)) / ((R)(k - j));
86  scratch[og] = (K(1.0) - a) * scratch[og-1];
87  bspline_help(k, x, scratch, j, ug + 1, og - 1, r);
88  a = (x + (R)(k - r - ug - 1)) / ((R)(k - j));
89  scratch[ug] = a * scratch[ug];
90  }
91 
92  for (og-- ; j <= g2; j++)
93  {
94  bspline_help(k, x, scratch, j, ug + 1, og, r);
95  a = (x + (R)(k - r - ug - 1)) / ((R)(k - j));
96  scratch[ug] = a * scratch[ug];
97  }
98 
99  for(; j < k; j++)
100  {
101  ug++;
102  bspline_help(k, x, scratch, j, ug, og, r);
103  }
104 
105  result_value = scratch[k-1];
106  }
107 
108  return(result_value);
109 }