NFFT  3.3.1
bspline.c
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 }