![]() |
NFFT
3.3.1
|
00001 /* 00002 * Copyright (c) 2002, 2016 Jens Keiner, Stefan Kunis, Daniel Potts 00003 * 00004 * This program is free software; you can redistribute it and/or modify it under 00005 * the terms of the GNU General Public License as published by the Free Software 00006 * Foundation; either version 2 of the License, or (at your option) any later 00007 * version. 00008 * 00009 * This program is distributed in the hope that it will be useful, but WITHOUT 00010 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 00011 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 00012 * details. 00013 * 00014 * You should have received a copy of the GNU General Public License along with 00015 * this program; if not, write to the Free Software Foundation, Inc., 51 00016 * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 00017 */ 00018 00019 #include "infft.h" 00020 00021 static inline void bspline_help(const INT k, const R x, R *scratch, const INT j, 00022 const INT ug, const INT og, const INT r) 00023 { 00024 INT i; /* Row index of the de Boor scheme */ 00025 INT idx; /* Index in scratch */ 00026 R a; /* Alpha of de Boor scheme */ 00027 00028 /* computation of one column */ 00029 for (i = og + r - k + 1, idx = og; idx >= ug; i--, idx--) 00030 { 00031 a = (x - (R)i) / ((R)(k - j)); 00032 scratch[idx] = (K(1.0) - a) * scratch[idx - 1] + a * scratch[idx]; 00033 } 00034 } 00035 00036 /* Evaluate the cardinal B-Spline B_{n-1} supported on [0,n]. */ 00037 R Y(bsplines)(const INT k, const R _x) 00038 { 00039 const R kk = (R)k; 00040 R result_value; 00041 INT r; 00042 INT g1, g2; /* boundaries */ 00043 INT j, idx, ug, og; /* indices */ 00044 R a; /* Alpha of de Boor scheme*/ 00045 R x = _x; 00046 R scratch[k]; 00047 00048 result_value = K(0.0); 00049 00050 if (K(0.0) < x && x < kk) 00051 { 00052 /* Exploit symmetry around k/2, maybe. */ 00053 if ( (kk - x) < x) 00054 { 00055 x = kk - x; 00056 } 00057 00058 r = (INT)LRINT(CEIL(x) - K(1.0)); 00059 00060 /* Do not use the explicit formula x^k / k! for first interval! De Boor's 00061 * algorithm is more accurate. See https://github.com/NFFT/nfft/issues/16. 00062 */ 00063 00064 for (idx = 0; idx < k; idx++) 00065 scratch[idx] = K(0.0); 00066 00067 scratch[k-r-1] = K(1.0); 00068 00069 /* Bounds of the algorithm. */ 00070 g1 = r; 00071 g2 = k - 1 - r; 00072 ug = g2; 00073 00074 /* g1 <= g2 */ 00075 00076 for (j = 1, og = g2 + 1; j <= g1; j++, og++) 00077 { 00078 a = (x + (R)(k - r - og - 1)) / ((R)(k - j)); 00079 scratch[og] = (K(1.0) - a) * scratch[og-1]; 00080 bspline_help(k, x, scratch, j, ug + 1, og - 1, r); 00081 a = (x + (R)(k - r - ug - 1)) / ((R)(k - j)); 00082 scratch[ug] = a * scratch[ug]; 00083 } 00084 00085 for (og-- ; j <= g2; j++) 00086 { 00087 bspline_help(k, x, scratch, j, ug + 1, og, r); 00088 a = (x + (R)(k - r - ug - 1)) / ((R)(k - j)); 00089 scratch[ug] = a * scratch[ug]; 00090 } 00091 00092 for(; j < k; j++) 00093 { 00094 ug++; 00095 bspline_help(k, x, scratch, j, ug, og, r); 00096 } 00097 00098 result_value = scratch[k-1]; 00099 } 00100 00101 return(result_value); 00102 }