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)
31 for (i = og + r - k + 1, idx = og; idx >= ug; i--, idx--)
33 a = (x - (R)i) / ((R)(k - j));
34 scratch[idx] = (K(1.0) - a) * scratch[idx - 1] + a * scratch[idx];
39 R Y(bspline)(
const INT k,
const R _x, R *scratch)
49 result_value = K(0.0);
51 if (K(0.0) < x && x < kk)
59 r = (INT)LRINT(CEIL(x) - K(1.0));
64 result_value = K(1.0);
65 for (j = 0; j < k - 1; j++)
66 result_value *= x/((R)(j+1));
71 for (idx = 0; idx < k; idx++)
72 scratch[idx] = K(0.0);
74 scratch[k-r-1] = K(1.0);
83 for (j = 1, og = g2 + 1; j <= g1; j++, og++)
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];
92 for (og-- ; j <= g2; j++)
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];
102 bspline_help(k, x, scratch, j, ug, og, r);
105 result_value = scratch[k-1];
108 return(result_value);