Bug Summary

File:src/nrrd/kernel.c
Location:line 3423, column 7
Description:Branch condition evaluates to a garbage value

Annotated Source Code

1/*
2 Teem: Tools to process and visualize scientific data and images .
3 Copyright (C) 2013, 2012, 2011, 2010, 2009 University of Chicago
4 Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann
5 Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah
6
7 This library is free software; you can redistribute it and/or
8 modify it under the terms of the GNU Lesser General Public License
9 (LGPL) as published by the Free Software Foundation; either
10 version 2.1 of the License, or (at your option) any later version.
11 The terms of redistributing and/or modifying this software also
12 include exceptions to the LGPL that facilitate static linking.
13
14 This library is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 Lesser General Public License for more details.
18
19 You should have received a copy of the GNU Lesser General Public License
20 along with this library; if not, write to Free Software Foundation, Inc.,
21 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
22*/
23
24#include "nrrd.h"
25
26/*
27** summary of information about how the kernel parameter vector is set:
28** Note that, annoyingly, nrrdKernelUsesScale (at end of this file)
29** has to be updated to record which kernels (including their derivatives)
30** use parm[0] for scale.
31
32 numParm parm[0] parm[1] parm[2]
33 nrrdKernelHann 2 scale cut-off
34 nrrdKernelBlackman 2 scale cut-off
35 nrrdKernelCatmullRom 0
36 nrrdKernelBSpline3 0
37nrrdKernelBSpline3ApproxInverse 0
38 nrrdKernelBSpline5 0
39nrrdKernelBSpline5ApproxInverse 0
40 nrrdKernelBSpline7 0
41nrrdKernelBSpline7ApproxInverse 0
42 nrrdKernelZero 1 scale
43 nrrdKernelBox 1 scale
44nrrdKernelCatmullRomSupportDebug 1 support
45 nrrdKernelBoxSupportDebug 1 support
46 nrrdKernelCos4SupportDebug 1 support
47 nrrdKernelCheap 1 scale
48nrrdKernelHermiteScaleSpaceFlag 0
49 nrrdKernelTent 1 scale
50 nrrdKernelForwDiff 1 scale
51 nrrdKernelCentDiff 1 scale
52 nrrdKernelBCCubic 3 scale B C
53 nrrdKernelAQuartic 2 scale A
54 nrrdKernelC3Quintic 0
55 nrrdKernelC4Hexic 0
56 nrrdKernelC4HexicApproxInverse 0
57 nrrdKernelC5Septic 0
58 nrrdKernelGaussian 2 sigma cut-off
59 nrrdKernelDiscreteGaussian 2 sigma cut-off
60 nrrdKernelTMF[][][] 1 a
61
62** Note that when parm[0] is named "scale", that parameter is optional,
63** and the default is 1.0, when given in string form
64** E.g. "tent" is understood as "tent:1",
65** but "gauss:4" isn't complete and won't parse; while "gauss:1,4" is good
66** See note above about nrrdKernelUsesScale (at end of this file)
67*/
68
69/* these functions replace what had been a lot of
70 identical functions for similar kernels */
71
72static double
73returnZero(const double *parm) {
74 AIR_UNUSED(parm)(void)(parm);
75 return 0.0;
76}
77
78static double
79returnOne(const double *parm) {
80 AIR_UNUSED(parm)(void)(parm);
81 return 1.0;
82}
83
84static double
85returnTwo(const double *parm) {
86 AIR_UNUSED(parm)(void)(parm);
87 return 2.0;
88}
89
90static double
91returnThree(const double *parm) {
92 AIR_UNUSED(parm)(void)(parm);
93 return 3.0;
94}
95
96static double
97returnFour(const double *parm) {
98 AIR_UNUSED(parm)(void)(parm);
99 return 4.0;
100}
101
102/* ------------------------------------------------------------ */
103
104/* learned: if you copy/paste the macros for these kernels into
105** other code, you *have* to make sure that the arguments for the
106** kernels that are supposed to be reals, are not passed as an
107** integral type (had this problem trying to re-use _BCCUBIC
108** with constant B="1" and C="0" and had trouble figuring out
109** why the kernel was give garbage results
110*/
111
112/* the "zero" kernel is here more as a template than for anything else
113 (as well as when you need the derivative of nrrdKernelForwDiff or
114 any other piece-wise constant kernels)
115 In particular the support method is pretty silly. */
116
117#define _ZERO(x)0 0
118
119static double
120_nrrdZeroSup(const double *parm) {
121 double S;
122
123 S = parm[0];
124 return S;
125}
126
127static double
128_nrrdZero1_d(double x, const double *parm) {
129 double S;
130
131 S = parm[0];
132 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x))/S;
133 return _ZERO(x)0/S;
134}
135
136static float
137_nrrdZero1_f(float x, const double *parm) {
138 float S;
139
140 S = AIR_CAST(float, parm[0])((float)(parm[0]));
141 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x))/S;
142 return _ZERO(x)0/S;
143}
144
145static void
146_nrrdZeroN_d(double *f, const double *x, size_t len, const double *parm) {
147 double S;
148 double t;
149 size_t i;
150
151 S = parm[0];
152 for (i=0; i<len; i++) {
153 t = x[i]; t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t))/S;
154 f[i] = _ZERO(t)0/S;
155 }
156}
157
158static void
159_nrrdZeroN_f(float *f, const float *x, size_t len, const double *parm) {
160 float t, S;
161 size_t i;
162
163 S = AIR_CAST(float, parm[0])((float)(parm[0]));
164 for (i=0; i<len; i++) {
165 t = x[i]; t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t))/S;
166 f[i] = _ZERO(t)0/S;
167 }
168}
169
170static NrrdKernel
171_nrrdKernelZero = {
172 "zero",
173 1, _nrrdZeroSup, returnZero,
174 _nrrdZero1_f, _nrrdZeroN_f, _nrrdZero1_d, _nrrdZeroN_d
175};
176NrrdKernel *const
177nrrdKernelZero = &_nrrdKernelZero;
178
179/* ------------------------------------------------------------ */
180
181#define _BOX(x)(x > 0.5 ? 0 : (x < 0.5 ? 1 : 0.5)) (x > 0.5 ? 0 : (x < 0.5 ? 1 : 0.5))
182
183static double
184_nrrdBoxSup(const double *parm) {
185 double S;
186
187 S = parm[0];
188 /* adding the 0.5 is to ensure that weights computed within the
189 support really do catch all the non-zero values */
190 return S/2 + 0.5;
191}
192
193static double
194_nrrdBox1_d(double x, const double *parm) {
195 double S;
196
197 S = parm[0];
198 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x))/S;
199 return _BOX(x)(x > 0.5 ? 0 : (x < 0.5 ? 1 : 0.5))/S;
200}
201
202static float
203_nrrdBox1_f(float x, const double *parm) {
204 float S;
205
206 S = AIR_CAST(float, parm[0])((float)(parm[0]));
207 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x))/S;
208 return AIR_CAST(float, _BOX(x)/S)((float)((x > 0.5 ? 0 : (x < 0.5 ? 1 : 0.5))/S));
209}
210
211static void
212_nrrdBoxN_d(double *f, const double *x, size_t len, const double *parm) {
213 double S;
214 double t;
215 size_t i;
216
217 S = parm[0];
218 for (i=0; i<len; i++) {
219 t = x[i]; t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t))/S;
220 f[i] = _BOX(t)(t > 0.5 ? 0 : (t < 0.5 ? 1 : 0.5))/S;
221 }
222}
223
224static void
225_nrrdBoxN_f(float *f, const float *x, size_t len, const double *parm) {
226 float t, S;
227 size_t i;
228
229 S = AIR_CAST(float, parm[0])((float)(parm[0]));
230 for (i=0; i<len; i++) {
231 t = x[i]; t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t))/S;
232 f[i] = AIR_CAST(float, _BOX(t)/S)((float)((t > 0.5 ? 0 : (t < 0.5 ? 1 : 0.5))/S));
233 }
234}
235
236static NrrdKernel
237_nrrdKernelBox = {
238 "box",
239 1, _nrrdBoxSup, returnOne,
240 _nrrdBox1_f, _nrrdBoxN_f, _nrrdBox1_d, _nrrdBoxN_d
241};
242NrrdKernel *const
243nrrdKernelBox = &_nrrdKernelBox;
244
245/* ------------------------------------------------------------ */
246
247static double
248_nrrdBoxSDSup(const double *parm) {
249
250 return parm[0];
251}
252
253static double
254_nrrdBoxSD1_d(double x, const double *parm) {
255 AIR_UNUSED(parm)(void)(parm);
256 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x));
257 return _BOX(x)(x > 0.5 ? 0 : (x < 0.5 ? 1 : 0.5));
258}
259
260static float
261_nrrdBoxSD1_f(float x, const double *parm) {
262 AIR_UNUSED(parm)(void)(parm);
263 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x));
264 return AIR_CAST(float, _BOX(x))((float)((x > 0.5 ? 0 : (x < 0.5 ? 1 : 0.5))));
265}
266
267static void
268_nrrdBoxSDN_d(double *f, const double *x, size_t len, const double *parm) {
269 size_t i;
270 AIR_UNUSED(parm)(void)(parm);
271 for (i=0; i<len; i++) {
272 double t;
273 t = AIR_ABS(x[i])((x[i]) > 0.0f ? (x[i]) : -(x[i]));
274 f[i] = _BOX(t)(t > 0.5 ? 0 : (t < 0.5 ? 1 : 0.5));
275 }
276}
277
278static void
279_nrrdBoxSDN_f(float *f, const float *x, size_t len, const double *parm) {
280 size_t i;
281 AIR_UNUSED(parm)(void)(parm);
282 for (i=0; i<len; i++) {
283 float t;
284 t = AIR_ABS(x[i])((x[i]) > 0.0f ? (x[i]) : -(x[i]));
285 f[i] = AIR_CAST(float, _BOX(t))((float)((t > 0.5 ? 0 : (t < 0.5 ? 1 : 0.5))));
286 }
287}
288
289static NrrdKernel
290_nrrdKernelBoxSupportDebug = {
291 "boxsup",
292 1, _nrrdBoxSDSup, returnOne,
293 _nrrdBoxSD1_f, _nrrdBoxSDN_f, _nrrdBoxSD1_d, _nrrdBoxSDN_d
294};
295NrrdKernel *const
296nrrdKernelBoxSupportDebug = &_nrrdKernelBoxSupportDebug;
297
298/* ------------------------------------------------------------ */
299
300#define COS4(x)(x > 0.5 ? 0.0 : cos(3.14159265358979323846*x)*cos(3.14159265358979323846
*x)*cos(3.14159265358979323846*x)*cos(3.14159265358979323846*
x)) (x > 0.5 \
301 ? 0.0 \
302 : cos(AIR_PI3.14159265358979323846*x)*cos(AIR_PI3.14159265358979323846*x)*cos(AIR_PI3.14159265358979323846*x)*cos(AIR_PI3.14159265358979323846*x))
303
304static double
305_nrrdCos4SDInt(const double *parm) {
306 AIR_UNUSED(parm)(void)(parm);
307 return 3.0/8.0;
308}
309
310static double
311_nrrdCos4SDSup(const double *parm) {
312
313 return parm[0];
314}
315
316static double
317_nrrdCos4SD1_d(double x, const double *parm) {
318 AIR_UNUSED(parm)(void)(parm);
319 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x));
320 return COS4(x)(x > 0.5 ? 0.0 : cos(3.14159265358979323846*x)*cos(3.14159265358979323846
*x)*cos(3.14159265358979323846*x)*cos(3.14159265358979323846*
x));
321}
322
323static float
324_nrrdCos4SD1_f(float x, const double *parm) {
325 AIR_UNUSED(parm)(void)(parm);
326 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x));
327 return AIR_CAST(float, COS4(x))((float)((x > 0.5 ? 0.0 : cos(3.14159265358979323846*x)*cos
(3.14159265358979323846*x)*cos(3.14159265358979323846*x)*cos(
3.14159265358979323846*x))));
328}
329
330static void
331_nrrdCos4SDN_d(double *f, const double *x, size_t len, const double *parm) {
332 size_t i;
333 AIR_UNUSED(parm)(void)(parm);
334 for (i=0; i<len; i++) {
335 double t;
336 t = AIR_ABS(x[i])((x[i]) > 0.0f ? (x[i]) : -(x[i]));
337 f[i] = COS4(t)(t > 0.5 ? 0.0 : cos(3.14159265358979323846*t)*cos(3.14159265358979323846
*t)*cos(3.14159265358979323846*t)*cos(3.14159265358979323846*
t));
338 }
339}
340
341static void
342_nrrdCos4SDN_f(float *f, const float *x, size_t len, const double *parm) {
343 size_t i;
344 AIR_UNUSED(parm)(void)(parm);
345 for (i=0; i<len; i++) {
346 float t;
347 t = AIR_ABS(x[i])((x[i]) > 0.0f ? (x[i]) : -(x[i]));
348 f[i] = AIR_CAST(float, COS4(t))((float)((t > 0.5 ? 0.0 : cos(3.14159265358979323846*t)*cos
(3.14159265358979323846*t)*cos(3.14159265358979323846*t)*cos(
3.14159265358979323846*t))));
349 }
350}
351
352static NrrdKernel
353_nrrdKernelCos4SupportDebug = {
354 "cos4sup",
355 1, _nrrdCos4SDSup, _nrrdCos4SDInt,
356 _nrrdCos4SD1_f, _nrrdCos4SDN_f, _nrrdCos4SD1_d, _nrrdCos4SDN_d
357};
358NrrdKernel *const
359nrrdKernelCos4SupportDebug = &_nrrdKernelCos4SupportDebug;
360
361/* ------------------------------------------------------------ */
362
363#define DCOS4(x)(x > 0.5 ? 0.0 : -4*3.14159265358979323846*(cos(3.14159265358979323846
*x)*cos(3.14159265358979323846*x)*cos(3.14159265358979323846*
x) *sin(3.14159265358979323846*x))) (x > 0.5 \
364 ? 0.0 \
365 : -4*AIR_PI3.14159265358979323846*(cos(AIR_PI3.14159265358979323846*x)*cos(AIR_PI3.14159265358979323846*x)*cos(AIR_PI3.14159265358979323846*x) \
366 *sin(AIR_PI3.14159265358979323846*x)))
367
368static double
369_nrrdDCos4SDSup(const double *parm) {
370 return parm[0];
371}
372
373static double
374_nrrdDCos4SD1_d(double x, const double *parm) {
375 int sgn;
376 AIR_UNUSED(parm)(void)(parm);
377 if (x < 0) { x = -x; sgn = -1; } else { sgn = 1; }
378 return sgn*DCOS4(x)(x > 0.5 ? 0.0 : -4*3.14159265358979323846*(cos(3.14159265358979323846
*x)*cos(3.14159265358979323846*x)*cos(3.14159265358979323846*
x) *sin(3.14159265358979323846*x)));
379}
380
381static float
382_nrrdDCos4SD1_f(float x, const double *parm) {
383 int sgn;
384 AIR_UNUSED(parm)(void)(parm);
385 if (x < 0) { x = -x; sgn = -1; } else { sgn = 1; }
386 return AIR_CAST(float, sgn*DCOS4(x))((float)(sgn*(x > 0.5 ? 0.0 : -4*3.14159265358979323846*(cos
(3.14159265358979323846*x)*cos(3.14159265358979323846*x)*cos(
3.14159265358979323846*x) *sin(3.14159265358979323846*x)))));
387}
388
389static void
390_nrrdDCos4SDN_d(double *f, const double *x, size_t len, const double *parm) {
391 int sgn;
392 double t;
393 size_t i;
394 AIR_UNUSED(parm)(void)(parm);
395 for (i=0; i<len; i++) {
396 t = x[i];
397 if (t < 0) { t = -t; sgn = -1; } else { sgn = 1; }
398 f[i] = sgn*DCOS4(t)(t > 0.5 ? 0.0 : -4*3.14159265358979323846*(cos(3.14159265358979323846
*t)*cos(3.14159265358979323846*t)*cos(3.14159265358979323846*
t) *sin(3.14159265358979323846*t)));
399 }
400}
401
402static void
403_nrrdDCos4SDN_f(float *f, const float *x, size_t len, const double *parm) {
404 int sgn;
405 float t;
406 size_t i;
407 AIR_UNUSED(parm)(void)(parm);
408 for (i=0; i<len; i++) {
409 t = x[i];
410 if (t < 0) { t = -t; sgn = -1; } else { sgn = 1; }
411 f[i] = AIR_CAST(float, sgn*DCOS4(t))((float)(sgn*(t > 0.5 ? 0.0 : -4*3.14159265358979323846*(cos
(3.14159265358979323846*t)*cos(3.14159265358979323846*t)*cos(
3.14159265358979323846*t) *sin(3.14159265358979323846*t)))));
412 }
413}
414
415static NrrdKernel
416_nrrdKernelCos4SupportDebugD = {
417 "cos4supD",
418 1, _nrrdDCos4SDSup, returnZero,
419 _nrrdDCos4SD1_f, _nrrdDCos4SDN_f, _nrrdDCos4SD1_d, _nrrdDCos4SDN_d
420};
421NrrdKernel *const
422nrrdKernelCos4SupportDebugD = &_nrrdKernelCos4SupportDebugD;
423
424/* ------------------------------------------------------------ */
425
426#define DDCOS4(x)(x > 0.5 ? 0.0 : -2*3.14159265358979323846*3.14159265358979323846
*(cos(3.14159265358979323846*2*x) + cos(3.14159265358979323846
*4*x))) (x > 0.5 \
427 ? 0.0 \
428 : -2*AIR_PI3.14159265358979323846*AIR_PI3.14159265358979323846*(cos(AIR_PI3.14159265358979323846*2*x) + cos(AIR_PI3.14159265358979323846*4*x)))
429
430static double
431_nrrdDDCos4SDSup(const double *parm) {
432 return parm[0];
433}
434
435static double
436_nrrdDDCos4SD1_d(double x, const double *parm) {
437 AIR_UNUSED(parm)(void)(parm);
438 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x));
439 return DDCOS4(x)(x > 0.5 ? 0.0 : -2*3.14159265358979323846*3.14159265358979323846
*(cos(3.14159265358979323846*2*x) + cos(3.14159265358979323846
*4*x)));
440}
441
442static float
443_nrrdDDCos4SD1_f(float x, const double *parm) {
444 AIR_UNUSED(parm)(void)(parm);
445 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x));
446 return AIR_CAST(float, DDCOS4(x))((float)((x > 0.5 ? 0.0 : -2*3.14159265358979323846*3.14159265358979323846
*(cos(3.14159265358979323846*2*x) + cos(3.14159265358979323846
*4*x)))));
447}
448
449static void
450_nrrdDDCos4SDN_d(double *f, const double *x, size_t len, const double *parm) {
451 size_t i;
452 AIR_UNUSED(parm)(void)(parm);
453 for (i=0; i<len; i++) {
454 double t;
455 t = AIR_ABS(x[i])((x[i]) > 0.0f ? (x[i]) : -(x[i]));
456 f[i] = DDCOS4(t)(t > 0.5 ? 0.0 : -2*3.14159265358979323846*3.14159265358979323846
*(cos(3.14159265358979323846*2*t) + cos(3.14159265358979323846
*4*t)));
457 }
458}
459
460static void
461_nrrdDDCos4SDN_f(float *f, const float *x, size_t len, const double *parm) {
462 size_t i;
463 AIR_UNUSED(parm)(void)(parm);
464 for (i=0; i<len; i++) {
465 float t;
466 t = AIR_ABS(x[i])((x[i]) > 0.0f ? (x[i]) : -(x[i]));
467 f[i] = AIR_CAST(float, DDCOS4(t))((float)((t > 0.5 ? 0.0 : -2*3.14159265358979323846*3.14159265358979323846
*(cos(3.14159265358979323846*2*t) + cos(3.14159265358979323846
*4*t)))));
468 }
469}
470
471static NrrdKernel
472_nrrdKernelCos4SupportDebugDD = {
473 "cos4supDD",
474 1, _nrrdDDCos4SDSup, returnZero,
475 _nrrdDDCos4SD1_f, _nrrdDDCos4SDN_f, _nrrdDDCos4SD1_d, _nrrdDDCos4SDN_d
476};
477NrrdKernel *const
478nrrdKernelCos4SupportDebugDD = &_nrrdKernelCos4SupportDebugDD;
479
480/* ------------------------------------------------------------ */
481
482#define DDDCOS4(x)(x > 0.5 ? 0.0 : 4*3.14159265358979323846*3.14159265358979323846
*3.14159265358979323846*(sin(2*3.14159265358979323846*x) + 2*
sin(4*3.14159265358979323846*x))) (x > 0.5 \
483 ? 0.0 \
484 : 4*AIR_PI3.14159265358979323846*AIR_PI3.14159265358979323846*AIR_PI3.14159265358979323846*(sin(2*AIR_PI3.14159265358979323846*x) + 2*sin(4*AIR_PI3.14159265358979323846*x)))
485
486static double
487_nrrdDDDCos4SDSup(const double *parm) {
488 return parm[0];
489}
490
491static double
492_nrrdDDDCos4SD1_d(double x, const double *parm) {
493 int sgn;
494 AIR_UNUSED(parm)(void)(parm);
495 if (x < 0) { x = -x; sgn = -1; } else { sgn = 1; }
496 return sgn*DDDCOS4(x)(x > 0.5 ? 0.0 : 4*3.14159265358979323846*3.14159265358979323846
*3.14159265358979323846*(sin(2*3.14159265358979323846*x) + 2*
sin(4*3.14159265358979323846*x)));
497}
498
499static float
500_nrrdDDDCos4SD1_f(float x, const double *parm) {
501 int sgn;
502 AIR_UNUSED(parm)(void)(parm);
503 if (x < 0) { x = -x; sgn = -1; } else { sgn = 1; }
504 return AIR_CAST(float, sgn*DDDCOS4(x))((float)(sgn*(x > 0.5 ? 0.0 : 4*3.14159265358979323846*3.14159265358979323846
*3.14159265358979323846*(sin(2*3.14159265358979323846*x) + 2*
sin(4*3.14159265358979323846*x)))));
505}
506
507static void
508_nrrdDDDCos4SDN_d(double *f, const double *x, size_t len, const double *parm) {
509 int sgn;
510 double t;
511 size_t i;
512 AIR_UNUSED(parm)(void)(parm);
513 for (i=0; i<len; i++) {
514 t = x[i];
515 if (t < 0) { t = -t; sgn = -1; } else { sgn = 1; }
516 f[i] = sgn*DDDCOS4(t)(t > 0.5 ? 0.0 : 4*3.14159265358979323846*3.14159265358979323846
*3.14159265358979323846*(sin(2*3.14159265358979323846*t) + 2*
sin(4*3.14159265358979323846*t)));
517 }
518}
519
520static void
521_nrrdDDDCos4SDN_f(float *f, const float *x, size_t len, const double *parm) {
522 int sgn;
523 float t;
524 size_t i;
525 AIR_UNUSED(parm)(void)(parm);
526 for (i=0; i<len; i++) {
527 t = x[i];
528 if (t < 0) { t = -t; sgn = -1; } else { sgn = 1; }
529 f[i] = AIR_CAST(float, sgn*DDDCOS4(t))((float)(sgn*(t > 0.5 ? 0.0 : 4*3.14159265358979323846*3.14159265358979323846
*3.14159265358979323846*(sin(2*3.14159265358979323846*t) + 2*
sin(4*3.14159265358979323846*t)))));
530 }
531}
532
533static NrrdKernel
534_nrrdKernelCos4SupportDebugDDD = {
535 "cos4supDDD",
536 1, _nrrdDDDCos4SDSup, returnZero,
537 _nrrdDDDCos4SD1_f, _nrrdDDDCos4SDN_f, _nrrdDDDCos4SD1_d, _nrrdDDDCos4SDN_d
538};
539NrrdKernel *const
540nrrdKernelCos4SupportDebugDDD = &_nrrdKernelCos4SupportDebugDDD;
541
542/* ------------------------------------------------------------ */
543
544/* The point here is that post-kernel-evaluation, we need to see
545 which sample is closest to the origin, and this is one way of
546 enabling that
547 SO: this kernel will not usefully report its integral or support! */
548#define _CHEAP(x)((x) > 0.0f ? (x) : -(x)) AIR_ABS(x)((x) > 0.0f ? (x) : -(x))
549
550static double
551_nrrdCheapSup(const double *parm) {
552 double S;
553
554 S = parm[0];
555 /* adding the 0.5 is to insure that weights computed within the
556 support really do catch all the non-zero values */
557 return S/2 + 0.5;
558}
559
560static double
561_nrrdCheap1_d(double x, const double *parm) {
562
563 return _CHEAP(x)((x) > 0.0f ? (x) : -(x))/parm[0];
564}
565
566static float
567_nrrdCheap1_f(float x, const double *parm) {
568
569 return AIR_CAST(float, _CHEAP(x)/parm[0])((float)(((x) > 0.0f ? (x) : -(x))/parm[0]));
570}
571
572static void
573_nrrdCheapN_d(double *f, const double *x, size_t len, const double *parm) {
574 double t;
575 size_t i;
576
577 for (i=0; i<len; i++) {
578 t = x[i];
579 f[i] = _CHEAP(t)((t) > 0.0f ? (t) : -(t))/parm[0];
580 }
581}
582
583static void
584_nrrdCheapN_f(float *f, const float *x, size_t len, const double *parm) {
585 float t;
586 size_t i;
587
588 for (i=0; i<len; i++) {
589 t = x[i];
590 f[i] = AIR_CAST(float, _CHEAP(t)/parm[0])((float)(((t) > 0.0f ? (t) : -(t))/parm[0]));
591 }
592}
593
594static NrrdKernel
595_nrrdKernelCheap = {
596 "cheap",
597 1, _nrrdCheapSup, returnOne,
598 _nrrdCheap1_f, _nrrdCheapN_f, _nrrdCheap1_d, _nrrdCheapN_d
599};
600NrrdKernel *const
601nrrdKernelCheap = &_nrrdKernelCheap;
602
603/* ------------------------------------------------------------ */
604
605#define _TENT(x)(x >= 1 ? 0 : 1 - x) (x >= 1 ? 0 : 1 - x)
606
607static double
608_nrrdTentSup(const double *parm) {
609 double S;
610
611 S = parm[0];
612 return S;
613}
614
615static double
616_nrrdTent1_d(double x, const double *parm) {
617 double S;
618
619 S = parm[0];
620 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x))/S;
621 return S ? _TENT(x)(x >= 1 ? 0 : 1 - x)/S : x == 0;
622}
623
624static float
625_nrrdTent1_f(float x, const double *parm) {
626 float S;
627
628 S = AIR_CAST(float, parm[0])((float)(parm[0]));
629 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x))/S;
630 return S ? _TENT(x)(x >= 1 ? 0 : 1 - x)/S : x == 0;
631}
632
633static void
634_nrrdTentN_d(double *f, const double *x, size_t len, const double *parm) {
635 double S;
636 double t;
637 size_t i;
638
639 S = parm[0];
640 for (i=0; i<len; i++) {
641 t = x[i]; t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t))/S;
642 f[i] = S ? _TENT(t)(t >= 1 ? 0 : 1 - t)/S : t == 0;
643 }
644}
645
646static void
647_nrrdTentN_f(float *f, const float *x, size_t len, const double *parm) {
648 float t, S;
649 size_t i;
650
651 S = AIR_CAST(float, parm[0])((float)(parm[0]));
652 for (i=0; i<len; i++) {
653 t = x[i]; t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t))/S;
654 f[i] = S ? _TENT(t)(t >= 1 ? 0 : 1 - t)/S : t == 0;
655 }
656}
657
658static NrrdKernel
659_nrrdKernelTent = {
660 "tent",
661 1, _nrrdTentSup, returnOne,
662 _nrrdTent1_f, _nrrdTentN_f, _nrrdTent1_d, _nrrdTentN_d
663};
664NrrdKernel *const
665nrrdKernelTent = &_nrrdKernelTent;
666
667/* ------------------------------------------------------------ */
668
669/*
670** NOTE: THERE IS NOT REALLY A HERMITE KERNEL (at least not yet,
671** because it takes both values and derivatives as arguments, which
672** the NrrdKernel currently can't handle). This isn't really a
673** kernel, its mostly a flag (hence the name), but it also has the
674** role of generating weights according to linear interpolation, which
675** is useful for the eventual spline evaluation.
676**
677** This hack is in sinister collusion with gage, to enable Hermite
678** interpolation for its stack reconstruction.
679*/
680
681static double
682_nrrdHermite1_d(double x, const double *parm) {
683 AIR_UNUSED(parm)(void)(parm);
684 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x));
685 return _TENT(x)(x >= 1 ? 0 : 1 - x);
686}
687
688static float
689_nrrdHermite1_f(float x, const double *parm) {
690 AIR_UNUSED(parm)(void)(parm);
691 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x));
692 return _TENT(x)(x >= 1 ? 0 : 1 - x);
693}
694
695static void
696_nrrdHermiteN_d(double *f, const double *x, size_t len, const double *parm) {
697 double t;
698 size_t i;
699 AIR_UNUSED(parm)(void)(parm);
700 for (i=0; i<len; i++) {
701 t = x[i]; t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t));
702 f[i] = _TENT(t)(t >= 1 ? 0 : 1 - t);
703 }
704}
705
706static void
707_nrrdHermiteN_f(float *f, const float *x, size_t len, const double *parm) {
708 float t;
709 size_t i;
710 AIR_UNUSED(parm)(void)(parm);
711 for (i=0; i<len; i++) {
712 t = x[i]; t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t));
713 f[i] = _TENT(t)(t >= 1 ? 0 : 1 - t);
714 }
715}
716
717/* HEY: should just re-use fields from nrrdKernelTent, instead
718 of creating new functions */
719static NrrdKernel
720_nrrdKernelHermiteScaleSpaceFlag = {
721 "hermiteSS",
722 0, returnOne, returnOne,
723 _nrrdHermite1_f, _nrrdHermiteN_f, _nrrdHermite1_d, _nrrdHermiteN_d
724};
725NrrdKernel *const
726nrrdKernelHermiteScaleSpaceFlag = &_nrrdKernelHermiteScaleSpaceFlag;
727
728/* ------------------------------------------------------------ */
729
730#define _FORDIF(x)(x < -1 ? 0.0f : (x < 0 ? 1.0f : (x < 1 ? -1.0f : 0.0f
))) (x < -1 ? 0.0f : \
731 (x < 0 ? 1.0f : \
732 (x < 1 ? -1.0f : 0.0f )))
733
734static double
735_nrrdFDSup(const double *parm) {
736 double S;
737
738 S = parm[0];
739 return S+0.0001; /* sigh */
740}
741
742static double
743_nrrdFD1_d(double x, const double *parm) {
744 double t, S;
745
746 S = parm[0];
747 t = x/S;
748 return _FORDIF(t)(t < -1 ? 0.0f : (t < 0 ? 1.0f : (t < 1 ? -1.0f : 0.0f
)))/(S*S);
749}
750
751static float
752_nrrdFD1_f(float x, const double *parm) {
753 float t, S;
754
755 S = AIR_CAST(float, parm[0])((float)(parm[0]));
756 t = x/S;
757 return _FORDIF(t)(t < -1 ? 0.0f : (t < 0 ? 1.0f : (t < 1 ? -1.0f : 0.0f
)))/(S*S);
758}
759
760static void
761_nrrdFDN_d(double *f, const double *x, size_t len, const double *parm) {
762 double t, S;
763 size_t i;
764
765 S = parm[0];
766 for (i=0; i<len; i++) {
767 t = x[i]/S;
768 f[i] = _FORDIF(t)(t < -1 ? 0.0f : (t < 0 ? 1.0f : (t < 1 ? -1.0f : 0.0f
)))/(S*S);
769 }
770}
771
772static void
773_nrrdFDN_f(float *f, const float *x, size_t len, const double *parm) {
774 float t, S;
775 size_t i;
776
777 S = AIR_CAST(float, parm[0])((float)(parm[0]));
778 for (i=0; i<len; i++) {
779 t = x[i]/S;
780 f[i] = _FORDIF(t)(t < -1 ? 0.0f : (t < 0 ? 1.0f : (t < 1 ? -1.0f : 0.0f
)))/(S*S);
781 }
782}
783
784static NrrdKernel
785_nrrdKernelFD = {
786 "fordif",
787 1, _nrrdFDSup, returnZero,
788 _nrrdFD1_f, _nrrdFDN_f, _nrrdFD1_d, _nrrdFDN_d
789};
790NrrdKernel *const
791nrrdKernelForwDiff = &_nrrdKernelFD;
792
793/* ------------------------------------------------------------ */
794
795#define _CENDIF(x)(x <= -2 ? 0 : (x <= -1 ? 0.5*x + 1 : (x <= 1 ? -0.5
*x : (x <= 2 ? 0.5*x - 1 : 0 )))) (x <= -2 ? 0 : \
796 (x <= -1 ? 0.5*x + 1 : \
797 (x <= 1 ? -0.5*x : \
798 (x <= 2 ? 0.5*x - 1 : 0 ))))
799
800static double
801_nrrdCDSup(const double *parm) {
802 double S;
803
804 S = parm[0];
805 return 2*S;
806}
807
808static double
809_nrrdCD1_d(double x, const double *parm) {
810 double S;
811
812 S = parm[0];
813 x /= S;
814 return _CENDIF(x)(x <= -2 ? 0 : (x <= -1 ? 0.5*x + 1 : (x <= 1 ? -0.5
*x : (x <= 2 ? 0.5*x - 1 : 0 ))))/(S*S);
815}
816
817static float
818_nrrdCD1_f(float x, const double *parm) {
819 float S;
820
821 S = AIR_CAST(float, parm[0])((float)(parm[0]));
822 x /= S;
823 return AIR_CAST(float, _CENDIF(x)/(S*S))((float)((x <= -2 ? 0 : (x <= -1 ? 0.5*x + 1 : (x <=
1 ? -0.5*x : (x <= 2 ? 0.5*x - 1 : 0 ))))/(S*S)));
824}
825
826static void
827_nrrdCDN_d(double *f, const double *x, size_t len, const double *parm) {
828 double S;
829 double t;
830 size_t i;
831
832 S = parm[0];
833 for (i=0; i<len; i++) {
834 t = x[i]/S;
835 f[i] = _CENDIF(t)(t <= -2 ? 0 : (t <= -1 ? 0.5*t + 1 : (t <= 1 ? -0.5
*t : (t <= 2 ? 0.5*t - 1 : 0 ))))/(S*S);
836 }
837}
838
839static void
840_nrrdCDN_f(float *f, const float *x, size_t len, const double *parm) {
841 float t, S;
842 size_t i;
843
844 S = AIR_CAST(float, parm[0])((float)(parm[0]));
845 for (i=0; i<len; i++) {
846 t = x[i]/S;
847 f[i] = AIR_CAST(float, _CENDIF(t)/(S*S))((float)((t <= -2 ? 0 : (t <= -1 ? 0.5*t + 1 : (t <=
1 ? -0.5*t : (t <= 2 ? 0.5*t - 1 : 0 ))))/(S*S)));
848 }
849}
850
851static NrrdKernel
852_nrrdKernelCD = {
853 "cendif",
854 1, _nrrdCDSup, returnZero,
855 _nrrdCD1_f, _nrrdCDN_f, _nrrdCD1_d, _nrrdCDN_d
856};
857NrrdKernel *const
858nrrdKernelCentDiff = &_nrrdKernelCD;
859
860/* ------------------------------------------------------------ */
861
862#define _BCCUBIC(x, B, C)(x >= 2.0 ? 0 : (x >= 1.0 ? (((-B/6 - C)*x + B + 5*C)*x
-2*B - 8*C)*x + 4*B/3 + 4*C : ((2 - 3*B/2 - C)*x - 3 + 2*B +
C)*x*x + 1 - B/3)) \
863 (x >= 2.0 ? 0 : \
864 (x >= 1.0 \
865 ? (((-B/6 - C)*x + B + 5*C)*x -2*B - 8*C)*x + 4*B/3 + 4*C \
866 : ((2 - 3*B/2 - C)*x - 3 + 2*B + C)*x*x + 1 - B/3))
867
868static double
869_nrrdBCSup(const double *parm) {
870 double S;
871
872 S = parm[0];
873 return 2*S;
874}
875
876static double
877_nrrdBC1_d(double x, const double *parm) {
878 double S;
879 double B, C;
880
881 S = parm[0]; B = parm[1]; C = parm[2];
882 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x))/S;
883 return _BCCUBIC(x, B, C)(x >= 2.0 ? 0 : (x >= 1.0 ? (((-B/6 - C)*x + B + 5*C)*x
-2*B - 8*C)*x + 4*B/3 + 4*C : ((2 - 3*B/2 - C)*x - 3 + 2*B +
C)*x*x + 1 - B/3))/S;
884}
885
886static float
887_nrrdBC1_f(float x, const double *parm) {
888 float B, C, S;
889
890 S = AIR_CAST(float, parm[0])((float)(parm[0]));
891 B = AIR_CAST(float, parm[1])((float)(parm[1]));
892 C = AIR_CAST(float, parm[2])((float)(parm[2]));
893 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x))/S;
894 return _BCCUBIC(x, B, C)(x >= 2.0 ? 0 : (x >= 1.0 ? (((-B/6 - C)*x + B + 5*C)*x
-2*B - 8*C)*x + 4*B/3 + 4*C : ((2 - 3*B/2 - C)*x - 3 + 2*B +
C)*x*x + 1 - B/3))/S;
895}
896
897static void
898_nrrdBCN_d(double *f, const double *x, size_t len, const double *parm) {
899 double S;
900 double t, B, C;
901 size_t i;
902
903 S = parm[0]; B = parm[1]; C = parm[2];
904 for (i=0; i<len; i++) {
905 t = x[i];
906 t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t))/S;
907 f[i] = _BCCUBIC(t, B, C)(t >= 2.0 ? 0 : (t >= 1.0 ? (((-B/6 - C)*t + B + 5*C)*t
-2*B - 8*C)*t + 4*B/3 + 4*C : ((2 - 3*B/2 - C)*t - 3 + 2*B +
C)*t*t + 1 - B/3))/S;
908 }
909}
910
911static void
912_nrrdBCN_f(float *f, const float *x, size_t len, const double *parm) {
913 float S, t, B, C;
914 size_t i;
915
916 S = AIR_CAST(float, parm[0])((float)(parm[0]));
917 B = AIR_CAST(float, parm[1])((float)(parm[1]));
918 C = AIR_CAST(float, parm[2])((float)(parm[2]));
919 for (i=0; i<len; i++) {
920 t = x[i];
921 t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t))/S;
922 f[i] = _BCCUBIC(t, B, C)(t >= 2.0 ? 0 : (t >= 1.0 ? (((-B/6 - C)*t + B + 5*C)*t
-2*B - 8*C)*t + 4*B/3 + 4*C : ((2 - 3*B/2 - C)*t - 3 + 2*B +
C)*t*t + 1 - B/3))/S;
923 }
924}
925
926static NrrdKernel
927_nrrdKernelBC = {
928 "BCcubic",
929 3, _nrrdBCSup, returnOne,
930 _nrrdBC1_f, _nrrdBCN_f, _nrrdBC1_d, _nrrdBCN_d
931};
932NrrdKernel *const
933nrrdKernelBCCubic = &_nrrdKernelBC;
934
935/* ------------------------------------------------------------ */
936
937#define _DBCCUBIC(x, B, C)(x >= 2.0 ? 0.0 : (x >= 1.0 ? ((-B/2 - 3*C)*x + 2*B + 10
*C)*x -2*B - 8*C : ((6 - 9*B/2 - 3*C)*x - 6 + 4*B + 2*C)*x)) \
938 (x >= 2.0 ? 0.0 : \
939 (x >= 1.0 \
940 ? ((-B/2 - 3*C)*x + 2*B + 10*C)*x -2*B - 8*C \
941 : ((6 - 9*B/2 - 3*C)*x - 6 + 4*B + 2*C)*x))
942
943static double
944_nrrdDBCSup(const double *parm) {
945 double S;
946
947 S = parm[0];
948 return 2*S;
949}
950
951static double
952_nrrdDBC1_d(double x, const double *parm) {
953 double S;
954 double B, C;
955 int sgn = 1;
956
957 S = parm[0]; B = parm[1]; C = parm[2];
958 if (x < 0) { x = -x; sgn = -1; }
959 x /= S;
960 return sgn*_DBCCUBIC(x, B, C)(x >= 2.0 ? 0.0 : (x >= 1.0 ? ((-B/2 - 3*C)*x + 2*B + 10
*C)*x -2*B - 8*C : ((6 - 9*B/2 - 3*C)*x - 6 + 4*B + 2*C)*x))/(S*S);
961}
962
963static float
964_nrrdDBC1_f(float x, const double *parm) {
965 float B, C, S;
966 int sgn = 1;
967
968 S = AIR_CAST(float, parm[0])((float)(parm[0]));
969 B = AIR_CAST(float, parm[1])((float)(parm[1]));
970 C = AIR_CAST(float, parm[2])((float)(parm[2]));
971 if (x < 0) { x = -x; sgn = -1; }
972 x /= S;
973 return AIR_CAST(float, sgn*_DBCCUBIC(x, B, C)/(S*S))((float)(sgn*(x >= 2.0 ? 0.0 : (x >= 1.0 ? ((-B/2 - 3*C
)*x + 2*B + 10*C)*x -2*B - 8*C : ((6 - 9*B/2 - 3*C)*x - 6 + 4
*B + 2*C)*x))/(S*S)));
974}
975
976static void
977_nrrdDBCN_d(double *f, const double *x, size_t len, const double *parm) {
978 double S;
979 double t, B, C;
980 size_t i;
981 int sgn;
982
983 S = parm[0]; B = parm[1]; C = parm[2];
984 for (i=0; i<len; i++) {
985 t = x[i]/S;
986 if (t < 0) { t = -t; sgn = -1; } else { sgn = 1; }
987 f[i] = sgn*_DBCCUBIC(t, B, C)(t >= 2.0 ? 0.0 : (t >= 1.0 ? ((-B/2 - 3*C)*t + 2*B + 10
*C)*t -2*B - 8*C : ((6 - 9*B/2 - 3*C)*t - 6 + 4*B + 2*C)*t))/(S*S);
988 }
989}
990
991static void
992_nrrdDBCN_f(float *f, const float *x, size_t len, const double *parm) {
993 float S, t, B, C;
994 int sgn;
995 size_t i;
996
997 S = AIR_CAST(float, parm[0])((float)(parm[0]));
998 B = AIR_CAST(float, parm[1])((float)(parm[1]));
999 C = AIR_CAST(float, parm[2])((float)(parm[2]));
1000 for (i=0; i<len; i++) {
1001 t = x[i]/S;
1002 if (t < 0) { t = -t; sgn = -1; } else { sgn = 1; }
1003 f[i] = AIR_CAST(float, sgn*_DBCCUBIC(t, B, C)/(S*S))((float)(sgn*(t >= 2.0 ? 0.0 : (t >= 1.0 ? ((-B/2 - 3*C
)*t + 2*B + 10*C)*t -2*B - 8*C : ((6 - 9*B/2 - 3*C)*t - 6 + 4
*B + 2*C)*t))/(S*S)));
1004 }
1005}
1006
1007static NrrdKernel
1008_nrrdKernelDBC = {
1009 "BCcubicD",
1010 3, _nrrdDBCSup, returnZero,
1011 _nrrdDBC1_f, _nrrdDBCN_f, _nrrdDBC1_d, _nrrdDBCN_d
1012};
1013NrrdKernel *const
1014nrrdKernelBCCubicD = &_nrrdKernelDBC;
1015
1016/* ------------------------------------------------------------ */
1017
1018#define _DDBCCUBIC(x, B, C)(x >= 2.0 ? 0 : (x >= 1.0 ? (-B - 6*C)*x + 2*B + 10*C :
(12 - 9*B - 6*C)*x - 6 + 4*B + 2*C )) \
1019 (x >= 2.0 ? 0 : \
1020 (x >= 1.0 \
1021 ? (-B - 6*C)*x + 2*B + 10*C \
1022 : (12 - 9*B - 6*C)*x - 6 + 4*B + 2*C ))
1023
1024static double
1025_nrrdDDBCSup(const double *parm) {
1026 double S;
1027
1028 S = parm[0];
1029 return 2*S;
1030}
1031
1032static double
1033_nrrdDDBC1_d(double x, const double *parm) {
1034 double S;
1035 double B, C;
1036
1037 S = parm[0]; B = parm[1]; C = parm[2];
1038 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x))/S;
1039 return _DDBCCUBIC(x, B, C)(x >= 2.0 ? 0 : (x >= 1.0 ? (-B - 6*C)*x + 2*B + 10*C :
(12 - 9*B - 6*C)*x - 6 + 4*B + 2*C ))/(S*S*S);
1040}
1041
1042static float
1043_nrrdDDBC1_f(float x, const double *parm) {
1044 float B, C, S;
1045
1046 S = AIR_CAST(float, parm[0])((float)(parm[0]));
1047 B = AIR_CAST(float, parm[1])((float)(parm[1]));
1048 C = AIR_CAST(float, parm[2])((float)(parm[2]));
1049 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x))/S;
1050 return _DDBCCUBIC(x, B, C)(x >= 2.0 ? 0 : (x >= 1.0 ? (-B - 6*C)*x + 2*B + 10*C :
(12 - 9*B - 6*C)*x - 6 + 4*B + 2*C ))/(S*S*S);
1051}
1052
1053static void
1054_nrrdDDBCN_d(double *f, const double *x, size_t len, const double *parm) {
1055 double S;
1056 double t, B, C;
1057 size_t i;
1058
1059 S = parm[0]; B = parm[1]; C = parm[2];
1060 for (i=0; i<len; i++) {
1061 t = x[i];
1062 t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t))/S;
1063 f[i] = _DDBCCUBIC(t, B, C)(t >= 2.0 ? 0 : (t >= 1.0 ? (-B - 6*C)*t + 2*B + 10*C :
(12 - 9*B - 6*C)*t - 6 + 4*B + 2*C ))/(S*S*S);
1064 }
1065}
1066
1067static void
1068_nrrdDDBCN_f(float *f, const float *x, size_t len, const double *parm) {
1069 float S, t, B, C;
1070 size_t i;
1071
1072 S = AIR_CAST(float, parm[0])((float)(parm[0]));
1073 B = AIR_CAST(float, parm[1])((float)(parm[1]));
1074 C = AIR_CAST(float, parm[2])((float)(parm[2]));
1075 for (i=0; i<len; i++) {
1076 t = x[i];
1077 t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t))/S;
1078 f[i] = _DDBCCUBIC(t, B, C)(t >= 2.0 ? 0 : (t >= 1.0 ? (-B - 6*C)*t + 2*B + 10*C :
(12 - 9*B - 6*C)*t - 6 + 4*B + 2*C ))/(S*S*S);
1079 }
1080}
1081
1082static NrrdKernel
1083_nrrdKernelDDBC = {
1084 "BCcubicDD",
1085 3, _nrrdDDBCSup, returnZero,
1086 _nrrdDDBC1_f, _nrrdDDBCN_f, _nrrdDDBC1_d, _nrrdDDBCN_d
1087};
1088NrrdKernel *const
1089nrrdKernelBCCubicDD = &_nrrdKernelDDBC;
1090
1091/* ------------------------------------------------------------ */
1092
1093/* if you've got the definition already, why not use it */
1094#define _CTMR(x)(x >= 2.0 ? 0 : (x >= 1.0 ? (((-0.0/6 - 0.5)*x + 0.0 + 5
*0.5)*x -2*0.0 - 8*0.5)*x + 4*0.0/3 + 4*0.5 : ((2 - 3*0.0/2 -
0.5)*x - 3 + 2*0.0 + 0.5)*x*x + 1 - 0.0/3)) _BCCUBIC(x, 0.0, 0.5)(x >= 2.0 ? 0 : (x >= 1.0 ? (((-0.0/6 - 0.5)*x + 0.0 + 5
*0.5)*x -2*0.0 - 8*0.5)*x + 4*0.0/3 + 4*0.5 : ((2 - 3*0.0/2 -
0.5)*x - 3 + 2*0.0 + 0.5)*x*x + 1 - 0.0/3))
1095
1096static double
1097_nrrdCTMR1_d(double x, const double *parm) {
1098 AIR_UNUSED(parm)(void)(parm);
1099 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x));
1100 return _CTMR(x)(x >= 2.0 ? 0 : (x >= 1.0 ? (((-0.0/6 - 0.5)*x + 0.0 + 5
*0.5)*x -2*0.0 - 8*0.5)*x + 4*0.0/3 + 4*0.5 : ((2 - 3*0.0/2 -
0.5)*x - 3 + 2*0.0 + 0.5)*x*x + 1 - 0.0/3));
1101}
1102
1103static float
1104_nrrdCTMR1_f(float x, const double *parm) {
1105 AIR_UNUSED(parm)(void)(parm);
1106 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x));
1107 return AIR_CAST(float, _CTMR(x))((float)((x >= 2.0 ? 0 : (x >= 1.0 ? (((-0.0/6 - 0.5)*x
+ 0.0 + 5*0.5)*x -2*0.0 - 8*0.5)*x + 4*0.0/3 + 4*0.5 : ((2 -
3*0.0/2 - 0.5)*x - 3 + 2*0.0 + 0.5)*x*x + 1 - 0.0/3))));
1108}
1109
1110static void
1111_nrrdCTMRN_d(double *f, const double *x, size_t len, const double *parm) {
1112 double t;
1113 size_t i;
1114 AIR_UNUSED(parm)(void)(parm);
1115 for (i=0; i<len; i++) {
1116 t = x[i];
1117 t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t));
1118 f[i] = _CTMR(t)(t >= 2.0 ? 0 : (t >= 1.0 ? (((-0.0/6 - 0.5)*t + 0.0 + 5
*0.5)*t -2*0.0 - 8*0.5)*t + 4*0.0/3 + 4*0.5 : ((2 - 3*0.0/2 -
0.5)*t - 3 + 2*0.0 + 0.5)*t*t + 1 - 0.0/3));
1119 }
1120}
1121
1122static void
1123_nrrdCTMRN_f(float *f, const float *x, size_t len, const double *parm) {
1124 float t;
1125 size_t i;
1126 AIR_UNUSED(parm)(void)(parm);
1127 for (i=0; i<len; i++) {
1128 t = x[i];
1129 t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t));
1130 f[i] = AIR_CAST(float, _CTMR(t))((float)((t >= 2.0 ? 0 : (t >= 1.0 ? (((-0.0/6 - 0.5)*t
+ 0.0 + 5*0.5)*t -2*0.0 - 8*0.5)*t + 4*0.0/3 + 4*0.5 : ((2 -
3*0.0/2 - 0.5)*t - 3 + 2*0.0 + 0.5)*t*t + 1 - 0.0/3))));
1131 }
1132}
1133
1134static NrrdKernel
1135_nrrdKernelCatmullRom = {
1136 "catmull-rom",
1137 0, returnTwo, returnOne,
1138 _nrrdCTMR1_f, _nrrdCTMRN_f, _nrrdCTMR1_d, _nrrdCTMRN_d
1139};
1140NrrdKernel *const
1141nrrdKernelCatmullRom = &_nrrdKernelCatmullRom;
1142
1143
1144static double
1145_nrrdCtmrSDSup(const double *parm) {
1146
1147 return AIR_MAX(2.0, parm[0])((2.0) > (parm[0]) ? (2.0) : (parm[0]));
1148}
1149
1150static NrrdKernel
1151_nrrdKernelCatmullRomSupportDebug = {
1152 "ctmrsup",
1153 1, _nrrdCtmrSDSup, returnOne,
1154 _nrrdCTMR1_f, _nrrdCTMRN_f, _nrrdCTMR1_d, _nrrdCTMRN_d
1155};
1156NrrdKernel *const
1157nrrdKernelCatmullRomSupportDebug = &_nrrdKernelCatmullRomSupportDebug;
1158
1159/* ------------------------------------------------------------ */
1160
1161/* if you've got the definition already, why not use it */
1162#define _DCTMR(x)(x >= 2.0 ? 0.0 : (x >= 1.0 ? ((-0.0/2 - 3*0.5)*x + 2*0.0
+ 10*0.5)*x -2*0.0 - 8*0.5 : ((6 - 9*0.0/2 - 3*0.5)*x - 6 + 4
*0.0 + 2*0.5)*x)) _DBCCUBIC(x, 0.0, 0.5)(x >= 2.0 ? 0.0 : (x >= 1.0 ? ((-0.0/2 - 3*0.5)*x + 2*0.0
+ 10*0.5)*x -2*0.0 - 8*0.5 : ((6 - 9*0.0/2 - 3*0.5)*x - 6 + 4
*0.0 + 2*0.5)*x))
1163
1164static double
1165_nrrdDCTMR1_d(double x, const double *parm) {
1166 int sgn;
1167 AIR_UNUSED(parm)(void)(parm);
1168 if (x < 0) { x = -x; sgn = -1; } else { sgn = 1; }
1169 return sgn*_DCTMR(x)(x >= 2.0 ? 0.0 : (x >= 1.0 ? ((-0.0/2 - 3*0.5)*x + 2*0.0
+ 10*0.5)*x -2*0.0 - 8*0.5 : ((6 - 9*0.0/2 - 3*0.5)*x - 6 + 4
*0.0 + 2*0.5)*x));
1170}
1171
1172static float
1173_nrrdDCTMR1_f(float x, const double *parm) {
1174 int sgn;
1175 AIR_UNUSED(parm)(void)(parm);
1176 if (x < 0) { x = -x; sgn = -1; } else { sgn = 1; }
1177 return AIR_CAST(float, sgn*_DCTMR(x))((float)(sgn*(x >= 2.0 ? 0.0 : (x >= 1.0 ? ((-0.0/2 - 3
*0.5)*x + 2*0.0 + 10*0.5)*x -2*0.0 - 8*0.5 : ((6 - 9*0.0/2 - 3
*0.5)*x - 6 + 4*0.0 + 2*0.5)*x))));
1178}
1179
1180static void
1181_nrrdDCTMRN_d(double *f, const double *x, size_t len, const double *parm) {
1182 int sgn;
1183 double t;
1184 size_t i;
1185 AIR_UNUSED(parm)(void)(parm);
1186 for (i=0; i<len; i++) {
1187 t = x[i];
1188 if (t < 0) { t = -t; sgn = -1; } else { sgn = 1; }
1189 f[i] = sgn*_DCTMR(t)(t >= 2.0 ? 0.0 : (t >= 1.0 ? ((-0.0/2 - 3*0.5)*t + 2*0.0
+ 10*0.5)*t -2*0.0 - 8*0.5 : ((6 - 9*0.0/2 - 3*0.5)*t - 6 + 4
*0.0 + 2*0.5)*t));
1190 }
1191}
1192
1193static void
1194_nrrdDCTMRN_f(float *f, const float *x, size_t len, const double *parm) {
1195 int sgn;
1196 float t;
1197 size_t i;
1198 AIR_UNUSED(parm)(void)(parm);
1199 for (i=0; i<len; i++) {
1200 t = x[i];
1201 if (t < 0) { t = -t; sgn = -1; } else { sgn = 1; }
1202 f[i] = AIR_CAST(float, sgn*_DCTMR(t))((float)(sgn*(t >= 2.0 ? 0.0 : (t >= 1.0 ? ((-0.0/2 - 3
*0.5)*t + 2*0.0 + 10*0.5)*t -2*0.0 - 8*0.5 : ((6 - 9*0.0/2 - 3
*0.5)*t - 6 + 4*0.0 + 2*0.5)*t))));
1203 }
1204}
1205
1206static NrrdKernel
1207_nrrdKernelCatmullRomD = {
1208 "catmull-romD",
1209 0, returnTwo, returnZero,
1210 _nrrdDCTMR1_f, _nrrdDCTMRN_f, _nrrdDCTMR1_d, _nrrdDCTMRN_d
1211};
1212NrrdKernel *const
1213nrrdKernelCatmullRomD = &_nrrdKernelCatmullRomD;
1214
1215static NrrdKernel
1216_nrrdKernelCatmullRomSupportDebugD = {
1217 "ctmrsupD",
1218 1, _nrrdCtmrSDSup, returnZero,
1219 _nrrdDCTMR1_f, _nrrdDCTMRN_f, _nrrdDCTMR1_d, _nrrdDCTMRN_d
1220};
1221NrrdKernel *const
1222nrrdKernelCatmullRomSupportDebugD = &_nrrdKernelCatmullRomSupportDebugD;
1223
1224/* ------------------------------------------------------------ */
1225
1226/* if you've got the definition already, why not use it */
1227#define _DDCTMR(x)(x >= 2.0 ? 0 : (x >= 1.0 ? (-0.0 - 6*0.5)*x + 2*0.0 + 10
*0.5 : (12 - 9*0.0 - 6*0.5)*x - 6 + 4*0.0 + 2*0.5 )) _DDBCCUBIC(x, 0.0, 0.5)(x >= 2.0 ? 0 : (x >= 1.0 ? (-0.0 - 6*0.5)*x + 2*0.0 + 10
*0.5 : (12 - 9*0.0 - 6*0.5)*x - 6 + 4*0.0 + 2*0.5 ))
1228
1229static double
1230_nrrdDDCTMR1_d(double x, const double *parm) {
1231 AIR_UNUSED(parm)(void)(parm);
1232 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x));
1233 return _DDCTMR(x)(x >= 2.0 ? 0 : (x >= 1.0 ? (-0.0 - 6*0.5)*x + 2*0.0 + 10
*0.5 : (12 - 9*0.0 - 6*0.5)*x - 6 + 4*0.0 + 2*0.5 ));
1234}
1235
1236static float
1237_nrrdDDCTMR1_f(float x, const double *parm) {
1238 AIR_UNUSED(parm)(void)(parm);
1239 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x));
1240 return AIR_CAST(float, _DDCTMR(x))((float)((x >= 2.0 ? 0 : (x >= 1.0 ? (-0.0 - 6*0.5)*x +
2*0.0 + 10*0.5 : (12 - 9*0.0 - 6*0.5)*x - 6 + 4*0.0 + 2*0.5 )
)));
1241}
1242
1243static void
1244_nrrdDDCTMRN_d(double *f, const double *x, size_t len, const double *parm) {
1245 double t;
1246 size_t i;
1247 AIR_UNUSED(parm)(void)(parm);
1248 for (i=0; i<len; i++) {
1249 t = x[i];
1250 t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t));
1251 f[i] = _DDCTMR(t)(t >= 2.0 ? 0 : (t >= 1.0 ? (-0.0 - 6*0.5)*t + 2*0.0 + 10
*0.5 : (12 - 9*0.0 - 6*0.5)*t - 6 + 4*0.0 + 2*0.5 ));
1252 }
1253}
1254
1255static void
1256_nrrdDDCTMRN_f(float *f, const float *x, size_t len, const double *parm) {
1257 float t;
1258 size_t i;
1259 AIR_UNUSED(parm)(void)(parm);
1260 for (i=0; i<len; i++) {
1261 t = x[i];
1262 t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t));
1263 f[i] = AIR_CAST(float, _DDCTMR(t))((float)((t >= 2.0 ? 0 : (t >= 1.0 ? (-0.0 - 6*0.5)*t +
2*0.0 + 10*0.5 : (12 - 9*0.0 - 6*0.5)*t - 6 + 4*0.0 + 2*0.5 )
)));
1264 }
1265}
1266
1267static NrrdKernel
1268_nrrdKernelCatmullRomDD = {
1269 "catmull-romDD",
1270 0, returnTwo, returnZero,
1271 _nrrdDDCTMR1_f, _nrrdDDCTMRN_f, _nrrdDDCTMR1_d, _nrrdDDCTMRN_d
1272};
1273NrrdKernel *const
1274nrrdKernelCatmullRomDD = &_nrrdKernelCatmullRomDD;
1275
1276static NrrdKernel
1277_nrrdKernelCatmullRomSupportDebugDD = {
1278 "ctmrsupDD",
1279 1, _nrrdCtmrSDSup, returnZero,
1280 _nrrdDDCTMR1_f, _nrrdDDCTMRN_f, _nrrdDDCTMR1_d, _nrrdDDCTMRN_d
1281};
1282NrrdKernel *const
1283nrrdKernelCatmullRomSupportDebugDD = &_nrrdKernelCatmullRomSupportDebugDD;
1284
1285/* ------------------------------------------------------------ */
1286
1287#define _AQUARTIC(x, A)(x >= 3.0 ? 0 : (x >= 2.0 ? A*(-54 + x*(81 + x*(-45 + x
*(11 - x)))) : (x >= 1.0 ? 4 - 6*A + x*(-10 + 25*A + x*(9 -
33*A + x*(-3.5 + 17*A + x*(0.5 - 3*A)))) : 1 + x*x*(-3 + 6*A
+ x*((2.5 - 10*A) + x*(-0.5 + 4*A)))))) \
1288 (x >= 3.0 ? 0 : \
1289 (x >= 2.0 \
1290 ? A*(-54 + x*(81 + x*(-45 + x*(11 - x)))) \
1291 : (x >= 1.0 \
1292 ? 4 - 6*A + x*(-10 + 25*A + x*(9 - 33*A \
1293 + x*(-3.5 + 17*A + x*(0.5 - 3*A)))) \
1294 : 1 + x*x*(-3 + 6*A + x*((2.5 - 10*A) + x*(-0.5 + 4*A))))))
1295
1296static double
1297_nrrdA4Sup(const double *parm) {
1298 double S;
1299
1300 S = parm[0];
1301 return 3*S;
1302}
1303
1304static double
1305_nrrdA41_d(double x, const double *parm) {
1306 double S;
1307 double A;
1308
1309 S = parm[0]; A = parm[1];
1310 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x))/S;
1311 return _AQUARTIC(x, A)(x >= 3.0 ? 0 : (x >= 2.0 ? A*(-54 + x*(81 + x*(-45 + x
*(11 - x)))) : (x >= 1.0 ? 4 - 6*A + x*(-10 + 25*A + x*(9 -
33*A + x*(-3.5 + 17*A + x*(0.5 - 3*A)))) : 1 + x*x*(-3 + 6*A
+ x*((2.5 - 10*A) + x*(-0.5 + 4*A))))))/S;
1312}
1313
1314static float
1315_nrrdA41_f(float x, const double *parm) {
1316 float A, S;
1317
1318 S = AIR_CAST(float, parm[0])((float)(parm[0])); A = AIR_CAST(float, parm[1])((float)(parm[1]));
1319 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x))/S;
1320 return AIR_CAST(float, _AQUARTIC(x, A)/S)((float)((x >= 3.0 ? 0 : (x >= 2.0 ? A*(-54 + x*(81 + x
*(-45 + x*(11 - x)))) : (x >= 1.0 ? 4 - 6*A + x*(-10 + 25*
A + x*(9 - 33*A + x*(-3.5 + 17*A + x*(0.5 - 3*A)))) : 1 + x*x
*(-3 + 6*A + x*((2.5 - 10*A) + x*(-0.5 + 4*A))))))/S));
1321}
1322
1323static void
1324_nrrdA4N_d(double *f, const double *x, size_t len, const double *parm) {
1325 double S;
1326 double t, A;
1327 size_t i;
1328
1329 S = parm[0]; A = parm[1];
1330 for (i=0; i<len; i++) {
1331 t = x[i];
1332 t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t))/S;
1333 f[i] = _AQUARTIC(t, A)(t >= 3.0 ? 0 : (t >= 2.0 ? A*(-54 + t*(81 + t*(-45 + t
*(11 - t)))) : (t >= 1.0 ? 4 - 6*A + t*(-10 + 25*A + t*(9 -
33*A + t*(-3.5 + 17*A + t*(0.5 - 3*A)))) : 1 + t*t*(-3 + 6*A
+ t*((2.5 - 10*A) + t*(-0.5 + 4*A))))))/S;
1334 }
1335}
1336
1337static void
1338_nrrdA4N_f(float *f, const float *x, size_t len, const double *parm) {
1339 float S, t, A;
1340 size_t i;
1341
1342 S = AIR_CAST(float, parm[0])((float)(parm[0])); A = AIR_CAST(float, parm[1])((float)(parm[1]));
1343 for (i=0; i<len; i++) {
1344 t = x[i];
1345 t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t))/S;
1346 f[i] = AIR_CAST(float, _AQUARTIC(t, A)/S)((float)((t >= 3.0 ? 0 : (t >= 2.0 ? A*(-54 + t*(81 + t
*(-45 + t*(11 - t)))) : (t >= 1.0 ? 4 - 6*A + t*(-10 + 25*
A + t*(9 - 33*A + t*(-3.5 + 17*A + t*(0.5 - 3*A)))) : 1 + t*t
*(-3 + 6*A + t*((2.5 - 10*A) + t*(-0.5 + 4*A))))))/S));
1347 }
1348}
1349
1350static NrrdKernel
1351_nrrdKernelA4 = {
1352 "Aquartic",
1353 2, _nrrdA4Sup, returnOne,
1354 _nrrdA41_f, _nrrdA4N_f, _nrrdA41_d, _nrrdA4N_d
1355};
1356NrrdKernel *const
1357nrrdKernelAQuartic = &_nrrdKernelA4;
1358
1359/* ------------------------------------------------------------ */
1360
1361#define _DAQUARTIC(x, A)(x >= 3.0 ? 0 : (x >= 2.0 ? A*(81 + x*(-90 + x*(33 - 4*
x))) : (x >= 1.0 ? -10 + 25*A + x*(18 - 66*A + x*(-10.5 + 51
*A + x*(2 - 12*A))) : x*(-6 + 12*A + x*(7.5 - 30*A + x*(-2 + 16
*A)))))) \
1362 (x >= 3.0 ? 0 : \
1363 (x >= 2.0 \
1364 ? A*(81 + x*(-90 + x*(33 - 4*x))) \
1365 : (x >= 1.0 \
1366 ? -10 + 25*A + x*(18 - 66*A + x*(-10.5 + 51*A + x*(2 - 12*A))) \
1367 : x*(-6 + 12*A + x*(7.5 - 30*A + x*(-2 + 16*A))))))
1368
1369static double
1370_nrrdDA4Sup(const double *parm) {
1371 double S;
1372
1373 S = parm[0];
1374 return 3*S;
1375}
1376
1377static double
1378_nrrdDA41_d(double x, const double *parm) {
1379 double S;
1380 double A;
1381 int sgn = 1;
1382
1383 S = parm[0]; A = parm[1];
1384 if (x < 0) { x = -x; sgn = -1; }
1385 x /= S;
1386 return sgn*_DAQUARTIC(x, A)(x >= 3.0 ? 0 : (x >= 2.0 ? A*(81 + x*(-90 + x*(33 - 4*
x))) : (x >= 1.0 ? -10 + 25*A + x*(18 - 66*A + x*(-10.5 + 51
*A + x*(2 - 12*A))) : x*(-6 + 12*A + x*(7.5 - 30*A + x*(-2 + 16
*A))))))/(S*S);
1387}
1388
1389static float
1390_nrrdDA41_f(float x, const double *parm) {
1391 float A, S;
1392 int sgn = 1;
1393
1394 S = AIR_CAST(float, parm[0])((float)(parm[0])); A = AIR_CAST(float, parm[1])((float)(parm[1]));
1395 if (x < 0) { x = -x; sgn = -1; }
1396 x /= S;
1397 return AIR_CAST(float, sgn*_DAQUARTIC(x, A)/(S*S))((float)(sgn*(x >= 3.0 ? 0 : (x >= 2.0 ? A*(81 + x*(-90
+ x*(33 - 4*x))) : (x >= 1.0 ? -10 + 25*A + x*(18 - 66*A +
x*(-10.5 + 51*A + x*(2 - 12*A))) : x*(-6 + 12*A + x*(7.5 - 30
*A + x*(-2 + 16*A))))))/(S*S)));
1398}
1399
1400static void
1401_nrrdDA4N_d(double *f, const double *x, size_t len, const double *parm) {
1402 double S;
1403 double t, A;
1404 size_t i;
1405 int sgn;
1406
1407 S = parm[0]; A = parm[1];
1408 for (i=0; i<len; i++) {
1409 t = x[i]/S;
1410 if (t < 0) { t = -t; sgn = -1; } else { sgn = 1; }
1411 f[i] = sgn*_DAQUARTIC(t, A)(t >= 3.0 ? 0 : (t >= 2.0 ? A*(81 + t*(-90 + t*(33 - 4*
t))) : (t >= 1.0 ? -10 + 25*A + t*(18 - 66*A + t*(-10.5 + 51
*A + t*(2 - 12*A))) : t*(-6 + 12*A + t*(7.5 - 30*A + t*(-2 + 16
*A))))))/(S*S);
1412 }
1413}
1414
1415static void
1416_nrrdDA4N_f(float *f, const float *x, size_t len, const double *parm) {
1417 float S, t, A;
1418 size_t i;
1419 int sgn;
1420
1421 S = AIR_CAST(float, parm[0])((float)(parm[0])); A = AIR_CAST(float, parm[1])((float)(parm[1]));
1422 for (i=0; i<len; i++) {
1423 t = x[i]/S;
1424 if (t < 0) { t = -t; sgn = -1; } else { sgn = 1; }
1425 f[i] = AIR_CAST(float, sgn*_DAQUARTIC(t, A)/(S*S))((float)(sgn*(t >= 3.0 ? 0 : (t >= 2.0 ? A*(81 + t*(-90
+ t*(33 - 4*t))) : (t >= 1.0 ? -10 + 25*A + t*(18 - 66*A +
t*(-10.5 + 51*A + t*(2 - 12*A))) : t*(-6 + 12*A + t*(7.5 - 30
*A + t*(-2 + 16*A))))))/(S*S)));
1426 }
1427}
1428
1429static NrrdKernel
1430_nrrdKernelDA4 = {
1431 "AquarticD",
1432 2, _nrrdDA4Sup, returnZero,
1433 _nrrdDA41_f, _nrrdDA4N_f, _nrrdDA41_d, _nrrdDA4N_d
1434};
1435NrrdKernel *const
1436nrrdKernelAQuarticD = &_nrrdKernelDA4;
1437
1438/* ------------------------------------------------------------ */
1439
1440#define _DDAQUARTIC(x, A)(x >= 3.0 ? 0 : (x >= 2.0 ? A*(-90 + x*(66 - x*12)) : (
x >= 1.0 ? 18 - 66*A + x*(-21 + 102*A + x*(6 - 36*A)) : -6
+ 12*A + x*(15 - 60*A + x*(-6 + 48*A))))) \
1441 (x >= 3.0 ? 0 : \
1442 (x >= 2.0 \
1443 ? A*(-90 + x*(66 - x*12)) \
1444 : (x >= 1.0 \
1445 ? 18 - 66*A + x*(-21 + 102*A + x*(6 - 36*A)) \
1446 : -6 + 12*A + x*(15 - 60*A + x*(-6 + 48*A)))))
1447
1448static double
1449_nrrdDDA4Sup(const double *parm) {
1450 double S;
1451
1452 S = parm[0];
1453 return 3*S;
1454}
1455
1456static double
1457_nrrdDDA41_d(double x, const double *parm) {
1458 double S;
1459 double A;
1460
1461 S = parm[0]; A = parm[1];
1462 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x))/S;
1463 return _DDAQUARTIC(x, A)(x >= 3.0 ? 0 : (x >= 2.0 ? A*(-90 + x*(66 - x*12)) : (
x >= 1.0 ? 18 - 66*A + x*(-21 + 102*A + x*(6 - 36*A)) : -6
+ 12*A + x*(15 - 60*A + x*(-6 + 48*A)))))/(S*S*S);
1464}
1465
1466static float
1467_nrrdDDA41_f(float x, const double *parm) {
1468 float S, A;
1469
1470 S = AIR_CAST(float, parm[0])((float)(parm[0])); A = AIR_CAST(float, parm[1])((float)(parm[1]));
1471 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x))/S;
1472 return _DDAQUARTIC(x, A)(x >= 3.0 ? 0 : (x >= 2.0 ? A*(-90 + x*(66 - x*12)) : (
x >= 1.0 ? 18 - 66*A + x*(-21 + 102*A + x*(6 - 36*A)) : -6
+ 12*A + x*(15 - 60*A + x*(-6 + 48*A)))))/(S*S*S);
1473}
1474
1475static void
1476_nrrdDDA4N_d(double *f, const double *x, size_t len, const double *parm) {
1477 double S;
1478 double t, A;
1479 size_t i;
1480
1481 S = parm[0]; A = parm[1];
1482 for (i=0; i<len; i++) {
1483 t = x[i];
1484 t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t))/S;
1485 f[i] = _DDAQUARTIC(t, A)(t >= 3.0 ? 0 : (t >= 2.0 ? A*(-90 + t*(66 - t*12)) : (
t >= 1.0 ? 18 - 66*A + t*(-21 + 102*A + t*(6 - 36*A)) : -6
+ 12*A + t*(15 - 60*A + t*(-6 + 48*A)))))/(S*S*S);
1486 }
1487}
1488
1489static void
1490_nrrdDDA4N_f(float *f, const float *x, size_t len, const double *parm) {
1491 float S, t, A;
1492 size_t i;
1493
1494 S = AIR_CAST(float, parm[0])((float)(parm[0])); A = AIR_CAST(float, parm[1])((float)(parm[1]));
1495 for (i=0; i<len; i++) {
1496 t = x[i];
1497 t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t))/S;
1498 f[i] = _DDAQUARTIC(t, A)(t >= 3.0 ? 0 : (t >= 2.0 ? A*(-90 + t*(66 - t*12)) : (
t >= 1.0 ? 18 - 66*A + t*(-21 + 102*A + t*(6 - 36*A)) : -6
+ 12*A + t*(15 - 60*A + t*(-6 + 48*A)))))/(S*S*S);
1499 }
1500}
1501
1502static NrrdKernel
1503_nrrdKernelDDA4 = {
1504 "AquarticDD",
1505 2, _nrrdDDA4Sup, returnZero,
1506 _nrrdDDA41_f, _nrrdDDA4N_f, _nrrdDDA41_d, _nrrdDDA4N_d
1507};
1508NrrdKernel *const
1509nrrdKernelAQuarticDD = &_nrrdKernelDDA4;
1510
1511/* ------------------------------------------------------------ */
1512
1513/*
1514** This is the unique, four-sample support, quintic, C^3 kernel, with 1st and
1515** 3rd derivatives zero at origin, which integrates to unity on [-2,2], which
1516** exactly reconstructs 0th and 1st order polynomials. Unfortunately it does
1517** NOT reconstruct 2nd order polynomials, so its not very useful. It worse
1518** than "one step down" from c4hexic (see below), since while its support and
1519** polynomial power is one less than c4hexic, it cannot reconstruct
1520** parabolas; c4hexic can reconstruct cubics.
1521**
1522** The same kernel is also available as nrrdKernelTMF w/ D,C,A = -1,3,2 ==
1523** TMF_dn_c3_2ef == "tmf:n,3,2" == nrrdKernelTMF[0][4][2], the only advantage
1524** here being that you have access to the first and second derivatives of
1525** this quintic kernel as nrrdKernelC3QuinticD and nrrdKernelC3QuinticDD.
1526**
1527** By the way, TMF_d0_c3_3ef == TMF_dn_c3_3ef == "tmf:n,3,3", which can (by
1528** definition) reconstruct parabolas, has four-sample support, and has
1529** piecewise polynomial order *seven*
1530*/
1531
1532#define _C3QUINTIC(x)(x >= 2.0 ? 0 : (x >= 1.0 ? 0.8 + x*x*(-2 + x*(2 + x*(-
0.75 + x*0.1))) : 0.7 + x*x*(-1 + x*x*(0.75 - x*0.3)))) \
1533 (x >= 2.0 ? 0 : \
1534 (x >= 1.0 \
1535 ? 0.8 + x*x*(-2 + x*(2 + x*(-0.75 + x*0.1))) \
1536 : 0.7 + x*x*(-1 + x*x*(0.75 - x*0.3))))
1537
1538static double
1539_c3quint1_d(double x, const double *parm) {
1540 AIR_UNUSED(parm)(void)(parm);
1541 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x));
1542 return _C3QUINTIC(x)(x >= 2.0 ? 0 : (x >= 1.0 ? 0.8 + x*x*(-2 + x*(2 + x*(-
0.75 + x*0.1))) : 0.7 + x*x*(-1 + x*x*(0.75 - x*0.3))));
1543}
1544
1545static float
1546_c3quint1_f(float x, const double *parm) {
1547 AIR_UNUSED(parm)(void)(parm);
1548 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x));
1549 return AIR_CAST(float, _C3QUINTIC(x))((float)((x >= 2.0 ? 0 : (x >= 1.0 ? 0.8 + x*x*(-2 + x*
(2 + x*(-0.75 + x*0.1))) : 0.7 + x*x*(-1 + x*x*(0.75 - x*0.3)
)))));
1550}
1551
1552static void
1553_c3quintN_d(double *f, const double *x, size_t len, const double *parm) {
1554 double t;
1555 size_t i;
1556 AIR_UNUSED(parm)(void)(parm);
1557 for (i=0; i<len; i++) {
1558 t = x[i];
1559 t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t));
1560 f[i] = _C3QUINTIC(t)(t >= 2.0 ? 0 : (t >= 1.0 ? 0.8 + t*t*(-2 + t*(2 + t*(-
0.75 + t*0.1))) : 0.7 + t*t*(-1 + t*t*(0.75 - t*0.3))));
1561 }
1562}
1563
1564static void
1565_c3quintN_f(float *f, const float *x, size_t len, const double *parm) {
1566 float t;
1567 size_t i;
1568 AIR_UNUSED(parm)(void)(parm);
1569 for (i=0; i<len; i++) {
1570 t = x[i];
1571 t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t));
1572 f[i] = AIR_CAST(float, _C3QUINTIC(t))((float)((t >= 2.0 ? 0 : (t >= 1.0 ? 0.8 + t*t*(-2 + t*
(2 + t*(-0.75 + t*0.1))) : 0.7 + t*t*(-1 + t*t*(0.75 - t*0.3)
)))));
1573 }
1574}
1575
1576static NrrdKernel
1577_c3quint = {
1578 "C3Quintic",
1579 0, returnTwo, returnOne,
1580 _c3quint1_f, _c3quintN_f, _c3quint1_d, _c3quintN_d
1581};
1582NrrdKernel *const
1583nrrdKernelC3Quintic = &_c3quint;
1584
1585/* ------------------------------------------------------------ */
1586
1587#define _DC3QUINTIC(x)(x >= 2.0 ? 0 : (x >= 1.0 ? x*(-4 + x*(6 + x*(-3 + x*0.5
))) : x*(-2 + x*x*(3 - x*1.5)))) \
1588 (x >= 2.0 ? 0 : \
1589 (x >= 1.0 \
1590 ? x*(-4 + x*(6 + x*(-3 + x*0.5))) \
1591 : x*(-2 + x*x*(3 - x*1.5))))
1592
1593static double
1594_Dc3quint1_d(double x, const double *parm) {
1595 int sgn = 1;
1596 AIR_UNUSED(parm)(void)(parm);
1597 if (x < 0) { x = -x; sgn = -1; }
1598 return sgn*_DC3QUINTIC(x)(x >= 2.0 ? 0 : (x >= 1.0 ? x*(-4 + x*(6 + x*(-3 + x*0.5
))) : x*(-2 + x*x*(3 - x*1.5))));
1599}
1600
1601static float
1602_Dc3quint1_f(float x, const double *parm) {
1603 int sgn = 1;
1604 AIR_UNUSED(parm)(void)(parm);
1605 if (x < 0) { x = -x; sgn = -1; }
1606 return AIR_CAST(float, sgn*_DC3QUINTIC(x))((float)(sgn*(x >= 2.0 ? 0 : (x >= 1.0 ? x*(-4 + x*(6 +
x*(-3 + x*0.5))) : x*(-2 + x*x*(3 - x*1.5))))));
1607}
1608
1609static void
1610_Dc3quintN_d(double *f, const double *x, size_t len, const double *parm) {
1611 double t;
1612 size_t i;
1613 int sgn;
1614 AIR_UNUSED(parm)(void)(parm);
1615 for (i=0; i<len; i++) {
1616 t = x[i];
1617 if (t < 0) { t = -t; sgn = -1; } else { sgn = 1; }
1618 f[i] = sgn*_DC3QUINTIC(t)(t >= 2.0 ? 0 : (t >= 1.0 ? t*(-4 + t*(6 + t*(-3 + t*0.5
))) : t*(-2 + t*t*(3 - t*1.5))));
1619 }
1620}
1621
1622static void
1623_Dc3quintN_f(float *f, const float *x, size_t len, const double *parm) {
1624 float t;
1625 size_t i;
1626 int sgn;
1627 AIR_UNUSED(parm)(void)(parm);
1628 for (i=0; i<len; i++) {
1629 t = x[i];
1630 if (t < 0) { t = -t; sgn = -1; } else { sgn = 1; }
1631 f[i] = AIR_CAST(float, sgn*_DC3QUINTIC(t))((float)(sgn*(t >= 2.0 ? 0 : (t >= 1.0 ? t*(-4 + t*(6 +
t*(-3 + t*0.5))) : t*(-2 + t*t*(3 - t*1.5))))));
1632 }
1633}
1634
1635static NrrdKernel
1636_nrrdKernelDC3Quintic = {
1637 "C3QuinticD",
1638 0, returnTwo, returnZero,
1639 _Dc3quint1_f, _Dc3quintN_f, _Dc3quint1_d, _Dc3quintN_d
1640};
1641NrrdKernel *const
1642nrrdKernelC3QuinticD = &_nrrdKernelDC3Quintic;
1643
1644/* ------------------------------------------------------------ */
1645
1646#define _DDC3QUINTIC(x)(x >= 2.0 ? 0 : (x >= 1.0 ? -4 + x*(12 + x*(-9 + x*2)) :
-2 + x*x*(9 - x*6))) \
1647 (x >= 2.0 ? 0 : \
1648 (x >= 1.0 \
1649 ? -4 + x*(12 + x*(-9 + x*2)) \
1650 : -2 + x*x*(9 - x*6)))
1651
1652static double
1653_DDc3quint1_d(double x, const double *parm) {
1654 AIR_UNUSED(parm)(void)(parm);
1655 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x));
1656 return _DDC3QUINTIC(x)(x >= 2.0 ? 0 : (x >= 1.0 ? -4 + x*(12 + x*(-9 + x*2)) :
-2 + x*x*(9 - x*6)));
1657}
1658
1659static float
1660_DDc3quint1_f(float x, const double *parm) {
1661 AIR_UNUSED(parm)(void)(parm);
1662 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x));
1663 return AIR_CAST(float, _DDC3QUINTIC(x))((float)((x >= 2.0 ? 0 : (x >= 1.0 ? -4 + x*(12 + x*(-9
+ x*2)) : -2 + x*x*(9 - x*6)))));
1664}
1665
1666static void
1667_DDc3quintN_d(double *f, const double *x, size_t len, const double *parm) {
1668 double t;
1669 size_t i;
1670 AIR_UNUSED(parm)(void)(parm);
1671 for (i=0; i<len; i++) {
1672 t = x[i];
1673 t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t));
1674 f[i] = _DDC3QUINTIC(t)(t >= 2.0 ? 0 : (t >= 1.0 ? -4 + t*(12 + t*(-9 + t*2)) :
-2 + t*t*(9 - t*6)));
1675 }
1676}
1677
1678static void
1679_DDc3quintN_f(float *f, const float *x, size_t len, const double *parm) {
1680 float t;
1681 size_t i;
1682 AIR_UNUSED(parm)(void)(parm);
1683 for (i=0; i<len; i++) {
1684 t = x[i];
1685 t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t));
1686 f[i] = AIR_CAST(float, _DDC3QUINTIC(t))((float)((t >= 2.0 ? 0 : (t >= 1.0 ? -4 + t*(12 + t*(-9
+ t*2)) : -2 + t*t*(9 - t*6)))));
1687 }
1688}
1689
1690static NrrdKernel
1691_DDc3quint = {
1692 "C3QuinticDD",
1693 0, returnTwo, returnZero,
1694 _DDc3quint1_f, _DDc3quintN_f, _DDc3quint1_d, _DDc3quintN_d
1695};
1696NrrdKernel *const
1697nrrdKernelC3QuinticDD = &_DDc3quint;
1698
1699/* ------------------------------------------------------------ */
1700
1701/*
1702** This is the unique, 6-sample support, hexic, C^4 kernel, with 1st and 3rd
1703** derivatives zero at origin, which integrates to unity on [-3,3], with 3rd
1704** order Taylor accuracy (errors start showing up when reconstructing 4th
1705** order polynomials) It doesn't interpolate, but its close, and it rings
1706** once.
1707**
1708** this is awfully close to, but not quite the same as, "tmf:n,3,4" ==
1709** TMF_dn_c3_4ef == nrrdKernelTMF[0][4][4], which is only C^3 smooth
1710*/
1711
1712#define _C4HEXIC(x)(x >= 3.0 ? 0 : (x >= 2.0 ? 1539.0/160.0 + x*(-189.0/8.0
+ x*(747.0/32.0 + x*(-12.0 + x*(109.0/32.0 + x*(-61.0/120.0 +
x/32.0))))) : (x >= 1.0 ? 3.0/160.0 + x*(35.0/8.0 + x*(-341.0
/32.0 + x*(10.0 + x*(-147.0/32.0 + x*(25.0/24.0 - x*3.0/32.0)
)))) : 69.0/80.0 + x*x*(-23.0/16.0 + x*x*(19.0/16.0 + x*(-7.0
/12.0 + x/16.0))) ))) \
1713 (x >= 3.0 \
1714 ? 0 \
1715 : (x >= 2.0 \
1716 ? 1539.0/160.0 + x*(-189.0/8.0 + x*(747.0/32.0 + x*(-12.0 + x*(109.0/32.0 + x*(-61.0/120.0 + x/32.0))))) \
1717 : (x >= 1.0 \
1718 ? 3.0/160.0 + x*(35.0/8.0 + x*(-341.0/32.0 + x*(10.0 + x*(-147.0/32.0 + x*(25.0/24.0 - x*3.0/32.0))))) \
1719 : 69.0/80.0 + x*x*(-23.0/16.0 + x*x*(19.0/16.0 + x*(-7.0/12.0 + x/16.0))) )))
1720
1721static double
1722_c4hex1_d(double x, const double *parm) {
1723 AIR_UNUSED(parm)(void)(parm);
1724 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x));
1725 return _C4HEXIC(x)(x >= 3.0 ? 0 : (x >= 2.0 ? 1539.0/160.0 + x*(-189.0/8.0
+ x*(747.0/32.0 + x*(-12.0 + x*(109.0/32.0 + x*(-61.0/120.0 +
x/32.0))))) : (x >= 1.0 ? 3.0/160.0 + x*(35.0/8.0 + x*(-341.0
/32.0 + x*(10.0 + x*(-147.0/32.0 + x*(25.0/24.0 - x*3.0/32.0)
)))) : 69.0/80.0 + x*x*(-23.0/16.0 + x*x*(19.0/16.0 + x*(-7.0
/12.0 + x/16.0))) )));
1726}
1727
1728static float
1729_c4hex1_f(float x, const double *parm) {
1730 AIR_UNUSED(parm)(void)(parm);
1731 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x));
1732 return AIR_CAST(float, _C4HEXIC(x))((float)((x >= 3.0 ? 0 : (x >= 2.0 ? 1539.0/160.0 + x*(
-189.0/8.0 + x*(747.0/32.0 + x*(-12.0 + x*(109.0/32.0 + x*(-61.0
/120.0 + x/32.0))))) : (x >= 1.0 ? 3.0/160.0 + x*(35.0/8.0
+ x*(-341.0/32.0 + x*(10.0 + x*(-147.0/32.0 + x*(25.0/24.0 -
x*3.0/32.0))))) : 69.0/80.0 + x*x*(-23.0/16.0 + x*x*(19.0/16.0
+ x*(-7.0/12.0 + x/16.0))) )))));
1733}
1734
1735static void
1736_c4hexN_d(double *f, const double *x, size_t len, const double *parm) {
1737 double t;
1738 size_t i;
1739 AIR_UNUSED(parm)(void)(parm);
1740 for (i=0; i<len; i++) {
1741 t = x[i];
1742 t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t));
1743 f[i] = _C4HEXIC(t)(t >= 3.0 ? 0 : (t >= 2.0 ? 1539.0/160.0 + t*(-189.0/8.0
+ t*(747.0/32.0 + t*(-12.0 + t*(109.0/32.0 + t*(-61.0/120.0 +
t/32.0))))) : (t >= 1.0 ? 3.0/160.0 + t*(35.0/8.0 + t*(-341.0
/32.0 + t*(10.0 + t*(-147.0/32.0 + t*(25.0/24.0 - t*3.0/32.0)
)))) : 69.0/80.0 + t*t*(-23.0/16.0 + t*t*(19.0/16.0 + t*(-7.0
/12.0 + t/16.0))) )));
1744 }
1745}
1746
1747static void
1748_c4hexN_f(float *f, const float *x, size_t len, const double *parm) {
1749 float t;
1750 size_t i;
1751 AIR_UNUSED(parm)(void)(parm);
1752 for (i=0; i<len; i++) {
1753 t = x[i];
1754 t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t));
1755 f[i] = AIR_CAST(float, _C4HEXIC(t))((float)((t >= 3.0 ? 0 : (t >= 2.0 ? 1539.0/160.0 + t*(
-189.0/8.0 + t*(747.0/32.0 + t*(-12.0 + t*(109.0/32.0 + t*(-61.0
/120.0 + t/32.0))))) : (t >= 1.0 ? 3.0/160.0 + t*(35.0/8.0
+ t*(-341.0/32.0 + t*(10.0 + t*(-147.0/32.0 + t*(25.0/24.0 -
t*3.0/32.0))))) : 69.0/80.0 + t*t*(-23.0/16.0 + t*t*(19.0/16.0
+ t*(-7.0/12.0 + t/16.0))) )))));
1756 }
1757}
1758
1759static NrrdKernel
1760_c4hex = {
1761 "C4Hexic",
1762 0, returnThree, returnOne,
1763 _c4hex1_f, _c4hexN_f, _c4hex1_d, _c4hexN_d
1764};
1765NrrdKernel *const
1766nrrdKernelC4Hexic = &_c4hex;
1767
1768/* ------------------------------------------------------------ */
1769
1770#define _DC4HEXIC(x)(x >= 3.0 ? 0 : (x >= 2.0 ? -189.0/8.0 + x*(747.0/16.0 +
x*(-36.0 + x*(109.0/8.0 + x*(-61.0/24.0 + x*(3.0/16.0))))) :
(x >= 1.0 ? 35.0/8.0 + x*(-341.0/16.0 + x*(30 + x*(-147.0
/8.0 + x*(125.0/24.0 + x*(-9.0/16.0))))) : x*(-23.0/8.0 + x*x
*(19.0/4.0 + x*(-35.0/12.0 + x*(3.0/8.0)))) ))) \
1771 (x >= 3.0 \
1772 ? 0 \
1773 : (x >= 2.0 \
1774 ? -189.0/8.0 + x*(747.0/16.0 + x*(-36.0 + x*(109.0/8.0 + x*(-61.0/24.0 + x*(3.0/16.0))))) \
1775 : (x >= 1.0 \
1776 ? 35.0/8.0 + x*(-341.0/16.0 + x*(30 + x*(-147.0/8.0 + x*(125.0/24.0 + x*(-9.0/16.0))))) \
1777 : x*(-23.0/8.0 + x*x*(19.0/4.0 + x*(-35.0/12.0 + x*(3.0/8.0)))) )))
1778
1779static double
1780_Dc4hex1_d(double x, const double *parm) {
1781 int sgn = 1;
1782 AIR_UNUSED(parm)(void)(parm);
1783 if (x < 0) { x = -x; sgn = -1; }
1784 return sgn*_DC4HEXIC(x)(x >= 3.0 ? 0 : (x >= 2.0 ? -189.0/8.0 + x*(747.0/16.0 +
x*(-36.0 + x*(109.0/8.0 + x*(-61.0/24.0 + x*(3.0/16.0))))) :
(x >= 1.0 ? 35.0/8.0 + x*(-341.0/16.0 + x*(30 + x*(-147.0
/8.0 + x*(125.0/24.0 + x*(-9.0/16.0))))) : x*(-23.0/8.0 + x*x
*(19.0/4.0 + x*(-35.0/12.0 + x*(3.0/8.0)))) )));
1785}
1786
1787static float
1788_Dc4hex1_f(float x, const double *parm) {
1789 int sgn = 1;
1790 AIR_UNUSED(parm)(void)(parm);
1791 if (x < 0) { x = -x; sgn = -1; }
1792 return AIR_CAST(float, sgn*_DC4HEXIC(x))((float)(sgn*(x >= 3.0 ? 0 : (x >= 2.0 ? -189.0/8.0 + x
*(747.0/16.0 + x*(-36.0 + x*(109.0/8.0 + x*(-61.0/24.0 + x*(3.0
/16.0))))) : (x >= 1.0 ? 35.0/8.0 + x*(-341.0/16.0 + x*(30
+ x*(-147.0/8.0 + x*(125.0/24.0 + x*(-9.0/16.0))))) : x*(-23.0
/8.0 + x*x*(19.0/4.0 + x*(-35.0/12.0 + x*(3.0/8.0)))) )))));
1793}
1794
1795static void
1796_Dc4hexN_d(double *f, const double *x, size_t len, const double *parm) {
1797 double t;
1798 size_t i;
1799 int sgn;
1800 AIR_UNUSED(parm)(void)(parm);
1801 for (i=0; i<len; i++) {
1802 t = x[i];
1803 if (t < 0) { t = -t; sgn = -1; } else { sgn = 1; }
1804 f[i] = sgn*_DC4HEXIC(t)(t >= 3.0 ? 0 : (t >= 2.0 ? -189.0/8.0 + t*(747.0/16.0 +
t*(-36.0 + t*(109.0/8.0 + t*(-61.0/24.0 + t*(3.0/16.0))))) :
(t >= 1.0 ? 35.0/8.0 + t*(-341.0/16.0 + t*(30 + t*(-147.0
/8.0 + t*(125.0/24.0 + t*(-9.0/16.0))))) : t*(-23.0/8.0 + t*t
*(19.0/4.0 + t*(-35.0/12.0 + t*(3.0/8.0)))) )));
1805 }
1806}
1807
1808static void
1809_Dc4hexN_f(float *f, const float *x, size_t len, const double *parm) {
1810 float t;
1811 size_t i;
1812 int sgn;
1813 AIR_UNUSED(parm)(void)(parm);
1814 for (i=0; i<len; i++) {
1815 t = x[i];
1816 if (t < 0) { t = -t; sgn = -1; } else { sgn = 1; }
1817 f[i] = AIR_CAST(float, sgn*_DC4HEXIC(t))((float)(sgn*(t >= 3.0 ? 0 : (t >= 2.0 ? -189.0/8.0 + t
*(747.0/16.0 + t*(-36.0 + t*(109.0/8.0 + t*(-61.0/24.0 + t*(3.0
/16.0))))) : (t >= 1.0 ? 35.0/8.0 + t*(-341.0/16.0 + t*(30
+ t*(-147.0/8.0 + t*(125.0/24.0 + t*(-9.0/16.0))))) : t*(-23.0
/8.0 + t*t*(19.0/4.0 + t*(-35.0/12.0 + t*(3.0/8.0)))) )))));
1818 }
1819}
1820
1821static NrrdKernel
1822_nrrdKernelDC4hexic = {
1823 "C4HexicD",
1824 0, returnThree, returnZero,
1825 _Dc4hex1_f, _Dc4hexN_f, _Dc4hex1_d, _Dc4hexN_d
1826};
1827NrrdKernel *const
1828nrrdKernelC4HexicD = &_nrrdKernelDC4hexic;
1829
1830/* ------------------------------------------------------------ */
1831
1832#define _DDC4HEXIC(x)(x >= 3.0 ? 0 : (x >= 2.0 ? 747.0/16.0 + x*(-72.0 + x*(
327.0/8.0 + x*(-61.0/6.0 + x*15.0/16.0))) : (x >= 1.0 ? -341.0
/16.0 + x*(60 + x*(-441.0/8.0 + x*(125.0/6.0 - x*45.0/16.0)))
: -23.0/8.0 + x*x*(57.0/4.0 + x*(-35.0/3.0 + x*(15.0/8.0))) )
)) \
1833 (x >= 3.0 \
1834 ? 0 \
1835 : (x >= 2.0 \
1836 ? 747.0/16.0 + x*(-72.0 + x*(327.0/8.0 + x*(-61.0/6.0 + x*15.0/16.0))) \
1837 : (x >= 1.0 \
1838 ? -341.0/16.0 + x*(60 + x*(-441.0/8.0 + x*(125.0/6.0 - x*45.0/16.0))) \
1839 : -23.0/8.0 + x*x*(57.0/4.0 + x*(-35.0/3.0 + x*(15.0/8.0))) )))
1840
1841static double
1842_DDc4hex1_d(double x, const double *parm) {
1843 AIR_UNUSED(parm)(void)(parm);
1844 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x));
1845 return _DDC4HEXIC(x)(x >= 3.0 ? 0 : (x >= 2.0 ? 747.0/16.0 + x*(-72.0 + x*(
327.0/8.0 + x*(-61.0/6.0 + x*15.0/16.0))) : (x >= 1.0 ? -341.0
/16.0 + x*(60 + x*(-441.0/8.0 + x*(125.0/6.0 - x*45.0/16.0)))
: -23.0/8.0 + x*x*(57.0/4.0 + x*(-35.0/3.0 + x*(15.0/8.0))) )
));
1846}
1847
1848static float
1849_DDc4hex1_f(float x, const double *parm) {
1850 AIR_UNUSED(parm)(void)(parm);
1851 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x));
1852 return AIR_CAST(float, _DDC4HEXIC(x))((float)((x >= 3.0 ? 0 : (x >= 2.0 ? 747.0/16.0 + x*(-72.0
+ x*(327.0/8.0 + x*(-61.0/6.0 + x*15.0/16.0))) : (x >= 1.0
? -341.0/16.0 + x*(60 + x*(-441.0/8.0 + x*(125.0/6.0 - x*45.0
/16.0))) : -23.0/8.0 + x*x*(57.0/4.0 + x*(-35.0/3.0 + x*(15.0
/8.0))) )))));
1853}
1854
1855static void
1856_DDc4hexN_d(double *f, const double *x, size_t len, const double *parm) {
1857 double t;
1858 size_t i;
1859 AIR_UNUSED(parm)(void)(parm);
1860 for (i=0; i<len; i++) {
1861 t = x[i];
1862 t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t));
1863 f[i] = _DDC4HEXIC(t)(t >= 3.0 ? 0 : (t >= 2.0 ? 747.0/16.0 + t*(-72.0 + t*(
327.0/8.0 + t*(-61.0/6.0 + t*15.0/16.0))) : (t >= 1.0 ? -341.0
/16.0 + t*(60 + t*(-441.0/8.0 + t*(125.0/6.0 - t*45.0/16.0)))
: -23.0/8.0 + t*t*(57.0/4.0 + t*(-35.0/3.0 + t*(15.0/8.0))) )
));
1864 }
1865}
1866
1867static void
1868_DDc4hexN_f(float *f, const float *x, size_t len, const double *parm) {
1869 float t;
1870 size_t i;
1871 AIR_UNUSED(parm)(void)(parm);
1872 for (i=0; i<len; i++) {
1873 t = x[i];
1874 t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t));
1875 f[i] = AIR_CAST(float, _DDC4HEXIC(t))((float)((t >= 3.0 ? 0 : (t >= 2.0 ? 747.0/16.0 + t*(-72.0
+ t*(327.0/8.0 + t*(-61.0/6.0 + t*15.0/16.0))) : (t >= 1.0
? -341.0/16.0 + t*(60 + t*(-441.0/8.0 + t*(125.0/6.0 - t*45.0
/16.0))) : -23.0/8.0 + t*t*(57.0/4.0 + t*(-35.0/3.0 + t*(15.0
/8.0))) )))));
1876 }
1877}
1878
1879static NrrdKernel
1880_DDc4hex = {
1881 "C4HexicDD",
1882 0, returnThree, returnZero,
1883 _DDc4hex1_f, _DDc4hexN_f, _DDc4hex1_d, _DDc4hexN_d
1884};
1885NrrdKernel *const
1886nrrdKernelC4HexicDD = &_DDc4hex;
1887
1888/* ------------------------------------------------------------ */
1889
1890#define _DDDC4HEXIC(x)(x >= 3.0 ? 0 : (x >= 2.0 ? -72.0 + x*(327.0/4.0 + x*(-
61.0/2.0 + 15*x/4)) : (x >= 1.0 ? 60 + x*(-441.0/4.0 + x*(
125.0/2.0 - 45*x/4)) : x*(57.0/2.0 + x*(-35 + 15*x/2)) ))) \
1891 (x >= 3.0 \
1892 ? 0 \
1893 : (x >= 2.0 \
1894 ? -72.0 + x*(327.0/4.0 + x*(-61.0/2.0 + 15*x/4)) \
1895 : (x >= 1.0 \
1896 ? 60 + x*(-441.0/4.0 + x*(125.0/2.0 - 45*x/4)) \
1897 : x*(57.0/2.0 + x*(-35 + 15*x/2)) \
1898 )))
1899
1900static double
1901_DDDc4hex1_d(double x, const double *parm) {
1902 int sgn = 1;
1903 AIR_UNUSED(parm)(void)(parm);
1904 if (x < 0) { x = -x; sgn = -1; }
1905 return sgn*_DDDC4HEXIC(x)(x >= 3.0 ? 0 : (x >= 2.0 ? -72.0 + x*(327.0/4.0 + x*(-
61.0/2.0 + 15*x/4)) : (x >= 1.0 ? 60 + x*(-441.0/4.0 + x*(
125.0/2.0 - 45*x/4)) : x*(57.0/2.0 + x*(-35 + 15*x/2)) )));
1906}
1907
1908static float
1909_DDDc4hex1_f(float x, const double *parm) {
1910 int sgn = 1;
1911 AIR_UNUSED(parm)(void)(parm);
1912 if (x < 0) { x = -x; sgn = -1; }
1913 return AIR_CAST(float, sgn*_DDDC4HEXIC(x))((float)(sgn*(x >= 3.0 ? 0 : (x >= 2.0 ? -72.0 + x*(327.0
/4.0 + x*(-61.0/2.0 + 15*x/4)) : (x >= 1.0 ? 60 + x*(-441.0
/4.0 + x*(125.0/2.0 - 45*x/4)) : x*(57.0/2.0 + x*(-35 + 15*x/
2)) )))));
1914}
1915
1916static void
1917_DDDc4hexN_d(double *f, const double *x, size_t len, const double *parm) {
1918 double t;
1919 size_t i;
1920 int sgn;
1921 AIR_UNUSED(parm)(void)(parm);
1922 for (i=0; i<len; i++) {
1923 t = x[i];
1924 if (t < 0) { t = -t; sgn = -1; } else { sgn = 1; }
1925 f[i] = sgn*_DDDC4HEXIC(t)(t >= 3.0 ? 0 : (t >= 2.0 ? -72.0 + t*(327.0/4.0 + t*(-
61.0/2.0 + 15*t/4)) : (t >= 1.0 ? 60 + t*(-441.0/4.0 + t*(
125.0/2.0 - 45*t/4)) : t*(57.0/2.0 + t*(-35 + 15*t/2)) )));
1926 }
1927}
1928
1929static void
1930_DDDc4hexN_f(float *f, const float *x, size_t len, const double *parm) {
1931 float t;
1932 size_t i;
1933 int sgn;
1934 AIR_UNUSED(parm)(void)(parm);
1935 for (i=0; i<len; i++) {
1936 t = x[i];
1937 if (t < 0) { t = -t; sgn = -1; } else { sgn = 1; }
1938 f[i] = AIR_CAST(float, sgn*_DDDC4HEXIC(t))((float)(sgn*(t >= 3.0 ? 0 : (t >= 2.0 ? -72.0 + t*(327.0
/4.0 + t*(-61.0/2.0 + 15*t/4)) : (t >= 1.0 ? 60 + t*(-441.0
/4.0 + t*(125.0/2.0 - 45*t/4)) : t*(57.0/2.0 + t*(-35 + 15*t/
2)) )))));
1939 }
1940}
1941
1942static NrrdKernel
1943_nrrdKernelDDDC4hexic = {
1944 "C4HexicDDD",
1945 0, returnThree, returnZero,
1946 _DDDc4hex1_f, _DDDc4hexN_f, _DDDc4hex1_d, _DDDc4hexN_d
1947};
1948NrrdKernel *const
1949nrrdKernelC4HexicDDD = &_nrrdKernelDDDC4hexic;
1950
1951
1952/* ------------- approximate inverse of c4h ------------------------- */
1953/* see comments around "_bspl3_ANI" in bsplKernel.c */
1954
1955static double
1956_c4hex_ANI_kvals[12] = {
1957 1.1906949847244948336,
1958 -0.13537708971729194940,
1959 0.047024535491780434571,
1960 -0.0088462060502312555095,
1961 0.0022443498051487024049,
1962 -0.00048639680369511914410,
1963 0.00011421848629250278186,
1964 -0.000025727759438407893986,
1965 5.9204264168395963454e-6,
1966 -1.3468552403175349134e-6,
1967 3.0637649910394681441e-7,
1968 -5.5762487950611026674e-8};
1969
1970static double
1971_c4hex_ANI_sup(const double *parm) {
1972 AIR_UNUSED(parm)(void)(parm);
1973 return 12.5;
1974}
1975
1976#define C4HEX_ANI(ret, tmp, x)tmp = ((unsigned int)(x+0.5)); if (tmp < 12) { ret = _c4hex_ANI_kvals
[tmp]; } else { ret = 0.0; } \
1977 tmp = AIR_CAST(unsigned int, x+0.5)((unsigned int)(x+0.5)); \
1978 if (tmp < 12) { \
1979 ret = _c4hex_ANI_kvals[tmp]; \
1980 } else { \
1981 ret = 0.0; \
1982 }
1983
1984static double
1985_c4hex_ANI_1d(double x, const double *parm) {
1986 double ax, r; unsigned int tmp;
1987 AIR_UNUSED(parm)(void)(parm);
1988 ax = AIR_ABS(x)((x) > 0.0f ? (x) : -(x));
1989 C4HEX_ANI(r, tmp, ax)tmp = ((unsigned int)(ax+0.5)); if (tmp < 12) { r = _c4hex_ANI_kvals
[tmp]; } else { r = 0.0; };
1990 return r;
1991}
1992
1993static float
1994_c4hex_ANI_1f(float x, const double *parm) {
1995 double ax, r; unsigned int tmp;
1996 AIR_UNUSED(parm)(void)(parm);
1997 ax = AIR_ABS(x)((x) > 0.0f ? (x) : -(x));
1998 C4HEX_ANI(r, tmp, ax)tmp = ((unsigned int)(ax+0.5)); if (tmp < 12) { r = _c4hex_ANI_kvals
[tmp]; } else { r = 0.0; };
1999 return AIR_CAST(float, r)((float)(r));
2000}
2001
2002static void
2003_c4hex_ANI_Nd(double *f, const double *x, size_t len, const double *parm) {
2004 double ax, r; unsigned int tmp;
2005 size_t i;
2006 AIR_UNUSED(parm)(void)(parm);
2007 for (i=0; i<len; i++) {
2008 ax = x[i]; ax = AIR_ABS(ax)((ax) > 0.0f ? (ax) : -(ax));
2009 C4HEX_ANI(r, tmp, ax)tmp = ((unsigned int)(ax+0.5)); if (tmp < 12) { r = _c4hex_ANI_kvals
[tmp]; } else { r = 0.0; };
2010 f[i] = r;
2011 }
2012}
2013
2014static void
2015_c4hex_ANI_Nf(float *f, const float *x, size_t len, const double *parm) {
2016 double ax, r; unsigned int tmp;
2017 size_t i;
2018 AIR_UNUSED(parm)(void)(parm);
2019 for (i=0; i<len; i++) {
2020 ax = x[i]; ax = AIR_ABS(ax)((ax) > 0.0f ? (ax) : -(ax));
2021 C4HEX_ANI(r, tmp, ax)tmp = ((unsigned int)(ax+0.5)); if (tmp < 12) { r = _c4hex_ANI_kvals
[tmp]; } else { r = 0.0; };
2022 f[i] = AIR_CAST(float, r)((float)(r));
2023 }
2024}
2025
2026static NrrdKernel
2027_nrrdKernelC4HexicApproxInverse = {
2028 "C4HexicAI", 0,
2029 _c4hex_ANI_sup, returnOne,
2030 _c4hex_ANI_1f, _c4hex_ANI_Nf,
2031 _c4hex_ANI_1d, _c4hex_ANI_Nd
2032};
2033NrrdKernel *const
2034nrrdKernelC4HexicApproxInverse = &_nrrdKernelC4HexicApproxInverse;
2035
2036/* ------------------------- c5septic ------------------------------ */
2037
2038/*
2039** This is the unique, 8-sample support, C^5 kernel, piecewise order-7
2040** with 4th order Taylor accuracy (errors start showing up when
2041** reconstructing 5th order polynomials). Coincidentally, it has zero
2042** 1st and 3rd deriv at the origin, and it integrates to unity on
2043** [-4,4]. It doesn't interpolate, but its close; it rings twice. */
2044
2045#define _C5SEPT0(x)(0.9379776601998824 + x*x*(-1.654320987654321 + x*x*(1.073045267489712
+ x*x*(-0.44997427983539096 + x*0.13978909465020575)))) (0.9379776601998824 + x*x*(-1.654320987654321 + x*x*(1.073045267489712 + x*x*(-0.44997427983539096 + x*0.13978909465020575))))
2046#define _C5SEPT1(x)(0.04651675485008818 + x*(-0.7377829218106996 + x*(0.9699074074074074
+ x*(0.18531378600823045 + x*(-0.7839506172839507 + x*(0.2357253086419753
+ x*(0.12021604938271604 - x*0.054552469135802466))))))) (0.04651675485008818 + x*(-0.7377829218106996 + x*(0.9699074074074074 + x*(0.18531378600823045 + x*(-0.7839506172839507 + x*(0.2357253086419753 + x*(0.12021604938271604 - x*0.054552469135802466)))))))
2047#define _C5SEPT2(x)(-0.01860670194003527 + x*(0.14022633744855967 + x*(-0.16296296296296298
+ x*(-0.09825102880658436 + x*(0.28858024691358025 + x*(-0.18858024691358025
+ x*(0.04405864197530864 - x*0.0013631687242798354))))))) (-0.01860670194003527 + x*(0.14022633744855967 + x*(-0.16296296296296298 + x*(-0.09825102880658436 + x*(0.28858024691358025 + x*(-0.18858024691358025 + x*(0.04405864197530864 - x*0.0013631687242798354)))))))
2048#define _C5SEPT3(x)(0.003101116990005879 + x*(-0.014223251028806585 + x*(0.02021604938271605
+ x*(0.003729423868312757 + x*(-0.0411522633744856 + x*(0.04714506172839506
+ x*(-0.023199588477366254 + x*0.004383450911228689))))))) (0.003101116990005879 + x*(-0.014223251028806585 + x*(0.02021604938271605 + x*(0.003729423868312757 + x*(-0.0411522633744856 + x*(0.04714506172839506 + x*(-0.023199588477366254 + x*0.004383450911228689)))))))
2049#define _C5SEPT(i, x)(0 == i ? (0.9379776601998824 + x*x*(-1.654320987654321 + x*x
*(1.073045267489712 + x*x*(-0.44997427983539096 + x*0.13978909465020575
)))) : (1 == i ? (0.04651675485008818 + x*(-0.7377829218106996
+ x*(0.9699074074074074 + x*(0.18531378600823045 + x*(-0.7839506172839507
+ x*(0.2357253086419753 + x*(0.12021604938271604 - x*0.054552469135802466
))))))) : (2 == i ? (-0.01860670194003527 + x*(0.14022633744855967
+ x*(-0.16296296296296298 + x*(-0.09825102880658436 + x*(0.28858024691358025
+ x*(-0.18858024691358025 + x*(0.04405864197530864 - x*0.0013631687242798354
))))))) : (3 == i ? (0.003101116990005879 + x*(-0.014223251028806585
+ x*(0.02021604938271605 + x*(0.003729423868312757 + x*(-0.0411522633744856
+ x*(0.04714506172839506 + x*(-0.023199588477366254 + x*0.004383450911228689
))))))) : 0.0)))) \
2050 (0 == i ? _C5SEPT0(x)(0.9379776601998824 + x*x*(-1.654320987654321 + x*x*(1.073045267489712
+ x*x*(-0.44997427983539096 + x*0.13978909465020575)))) \
2051 : (1 == i ? _C5SEPT1(x)(0.04651675485008818 + x*(-0.7377829218106996 + x*(0.9699074074074074
+ x*(0.18531378600823045 + x*(-0.7839506172839507 + x*(0.2357253086419753
+ x*(0.12021604938271604 - x*0.054552469135802466))))))) \
2052 : (2 == i ? _C5SEPT2(x)(-0.01860670194003527 + x*(0.14022633744855967 + x*(-0.16296296296296298
+ x*(-0.09825102880658436 + x*(0.28858024691358025 + x*(-0.18858024691358025
+ x*(0.04405864197530864 - x*0.0013631687242798354))))))) \
2053 : (3 == i ? _C5SEPT3(x)(0.003101116990005879 + x*(-0.014223251028806585 + x*(0.02021604938271605
+ x*(0.003729423868312757 + x*(-0.0411522633744856 + x*(0.04714506172839506
+ x*(-0.023199588477366254 + x*0.004383450911228689))))))) \
2054 : 0.0))))
2055
2056static double
2057_c5sept1_d(double x, const double *parm) {
2058 unsigned int xi;
2059 AIR_UNUSED(parm)(void)(parm);
2060 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x));
2061 xi = AIR_CAST(unsigned int, x)((unsigned int)(x));
2062 x -= xi;
2063 return _C5SEPT(xi, x)(0 == xi ? (0.9379776601998824 + x*x*(-1.654320987654321 + x*
x*(1.073045267489712 + x*x*(-0.44997427983539096 + x*0.13978909465020575
)))) : (1 == xi ? (0.04651675485008818 + x*(-0.7377829218106996
+ x*(0.9699074074074074 + x*(0.18531378600823045 + x*(-0.7839506172839507
+ x*(0.2357253086419753 + x*(0.12021604938271604 - x*0.054552469135802466
))))))) : (2 == xi ? (-0.01860670194003527 + x*(0.14022633744855967
+ x*(-0.16296296296296298 + x*(-0.09825102880658436 + x*(0.28858024691358025
+ x*(-0.18858024691358025 + x*(0.04405864197530864 - x*0.0013631687242798354
))))))) : (3 == xi ? (0.003101116990005879 + x*(-0.014223251028806585
+ x*(0.02021604938271605 + x*(0.003729423868312757 + x*(-0.0411522633744856
+ x*(0.04714506172839506 + x*(-0.023199588477366254 + x*0.004383450911228689
))))))) : 0.0))));
2064}
2065
2066static float
2067_c5sept1_f(float x, const double *parm) {
2068 unsigned int xi;
2069 AIR_UNUSED(parm)(void)(parm);
2070 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x));
2071 xi = AIR_CAST(unsigned int, x)((unsigned int)(x));
2072 x -= AIR_CAST(float, xi)((float)(xi));
2073 return AIR_CAST(float, _C5SEPT(xi, x))((float)((0 == xi ? (0.9379776601998824 + x*x*(-1.654320987654321
+ x*x*(1.073045267489712 + x*x*(-0.44997427983539096 + x*0.13978909465020575
)))) : (1 == xi ? (0.04651675485008818 + x*(-0.7377829218106996
+ x*(0.9699074074074074 + x*(0.18531378600823045 + x*(-0.7839506172839507
+ x*(0.2357253086419753 + x*(0.12021604938271604 - x*0.054552469135802466
))))))) : (2 == xi ? (-0.01860670194003527 + x*(0.14022633744855967
+ x*(-0.16296296296296298 + x*(-0.09825102880658436 + x*(0.28858024691358025
+ x*(-0.18858024691358025 + x*(0.04405864197530864 - x*0.0013631687242798354
))))))) : (3 == xi ? (0.003101116990005879 + x*(-0.014223251028806585
+ x*(0.02021604938271605 + x*(0.003729423868312757 + x*(-0.0411522633744856
+ x*(0.04714506172839506 + x*(-0.023199588477366254 + x*0.004383450911228689
))))))) : 0.0))))));
2074}
2075
2076static void
2077_c5septN_d(double *f, const double *x, size_t len, const double *parm) {
2078 unsigned int ti;
2079 double t;
2080 size_t i;
2081 AIR_UNUSED(parm)(void)(parm);
2082 for (i=0; i<len; i++) {
2083 t = x[i];
2084 t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t));
2085 ti = AIR_CAST(unsigned int, t)((unsigned int)(t));
2086 t -= ti;
2087 f[i] = _C5SEPT(ti, t)(0 == ti ? (0.9379776601998824 + t*t*(-1.654320987654321 + t*
t*(1.073045267489712 + t*t*(-0.44997427983539096 + t*0.13978909465020575
)))) : (1 == ti ? (0.04651675485008818 + t*(-0.7377829218106996
+ t*(0.9699074074074074 + t*(0.18531378600823045 + t*(-0.7839506172839507
+ t*(0.2357253086419753 + t*(0.12021604938271604 - t*0.054552469135802466
))))))) : (2 == ti ? (-0.01860670194003527 + t*(0.14022633744855967
+ t*(-0.16296296296296298 + t*(-0.09825102880658436 + t*(0.28858024691358025
+ t*(-0.18858024691358025 + t*(0.04405864197530864 - t*0.0013631687242798354
))))))) : (3 == ti ? (0.003101116990005879 + t*(-0.014223251028806585
+ t*(0.02021604938271605 + t*(0.003729423868312757 + t*(-0.0411522633744856
+ t*(0.04714506172839506 + t*(-0.023199588477366254 + t*0.004383450911228689
))))))) : 0.0))));
2088 }
2089}
2090
2091static void
2092_c5septN_f(float *f, const float *x, size_t len, const double *parm) {
2093 unsigned int ti;
2094 float t;
2095 size_t i;
2096 AIR_UNUSED(parm)(void)(parm);
2097 for (i=0; i<len; i++) {
2098 t = x[i];
2099 t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t));
2100 ti = AIR_CAST(unsigned int, t)((unsigned int)(t));
2101 t -= AIR_CAST(float, ti)((float)(ti));
2102 f[i] = AIR_CAST(float, _C5SEPT(ti, t))((float)((0 == ti ? (0.9379776601998824 + t*t*(-1.654320987654321
+ t*t*(1.073045267489712 + t*t*(-0.44997427983539096 + t*0.13978909465020575
)))) : (1 == ti ? (0.04651675485008818 + t*(-0.7377829218106996
+ t*(0.9699074074074074 + t*(0.18531378600823045 + t*(-0.7839506172839507
+ t*(0.2357253086419753 + t*(0.12021604938271604 - t*0.054552469135802466
))))))) : (2 == ti ? (-0.01860670194003527 + t*(0.14022633744855967
+ t*(-0.16296296296296298 + t*(-0.09825102880658436 + t*(0.28858024691358025
+ t*(-0.18858024691358025 + t*(0.04405864197530864 - t*0.0013631687242798354
))))))) : (3 == ti ? (0.003101116990005879 + t*(-0.014223251028806585
+ t*(0.02021604938271605 + t*(0.003729423868312757 + t*(-0.0411522633744856
+ t*(0.04714506172839506 + t*(-0.023199588477366254 + t*0.004383450911228689
))))))) : 0.0))))));
2103 }
2104}
2105
2106static NrrdKernel
2107_c5sept = {
2108 "C5Septic",
2109 0, returnFour, returnOne,
2110 _c5sept1_f, _c5septN_f, _c5sept1_d, _c5septN_d
2111};
2112NrrdKernel *const
2113nrrdKernelC5Septic = &_c5sept;
2114
2115#define _DC5SEPT0(x)(x*(-3.308641975308642 + x*x*(4.292181069958848 + x*x*(-2.6998456790123457
+ x*0.9785236625514403)))) (x*(-3.308641975308642 + x*x*(4.292181069958848 + x*x*(-2.6998456790123457 + x*0.9785236625514403))))
2116#define _DC5SEPT1(x)(-0.7377829218106996 + x*(1.9398148148148149 + x*(0.5559413580246914
+ x*(-3.1358024691358026 + x*(1.1786265432098766 + x*(0.7212962962962963
- x*0.3818672839506173)))))) (-0.7377829218106996 + x*(1.9398148148148149 + x*(0.5559413580246914 + x*(-3.1358024691358026 + x*(1.1786265432098766 + x*(0.7212962962962963 - x*0.3818672839506173))))))
2117#define _DC5SEPT2(x)(0.14022633744855967 + x*(-0.32592592592592595 + x*(-0.29475308641975306
+ x*(1.154320987654321 + x*(-0.9429012345679012 + x*(0.26435185185185184
- x*0.009542181069958848)))))) (0.14022633744855967 + x*(-0.32592592592592595 + x*(-0.29475308641975306 + x*(1.154320987654321 + x*(-0.9429012345679012 + x*(0.26435185185185184 - x*0.009542181069958848))))))
2118#define _DC5SEPT3(x)(-0.014223251028806585 + x*(0.0404320987654321 + x*(0.011188271604938271
+ x*(-0.1646090534979424 + x*(0.2357253086419753 + x*(-0.13919753086419753
+ x*0.03068415637860082)))))) (-0.014223251028806585 + x*(0.0404320987654321 + x*(0.011188271604938271 + x*(-0.1646090534979424 + x*(0.2357253086419753 + x*(-0.13919753086419753 + x*0.03068415637860082))))))
2119#define _DC5SEPT(i, x)(0 == i ? (x*(-3.308641975308642 + x*x*(4.292181069958848 + x
*x*(-2.6998456790123457 + x*0.9785236625514403)))) : (1 == i ?
(-0.7377829218106996 + x*(1.9398148148148149 + x*(0.5559413580246914
+ x*(-3.1358024691358026 + x*(1.1786265432098766 + x*(0.7212962962962963
- x*0.3818672839506173)))))) : (2 == i ? (0.14022633744855967
+ x*(-0.32592592592592595 + x*(-0.29475308641975306 + x*(1.154320987654321
+ x*(-0.9429012345679012 + x*(0.26435185185185184 - x*0.009542181069958848
)))))) : (3 == i ? (-0.014223251028806585 + x*(0.0404320987654321
+ x*(0.011188271604938271 + x*(-0.1646090534979424 + x*(0.2357253086419753
+ x*(-0.13919753086419753 + x*0.03068415637860082)))))) : 0.0
)))) \
2120 (0 == i ? _DC5SEPT0(x)(x*(-3.308641975308642 + x*x*(4.292181069958848 + x*x*(-2.6998456790123457
+ x*0.9785236625514403)))) \
2121 : (1 == i ? _DC5SEPT1(x)(-0.7377829218106996 + x*(1.9398148148148149 + x*(0.5559413580246914
+ x*(-3.1358024691358026 + x*(1.1786265432098766 + x*(0.7212962962962963
- x*0.3818672839506173)))))) \
2122 : (2 == i ? _DC5SEPT2(x)(0.14022633744855967 + x*(-0.32592592592592595 + x*(-0.29475308641975306
+ x*(1.154320987654321 + x*(-0.9429012345679012 + x*(0.26435185185185184
- x*0.009542181069958848)))))) \
2123 : (3 == i ? _DC5SEPT3(x)(-0.014223251028806585 + x*(0.0404320987654321 + x*(0.011188271604938271
+ x*(-0.1646090534979424 + x*(0.2357253086419753 + x*(-0.13919753086419753
+ x*0.03068415637860082)))))) \
2124 : 0.0))))
2125
2126static double
2127_dc5sept1_d(double x, const double *parm) {
2128 unsigned int xi;
2129 int sgn = 1;
2130 AIR_UNUSED(parm)(void)(parm);
2131 if (x < 0) { x = -x; sgn = -1; }
2132 xi = AIR_CAST(unsigned int, x)((unsigned int)(x));
2133 x -= xi;
2134 return sgn*_DC5SEPT(xi, x)(0 == xi ? (x*(-3.308641975308642 + x*x*(4.292181069958848 + x
*x*(-2.6998456790123457 + x*0.9785236625514403)))) : (1 == xi
? (-0.7377829218106996 + x*(1.9398148148148149 + x*(0.5559413580246914
+ x*(-3.1358024691358026 + x*(1.1786265432098766 + x*(0.7212962962962963
- x*0.3818672839506173)))))) : (2 == xi ? (0.14022633744855967
+ x*(-0.32592592592592595 + x*(-0.29475308641975306 + x*(1.154320987654321
+ x*(-0.9429012345679012 + x*(0.26435185185185184 - x*0.009542181069958848
)))))) : (3 == xi ? (-0.014223251028806585 + x*(0.0404320987654321
+ x*(0.011188271604938271 + x*(-0.1646090534979424 + x*(0.2357253086419753
+ x*(-0.13919753086419753 + x*0.03068415637860082)))))) : 0.0
))));
2135}
2136
2137static float
2138_dc5sept1_f(float x, const double *parm) {
2139 unsigned int xi;
2140 int sgn = 1;
2141 AIR_UNUSED(parm)(void)(parm);
2142 if (x < 0) { x = -x; sgn = -1; }
2143 xi = AIR_CAST(unsigned int, x)((unsigned int)(x));
2144 x -= AIR_CAST(float, xi)((float)(xi));
2145 return AIR_CAST(float, sgn*_DC5SEPT(xi, x))((float)(sgn*(0 == xi ? (x*(-3.308641975308642 + x*x*(4.292181069958848
+ x*x*(-2.6998456790123457 + x*0.9785236625514403)))) : (1 ==
xi ? (-0.7377829218106996 + x*(1.9398148148148149 + x*(0.5559413580246914
+ x*(-3.1358024691358026 + x*(1.1786265432098766 + x*(0.7212962962962963
- x*0.3818672839506173)))))) : (2 == xi ? (0.14022633744855967
+ x*(-0.32592592592592595 + x*(-0.29475308641975306 + x*(1.154320987654321
+ x*(-0.9429012345679012 + x*(0.26435185185185184 - x*0.009542181069958848
)))))) : (3 == xi ? (-0.014223251028806585 + x*(0.0404320987654321
+ x*(0.011188271604938271 + x*(-0.1646090534979424 + x*(0.2357253086419753
+ x*(-0.13919753086419753 + x*0.03068415637860082)))))) : 0.0
))))));
2146}
2147
2148static void
2149_dc5septN_d(double *f, const double *x, size_t len, const double *parm) {
2150 unsigned int ti;
2151 double t;
2152 size_t i;
2153 int sgn;
2154 AIR_UNUSED(parm)(void)(parm);
2155 for (i=0; i<len; i++) {
2156 t = x[i];
2157 if (t < 0) { t = -t; sgn = -1; } else { sgn = 1; }
2158 ti = AIR_CAST(unsigned int, t)((unsigned int)(t));
2159 t -= ti;
2160 f[i] = sgn*_DC5SEPT(ti, t)(0 == ti ? (t*(-3.308641975308642 + t*t*(4.292181069958848 + t
*t*(-2.6998456790123457 + t*0.9785236625514403)))) : (1 == ti
? (-0.7377829218106996 + t*(1.9398148148148149 + t*(0.5559413580246914
+ t*(-3.1358024691358026 + t*(1.1786265432098766 + t*(0.7212962962962963
- t*0.3818672839506173)))))) : (2 == ti ? (0.14022633744855967
+ t*(-0.32592592592592595 + t*(-0.29475308641975306 + t*(1.154320987654321
+ t*(-0.9429012345679012 + t*(0.26435185185185184 - t*0.009542181069958848
)))))) : (3 == ti ? (-0.014223251028806585 + t*(0.0404320987654321
+ t*(0.011188271604938271 + t*(-0.1646090534979424 + t*(0.2357253086419753
+ t*(-0.13919753086419753 + t*0.03068415637860082)))))) : 0.0
))));
2161 }
2162}
2163
2164static void
2165_dc5septN_f(float *f, const float *x, size_t len, const double *parm) {
2166 unsigned int ti;
2167 float t;
2168 size_t i;
2169 int sgn = 1;
2170 AIR_UNUSED(parm)(void)(parm);
2171 for (i=0; i<len; i++) {
2172 t = x[i];
2173 if (t < 0) { t = -t; sgn = -1; } else { sgn = 1; }
2174 ti = AIR_CAST(unsigned int, t)((unsigned int)(t));
2175 t -= AIR_CAST(float, ti)((float)(ti));
2176 f[i] = AIR_CAST(float, sgn*_DC5SEPT(ti, t))((float)(sgn*(0 == ti ? (t*(-3.308641975308642 + t*t*(4.292181069958848
+ t*t*(-2.6998456790123457 + t*0.9785236625514403)))) : (1 ==
ti ? (-0.7377829218106996 + t*(1.9398148148148149 + t*(0.5559413580246914
+ t*(-3.1358024691358026 + t*(1.1786265432098766 + t*(0.7212962962962963
- t*0.3818672839506173)))))) : (2 == ti ? (0.14022633744855967
+ t*(-0.32592592592592595 + t*(-0.29475308641975306 + t*(1.154320987654321
+ t*(-0.9429012345679012 + t*(0.26435185185185184 - t*0.009542181069958848
)))))) : (3 == ti ? (-0.014223251028806585 + t*(0.0404320987654321
+ t*(0.011188271604938271 + t*(-0.1646090534979424 + t*(0.2357253086419753
+ t*(-0.13919753086419753 + t*0.03068415637860082)))))) : 0.0
))))));
2177 }
2178}
2179
2180static NrrdKernel
2181_dc5sept = {
2182 "C5SepticD",
2183 0, returnFour, returnZero,
2184 _dc5sept1_f, _dc5septN_f, _dc5sept1_d, _dc5septN_d
2185};
2186NrrdKernel *const
2187nrrdKernelC5SepticD = &_dc5sept;
2188
2189#define _DDC5SEPT0(x)(-3.308641975308642 + x*x*(12.876543209876543 + x*x*(-13.499228395061728
+ x*5.871141975308642))) (-3.308641975308642 + x*x*(12.876543209876543 + x*x*(-13.499228395061728 + x*5.871141975308642)))
2190#define _DDC5SEPT1(x)(1.9398148148148149 + x*(1.1118827160493827 + x*(-9.407407407407407
+ x*(4.714506172839506 + x*(3.6064814814814814 - x*2.2912037037037036
))))) (1.9398148148148149 + x*(1.1118827160493827 + x*(-9.407407407407407 + x*(4.714506172839506 + x*(3.6064814814814814 - x*2.2912037037037036)))))
2191#define _DDC5SEPT2(x)(-0.32592592592592595 + x*(-0.5895061728395061 + x*(3.462962962962963
+ x*(-3.771604938271605 + x*(1.3217592592592593 - x*0.05725308641975309
))))) (-0.32592592592592595 + x*(-0.5895061728395061 + x*(3.462962962962963 + x*(-3.771604938271605 + x*(1.3217592592592593 - x*0.05725308641975309)))))
2192#define _DDC5SEPT3(x)(0.0404320987654321 + x*(0.022376543209876542 + x*(-0.49382716049382713
+ x*(0.9429012345679012 + x*(-0.6959876543209876 + x*0.18410493827160493
))))) (0.0404320987654321 + x*(0.022376543209876542 + x*(-0.49382716049382713 + x*(0.9429012345679012 + x*(-0.6959876543209876 + x*0.18410493827160493)))))
2193#define _DDC5SEPT(i, x)(0 == i ? (-3.308641975308642 + x*x*(12.876543209876543 + x*x
*(-13.499228395061728 + x*5.871141975308642))) : (1 == i ? (1.9398148148148149
+ x*(1.1118827160493827 + x*(-9.407407407407407 + x*(4.714506172839506
+ x*(3.6064814814814814 - x*2.2912037037037036))))) : (2 == i
? (-0.32592592592592595 + x*(-0.5895061728395061 + x*(3.462962962962963
+ x*(-3.771604938271605 + x*(1.3217592592592593 - x*0.05725308641975309
))))) : (3 == i ? (0.0404320987654321 + x*(0.022376543209876542
+ x*(-0.49382716049382713 + x*(0.9429012345679012 + x*(-0.6959876543209876
+ x*0.18410493827160493))))) : 0.0)))) \
2194 (0 == i ? _DDC5SEPT0(x)(-3.308641975308642 + x*x*(12.876543209876543 + x*x*(-13.499228395061728
+ x*5.871141975308642))) \
2195 : (1 == i ? _DDC5SEPT1(x)(1.9398148148148149 + x*(1.1118827160493827 + x*(-9.407407407407407
+ x*(4.714506172839506 + x*(3.6064814814814814 - x*2.2912037037037036
))))) \
2196 : (2 == i ? _DDC5SEPT2(x)(-0.32592592592592595 + x*(-0.5895061728395061 + x*(3.462962962962963
+ x*(-3.771604938271605 + x*(1.3217592592592593 - x*0.05725308641975309
))))) \
2197 : (3 == i ? _DDC5SEPT3(x)(0.0404320987654321 + x*(0.022376543209876542 + x*(-0.49382716049382713
+ x*(0.9429012345679012 + x*(-0.6959876543209876 + x*0.18410493827160493
))))) \
2198 : 0.0))))
2199
2200static double
2201_ddc5sept1_d(double x, const double *parm) {
2202 unsigned int xi;
2203 AIR_UNUSED(parm)(void)(parm);
2204 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x));
2205 xi = AIR_CAST(unsigned int, x)((unsigned int)(x));
2206 x -= xi;
2207 return _DDC5SEPT(xi, x)(0 == xi ? (-3.308641975308642 + x*x*(12.876543209876543 + x*
x*(-13.499228395061728 + x*5.871141975308642))) : (1 == xi ? (
1.9398148148148149 + x*(1.1118827160493827 + x*(-9.407407407407407
+ x*(4.714506172839506 + x*(3.6064814814814814 - x*2.2912037037037036
))))) : (2 == xi ? (-0.32592592592592595 + x*(-0.5895061728395061
+ x*(3.462962962962963 + x*(-3.771604938271605 + x*(1.3217592592592593
- x*0.05725308641975309))))) : (3 == xi ? (0.0404320987654321
+ x*(0.022376543209876542 + x*(-0.49382716049382713 + x*(0.9429012345679012
+ x*(-0.6959876543209876 + x*0.18410493827160493))))) : 0.0)
)));
2208}
2209
2210static float
2211_ddc5sept1_f(float x, const double *parm) {
2212 unsigned int xi;
2213 AIR_UNUSED(parm)(void)(parm);
2214 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x));
2215 xi = AIR_CAST(unsigned int, x)((unsigned int)(x));
2216 x -= AIR_CAST(float, xi)((float)(xi));
2217 return AIR_CAST(float, _DDC5SEPT(xi, x))((float)((0 == xi ? (-3.308641975308642 + x*x*(12.876543209876543
+ x*x*(-13.499228395061728 + x*5.871141975308642))) : (1 == xi
? (1.9398148148148149 + x*(1.1118827160493827 + x*(-9.407407407407407
+ x*(4.714506172839506 + x*(3.6064814814814814 - x*2.2912037037037036
))))) : (2 == xi ? (-0.32592592592592595 + x*(-0.5895061728395061
+ x*(3.462962962962963 + x*(-3.771604938271605 + x*(1.3217592592592593
- x*0.05725308641975309))))) : (3 == xi ? (0.0404320987654321
+ x*(0.022376543209876542 + x*(-0.49382716049382713 + x*(0.9429012345679012
+ x*(-0.6959876543209876 + x*0.18410493827160493))))) : 0.0)
)))));
2218}
2219
2220static void
2221_ddc5septN_d(double *f, const double *x, size_t len, const double *parm) {
2222 unsigned int ti;
2223 double t;
2224 size_t i;
2225 AIR_UNUSED(parm)(void)(parm);
2226 for (i=0; i<len; i++) {
2227 t = x[i];
2228 t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t));
2229 ti = AIR_CAST(unsigned int, t)((unsigned int)(t));
2230 t -= ti;
2231 f[i] = _DDC5SEPT(ti, t)(0 == ti ? (-3.308641975308642 + t*t*(12.876543209876543 + t*
t*(-13.499228395061728 + t*5.871141975308642))) : (1 == ti ? (
1.9398148148148149 + t*(1.1118827160493827 + t*(-9.407407407407407
+ t*(4.714506172839506 + t*(3.6064814814814814 - t*2.2912037037037036
))))) : (2 == ti ? (-0.32592592592592595 + t*(-0.5895061728395061
+ t*(3.462962962962963 + t*(-3.771604938271605 + t*(1.3217592592592593
- t*0.05725308641975309))))) : (3 == ti ? (0.0404320987654321
+ t*(0.022376543209876542 + t*(-0.49382716049382713 + t*(0.9429012345679012
+ t*(-0.6959876543209876 + t*0.18410493827160493))))) : 0.0)
)));
2232 }
2233}
2234
2235static void
2236_ddc5septN_f(float *f, const float *x, size_t len, const double *parm) {
2237 unsigned int ti;
2238 float t;
2239 size_t i;
2240 AIR_UNUSED(parm)(void)(parm);
2241 for (i=0; i<len; i++) {
2242 t = x[i];
2243 t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t));
2244 ti = AIR_CAST(unsigned int, t)((unsigned int)(t));
2245 t -= AIR_CAST(float, ti)((float)(ti));
2246 f[i] = AIR_CAST(float, _DDC5SEPT(ti, t))((float)((0 == ti ? (-3.308641975308642 + t*t*(12.876543209876543
+ t*t*(-13.499228395061728 + t*5.871141975308642))) : (1 == ti
? (1.9398148148148149 + t*(1.1118827160493827 + t*(-9.407407407407407
+ t*(4.714506172839506 + t*(3.6064814814814814 - t*2.2912037037037036
))))) : (2 == ti ? (-0.32592592592592595 + t*(-0.5895061728395061
+ t*(3.462962962962963 + t*(-3.771604938271605 + t*(1.3217592592592593
- t*0.05725308641975309))))) : (3 == ti ? (0.0404320987654321
+ t*(0.022376543209876542 + t*(-0.49382716049382713 + t*(0.9429012345679012
+ t*(-0.6959876543209876 + t*0.18410493827160493))))) : 0.0)
)))));
2247 }
2248}
2249
2250static NrrdKernel
2251_ddc5sept = {
2252 "C5SepticDD",
2253 0, returnFour, returnZero,
2254 _ddc5sept1_f, _ddc5septN_f, _ddc5sept1_d, _ddc5septN_d
2255};
2256NrrdKernel *const
2257nrrdKernelC5SepticDD = &_ddc5sept;
2258
2259#define _DDDC5SEPT0(x)(x*(25.75308641975309 + x*x*(-53.99691358024691 + x*29.35570987654321
))) (x*(25.75308641975309 + x*x*(-53.99691358024691 + x*29.35570987654321)))
2260#define _DDDC5SEPT1(x)(1.111882716049383 + x*(-18.81481481481481 + x*(14.14351851851852
+ x*(14.42592592592593 - x*11.45601851851852)))) (1.111882716049383 + x*(-18.81481481481481 + x*(14.14351851851852 + x*(14.42592592592593 - x*11.45601851851852))))
2261#define _DDDC5SEPT2(x)(-0.5895061728395062 + x*(6.925925925925926 + x*(-11.31481481481481
+ x*(5.287037037037037 - x*0.2862654320987654)))) (-0.5895061728395062 + x*(6.925925925925926 + x*(-11.31481481481481 + x*(5.287037037037037 - x*0.2862654320987654))))
2262#define _DDDC5SEPT3(x)(0.02237654320987654 + x*(-0.9876543209876543 + x*(2.828703703703704
+ x*(-2.783950617283951 + x*0.9205246913580247)))) (0.02237654320987654 + x*(-0.9876543209876543 + x*(2.828703703703704 + x*(-2.783950617283951 + x*0.9205246913580247))))
2263#define _DDDC5SEPT(i, x)(0 == i ? (x*(25.75308641975309 + x*x*(-53.99691358024691 + x
*29.35570987654321))) : (1 == i ? (1.111882716049383 + x*(-18.81481481481481
+ x*(14.14351851851852 + x*(14.42592592592593 - x*11.45601851851852
)))) : (2 == i ? (-0.5895061728395062 + x*(6.925925925925926 +
x*(-11.31481481481481 + x*(5.287037037037037 - x*0.2862654320987654
)))) : (3 == i ? (0.02237654320987654 + x*(-0.9876543209876543
+ x*(2.828703703703704 + x*(-2.783950617283951 + x*0.9205246913580247
)))) : 0.0)))) \
2264 (0 == i ? _DDDC5SEPT0(x)(x*(25.75308641975309 + x*x*(-53.99691358024691 + x*29.35570987654321
))) \
2265 : (1 == i ? _DDDC5SEPT1(x)(1.111882716049383 + x*(-18.81481481481481 + x*(14.14351851851852
+ x*(14.42592592592593 - x*11.45601851851852)))) \
2266 : (2 == i ? _DDDC5SEPT2(x)(-0.5895061728395062 + x*(6.925925925925926 + x*(-11.31481481481481
+ x*(5.287037037037037 - x*0.2862654320987654)))) \
2267 : (3 == i ? _DDDC5SEPT3(x)(0.02237654320987654 + x*(-0.9876543209876543 + x*(2.828703703703704
+ x*(-2.783950617283951 + x*0.9205246913580247)))) \
2268 : 0.0))))
2269
2270static double
2271_dddc5sept1_d(double x, const double *parm) {
2272 unsigned int xi;
2273 int sgn = 1;
2274 AIR_UNUSED(parm)(void)(parm);
2275 if (x < 0) { x = -x; sgn = -1; }
2276 xi = AIR_CAST(unsigned int, x)((unsigned int)(x));
2277 x -= xi;
2278 return sgn*_DDDC5SEPT(xi, x)(0 == xi ? (x*(25.75308641975309 + x*x*(-53.99691358024691 + x
*29.35570987654321))) : (1 == xi ? (1.111882716049383 + x*(-18.81481481481481
+ x*(14.14351851851852 + x*(14.42592592592593 - x*11.45601851851852
)))) : (2 == xi ? (-0.5895061728395062 + x*(6.925925925925926
+ x*(-11.31481481481481 + x*(5.287037037037037 - x*0.2862654320987654
)))) : (3 == xi ? (0.02237654320987654 + x*(-0.9876543209876543
+ x*(2.828703703703704 + x*(-2.783950617283951 + x*0.9205246913580247
)))) : 0.0))));
2279}
2280
2281static float
2282_dddc5sept1_f(float x, const double *parm) {
2283 unsigned int xi;
2284 int sgn = 1;
2285 AIR_UNUSED(parm)(void)(parm);
2286 if (x < 0) { x = -x; sgn = -1; }
2287 xi = AIR_CAST(unsigned int, x)((unsigned int)(x));
2288 x -= AIR_CAST(float, xi)((float)(xi));
2289 return AIR_CAST(float, sgn*_DDDC5SEPT(xi, x))((float)(sgn*(0 == xi ? (x*(25.75308641975309 + x*x*(-53.99691358024691
+ x*29.35570987654321))) : (1 == xi ? (1.111882716049383 + x
*(-18.81481481481481 + x*(14.14351851851852 + x*(14.42592592592593
- x*11.45601851851852)))) : (2 == xi ? (-0.5895061728395062 +
x*(6.925925925925926 + x*(-11.31481481481481 + x*(5.287037037037037
- x*0.2862654320987654)))) : (3 == xi ? (0.02237654320987654
+ x*(-0.9876543209876543 + x*(2.828703703703704 + x*(-2.783950617283951
+ x*0.9205246913580247)))) : 0.0))))));
2290}
2291
2292static void
2293_dddc5septN_d(double *f, const double *x, size_t len, const double *parm) {
2294 unsigned int ti;
2295 double t;
2296 size_t i;
2297 int sgn;
2298 AIR_UNUSED(parm)(void)(parm);
2299 for (i=0; i<len; i++) {
2300 t = x[i];
2301 if (t < 0) { t = -t; sgn = -1; } else { sgn = 1; }
2302 ti = AIR_CAST(unsigned int, t)((unsigned int)(t));
2303 t -= ti;
2304 f[i] = sgn*_DDDC5SEPT(ti, t)(0 == ti ? (t*(25.75308641975309 + t*t*(-53.99691358024691 + t
*29.35570987654321))) : (1 == ti ? (1.111882716049383 + t*(-18.81481481481481
+ t*(14.14351851851852 + t*(14.42592592592593 - t*11.45601851851852
)))) : (2 == ti ? (-0.5895061728395062 + t*(6.925925925925926
+ t*(-11.31481481481481 + t*(5.287037037037037 - t*0.2862654320987654
)))) : (3 == ti ? (0.02237654320987654 + t*(-0.9876543209876543
+ t*(2.828703703703704 + t*(-2.783950617283951 + t*0.9205246913580247
)))) : 0.0))));
2305 }
2306}
2307
2308static void
2309_dddc5septN_f(float *f, const float *x, size_t len, const double *parm) {
2310 unsigned int ti;
2311 float t;
2312 size_t i;
2313 int sgn = 1;
2314 AIR_UNUSED(parm)(void)(parm);
2315 for (i=0; i<len; i++) {
2316 t = x[i];
2317 if (t < 0) { t = -t; sgn = -1; } else { sgn = 1; }
2318 ti = AIR_CAST(unsigned int, t)((unsigned int)(t));
2319 t -= AIR_CAST(float, ti)((float)(ti));
2320 f[i] = AIR_CAST(float, sgn*_DDDC5SEPT(ti, t))((float)(sgn*(0 == ti ? (t*(25.75308641975309 + t*t*(-53.99691358024691
+ t*29.35570987654321))) : (1 == ti ? (1.111882716049383 + t
*(-18.81481481481481 + t*(14.14351851851852 + t*(14.42592592592593
- t*11.45601851851852)))) : (2 == ti ? (-0.5895061728395062 +
t*(6.925925925925926 + t*(-11.31481481481481 + t*(5.287037037037037
- t*0.2862654320987654)))) : (3 == ti ? (0.02237654320987654
+ t*(-0.9876543209876543 + t*(2.828703703703704 + t*(-2.783950617283951
+ t*0.9205246913580247)))) : 0.0))))));
2321 }
2322}
2323
2324static NrrdKernel
2325_dddc5sept = {
2326 "C5SepticDDD",
2327 0, returnFour, returnZero,
2328 _dddc5sept1_f, _dddc5septN_f, _dddc5sept1_d, _dddc5septN_d
2329};
2330NrrdKernel *const
2331nrrdKernelC5SepticDDD = &_dddc5sept;
2332
2333/* note that this implies a much more accurate inverse than is given
2334 for the splines or other kernels; this is a consequence of GLK
2335 re-purposing the Mathematica expressions for the bpsln7 inverse,
2336 which is unfortunate: c5septic is nearly interpolating, so far
2337 fewer terms would suffice */
2338#define C5SEPT_AI_LEN26 26
2339static double
2340_c5sept_ANI_kvals[C5SEPT_AI_LEN26] = {
2341 1.072662863909143543451743,
2342 -0.05572032001443187610952953,
2343 0.02453993146603215267432700,
2344 -0.005922375635530229254855750,
2345 0.0009781882769025851448681918,
2346 -0.0002499281491108793120774480,
2347 0.00005196973116762945530292666,
2348 -0.00001090007030248955413371701,
2349 2.425581976693179040189747e-6,
2350 -5.143328756144314306529358e-7,
2351 1.109572595055083858393193e-7,
2352 -2.400323559797703961855318e-8,
2353 5.151959978856239469272136e-9,
2354 -1.111431289951609447815300e-9,
2355 2.394624249806782312051293e-10,
2356 -5.155654818408366273965675e-11,
2357 1.111091879440261302739584e-11,
2358 -2.393303168472914213194646e-12,
2359 5.155527554392687415035171e-13,
2360 -1.110707782883835600110634e-13,
2361 2.392644268109899456937863e-14,
2362 -5.154377464414088544322752e-15,
2363 1.110376957658466603291262e-15,
2364 -2.391459139266885907929865e-16,
2365 5.137528165538909945741180e-17,
2366 -9.024576392408067896802608e-18};
2367
2368static double
2369_c5sept_ANI_sup(const double *parm) {
2370 AIR_UNUSED(parm)(void)(parm);
2371 return C5SEPT_AI_LEN26 + 0.5;
2372}
2373
2374static double
2375_c5sept_ANI_int(const double *parm) {
2376 AIR_UNUSED(parm)(void)(parm);
2377 return 1.0;
2378}
2379
2380#define C5SEPT_ANI(ret, tmp, x)tmp = ((unsigned int)(x+0.5)); if (tmp < 26) { ret = _c5sept_ANI_kvals
[tmp]; } else { ret = 0.0; } \
2381 tmp = AIR_CAST(unsigned int, x+0.5)((unsigned int)(x+0.5)); \
2382 if (tmp < C5SEPT_AI_LEN26) { \
2383 ret = _c5sept_ANI_kvals[tmp]; \
2384 } else { \
2385 ret = 0.0; \
2386 }
2387
2388static double
2389_c5sept_ANI_1d(double x, const double *parm) {
2390 double ax, r; unsigned int tmp;
2391 AIR_UNUSED(parm)(void)(parm);
2392
2393 ax = AIR_ABS(x)((x) > 0.0f ? (x) : -(x));
2394 C5SEPT_ANI(r, tmp, ax)tmp = ((unsigned int)(ax+0.5)); if (tmp < 26) { r = _c5sept_ANI_kvals
[tmp]; } else { r = 0.0; };
2395 return r;
2396}
2397
2398static float
2399_c5sept_ANI_1f(float x, const double *parm) {
2400 double ax, r; unsigned int tmp;
2401 AIR_UNUSED(parm)(void)(parm);
2402
2403 ax = AIR_ABS(x)((x) > 0.0f ? (x) : -(x));
2404 C5SEPT_ANI(r, tmp, ax)tmp = ((unsigned int)(ax+0.5)); if (tmp < 26) { r = _c5sept_ANI_kvals
[tmp]; } else { r = 0.0; };
2405 return AIR_CAST(float, r)((float)(r));
2406}
2407
2408static void
2409_c5sept_ANI_Nd(double *f, const double *x, size_t len, const double *parm) {
2410 double ax, r; unsigned int tmp;
2411 size_t i;
2412 AIR_UNUSED(parm)(void)(parm);
2413
2414 for (i=0; i<len; i++) {
2415 ax = x[i]; ax = AIR_ABS(ax)((ax) > 0.0f ? (ax) : -(ax));
2416 C5SEPT_ANI(r, tmp, ax)tmp = ((unsigned int)(ax+0.5)); if (tmp < 26) { r = _c5sept_ANI_kvals
[tmp]; } else { r = 0.0; };
2417 f[i] = r;
2418 }
2419}
2420
2421static void
2422_c5sept_ANI_Nf(float *f, const float *x, size_t len, const double *parm) {
2423 double ax, r; unsigned int tmp;
2424 size_t i;
2425 AIR_UNUSED(parm)(void)(parm);
2426
2427 for (i=0; i<len; i++) {
2428 ax = x[i]; ax = AIR_ABS(ax)((ax) > 0.0f ? (ax) : -(ax));
2429 C5SEPT_ANI(r, tmp, ax)tmp = ((unsigned int)(ax+0.5)); if (tmp < 26) { r = _c5sept_ANI_kvals
[tmp]; } else { r = 0.0; };
2430 f[i] = AIR_CAST(float, r)((float)(r));
2431 }
2432}
2433
2434static NrrdKernel
2435_nrrdKernelC5SepticApproxInverse = {
2436 "C5SepticAI", 0,
2437 _c5sept_ANI_sup, _c5sept_ANI_int,
2438 _c5sept_ANI_1f, _c5sept_ANI_Nf,
2439 _c5sept_ANI_1d, _c5sept_ANI_Nd
2440};
2441NrrdKernel *const
2442nrrdKernelC5SepticApproxInverse = &_nrrdKernelC5SepticApproxInverse;
2443
2444/* ------------------------------------------------------------ */
2445
2446#define _GAUSS(x, sig, cut)( x >= sig*cut ? 0 : exp(-x*x/(2.0*sig*sig))/(sig*2.50662827463100050241
)) ( \
2447 x >= sig*cut ? 0 \
2448 : exp(-x*x/(2.0*sig*sig))/(sig*2.50662827463100050241))
2449
2450static double
2451_nrrdGInt(const double *parm) {
2452 double cut;
2453
2454 cut = parm[1];
2455 return airErf(cut/sqrt(2.0));
2456}
2457
2458static double
2459_nrrdGSup(const double *parm) {
2460 double sig, cut;
2461
2462 sig = parm[0];
2463 cut = parm[1];
2464 return sig*cut;
2465}
2466
2467static double
2468_nrrdG1_d(double x, const double *parm) {
2469 double sig, cut;
2470
2471 sig = parm[0];
2472 cut = parm[1];
2473 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x));
2474 return _GAUSS(x, sig, cut)( x >= sig*cut ? 0 : exp(-x*x/(2.0*sig*sig))/(sig*2.50662827463100050241
));
2475}
2476
2477static float
2478_nrrdG1_f(float x, const double *parm) {
2479 float sig, cut;
2480
2481 sig = AIR_CAST(float, parm[0])((float)(parm[0]));
2482 cut = AIR_CAST(float, parm[1])((float)(parm[1]));
2483 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x));
2484 return AIR_CAST(float, _GAUSS(x, sig, cut))((float)(( x >= sig*cut ? 0 : exp(-x*x/(2.0*sig*sig))/(sig
*2.50662827463100050241))));
2485}
2486
2487static void
2488_nrrdGN_d(double *f, const double *x, size_t len, const double *parm) {
2489 double sig, cut, t;
2490 size_t i;
2491
2492 sig = parm[0];
2493 cut = parm[1];
2494 for (i=0; i<len; i++) {
2495 t = x[i];
2496 t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t));
2497 f[i] = _GAUSS(t, sig, cut)( t >= sig*cut ? 0 : exp(-t*t/(2.0*sig*sig))/(sig*2.50662827463100050241
));
2498 }
2499}
2500
2501static void
2502_nrrdGN_f(float *f, const float *x, size_t len, const double *parm) {
2503 float sig, cut, t;
2504 size_t i;
2505
2506 sig = AIR_CAST(float, parm[0])((float)(parm[0]));
2507 cut = AIR_CAST(float, parm[1])((float)(parm[1]));
2508 for (i=0; i<len; i++) {
2509 t = x[i];
2510 t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t));
2511 f[i] = AIR_CAST(float, _GAUSS(t, sig, cut))((float)(( t >= sig*cut ? 0 : exp(-t*t/(2.0*sig*sig))/(sig
*2.50662827463100050241))));
2512 }
2513}
2514
2515static NrrdKernel
2516_nrrdKernelG = {
2517 "gauss",
2518 2, _nrrdGSup, _nrrdGInt,
2519 _nrrdG1_f, _nrrdGN_f, _nrrdG1_d, _nrrdGN_d
2520};
2521NrrdKernel *const
2522nrrdKernelGaussian = &_nrrdKernelG;
2523
2524/* ------------------------------------------------------------ */
2525
2526#define _DISCRETEGAUSS(xx, sig, abscut)(sig > 0 ? (xx > abscut ? 0 : airBesselInExpScaled(((int
)(xx + 0.5)), sig*sig)) : xx <= 0.5) \
2527 (sig > 0 \
2528 ? (xx > abscut \
2529 ? 0 \
2530 : airBesselInExpScaled(AIR_CAST(int, xx + 0.5)((int)(xx + 0.5)), sig*sig)) \
2531 : xx <= 0.5)
2532
2533/* the last line used to be AIR_MAX(0.5, (ret)), but the problem was
2534 that even for reasonable cutoffs, like sigma=6, there would be a
2535 sudden change in the non-zero values at the edge of the kernel with
2536 slowly increasing sigma. The real solution is to have a smarter way
2537 of determining where to cut-off this particular kernel, the
2538 discrete gaussian, when the values are at the low level one would
2539 expect with "sigma=6" when talking about a continuous Gaussian */
2540#define _DGABSCUT(ret, sig, cut)(ret) = 0.5 + ceil((sig)*(cut)); (ret) = ((2.5) > ((ret)) ?
(2.5) : ((ret))) \
2541 (ret) = 0.5 + ceil((sig)*(cut)); \
2542 (ret) = AIR_MAX(2.5, (ret))((2.5) > ((ret)) ? (2.5) : ((ret)))
2543
2544static double
2545_nrrdDiscGaussianSup(const double *parm) {
2546 double ret;
2547
2548 _DGABSCUT(ret, parm[0], parm[1])(ret) = 0.5 + ceil((parm[0])*(parm[1])); (ret) = ((2.5) > (
(ret)) ? (2.5) : ((ret)));
2549 return ret;
2550}
2551
2552/* the functions are out of their usual definition order because we
2553 use the function evaluation to determine the integral, rather than
2554 trying to do it analytically */
2555
2556static double
2557_nrrdDiscGaussian1_d(double xx, const double *parm) {
2558 double sig, cut;
2559
2560 sig = parm[0];
2561 _DGABSCUT(cut, sig, parm[1])(cut) = 0.5 + ceil((sig)*(parm[1])); (cut) = ((2.5) > ((cut
)) ? (2.5) : ((cut)));
2562 xx = AIR_ABS(xx)((xx) > 0.0f ? (xx) : -(xx));
2563 return _DISCRETEGAUSS(xx, sig, cut)(sig > 0 ? (xx > cut ? 0 : airBesselInExpScaled(((int)(
xx + 0.5)), sig*sig)) : xx <= 0.5);
2564}
2565
2566static double
2567_nrrdDiscGaussianInt(const double *parm) {
2568 double sum, cut;
2569 int ii, supp;
2570
2571 _DGABSCUT(cut, parm[0], parm[1])(cut) = 0.5 + ceil((parm[0])*(parm[1])); (cut) = ((2.5) > (
(cut)) ? (2.5) : ((cut)));
2572 supp = AIR_CAST(int, cut)((int)(cut));
2573 sum = 0.0;
2574 for (ii=-supp; ii<=supp; ii++) {
2575 sum += _nrrdDiscGaussian1_d(ii, parm);
2576 }
2577 return sum;
2578}
2579
2580static float
2581_nrrdDiscGaussian1_f(float xx, const double *parm) {
2582 double sig, cut;
2583
2584 sig = parm[0];
2585 _DGABSCUT(cut, sig, parm[1])(cut) = 0.5 + ceil((sig)*(parm[1])); (cut) = ((2.5) > ((cut
)) ? (2.5) : ((cut)));
2586 xx = AIR_ABS(xx)((xx) > 0.0f ? (xx) : -(xx));
2587 return AIR_CAST(float, _DISCRETEGAUSS(xx, sig, cut))((float)((sig > 0 ? (xx > cut ? 0 : airBesselInExpScaled
(((int)(xx + 0.5)), sig*sig)) : xx <= 0.5)));
2588}
2589
2590static void
2591_nrrdDiscGaussianN_d(double *f, const double *x,
2592 size_t len, const double *parm) {
2593 double sig, cut, tt;
2594 size_t ii;
2595
2596 sig = parm[0];
2597 _DGABSCUT(cut, sig, parm[1])(cut) = 0.5 + ceil((sig)*(parm[1])); (cut) = ((2.5) > ((cut
)) ? (2.5) : ((cut)));
2598 for (ii=0; ii<len; ii++) {
2599 tt = AIR_ABS(x[ii])((x[ii]) > 0.0f ? (x[ii]) : -(x[ii]));
2600 f[ii] = _DISCRETEGAUSS(tt, sig, cut)(sig > 0 ? (tt > cut ? 0 : airBesselInExpScaled(((int)(
tt + 0.5)), sig*sig)) : tt <= 0.5);
2601 }
2602}
2603
2604static void
2605_nrrdDiscGaussianN_f(float *f, const float *x, size_t len, const double *parm) {
2606 double sig, cut, tt;
2607 size_t ii;
2608
2609 sig = parm[0];
2610 _DGABSCUT(cut, sig, parm[1])(cut) = 0.5 + ceil((sig)*(parm[1])); (cut) = ((2.5) > ((cut
)) ? (2.5) : ((cut)));
2611 for (ii=0; ii<len; ii++) {
2612 tt = AIR_ABS(x[ii])((x[ii]) > 0.0f ? (x[ii]) : -(x[ii]));
2613 f[ii] = AIR_CAST(float, _DISCRETEGAUSS(tt, sig, cut))((float)((sig > 0 ? (tt > cut ? 0 : airBesselInExpScaled
(((int)(tt + 0.5)), sig*sig)) : tt <= 0.5)));
2614 }
2615}
2616
2617static NrrdKernel
2618_nrrdKernelDiscreteGaussian = {
2619 "discretegauss", 2,
2620 _nrrdDiscGaussianSup, _nrrdDiscGaussianInt,
2621 _nrrdDiscGaussian1_f, _nrrdDiscGaussianN_f,
2622 _nrrdDiscGaussian1_d, _nrrdDiscGaussianN_d
2623};
2624NrrdKernel *const
2625nrrdKernelDiscreteGaussian = &_nrrdKernelDiscreteGaussian;
2626
2627/* The current implementation of nrrdKernelDiscreteGaussian, with the current
2628 implementation of airBesselInExpScaled, has problems for large
2629 sigma. Until those are fixed, this is a suggested limit on how big sigma
2630 (kparm[0]) can be and still have an accurate blurring kernel. This
2631 problem can be avoided completely by doing the blurring in frequency
2632 space, which is implemented in teem/src/gage/stackBlur.c's
2633 _stackBlurDiscreteGaussFFT(), although that has its own problem: in places
2634 where a signal really should zero, the FFT can produce some very
2635 low-amplitude noise (and hence new extrema) */
2636const double nrrdKernelDiscreteGaussianGoodSigmaMax = 6.0;
2637
2638/* ------------------------------------------------------------ */
2639
2640#define _DGAUSS(x, sig, cut)( x >= sig*cut ? 0 : -exp(-x*x/(2.0*sig*sig))*x/(sig*sig*sig
*2.50662827463100050241)) ( \
2641 x >= sig*cut ? 0 \
2642 : -exp(-x*x/(2.0*sig*sig))*x/(sig*sig*sig*2.50662827463100050241))
2643
2644static double
2645_nrrdDGInt(const double *parm) {
2646 AIR_UNUSED(parm)(void)(parm);
2647 return 0;
2648}
2649
2650static double
2651_nrrdDGSup(const double *parm) {
2652 double sig, cut;
2653
2654 sig = parm[0];
2655 cut = parm[1];
2656 return sig*cut;
2657}
2658
2659static double
2660_nrrdDG1_d(double x, const double *parm) {
2661 double sig, cut;
2662 int sgn = 1;
2663
2664 sig = parm[0];
2665 cut = parm[1];
2666 if (x < 0) { x = -x; sgn = -1; }
2667 return sgn*_DGAUSS(x, sig, cut)( x >= sig*cut ? 0 : -exp(-x*x/(2.0*sig*sig))*x/(sig*sig*sig
*2.50662827463100050241));
2668}
2669
2670static float
2671_nrrdDG1_f(float x, const double *parm) {
2672 float sig, cut;
2673 int sgn = 1;
2674
2675 sig = AIR_CAST(float, parm[0])((float)(parm[0]));
2676 cut = AIR_CAST(float, parm[1])((float)(parm[1]));
2677 if (x < 0) { x = -x; sgn = -1; }
2678 return AIR_CAST(float, sgn*_DGAUSS(x, sig, cut))((float)(sgn*( x >= sig*cut ? 0 : -exp(-x*x/(2.0*sig*sig))
*x/(sig*sig*sig*2.50662827463100050241))));
2679}
2680
2681static void
2682_nrrdDGN_d(double *f, const double *x, size_t len, const double *parm) {
2683 double sig, cut, t;
2684 size_t i;
2685 int sgn;
2686
2687 sig = parm[0];
2688 cut = parm[1];
2689 for (i=0; i<len; i++) {
2690 t = x[i];
2691 if (t < 0) { t = -t; sgn = -1; } else { sgn = 1; }
2692 f[i] = sgn*_DGAUSS(t, sig, cut)( t >= sig*cut ? 0 : -exp(-t*t/(2.0*sig*sig))*t/(sig*sig*sig
*2.50662827463100050241));
2693 }
2694}
2695
2696static void
2697_nrrdDGN_f(float *f, const float *x, size_t len, const double *parm) {
2698 float sig, cut, t;
2699 size_t i;
2700 int sgn;
2701
2702 sig = AIR_CAST(float, parm[0])((float)(parm[0]));
2703 cut = AIR_CAST(float, parm[1])((float)(parm[1]));
2704 for (i=0; i<len; i++) {
2705 t = x[i];
2706 if (t < 0) { t = -t; sgn = -1; } else { sgn = 1; }
2707 f[i] = AIR_CAST(float, sgn*_DGAUSS(t, sig, cut))((float)(sgn*( t >= sig*cut ? 0 : -exp(-t*t/(2.0*sig*sig))
*t/(sig*sig*sig*2.50662827463100050241))));
2708 }
2709}
2710
2711static NrrdKernel
2712_nrrdKernelDG = {
2713 "gaussD",
2714 2, _nrrdDGSup, _nrrdDGInt,
2715 _nrrdDG1_f, _nrrdDGN_f, _nrrdDG1_d, _nrrdDGN_d
2716};
2717NrrdKernel *const
2718nrrdKernelGaussianD = &_nrrdKernelDG;
2719
2720/* ------------------------------------------------------------ */
2721
2722#define _DDGAUSS(x, sig, cut)( x >= sig*cut ? 0 : exp(-x*x/(2.0*sig*sig))*(x*x-sig*sig)
/ (sig*sig*sig*sig*sig*2.50662827463100050241)) ( \
2723 x >= sig*cut ? 0 \
2724 : exp(-x*x/(2.0*sig*sig))*(x*x-sig*sig) / \
2725 (sig*sig*sig*sig*sig*2.50662827463100050241))
2726
2727static double
2728_nrrdDDGInt(const double *parm) {
2729 double sig, cut;
2730
2731 sig = parm[0];
2732 cut = parm[1];
2733 return -0.79788456080286535587*cut*exp(-cut*cut/2)/(sig*sig);
2734}
2735
2736static double
2737_nrrdDDGSup(const double *parm) {
2738 double sig, cut;
2739
2740 sig = parm[0];
2741 cut = parm[1];
2742 return sig*cut;
2743}
2744
2745static double
2746_nrrdDDG1_d(double x, const double *parm) {
2747 double sig, cut;
2748
2749 sig = parm[0];
2750 cut = parm[1];
2751 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x));
2752 return _DDGAUSS(x, sig, cut)( x >= sig*cut ? 0 : exp(-x*x/(2.0*sig*sig))*(x*x-sig*sig)
/ (sig*sig*sig*sig*sig*2.50662827463100050241));
2753}
2754
2755static float
2756_nrrdDDG1_f(float x, const double *parm) {
2757 float sig, cut;
2758
2759 sig = AIR_CAST(float, parm[0])((float)(parm[0]));
2760 cut = AIR_CAST(float, parm[1])((float)(parm[1]));
2761 x = AIR_ABS(x)((x) > 0.0f ? (x) : -(x));
2762 return AIR_CAST(float, _DDGAUSS(x, sig, cut))((float)(( x >= sig*cut ? 0 : exp(-x*x/(2.0*sig*sig))*(x*x
-sig*sig) / (sig*sig*sig*sig*sig*2.50662827463100050241))));
2763}
2764
2765static void
2766_nrrdDDGN_d(double *f, const double *x, size_t len, const double *parm) {
2767 double sig, cut, t;
2768 size_t i;
2769
2770 sig = parm[0];
2771 cut = parm[1];
2772 for (i=0; i<len; i++) {
2773 t = x[i];
2774 t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t));
2775 f[i] = _DDGAUSS(t, sig, cut)( t >= sig*cut ? 0 : exp(-t*t/(2.0*sig*sig))*(t*t-sig*sig)
/ (sig*sig*sig*sig*sig*2.50662827463100050241));
2776 }
2777}
2778
2779static void
2780_nrrdDDGN_f(float *f, const float *x, size_t len, const double *parm) {
2781 float sig, cut, t;
2782 size_t i;
2783
2784 sig = AIR_CAST(float, parm[0])((float)(parm[0]));
2785 cut = AIR_CAST(float, parm[1])((float)(parm[1]));
2786 for (i=0; i<len; i++) {
2787 t = x[i];
2788 t = AIR_ABS(t)((t) > 0.0f ? (t) : -(t));
2789 f[i] = AIR_CAST(float, _DDGAUSS(t, sig, cut))((float)(( t >= sig*cut ? 0 : exp(-t*t/(2.0*sig*sig))*(t*t
-sig*sig) / (sig*sig*sig*sig*sig*2.50662827463100050241))));
2790 }
2791}
2792
2793static NrrdKernel
2794_nrrdKernelDDG = {
2795 "gaussDD",
2796 2, _nrrdDDGSup, _nrrdDDGInt,
2797 _nrrdDDG1_f, _nrrdDDGN_f, _nrrdDDG1_d, _nrrdDDGN_d
2798};
2799NrrdKernel *const
2800nrrdKernelGaussianDD = &_nrrdKernelDDG;
2801
2802
2803/* ------------------------------------------------------------ */
2804
2805NrrdKernel *
2806_nrrdKernelStrToKern(char *str) {
2807
2808 if (!strcmp("zero", str)) return nrrdKernelZero;
2809 if (!strcmp("box", str)) return nrrdKernelBox;
2810 if (!strcmp("boxsup", str)) return nrrdKernelBoxSupportDebug;
2811 if (!strcmp("cos4sup", str)) return nrrdKernelCos4SupportDebug;
2812 if (!strcmp("cos4supd", str)) return nrrdKernelCos4SupportDebugD;
2813 if (!strcmp("cos4supdd", str)) return nrrdKernelCos4SupportDebugDD;
2814 if (!strcmp("cos4supddd", str)) return nrrdKernelCos4SupportDebugDDD;
2815 if (!strcmp("cheap", str)) return nrrdKernelCheap;
2816 if (!strcmp("hermiteflag", str)) return nrrdKernelHermiteScaleSpaceFlag;
2817 if (!strcmp("hermite", str)) return nrrdKernelHermiteScaleSpaceFlag;
2818 if (!strcmp("hermitess", str)) return nrrdKernelHermiteScaleSpaceFlag;
2819 if (!strcmp("herm", str)) return nrrdKernelHermiteScaleSpaceFlag;
2820 if (!strcmp("tent", str)) return nrrdKernelTent;
2821 if (!strcmp("tentd", str)) return nrrdKernelForwDiff;
2822 if (!strcmp("forwdiff", str)) return nrrdKernelForwDiff;
2823 if (!strcmp("fordif", str)) return nrrdKernelForwDiff;
2824 if (!strcmp("centdiff", str)) return nrrdKernelCentDiff;
2825 if (!strcmp("cendif", str)) return nrrdKernelCentDiff;
2826 if (!strcmp("bccubic", str)) return nrrdKernelBCCubic;
2827 if (!strcmp("cubic", str)) return nrrdKernelBCCubic;
2828 if (!strcmp("bccubicd", str)) return nrrdKernelBCCubicD;
2829 if (!strcmp("cubicd", str)) return nrrdKernelBCCubicD;
2830 if (!strcmp("bccubicdd", str)) return nrrdKernelBCCubicDD;
2831 if (!strcmp("cubicdd", str)) return nrrdKernelBCCubicDD;
2832 if (!strcmp("ctmr", str)) return nrrdKernelCatmullRom;
2833 if (!strcmp("catmull-rom", str)) return nrrdKernelCatmullRom;
2834 if (!strcmp("ctmrd", str)) return nrrdKernelCatmullRomD;
2835 if (!strcmp("catmull-romd", str)) return nrrdKernelCatmullRomD;
2836 if (!strcmp("ctmrdd", str)) return nrrdKernelCatmullRomDD;
2837 if (!strcmp("catmull-romdd", str)) return nrrdKernelCatmullRomDD;
2838 if (!strcmp("ctmrsup", str)) return nrrdKernelCatmullRomSupportDebug;
2839 if (!strcmp("ctmrsupd", str)) return nrrdKernelCatmullRomSupportDebugD;
2840 if (!strcmp("ctmrsupdd", str)) return nrrdKernelCatmullRomSupportDebugDD;
2841 if (!strcmp("aquartic", str)) return nrrdKernelAQuartic;
2842 if (!strcmp("quartic", str)) return nrrdKernelAQuartic;
2843 if (!strcmp("aquarticd", str)) return nrrdKernelAQuarticD;
2844 if (!strcmp("quarticd", str)) return nrrdKernelAQuarticD;
2845 if (!strcmp("aquarticdd", str)) return nrrdKernelAQuarticDD;
2846 if (!strcmp("quarticdd", str)) return nrrdKernelAQuarticDD;
2847 if (!strcmp("c3quintic", str)) return nrrdKernelC3Quintic;
2848 if (!strcmp("c3q", str)) return nrrdKernelC3Quintic;
2849 if (!strcmp("c3quinticd", str)) return nrrdKernelC3QuinticD;
2850 if (!strcmp("c3qd", str)) return nrrdKernelC3QuinticD;
2851 if (!strcmp("c3quinticdd", str)) return nrrdKernelC3QuinticDD;
2852 if (!strcmp("c3qdd", str)) return nrrdKernelC3QuinticDD;
2853 if (!strcmp("c4hexic", str)) return nrrdKernelC4Hexic;
2854 if (!strcmp("c4hai", str)) return nrrdKernelC4HexicApproxInverse;
2855 if (!strcmp("c4hexicai", str)) return nrrdKernelC4HexicApproxInverse;
2856 if (!strcmp("c4h", str)) return nrrdKernelC4Hexic;
2857 if (!strcmp("c4hexicd", str)) return nrrdKernelC4HexicD;
2858 if (!strcmp("c4hd", str)) return nrrdKernelC4HexicD;
2859 if (!strcmp("c4hexicdd", str)) return nrrdKernelC4HexicDD;
2860 if (!strcmp("c4hdd", str)) return nrrdKernelC4HexicDD;
2861 if (!strcmp("c4hexicddd", str)) return nrrdKernelC4HexicDDD;
2862 if (!strcmp("c4hddd", str)) return nrrdKernelC4HexicDDD;
2863 if (!strcmp("c5septic", str)) return nrrdKernelC5Septic;
2864 if (!strcmp("c5septicd", str)) return nrrdKernelC5SepticD;
2865 if (!strcmp("c5septicdd", str)) return nrrdKernelC5SepticDD;
2866 if (!strcmp("c5septicddd", str))return nrrdKernelC5SepticDDD;
2867 if (!strcmp("c5septicai", str)) return nrrdKernelC5SepticApproxInverse;
2868 if (!strcmp("gaussian", str)) return nrrdKernelGaussian;
2869 if (!strcmp("gauss", str)) return nrrdKernelGaussian;
2870 if (!strcmp("gaussiand", str)) return nrrdKernelGaussianD;
2871 if (!strcmp("gaussd", str)) return nrrdKernelGaussianD;
2872 if (!strcmp("gd", str)) return nrrdKernelGaussianD;
2873 if (!strcmp("gaussiandd", str)) return nrrdKernelGaussianDD;
2874 if (!strcmp("gaussdd", str)) return nrrdKernelGaussianDD;
2875 if (!strcmp("gdd", str)) return nrrdKernelGaussianDD;
2876 if (!strcmp("ds", str)) return nrrdKernelDiscreteGaussian;
2877 if (!strcmp("dscale", str)) return nrrdKernelDiscreteGaussian;
2878 if (!strcmp("dg", str)) return nrrdKernelDiscreteGaussian;
2879 if (!strcmp("dgauss", str)) return nrrdKernelDiscreteGaussian;
2880 if (!strcmp("dgaussian", str)) return nrrdKernelDiscreteGaussian;
2881 if (!strcmp("discretegauss", str)) return nrrdKernelDiscreteGaussian;
2882 if (!strcmp("hann", str)) return nrrdKernelHann;
2883 if (!strcmp("hannd", str)) return nrrdKernelHannD;
2884 if (!strcmp("hanndd", str)) return nrrdKernelHannDD;
2885 if (!strcmp("bkmn", str)) return nrrdKernelBlackman;
2886 if (!strcmp("black", str)) return nrrdKernelBlackman;
2887 if (!strcmp("blackman", str)) return nrrdKernelBlackman;
2888 if (!strcmp("bkmnd", str)) return nrrdKernelBlackmanD;
2889 if (!strcmp("blackd", str)) return nrrdKernelBlackmanD;
2890 if (!strcmp("blackmand", str)) return nrrdKernelBlackmanD;
2891 if (!strcmp("bkmndd", str)) return nrrdKernelBlackmanDD;
2892 if (!strcmp("blackdd", str)) return nrrdKernelBlackmanDD;
2893 if (!strcmp("blackmandd", str)) return nrrdKernelBlackmanDD;
2894 if (!strcmp("bspl1", str)) return nrrdKernelBSpline1;
2895 if (!strcmp("bspln1", str)) return nrrdKernelBSpline1;
2896 if (!strcmp("bspl1d", str)) return nrrdKernelBSpline1D;
2897 if (!strcmp("bspln1d", str)) return nrrdKernelBSpline1D;
2898 if (!strcmp("bspl2", str)) return nrrdKernelBSpline2;
2899 if (!strcmp("bspln2", str)) return nrrdKernelBSpline2;
2900 if (!strcmp("bspl2d", str)) return nrrdKernelBSpline2D;
2901 if (!strcmp("bspln2d", str)) return nrrdKernelBSpline2D;
2902 if (!strcmp("bspl2dd", str)) return nrrdKernelBSpline2DD;
2903 if (!strcmp("bspln2dd", str)) return nrrdKernelBSpline2DD;
2904 if (!strcmp("bspl3", str)) return nrrdKernelBSpline3;
2905 if (!strcmp("bspln3", str)) return nrrdKernelBSpline3;
2906 if (!strcmp("bspl3ai", str)) return nrrdKernelBSpline3ApproxInverse;
2907 if (!strcmp("bspln3ai", str)) return nrrdKernelBSpline3ApproxInverse;
2908 if (!strcmp("bspl3d", str)) return nrrdKernelBSpline3D;
2909 if (!strcmp("bspln3d", str)) return nrrdKernelBSpline3D;
2910 if (!strcmp("bspl3dd", str)) return nrrdKernelBSpline3DD;
2911 if (!strcmp("bspln3dd", str)) return nrrdKernelBSpline3DD;
2912 if (!strcmp("bspl3ddd", str)) return nrrdKernelBSpline3DDD;
2913 if (!strcmp("bspln3ddd", str)) return nrrdKernelBSpline3DDD;
2914 if (!strcmp("bspl4", str)) return nrrdKernelBSpline4;
2915 if (!strcmp("bspln4", str)) return nrrdKernelBSpline4;
2916 if (!strcmp("bspl4d", str)) return nrrdKernelBSpline4D;
2917 if (!strcmp("bspln4d", str)) return nrrdKernelBSpline4D;
2918 if (!strcmp("bspl4dd", str)) return nrrdKernelBSpline4DD;
2919 if (!strcmp("bspln4dd", str)) return nrrdKernelBSpline4DD;
2920 if (!strcmp("bspl4ddd", str)) return nrrdKernelBSpline4DDD;
2921 if (!strcmp("bspln4ddd", str)) return nrrdKernelBSpline4DDD;
2922 if (!strcmp("bspl5", str)) return nrrdKernelBSpline5;
2923 if (!strcmp("bspln5", str)) return nrrdKernelBSpline5;
2924 if (!strcmp("bspl5ai", str)) return nrrdKernelBSpline5ApproxInverse;
2925 if (!strcmp("bspln5ai", str)) return nrrdKernelBSpline5ApproxInverse;
2926 if (!strcmp("bspl5d", str)) return nrrdKernelBSpline5D;
2927 if (!strcmp("bspln5d", str)) return nrrdKernelBSpline5D;
2928 if (!strcmp("bspl5dd", str)) return nrrdKernelBSpline5DD;
2929 if (!strcmp("bspln5dd", str)) return nrrdKernelBSpline5DD;
2930 if (!strcmp("bspl5ddd", str)) return nrrdKernelBSpline5DDD;
2931 if (!strcmp("bspln5ddd", str)) return nrrdKernelBSpline5DDD;
2932 if (!strcmp("bspl6", str)) return nrrdKernelBSpline6;
2933 if (!strcmp("bspln6", str)) return nrrdKernelBSpline6;
2934 if (!strcmp("bspl6d", str)) return nrrdKernelBSpline6D;
2935 if (!strcmp("bspln6d", str)) return nrrdKernelBSpline6D;
2936 if (!strcmp("bspl6dd", str)) return nrrdKernelBSpline6DD;
2937 if (!strcmp("bspln6dd", str)) return nrrdKernelBSpline6DD;
2938 if (!strcmp("bspl6ddd", str)) return nrrdKernelBSpline6DDD;
2939 if (!strcmp("bspln6ddd", str)) return nrrdKernelBSpline6DDD;
2940 if (!strcmp("bspl7", str)) return nrrdKernelBSpline7;
2941 if (!strcmp("bspln7", str)) return nrrdKernelBSpline7;
2942 if (!strcmp("bspl7ai", str)) return nrrdKernelBSpline7ApproxInverse;
2943 if (!strcmp("bspln7ai", str)) return nrrdKernelBSpline7ApproxInverse;
2944 if (!strcmp("bspl7d", str)) return nrrdKernelBSpline7D;
2945 if (!strcmp("bspln7d", str)) return nrrdKernelBSpline7D;
2946 if (!strcmp("bspl7dd", str)) return nrrdKernelBSpline7DD;
2947 if (!strcmp("bspln7dd", str)) return nrrdKernelBSpline7DD;
2948 if (!strcmp("bspl7ddd", str)) return nrrdKernelBSpline7DDD;
2949 if (!strcmp("bspln7ddd", str)) return nrrdKernelBSpline7DDD;
2950 return NULL((void*)0);
2951}
2952
2953/* this returns a number between -1 and max;
2954 it does NOT do the increment-by-one;
2955 it does NOT do range checking */
2956int
2957_nrrdKernelParseTMFInt(int *val, char *str) {
2958 static const char me[]="nrrdKernelParseTMFInt";
2959
2960 if (!strcmp("n", str)) {
2961 *val = -1;
2962 } else {
2963 if (1 != sscanf(str, "%d", val)) {
2964 biffAddf(NRRDnrrdBiffKey, "%s: couldn't parse \"%s\" as int", me, str);
2965 return 1;
2966 }
2967 }
2968 return 0;
2969}
2970
2971int
2972nrrdKernelParse(const NrrdKernel **kernelP,
2973 double *parm, const char *_str) {
2974 static const char me[]="nrrdKernelParse";
2975 char str[AIR_STRLEN_HUGE(1024+1)],
2976 kstr[AIR_STRLEN_MED(256+1)], *_pstr=NULL((void*)0), *pstr,
2977 *tmfStr[4] = {NULL((void*)0), NULL((void*)0), NULL((void*)0), NULL((void*)0)};
2978 int tmfD, tmfC, tmfA;
2979 unsigned int jj, haveParm, needParm;
2980 airArray *mop;
2981
2982 if (!(kernelP && parm && _str)) {
2983 biffAddf(NRRDnrrdBiffKey, "%s: got NULL pointer", me);
2984 return 1;
2985 }
2986
2987 /* [jorik] (if i understood this correctly) parm is always of length
2988 ** NRRD_KERNEL_PARMS_NUM. We have to clear all parameters here, since
2989 ** nrrdKernelSet copies all arguments into its own array later, and
2990 ** copying uninitialised memory is bad (it traps my memory debugger).
2991 */
2992 for (jj=0; jj<NRRD_KERNEL_PARMS_NUM8; jj++) {
2993 parm[jj] = 0;
2994 }
2995
2996 airStrcpy(str, AIR_STRLEN_HUGE(1024+1), _str);
2997 strcpy(kstr, "")__builtin___strcpy_chk (kstr, "", __builtin_object_size (kstr
, 2 > 1 ? 1 : 0));
2998 pstr = NULL((void*)0);
2999 pstr = strchr(str, ':');
3000 if (pstr) {
3001 *pstr = '\0';
3002 _pstr = ++pstr;
3003 }
3004 strcpy(kstr, str)__builtin___strcpy_chk (kstr, str, __builtin_object_size (kstr
, 2 > 1 ? 1 : 0));
3005 airToLower(kstr);
3006 mop = airMopNew();
3007 /* first see if its a TMF, then try parsing it as the other stuff */
3008 if (kstr == strstr(kstr, "tmf")) {
3009 if (4 == airParseStrS(tmfStr, pstr, ",", 4)) {
3010 airMopAdd(mop, tmfStr[0], airFree, airMopAlways);
3011 airMopAdd(mop, tmfStr[1], airFree, airMopAlways);
3012 airMopAdd(mop, tmfStr[2], airFree, airMopAlways);
3013 airMopAdd(mop, tmfStr[3], airFree, airMopAlways);
3014 /* a TMF with a parameter: D,C,A,a */
3015 if (1 != airSingleSscanf(tmfStr[3], "%lg", parm)) {
3016 biffAddf(NRRDnrrdBiffKey, "%s: couldn't parse TMF parameter \"%s\" as double",
3017 me, tmfStr[3]);
3018 airMopError(mop); return 1;
3019 }
3020 } else if (3 == airParseStrS(tmfStr, pstr, ",", 3)) {
3021 airMopAdd(mop, tmfStr[0], airFree, airMopAlways);
3022 airMopAdd(mop, tmfStr[1], airFree, airMopAlways);
3023 airMopAdd(mop, tmfStr[2], airFree, airMopAlways);
3024 /* a TMF without a parameter: D,C,A ==> a=0.0 */
3025 parm[0] = 0.0;
3026 } else {
3027 biffAddf(NRRDnrrdBiffKey, "%s: TMF kernels require 3 arguments D, C, A "
3028 "in the form tmf:D,C,A", me);
3029 airMopError(mop); return 1;
3030 }
3031 if (_nrrdKernelParseTMFInt(&tmfD, tmfStr[0])
3032 || _nrrdKernelParseTMFInt(&tmfC, tmfStr[1])
3033 || _nrrdKernelParseTMFInt(&tmfA, tmfStr[2])) {
3034 biffAddf(NRRDnrrdBiffKey, "%s: problem parsing \"%s,%s,%s\" as D,C,A "
3035 "for TMF kernel", me, tmfStr[0], tmfStr[1], tmfStr[2]);
3036 airMopError(mop); return 1;
3037 }
3038 if (!AIR_IN_CL(-1, tmfD, (int)nrrdKernelTMF_maxD)((-1) <= (tmfD) && (tmfD) <= ((int)nrrdKernelTMF_maxD
))) {
3039 biffAddf(NRRDnrrdBiffKey, "%s: derivative value %d outside range [-1,%d]",
3040 me, tmfD, nrrdKernelTMF_maxD);
3041 airMopError(mop); return 1;
3042 }
3043 if (!AIR_IN_CL(-1, tmfC, (int)nrrdKernelTMF_maxC)((-1) <= (tmfC) && (tmfC) <= ((int)nrrdKernelTMF_maxC
))) {
3044 biffAddf(NRRDnrrdBiffKey, "%s: continuity value %d outside range [-1,%d]",
3045 me, tmfC, nrrdKernelTMF_maxC);
3046 airMopError(mop); return 1;
3047 }
3048 if (!AIR_IN_CL(1, tmfA, (int)nrrdKernelTMF_maxA)((1) <= (tmfA) && (tmfA) <= ((int)nrrdKernelTMF_maxA
))) {
3049 biffAddf(NRRDnrrdBiffKey, "%s: accuracy value %d outside range [1,%d]",
3050 me, tmfA, nrrdKernelTMF_maxA);
3051 airMopError(mop); return 1;
3052 }
3053 /*
3054 fprintf(stderr, "!%s: D,C,A = %d,%d,%d --> %d,%d,%d\n", me,
3055 tmfD, tmfC, tmfA, tmfD+1, tmfC+1, tmfA);
3056 */
3057 *kernelP = nrrdKernelTMF[tmfD+1][tmfC+1][tmfA];
3058 } else {
3059 /* its not a TMF */
3060 if (!(*kernelP = _nrrdKernelStrToKern(kstr))) {
3061 biffAddf(NRRDnrrdBiffKey, "%s: kernel \"%s\" not recognized", me, kstr);
3062 airMopError(mop); return 1;
3063 }
3064 if ((*kernelP)->numParm > NRRD_KERNEL_PARMS_NUM8) {
3065 biffAddf(NRRDnrrdBiffKey, "%s: kernel \"%s\" requests %d parameters > max %d",
3066 me, kstr, (*kernelP)->numParm, NRRD_KERNEL_PARMS_NUM8);
3067 airMopError(mop); return 1;
3068 }
3069 if (*kernelP == nrrdKernelGaussian ||
3070 *kernelP == nrrdKernelGaussianD ||
3071 *kernelP == nrrdKernelGaussianDD ||
3072 *kernelP == nrrdKernelDiscreteGaussian ||
3073 *kernelP == nrrdKernelBoxSupportDebug ||
3074 *kernelP == nrrdKernelCos4SupportDebug ||
3075 *kernelP == nrrdKernelCos4SupportDebugD ||
3076 *kernelP == nrrdKernelCos4SupportDebugDD ||
3077 *kernelP == nrrdKernelCos4SupportDebugDDD) {
3078 /* for these kernels, we need all the parameters given explicitly */
3079 needParm = (*kernelP)->numParm;
3080 } else {
3081 /* For everything else (note that TMF kernels are handled
3082 separately), we can make do with one less than the required,
3083 by using the default spacing */
3084 needParm = ((*kernelP)->numParm > 0
3085 ? (*kernelP)->numParm - 1
3086 : 0);
3087 }
3088 if (needParm > 0 && !pstr) {
3089 biffAddf(NRRDnrrdBiffKey, "%s: didn't get any of %d required doubles after "
3090 "colon in \"%s\"",
3091 me, needParm, kstr);
3092 airMopError(mop); return 1;
3093 }
3094 for (haveParm=0; haveParm<(*kernelP)->numParm; haveParm++) {
3095 if (!pstr)
3096 break;
3097 if (1 != airSingleSscanf(pstr, "%lg", parm+haveParm)) {
3098 biffAddf(NRRDnrrdBiffKey, "%s: trouble parsing \"%s\" as double (in \"%s\")",
3099 me, _pstr, _str);
3100 airMopError(mop); return 1;
3101 }
3102 if ((pstr = strchr(pstr, ','))) {
3103 pstr++;
3104 if (!*pstr) {
3105 biffAddf(NRRDnrrdBiffKey, "%s: nothing after last comma in \"%s\" (in \"%s\")",
3106 me, _pstr, _str);
3107 airMopError(mop); return 1;
3108 }
3109 }
3110 }
3111 /* haveParm is now the number of parameters that were parsed. */
3112 if (haveParm < needParm) {
3113 biffAddf(NRRDnrrdBiffKey, "%s: parsed only %d of %d required doubles "
3114 "from \"%s\" (in \"%s\")",
3115 me, haveParm, needParm, _pstr, _str);
3116 airMopError(mop); return 1;
3117 } else if (haveParm == needParm &&
3118 needParm == (*kernelP)->numParm-1) {
3119 /* shift up parsed values, and set parm[0] to default */
3120 for (jj=haveParm; jj>=1; jj--) {
3121 parm[jj] = parm[jj-1];
3122 }
3123 parm[0] = nrrdDefaultKernelParm0;
3124 } else {
3125 if (pstr) {
3126 biffAddf(NRRDnrrdBiffKey, "%s: \"%s\" (in \"%s\") has more than %d doubles",
3127 me, _pstr, _str, (*kernelP)->numParm);
3128 airMopError(mop); return 1;
3129 }
3130 }
3131 }
3132 /*
3133 fprintf(stderr, "%s: %g %g %g %g %g\n", me,
3134 parm[0], parm[1], parm[2], parm[3], parm[4]);
3135 */
3136 airMopOkay(mop);
3137 return 0;
3138}
3139
3140int
3141nrrdKernelSpecParse(NrrdKernelSpec *ksp, const char *str) {
3142 static const char me[]="nrrdKernelSpecParse";
3143 const NrrdKernel *kern;
3144 double kparm[NRRD_KERNEL_PARMS_NUM8];
3145
3146 if (!( ksp && str )) {
3147 biffAddf(NRRDnrrdBiffKey, "%s: got NULL pointer", me);
3148 return 1;
3149 }
3150 if (nrrdKernelParse(&kern, kparm, str)) {
3151 biffAddf(NRRDnrrdBiffKey, "%s: ", me);
3152 return 1;
3153 }
3154 nrrdKernelSpecSet(ksp, kern, kparm);
3155 return 0;
3156}
3157
3158/*
3159** note that the given string has to be allocated for a certain size
3160** which is plenty big
3161*/
3162int
3163nrrdKernelSpecSprint(char str[AIR_STRLEN_LARGE(512+1)], const NrrdKernelSpec *ksp) {
3164 static const char me[]="nrrdKernelSpecSprint";
3165 unsigned int warnLen = AIR_STRLEN_LARGE(512+1)/3;
3166 char stmp[AIR_STRLEN_LARGE(512+1)];
3167
3168 if (!( str && ksp )) {
3169 biffAddf(NRRDnrrdBiffKey, "%s: got NULL pointer", me);
3170 return 1;
3171 }
3172 if (strlen(ksp->kernel->name) > warnLen) {
3173 biffAddf(NRRDnrrdBiffKey, "%s: kernel name (len %s) might lead to overflow", me,
3174 airSprintSize_t(stmp, strlen(ksp->kernel->name)));
3175 return 1;
3176 }
3177 if (strstr(ksp->kernel->name, "TMF")) {
3178 /* these are handled differently; the identification of the
3179 kernel is actually packaged as kernel parameters */
3180 if (!(ksp->kernel->name == strstr(ksp->kernel->name, "TMF"))) {
3181 biffAddf(NRRDnrrdBiffKey, "%s: TMF kernel name %s didn't start with TMF",
3182 me, ksp->kernel->name);
3183 return 1;
3184 }
3185 /* 0123456789012 */
3186 /* TMF_dX_cX_Xef */
3187 if (!( 13 == strlen(ksp->kernel->name)
3188 && '_' == ksp->kernel->name[3]
3189 && '_' == ksp->kernel->name[6]
3190 && '_' == ksp->kernel->name[9] )) {
3191 biffAddf(NRRDnrrdBiffKey, "%s: sorry, expected strlen(%s) = 13 with 3 _s",
3192 me, ksp->kernel->name);
3193 return 1;
3194 }
3195 sprintf(str, "tmf:%c,%c,%c",__builtin___sprintf_chk (str, 0, __builtin_object_size (str, 2
> 1 ? 1 : 0), "tmf:%c,%c,%c", ksp->kernel->name[5],
ksp->kernel->name[8], ksp->kernel->name[10])
3196 ksp->kernel->name[5],__builtin___sprintf_chk (str, 0, __builtin_object_size (str, 2
> 1 ? 1 : 0), "tmf:%c,%c,%c", ksp->kernel->name[5],
ksp->kernel->name[8], ksp->kernel->name[10])
3197 ksp->kernel->name[8],__builtin___sprintf_chk (str, 0, __builtin_object_size (str, 2
> 1 ? 1 : 0), "tmf:%c,%c,%c", ksp->kernel->name[5],
ksp->kernel->name[8], ksp->kernel->name[10])
3198 ksp->kernel->name[10])__builtin___sprintf_chk (str, 0, __builtin_object_size (str, 2
> 1 ? 1 : 0), "tmf:%c,%c,%c", ksp->kernel->name[5],
ksp->kernel->name[8], ksp->kernel->name[10]);
3199 /* see if the single parm should be added on */
3200 if (0.0 != ksp->parm[0]) {
3201 sprintf(stmp, ",%.17g", ksp->parm[0])__builtin___sprintf_chk (stmp, 0, __builtin_object_size (stmp
, 2 > 1 ? 1 : 0), ",%.17g", ksp->parm[0]);
3202 strcat(str, stmp)__builtin___strcat_chk (str, stmp, __builtin_object_size (str
, 2 > 1 ? 1 : 0));
3203 }
3204 } else {
3205 strcpy(str, ksp->kernel->name)__builtin___strcpy_chk (str, ksp->kernel->name, __builtin_object_size
(str, 2 > 1 ? 1 : 0));
3206 if (ksp->kernel->numParm) {
3207 unsigned int pi;
3208 for (pi=0; pi<ksp->kernel->numParm; pi++) {
3209 sprintf(stmp, "%c%.17g", (!pi ? ':' : ','), ksp->parm[pi])__builtin___sprintf_chk (stmp, 0, __builtin_object_size (stmp
, 2 > 1 ? 1 : 0), "%c%.17g", (!pi ? ':' : ','), ksp->parm
[pi]);
3210 if (strlen(str) + strlen(stmp) > warnLen) {
3211 biffAddf(NRRDnrrdBiffKey, "%s: kernel parm %u could overflow", me, pi);
3212 return 1;
3213 }
3214 strcat(str, stmp)__builtin___strcat_chk (str, stmp, __builtin_object_size (str
, 2 > 1 ? 1 : 0));
3215 }
3216 }
3217 }
3218 return 0;
3219}
3220
3221int
3222nrrdKernelSprint(char str[AIR_STRLEN_LARGE(512+1)], const NrrdKernel *kernel,
3223 const double kparm[NRRD_KERNEL_PARMS_NUM8]) {
3224 static const char me[]="nrrdKernelSprint";
3225 NrrdKernelSpec ksp;
3226
3227 nrrdKernelSpecSet(&ksp, kernel, kparm);
3228 if (nrrdKernelSpecSprint(str, &ksp)) {
3229 biffAddf(NRRDnrrdBiffKey, "%s: trouble", me);
3230 return 1;
3231 }
3232 return 0;
3233}
3234
3235/*
3236** This DOES make an effort to set *differ based on "ordering" (-1 or +1)
3237** but HEY that's very contrived; why bother?
3238*/
3239int
3240nrrdKernelCompare(const NrrdKernel *kernA,
3241 const double parmA[NRRD_KERNEL_PARMS_NUM8],
3242 const NrrdKernel *kernB,
3243 const double parmB[NRRD_KERNEL_PARMS_NUM8],
3244 int *differ, char explain[AIR_STRLEN_LARGE(512+1)]) {
3245 static const char me[]="nrrdKernelCompare";
3246 unsigned int pnum, pidx;
3247
3248 if (!(kernA && kernB && differ)) {
11
←
Taking false branch
→
3249 biffAddf(NRRDnrrdBiffKey, "%s: got NULL pointer (%p, %p, or %p)", me,
3250 AIR_CVOIDP(kernA)((const void *)(kernA)), AIR_CVOIDP(kernB)((const void *)(kernB)), AIR_VOIDP(differ)((void *)(differ)));
3251 return 1;
3252 }
3253 if (kernA != kernB) {
12
←
Assuming 'kernA' is equal to 'kernB'
→
13
←
Taking false branch
→
3254 *differ = kernA < kernB ? -1 : 1;
3255 if (explain) {
3256 sprintf(explain, "kernA %s kernB", *differ < 0 ? "<" : ">")__builtin___sprintf_chk (explain, 0, __builtin_object_size (explain
, 2 > 1 ? 1 : 0), "kernA %s kernB", *differ < 0 ? "<"
: ">");
3257 }
3258 return 0;
3259 }
3260 pnum = kernA->numParm;
3261 if (!pnum) {
14
←
Assuming 'pnum' is not equal to 0
→
15
←
Taking false branch
→
3262 /* actually, we're done: no parms and kernels are equal */
3263 *differ = 0;
3264 return 0;
3265 }
3266 if (!(parmA && parmB)) {
16
←
Taking true branch
→
3267 biffAddf(NRRDnrrdBiffKey, "%s: kernel %s needs %u parms but got NULL parm vectors",
3268 me, kernA->name ? kernA->name : "(unnamed)", pnum);
17
←
'?' condition is true
→
3269 return 0;
3270 }
3271 for (pidx=0; pidx<pnum; pidx++) {
3272 if (parmA[pidx] != parmB[pidx]) {
3273 *differ = parmA[pidx] < parmB[pidx] ? -1 : 1;
3274 if (explain) {
3275 sprintf(explain, "parmA[%u]=%f %s parmB[%u]=%f", pidx, parmA[pidx],__builtin___sprintf_chk (explain, 0, __builtin_object_size (explain
, 2 > 1 ? 1 : 0), "parmA[%u]=%f %s parmB[%u]=%f", pidx, parmA
[pidx], *differ < 0 ? "<" : ">", pidx, parmB[pidx])
3276 *differ < 0 ? "<" : ">", pidx, parmB[pidx])__builtin___sprintf_chk (explain, 0, __builtin_object_size (explain
, 2 > 1 ? 1 : 0), "parmA[%u]=%f %s parmB[%u]=%f", pidx, parmA
[pidx], *differ < 0 ? "<" : ">", pidx, parmB[pidx]);
3277 }
3278 return 0;
3279 }
3280 }
3281
3282 /* so far nothing unequal */
3283 *differ = 0;
3284 return 0;
3285}
3286
3287/*
3288** This DOES NOT make an effort to set *differ based on "ordering";
3289*/
3290int
3291nrrdKernelSpecCompare(const NrrdKernelSpec *aa,
3292 const NrrdKernelSpec *bb,
3293 int *differ, char explain[AIR_STRLEN_LARGE(512+1)]) {
3294 static const char me[]="nrrdKernelSpecCompare";
3295 char subexplain[AIR_STRLEN_LARGE(512+1)];
3296
3297 if (!( differ )) {
3298 biffAddf(NRRDnrrdBiffKey, "%s: got NULL differ", me);
3299 return 1;
3300 }
3301 if (!!aa != !!bb) {
3302 if (explain) {
3303 sprintf(explain, "different NULL-ities of kspec itself %s != %s",__builtin___sprintf_chk (explain, 0, __builtin_object_size (explain
, 2 > 1 ? 1 : 0), "different NULL-ities of kspec itself %s != %s"
, aa ? "non-NULL" : "NULL", bb ? "non-NULL" : "NULL")
3304 aa ? "non-NULL" : "NULL",__builtin___sprintf_chk (explain, 0, __builtin_object_size (explain
, 2 > 1 ? 1 : 0), "different NULL-ities of kspec itself %s != %s"
, aa ? "non-NULL" : "NULL", bb ? "non-NULL" : "NULL")
3305 bb ? "non-NULL" : "NULL")__builtin___sprintf_chk (explain, 0, __builtin_object_size (explain
, 2 > 1 ? 1 : 0), "different NULL-ities of kspec itself %s != %s"
, aa ? "non-NULL" : "NULL", bb ? "non-NULL" : "NULL");
3306 }
3307 *differ = 1; return 0;
3308 }
3309 if (!aa) {
3310 /* got two NULL kernel specs ==> equal */
3311 *differ = 0; return 0;
3312 }
3313 if (!!aa->kernel != !!bb->kernel) {
3314 if (explain) {
3315 sprintf(explain, "different NULL-ities of kspec->kernel %s != %s",__builtin___sprintf_chk (explain, 0, __builtin_object_size (explain
, 2 > 1 ? 1 : 0), "different NULL-ities of kspec->kernel %s != %s"
, aa->kernel ? "non-NULL" : "NULL", bb->kernel ? "non-NULL"
: "NULL")
3316 aa->kernel ? "non-NULL" : "NULL",__builtin___sprintf_chk (explain, 0, __builtin_object_size (explain
, 2 > 1 ? 1 : 0), "different NULL-ities of kspec->kernel %s != %s"
, aa->kernel ? "non-NULL" : "NULL", bb->kernel ? "non-NULL"
: "NULL")
3317 bb->kernel ? "non-NULL" : "NULL")__builtin___sprintf_chk (explain, 0, __builtin_object_size (explain
, 2 > 1 ? 1 : 0), "different NULL-ities of kspec->kernel %s != %s"
, aa->kernel ? "non-NULL" : "NULL", bb->kernel ? "non-NULL"
: "NULL");
3318 }
3319 *differ = 1; return 0;
3320 }
3321 if (!aa->kernel) {
3322 /* both kernels NULL, can't do anything informative with parms */
3323 *differ = 0; return 0;
3324 }
3325 if (nrrdKernelCompare(aa->kernel, aa->parm,
3326 bb->kernel, bb->parm,
3327 differ, subexplain)) {
3328 biffAddf(NRRDnrrdBiffKey, "%s: trouble comparing kernels", me);
3329 return 1;
3330 }
3331 if (*differ) {
3332 if (explain) {
3333 sprintf(explain, "kern/parm pairs differ: %s", subexplain)__builtin___sprintf_chk (explain, 0, __builtin_object_size (explain
, 2 > 1 ? 1 : 0), "kern/parm pairs differ: %s", subexplain
);
3334 }
3335 *differ = 1; /* losing ordering info (of dubious value) */
3336 return 0;
3337 }
3338 *differ = 0;
3339 return 0;
3340}
3341
3342/*
3343******** nrrdKernelCheck
3344**
3345** Makes sure a given kernel is behaving as expected
3346**
3347** Tests:
3348** nrrdKernelSprint
3349** nrrdKernelParse
3350** nrrdKernelCompare
3351** nrrdKernelSpecNew
3352** nrrdKernelSpecNix
3353** nrrdKernelSpecSet
3354** nrrdKernelSpecSprint
3355** nrrdKernelSpecParse
3356** and also exercises all the ways of evaluating the kernel and
3357** makes sure they all agree, and agree with the integral kernel, if given
3358*/
3359int
3360nrrdKernelCheck(const NrrdKernel *kern,
3361 const double parm[NRRD_KERNEL_PARMS_NUM8],
3362 size_t evalNum, double epsilon,
3363 unsigned int diffOkEvalMax,
3364 unsigned int diffOkIntglMax,
3365 const NrrdKernel *ikern,
3366 const double iparm[NRRD_KERNEL_PARMS_NUM8]) {
3367 const NrrdKernel *parsedkern;
3368 double parsedparm[NRRD_KERNEL_PARMS_NUM8], supp, integral;
3369 static const char me[]="nrrdKernelCheck";
3370 char kstr[AIR_STRLEN_LARGE(512+1)], kspstr[AIR_STRLEN_LARGE(512+1)],
3371 explain[AIR_STRLEN_LARGE(512+1)], stmp[AIR_STRLEN_SMALL(128+1)];
3372 int differ;
1
'differ' declared without an initial value
→
3373 size_t evalIdx;
3374 double *dom_d, *ran_d, wee;
3375 float *dom_f, *ran_f;
3376 unsigned int diffOkEvalNum, diffOkIntglNum;
3377 NrrdKernelSpec *kspA, *kspB;
3378 airArray *mop;
3379
3380 mop = airMopNew();
3381 if (!kern) {
2
←
Assuming 'kern' is non-null
→
3
←
Taking false branch
→
3382 biffAddf(NRRDnrrdBiffKey, "%s: got NULL kernel", me);
3383 airMopError(mop); return 1;
3384 }
3385 if (!(evalNum > 20)) {
4
←
Assuming 'evalNum' is > 20
→
5
←
Taking false branch
→
3386 biffAddf(NRRDnrrdBiffKey, "%s: need evalNum > 20", me);
3387 airMopError(mop); return 1;
3388 }
3389 if (!(kern->support && kern->integral
6
←
Taking false branch
→
3390 && kern->eval1_f && kern->evalN_f
3391 && kern->eval1_d && kern->evalN_d)) {
3392 biffAddf(NRRDnrrdBiffKey, "%s: kernel has NULL fields (%d,%d,%d,%d,%d,%d)", me,
3393 !!(kern->support), !!(kern->integral),
3394 !!(kern->eval1_f), !!(kern->evalN_f),
3395 !!(kern->eval1_d), !!(kern->evalN_d));
3396 airMopError(mop); return 1;
3397 }
3398 kspA = nrrdKernelSpecNew();
3399 airMopAdd(mop, kspA, (airMopper)nrrdKernelSpecNix, airMopAlways);
3400 kspB = nrrdKernelSpecNew();
3401 airMopAdd(mop, kspB, (airMopper)nrrdKernelSpecNix, airMopAlways);
3402 nrrdKernelSpecSet(kspA, kern, parm);
3403 if (nrrdKernelSprint(kstr, kern, parm)
7
←
Taking false branch
→
3404 || nrrdKernelSpecSprint(kspstr, kspA)) {
3405 biffAddf(NRRDnrrdBiffKey, "%s: trouble", me);
3406 airMopError(mop); return 1;
3407 }
3408 if (strcmp(kstr, kspstr)) {
8
←
Taking false branch
→
3409 biffAddf(NRRDnrrdBiffKey, "%s: sprinted kernel |%s| != kspec |%s|", me, kstr, kspstr);
3410 airMopError(mop); return 1;
3411 }
3412 if (nrrdKernelParse(&parsedkern, parsedparm, kstr)
9
←
Taking false branch
→
3413 || nrrdKernelSpecParse(kspB, kstr)) {
3414 biffAddf(NRRDnrrdBiffKey, "%s: trouble parsing |%s| back to kern/parm pair or kspec",
3415 me, kstr);
3416 airMopError(mop); return 1;
3417 }
3418 if (nrrdKernelCompare(kern, parm, parsedkern, parsedparm,
10
←
Calling 'nrrdKernelCompare'
→
18
←
Returning from 'nrrdKernelCompare'
→
19
←
Taking false branch
→
3419 &differ, explain)) {
3420 biffAddf(NRRDnrrdBiffKey, "%s: trouble comparing kern/parm pairs", me);
3421 airMopError(mop); return 1;
3422 }
3423 if (differ) {
20
←
Branch condition evaluates to a garbage value
3424 biffAddf(NRRDnrrdBiffKey, "%s: given and re-parsed kernels differ: %s", me, explain);
3425 airMopError(mop); return 1;
3426 }
3427 if (nrrdKernelCompare(kspA->kernel, kspA->parm,
3428 kspB->kernel, kspB->parm,
3429 &differ, explain)) {
3430 biffAddf(NRRDnrrdBiffKey, "%s: trouble comparing kspecs", me);
3431 airMopError(mop); return 1;
3432 }
3433 if (differ) {
3434 biffAddf(NRRDnrrdBiffKey, "%s: given and re-parsed kspecs differ: %s", me, explain);
3435 airMopError(mop); return 1;
3436 }
3437
3438 supp = kern->support(parm);
3439 /* wee is the step between evaluation points */
3440 wee = 2*supp/AIR_CAST(double, evalNum)((double)(evalNum));
3441 if ( (kern->eval1_d)(supp+wee/1000, parm) ||
3442 (kern->eval1_d)(supp+wee, parm) ||
3443 (kern->eval1_d)(supp+10*wee, parm) ||
3444 (kern->eval1_d)(-supp-wee/1000, parm) ||
3445 (kern->eval1_d)(-supp-wee, parm) ||
3446 (kern->eval1_d)(-supp-10*wee, parm) ) {
3447 if (nrrdKernelCheap != kern) {
3448 /* the "cheap" kernel alone gets a pass on reporting its support */
3449 biffAddf(NRRDnrrdBiffKey, "%s: kern %s is non-zero outside support %g",
3450 me, kstr, supp);
3451 airMopError(mop); return 1;
3452 }
3453 }
3454 /* allocate domain and range for both float and double */
3455 dom_d = AIR_CALLOC(evalNum, double)(double*)(calloc((evalNum), sizeof(double)));
3456 airMopAdd(mop, dom_d, airFree, airMopAlways);
3457 ran_d = AIR_CALLOC(evalNum, double)(double*)(calloc((evalNum), sizeof(double)));
3458 airMopAdd(mop, ran_d, airFree, airMopAlways);
3459 dom_f = AIR_CALLOC(evalNum, float)(float*)(calloc((evalNum), sizeof(float)));
3460 airMopAdd(mop, dom_f, airFree, airMopAlways);
3461 ran_f = AIR_CALLOC(evalNum, float)(float*)(calloc((evalNum), sizeof(float)));
3462 airMopAdd(mop, ran_f, airFree, airMopAlways);
3463 if (!( dom_d && ran_d && dom_f && ran_f )) {
3464 biffAddf(NRRDnrrdBiffKey, "%s: couldn't alloc buffers for %s values for %s",
3465 me, airSprintSize_t(stmp, evalNum), kstr);
3466 airMopError(mop); return 1;
3467 }
3468 for (evalIdx=0; evalIdx<evalNum; evalIdx++) {
3469 dom_d[evalIdx] = AIR_AFFINE(-0.5, evalIdx,( ((double)(supp)-(-supp))*((double)(evalIdx)-(-0.5)) / ((double
)(((double)(evalNum))-0.5)-(-0.5)) + (-supp))
3470 AIR_CAST(double, evalNum)-0.5,( ((double)(supp)-(-supp))*((double)(evalIdx)-(-0.5)) / ((double
)(((double)(evalNum))-0.5)-(-0.5)) + (-supp))
3471 -supp, supp)( ((double)(supp)-(-supp))*((double)(evalIdx)-(-0.5)) / ((double
)(((double)(evalNum))-0.5)-(-0.5)) + (-supp));
3472 dom_f[evalIdx] = AIR_CAST(float, dom_d[evalIdx])((float)(dom_d[evalIdx]));
3473 }
3474 /* do the vector evaluations */
3475 kern->evalN_f(ran_f, dom_f, evalNum, parm);
3476 kern->evalN_d(ran_d, dom_d, evalNum, parm);
3477 /*
3478 for (evalIdx=0; evalIdx<evalNum; evalIdx++) {
3479 fprintf(stderr, "%u %g --> %g\n", AIR_CAST(unsigned int, evalIdx),
3480 dom_d[evalIdx], ran_d[evalIdx]);
3481 }
3482 */
3483 /* compare evaluations (and maybe derivatives) and numerically compute
3484 integral */
3485 diffOkEvalNum = 0;
3486 diffOkIntglNum = 0;
3487 integral = 0.0;
3488 for (evalIdx=0; evalIdx<evalNum; evalIdx++) {
3489 double single_f, single_d;
3490 single_f = kern->eval1_f(dom_f[evalIdx], parm);
3491 single_d = kern->eval1_d(dom_d[evalIdx], parm);
3492 integral += single_d;
3493 /* single float vs vector float */
3494 if (nrrdKernelForwDiff == kern
3495 || nrrdKernelBCCubic == kern
3496 || nrrdKernelBCCubicDD == kern
3497 || nrrdKernelAQuarticDD == kern) {
3498 /* HEY this is crazy: need a special epsilon for these kernels;
3499 WHY WHY do these kernels evaluate to different things in the
3500 single versus the vector case? */
3501 float specEps;
3502 if (nrrdKernelForwDiff == kern) {
3503 specEps = 5e-9f;
3504 } else if (nrrdKernelBCCubic == kern) {
3505 specEps = 5e-8f;
3506 } else if (nrrdKernelBCCubicDD == kern) {
3507 specEps = 5e-8f;
3508 } else if (nrrdKernelAQuarticDD == kern) {
3509 specEps = 5e-8f;
3510 } else {
3511 specEps = 0.0;
3512 }
3513 if (fabs(single_f - ran_f[evalIdx]) > specEps) {
3514 biffAddf(NRRDnrrdBiffKey, "%s: %s (eval1_f(%.17g)=%.17g) != "
3515 "(evalN_f(%.17g)=%.17g) by %.17g > %.17g",
3516 me, kstr, dom_f[evalIdx], single_f,
3517 dom_f[evalIdx], ran_f[evalIdx],
3518 fabs(single_f - ran_f[evalIdx]), specEps);
3519 airMopError(mop); return 1;
3520 }
3521 } else {
3522 if (single_f != ran_f[evalIdx]) {
3523 biffAddf(NRRDnrrdBiffKey, "%s: %s (eval1_f(%.17g)=%.17g) != "
3524 "(evalN_f(%.17g)=%.17g)",
3525 me, kstr, dom_f[evalIdx], single_f,
3526 dom_f[evalIdx], ran_f[evalIdx]);
3527 airMopError(mop); return 1;
3528 }
3529 }
3530 /* single double vs vector double */
3531 if (single_d != ran_d[evalIdx]) {
3532 biffAddf(NRRDnrrdBiffKey, "%s: %s (eval1_d(%.17g)=%.17g) != (evalN_d(%.17g)=%.17g)",
3533 me, kstr, dom_d[evalIdx], single_d,
3534 dom_d[evalIdx], ran_d[evalIdx]);
3535 airMopError(mop); return 1;
3536 }
3537 /* single float vs single double */
3538 if (fabs(single_f - single_d) > epsilon) {
3539 diffOkEvalNum++;
3540 if (diffOkEvalNum > diffOkEvalMax) {
3541 biffAddf(NRRDnrrdBiffKey,
3542 "%s: %s |eval1_f(%.17g)=%.17g) - (eval1_d(%.17g)=%.17g)|"
3543 " %.17g > epsilon %.17g too many times (%u > %u)", me,
3544 kstr, dom_f[evalIdx], single_f, dom_d[evalIdx], single_d,
3545 fabs(single_f - single_d), epsilon,
3546 diffOkEvalNum, diffOkEvalMax);
3547 airMopError(mop); return 1;
3548 }
3549 }
3550 /* check whether we're the derivative of ikern */
3551 if (ikern) {
3552 double forw, back, ndrv;
3553 forw = ikern->eval1_d(dom_d[evalIdx] + wee/2, iparm);
3554 back = ikern->eval1_d(dom_d[evalIdx] - wee/2, iparm);
3555 ndrv = (forw - back)/wee;
3556 if (fabs(ndrv - single_d) > epsilon) {
3557 diffOkIntglNum++;
3558 if (diffOkIntglNum > diffOkIntglMax) {
3559 biffAddf(NRRDnrrdBiffKey, "%s: %s(%.17g) |num deriv(%s) %.17g - %.17g| "
3560 "%.17g > %.17g too many times (%u > %u)",
3561 me, kstr, dom_d[evalIdx], ikern->name, ndrv, single_d,
3562 fabs(ndrv - single_d), epsilon,
3563 diffOkIntglNum, diffOkIntglMax);
3564 airMopError(mop); return 1;
3565 }
3566 }
3567 }
3568 }
3569 integral *= 2*supp/(AIR_CAST(double, evalNum)((double)(evalNum)));
3570 /* the "cheap" kernel alone gets a pass on reporting its integral */
3571 if (nrrdKernelCheap != kern) {
3572 double hackeps=10;
3573 /* hackeps is clearly a hack to permit the integral to have greater
3574 error than any single evaluation; there must be a more principled
3575 way to set this */
3576 if (fabs(integral - kern->integral(parm)) > hackeps*epsilon) {
3577 biffAddf(NRRDnrrdBiffKey, "%s: %s |numerical integral %.17g - claimed %.17g| "
3578 "%.17g > %.17g", me, kstr, integral, kern->integral(parm),
3579 fabs(integral - kern->integral(parm)), hackeps*epsilon);
3580 airMopError(mop); return 1;
3581 }
3582 }
3583
3584 /* HEY check being derivative of ikern/iparm */
3585 AIR_UNUSED(ikern)(void)(ikern);
3586 AIR_UNUSED(iparm)(void)(iparm);
3587
3588 airMopOkay(mop);
3589 return 0;
3590}
3591
3592int
3593nrrdKernelParm0IsScale(const NrrdKernel *kern) {
3594 int ret;
3595
3596 if (!kern) {
3597 ret = 0;
3598 } else if (nrrdKernelHann == kern ||
3599 nrrdKernelHannD == kern ||
3600 nrrdKernelHannDD == kern ||
3601 nrrdKernelBlackman == kern ||
3602 nrrdKernelBlackmanD == kern ||
3603 nrrdKernelBlackmanDD == kern ||
3604 nrrdKernelZero == kern ||
3605 nrrdKernelBox == kern ||
3606 nrrdKernelCheap == kern ||
3607 nrrdKernelTent == kern ||
3608 nrrdKernelForwDiff == kern ||
3609 nrrdKernelCentDiff == kern ||
3610 nrrdKernelBCCubic == kern ||
3611 nrrdKernelBCCubicD == kern ||
3612 nrrdKernelBCCubicDD == kern ||
3613 nrrdKernelAQuartic == kern ||
3614 nrrdKernelAQuarticD == kern ||
3615 nrrdKernelAQuarticDD == kern ||
3616 nrrdKernelGaussian == kern ||
3617 nrrdKernelGaussianD == kern ||
3618 nrrdKernelGaussianDD == kern ||
3619 nrrdKernelDiscreteGaussian == kern) {
3620 ret = 1;
3621 } else {
3622 ret = 0;
3623 }
3624 return ret;
3625}