/Users/scm/git/github/teem/src/nrrd/kernel.c

Bug Summary

File:	src/nrrd/kernel.c
Location:	line 165, column 15
Description:	Value stored to 't' is never read

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
37	nrrdKernelBSpline3ApproxInverse 0
38	nrrdKernelBSpline5 0
39	nrrdKernelBSpline5ApproxInverse 0
40	nrrdKernelBSpline7 0
41	nrrdKernelBSpline7ApproxInverse 0
42	nrrdKernelZero 1 scale
43	nrrdKernelBox 1 scale
44	nrrdKernelCatmullRomSupportDebug 1 support
45	nrrdKernelBoxSupportDebug 1 support
46	nrrdKernelCos4SupportDebug 1 support
47	nrrdKernelCheap 1 scale
48	nrrdKernelHermiteScaleSpaceFlag 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
72	static double
73	returnZero(const double *parm) {
74	AIR_UNUSED(parm)(void)(parm);
75	return 0.0;
76	}
77
78	static double
79	returnOne(const double *parm) {
80	AIR_UNUSED(parm)(void)(parm);
81	return 1.0;
82	}
83
84	static double
85	returnTwo(const double *parm) {
86	AIR_UNUSED(parm)(void)(parm);
87	return 2.0;
88	}
89
90	static double
91	returnThree(const double *parm) {
92	AIR_UNUSED(parm)(void)(parm);
93	return 3.0;
94	}
95
96	static double
97	returnFour(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
119	static double
120	_nrrdZeroSup(const double *parm) {
121	double S;
122
123	S = parm[0];
124	return S;
125	}
126
127	static 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
136	static 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
145	static 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
158	static 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;
	Value stored to 't' is never read
166	f[i] = _ZERO(t)0/S;
167	}
168	}
169
170	static NrrdKernel
171	_nrrdKernelZero = {
172	"zero",
173	1, _nrrdZeroSup, returnZero,
174	_nrrdZero1_f, _nrrdZeroN_f, _nrrdZero1_d, _nrrdZeroN_d
175	};
176	NrrdKernel *const
177	nrrdKernelZero = &_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
183	static 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
193	static 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
202	static 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
211	static 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
224	static 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
236	static NrrdKernel
237	_nrrdKernelBox = {
238	"box",
239	1, _nrrdBoxSup, returnOne,
240	_nrrdBox1_f, _nrrdBoxN_f, _nrrdBox1_d, _nrrdBoxN_d
241	};
242	NrrdKernel *const
243	nrrdKernelBox = &_nrrdKernelBox;
244
245	/* ------------------------------------------------------------ */
246
247	static double
248	_nrrdBoxSDSup(const double *parm) {
249
250	return parm[0];
251	}
252
253	static 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
260	static 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
267	static 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
278	static 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
289	static NrrdKernel
290	_nrrdKernelBoxSupportDebug = {
291	"boxsup",
292	1, _nrrdBoxSDSup, returnOne,
293	_nrrdBoxSD1_f, _nrrdBoxSDN_f, _nrrdBoxSD1_d, _nrrdBoxSDN_d
294	};
295	NrrdKernel *const
296	nrrdKernelBoxSupportDebug = &_nrrdKernelBoxSupportDebug;
297
298	/* ------------------------------------------------------------ */
299
300	#define COS4(x)(x > 0.5 ? 0.0 : cos(3.14159265358979323846x)cos(3.14159265358979323846 x)cos(3.14159265358979323846x)cos(3.14159265358979323846* x)) (x > 0.5 \
301	? 0.0 \
302	: cos(AIR_PI3.14159265358979323846x)cos(AIR_PI3.14159265358979323846x)cos(AIR_PI3.14159265358979323846x)cos(AIR_PI3.14159265358979323846*x))
303
304	static double
305	_nrrdCos4SDInt(const double *parm) {
306	AIR_UNUSED(parm)(void)(parm);
307	return 3.0/8.0;
308	}
309
310	static double
311	_nrrdCos4SDSup(const double *parm) {
312
313	return parm[0];
314	}
315
316	static 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.14159265358979323846x)cos(3.14159265358979323846 x)cos(3.14159265358979323846x)cos(3.14159265358979323846* x));
321	}
322
323	static 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.14159265358979323846x)cos (3.14159265358979323846x)cos(3.14159265358979323846x)cos( 3.14159265358979323846*x))));
328	}
329
330	static 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.14159265358979323846t)cos(3.14159265358979323846 t)cos(3.14159265358979323846t)cos(3.14159265358979323846* t));
338	}
339	}
340
341	static 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.14159265358979323846t)cos (3.14159265358979323846t)cos(3.14159265358979323846t)cos( 3.14159265358979323846*t))));
349	}
350	}
351
352	static NrrdKernel
353	_nrrdKernelCos4SupportDebug = {
354	"cos4sup",
355	1, _nrrdCos4SDSup, _nrrdCos4SDInt,
356	_nrrdCos4SD1_f, _nrrdCos4SDN_f, _nrrdCos4SD1_d, _nrrdCos4SDN_d
357	};
358	NrrdKernel *const
359	nrrdKernelCos4SupportDebug = &_nrrdKernelCos4SupportDebug;
360
361	/* ------------------------------------------------------------ */
362
363	#define DCOS4(x)(x > 0.5 ? 0.0 : -43.14159265358979323846(cos(3.14159265358979323846 x)cos(3.14159265358979323846x)cos(3.14159265358979323846* x) sin(3.14159265358979323846x))) (x > 0.5 \
364	? 0.0 \
365	: -4AIR_PI3.14159265358979323846(cos(AIR_PI3.14159265358979323846x)cos(AIR_PI3.14159265358979323846x)cos(AIR_PI3.14159265358979323846*x) \
366	sin(AIR_PI3.14159265358979323846x)))
367
368	static double
369	_nrrdDCos4SDSup(const double *parm) {
370	return parm[0];
371	}
372
373	static 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 sgnDCOS4(x)(x > 0.5 ? 0.0 : -43.14159265358979323846(cos(3.14159265358979323846 x)cos(3.14159265358979323846x)cos(3.14159265358979323846 x) sin(3.14159265358979323846x)));
379	}
380
381	static 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, sgnDCOS4(x))((float)(sgn(x > 0.5 ? 0.0 : -43.14159265358979323846(cos (3.14159265358979323846x)cos(3.14159265358979323846x)cos( 3.14159265358979323846x) sin(3.14159265358979323846*x)))));
387	}
388
389	static 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] = sgnDCOS4(t)(t > 0.5 ? 0.0 : -43.14159265358979323846(cos(3.14159265358979323846 t)cos(3.14159265358979323846t)cos(3.14159265358979323846 t) sin(3.14159265358979323846t)));
399	}
400	}
401
402	static 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, sgnDCOS4(t))((float)(sgn(t > 0.5 ? 0.0 : -43.14159265358979323846(cos (3.14159265358979323846t)cos(3.14159265358979323846t)cos( 3.14159265358979323846t) sin(3.14159265358979323846*t)))));
412	}
413	}
414
415	static NrrdKernel
416	_nrrdKernelCos4SupportDebugD = {
417	"cos4supD",
418	1, _nrrdDCos4SDSup, returnZero,
419	_nrrdDCos4SD1_f, _nrrdDCos4SDN_f, _nrrdDCos4SD1_d, _nrrdDCos4SDN_d
420	};
421	NrrdKernel *const
422	nrrdKernelCos4SupportDebugD = &_nrrdKernelCos4SupportDebugD;
423
424	/* ------------------------------------------------------------ */
425
426	#define DDCOS4(x)(x > 0.5 ? 0.0 : -23.141592653589793238463.14159265358979323846 (cos(3.141592653589793238462x) + cos(3.14159265358979323846 4*x))) (x > 0.5 \
427	? 0.0 \
428	: -2AIR_PI3.14159265358979323846AIR_PI3.14159265358979323846(cos(AIR_PI3.141592653589793238462x) + cos(AIR_PI3.141592653589793238464*x)))
429
430	static double
431	_nrrdDDCos4SDSup(const double *parm) {
432	return parm[0];
433	}
434
435	static 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 : -23.141592653589793238463.14159265358979323846 (cos(3.141592653589793238462x) + cos(3.14159265358979323846 4*x)));
440	}
441
442	static 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 : -23.141592653589793238463.14159265358979323846 (cos(3.141592653589793238462x) + cos(3.14159265358979323846 4*x)))));
447	}
448
449	static 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 : -23.141592653589793238463.14159265358979323846 (cos(3.141592653589793238462t) + cos(3.14159265358979323846 4*t)));
457	}
458	}
459
460	static 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 : -23.141592653589793238463.14159265358979323846 (cos(3.141592653589793238462t) + cos(3.14159265358979323846 4*t)))));
468	}
469	}
470
471	static NrrdKernel
472	_nrrdKernelCos4SupportDebugDD = {
473	"cos4supDD",
474	1, _nrrdDDCos4SDSup, returnZero,
475	_nrrdDDCos4SD1_f, _nrrdDDCos4SDN_f, _nrrdDDCos4SD1_d, _nrrdDDCos4SDN_d
476	};
477	NrrdKernel *const
478	nrrdKernelCos4SupportDebugDD = &_nrrdKernelCos4SupportDebugDD;
479
480	/* ------------------------------------------------------------ */
481
482	#define DDDCOS4(x)(x > 0.5 ? 0.0 : 43.141592653589793238463.14159265358979323846 3.14159265358979323846(sin(23.14159265358979323846x) + 2* sin(43.14159265358979323846x))) (x > 0.5 \
483	? 0.0 \
484	: 4AIR_PI3.14159265358979323846AIR_PI3.14159265358979323846AIR_PI3.14159265358979323846(sin(2AIR_PI3.14159265358979323846x) + 2sin(4AIR_PI3.14159265358979323846*x)))
485
486	static double
487	_nrrdDDDCos4SDSup(const double *parm) {
488	return parm[0];
489	}
490
491	static 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 sgnDDDCOS4(x)(x > 0.5 ? 0.0 : 43.141592653589793238463.14159265358979323846 3.14159265358979323846(sin(23.14159265358979323846x) + 2 sin(43.14159265358979323846x)));
497	}
498
499	static 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, sgnDDDCOS4(x))((float)(sgn(x > 0.5 ? 0.0 : 43.141592653589793238463.14159265358979323846 3.14159265358979323846(sin(23.14159265358979323846x) + 2* sin(43.14159265358979323846x)))));
505	}
506
507	static 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] = sgnDDDCOS4(t)(t > 0.5 ? 0.0 : 43.141592653589793238463.14159265358979323846 3.14159265358979323846(sin(23.14159265358979323846t) + 2 sin(43.14159265358979323846t)));
517	}
518	}
519
520	static 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, sgnDDDCOS4(t))((float)(sgn(t > 0.5 ? 0.0 : 43.141592653589793238463.14159265358979323846 3.14159265358979323846(sin(23.14159265358979323846t) + 2* sin(43.14159265358979323846t)))));
530	}
531	}
532
533	static NrrdKernel
534	_nrrdKernelCos4SupportDebugDDD = {
535	"cos4supDDD",
536	1, _nrrdDDDCos4SDSup, returnZero,
537	_nrrdDDDCos4SD1_f, _nrrdDDDCos4SDN_f, _nrrdDDDCos4SD1_d, _nrrdDDDCos4SDN_d
538	};
539	NrrdKernel *const
540	nrrdKernelCos4SupportDebugDDD = &_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
550	static 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
560	static double
561	_nrrdCheap1_d(double x, const double *parm) {
562
563	return _CHEAP(x)((x) > 0.0f ? (x) : -(x))/parm[0];
564	}
565
566	static 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
572	static 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
583	static 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
594	static NrrdKernel
595	_nrrdKernelCheap = {
596	"cheap",
597	1, _nrrdCheapSup, returnOne,
598	_nrrdCheap1_f, _nrrdCheapN_f, _nrrdCheap1_d, _nrrdCheapN_d
599	};
600	NrrdKernel *const
601	nrrdKernelCheap = &_nrrdKernelCheap;
602
603	/* ------------------------------------------------------------ */
604
605	#define _TENT(x)(x >= 1 ? 0 : 1 - x) (x >= 1 ? 0 : 1 - x)
606
607	static double
608	_nrrdTentSup(const double *parm) {
609	double S;
610
611	S = parm[0];
612	return S;
613	}
614
615	static 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
624	static 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
633	static 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
646	static 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
658	static NrrdKernel
659	_nrrdKernelTent = {
660	"tent",
661	1, _nrrdTentSup, returnOne,
662	_nrrdTent1_f, _nrrdTentN_f, _nrrdTent1_d, _nrrdTentN_d
663	};
664	NrrdKernel *const
665	nrrdKernelTent = &_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
681	static 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
688	static 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
695	static 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
706	static 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 */
719	static NrrdKernel
720	_nrrdKernelHermiteScaleSpaceFlag = {
721	"hermiteSS",
722	0, returnOne, returnOne,
723	_nrrdHermite1_f, _nrrdHermiteN_f, _nrrdHermite1_d, _nrrdHermiteN_d
724	};
725	NrrdKernel *const
726	nrrdKernelHermiteScaleSpaceFlag = &_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
734	static double
735	_nrrdFDSup(const double *parm) {
736	double S;
737
738	S = parm[0];
739	return S+0.0001; /* sigh */
740	}
741
742	static 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
751	static 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
760	static 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
772	static 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
784	static NrrdKernel
785	_nrrdKernelFD = {
786	"fordif",
787	1, _nrrdFDSup, returnZero,
788	_nrrdFD1_f, _nrrdFDN_f, _nrrdFD1_d, _nrrdFDN_d
789	};
790	NrrdKernel *const
791	nrrdKernelForwDiff = &_nrrdKernelFD;
792
793	/* ------------------------------------------------------------ */
794
795	#define _CENDIF(x)(x <= -2 ? 0 : (x <= -1 ? 0.5x + 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
800	static double
801	_nrrdCDSup(const double *parm) {
802	double S;
803
804	S = parm[0];
805	return 2*S;
806	}
807
808	static 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.5x + 1 : (x <= 1 ? -0.5 x : (x <= 2 ? 0.5x - 1 : 0 ))))/(SS);
815	}
816
817	static 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)/(SS))((float)((x <= -2 ? 0 : (x <= -1 ? 0.5x + 1 : (x <= 1 ? -0.5x : (x <= 2 ? 0.5x - 1 : 0 ))))/(S*S)));
824	}
825
826	static 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.5t + 1 : (t <= 1 ? -0.5 t : (t <= 2 ? 0.5t - 1 : 0 ))))/(SS);
836	}
837	}
838
839	static 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)/(SS))((float)((t <= -2 ? 0 : (t <= -1 ? 0.5t + 1 : (t <= 1 ? -0.5t : (t <= 2 ? 0.5t - 1 : 0 ))))/(S*S)));
848	}
849	}
850
851	static NrrdKernel
852	_nrrdKernelCD = {
853	"cendif",
854	1, _nrrdCDSup, returnZero,
855	_nrrdCD1_f, _nrrdCDN_f, _nrrdCD1_d, _nrrdCDN_d
856	};
857	NrrdKernel *const
858	nrrdKernelCentDiff = &_nrrdKernelCD;
859
860	/* ------------------------------------------------------------ */
861
862	#define _BCCUBIC(x, B, C)(x >= 2.0 ? 0 : (x >= 1.0 ? (((-B/6 - C)x + B + 5C)x -2B - 8C)x + 4B/3 + 4C : ((2 - 3B/2 - C)x - 3 + 2B + C)x*x + 1 - B/3)) \
863	(x >= 2.0 ? 0 : \
864	(x >= 1.0 \
865	? (((-B/6 - C)x + B + 5C)x -2B - 8C)x + 4B/3 + 4C \
866	: ((2 - 3B/2 - C)x - 3 + 2B + C)x*x + 1 - B/3))
867
868	static double
869	_nrrdBCSup(const double *parm) {
870	double S;
871
872	S = parm[0];
873	return 2*S;
874	}
875
876	static 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 + 5C)x -2B - 8C)x + 4B/3 + 4C : ((2 - 3B/2 - C)x - 3 + 2B + C)x*x + 1 - B/3))/S;
884	}
885
886	static 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 + 5C)x -2B - 8C)x + 4B/3 + 4C : ((2 - 3B/2 - C)x - 3 + 2B + C)x*x + 1 - B/3))/S;
895	}
896
897	static 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 + 5C)t -2B - 8C)t + 4B/3 + 4C : ((2 - 3B/2 - C)t - 3 + 2B + C)t*t + 1 - B/3))/S;
908	}
909	}
910
911	static 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 + 5C)t -2B - 8C)t + 4B/3 + 4C : ((2 - 3B/2 - C)t - 3 + 2B + C)t*t + 1 - B/3))/S;
923	}
924	}
925
926	static NrrdKernel
927	_nrrdKernelBC = {
928	"BCcubic",
929	3, _nrrdBCSup, returnOne,
930	_nrrdBC1_f, _nrrdBCN_f, _nrrdBC1_d, _nrrdBCN_d
931	};
932	NrrdKernel *const
933	nrrdKernelBCCubic = &_nrrdKernelBC;
934
935	/* ------------------------------------------------------------ */
936
937	#define _DBCCUBIC(x, B, C)(x >= 2.0 ? 0.0 : (x >= 1.0 ? ((-B/2 - 3C)x + 2B + 10 C)x -2B - 8C : ((6 - 9B/2 - 3C)x - 6 + 4B + 2C)*x)) \
938	(x >= 2.0 ? 0.0 : \
939	(x >= 1.0 \
940	? ((-B/2 - 3C)x + 2B + 10C)x -2B - 8*C \
941	: ((6 - 9B/2 - 3C)x - 6 + 4B + 2C)x))
942
943	static double
944	_nrrdDBCSup(const double *parm) {
945	double S;
946
947	S = parm[0];
948	return 2*S;
949	}
950
951	static 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 - 3C)x + 2B + 10 C)x -2B - 8C : ((6 - 9B/2 - 3C)x - 6 + 4B + 2C)x))/(S*S);
961	}
962
963	static 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)/(SS))((float)(sgn(x >= 2.0 ? 0.0 : (x >= 1.0 ? ((-B/2 - 3C )x + 2B + 10C)x -2B - 8C : ((6 - 9B/2 - 3C)x - 6 + 4 B + 2C)x))/(S*S)));
974	}
975
976	static 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 - 3C)t + 2B + 10 C)t -2B - 8C : ((6 - 9B/2 - 3C)t - 6 + 4B + 2C)t))/(S*S);
988	}
989	}
990
991	static 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)/(SS))((float)(sgn(t >= 2.0 ? 0.0 : (t >= 1.0 ? ((-B/2 - 3C )t + 2B + 10C)t -2B - 8C : ((6 - 9B/2 - 3C)t - 6 + 4 B + 2C)t))/(S*S)));
1004	}
1005	}
1006
1007	static NrrdKernel
1008	_nrrdKernelDBC = {
1009	"BCcubicD",
1010	3, _nrrdDBCSup, returnZero,
1011	_nrrdDBC1_f, _nrrdDBCN_f, _nrrdDBC1_d, _nrrdDBCN_d
1012	};
1013	NrrdKernel *const
1014	nrrdKernelBCCubicD = &_nrrdKernelDBC;
1015
1016	/* ------------------------------------------------------------ */
1017
1018	#define _DDBCCUBIC(x, B, C)(x >= 2.0 ? 0 : (x >= 1.0 ? (-B - 6C)x + 2B + 10C : (12 - 9B - 6C)x - 6 + 4B + 2*C )) \
1019	(x >= 2.0 ? 0 : \
1020	(x >= 1.0 \
1021	? (-B - 6C)x + 2B + 10C \
1022	: (12 - 9B - 6C)x - 6 + 4B + 2*C ))
1023
1024	static double
1025	_nrrdDDBCSup(const double *parm) {
1026	double S;
1027
1028	S = parm[0];
1029	return 2*S;
1030	}
1031
1032	static 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 - 6C)x + 2B + 10C : (12 - 9B - 6C)x - 6 + 4B + 2C ))/(SS*S);
1040	}
1041
1042	static 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 - 6C)x + 2B + 10C : (12 - 9B - 6C)x - 6 + 4B + 2C ))/(SS*S);
1051	}
1052
1053	static 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 - 6C)t + 2B + 10C : (12 - 9B - 6C)t - 6 + 4B + 2C ))/(SS*S);
1064	}
1065	}
1066
1067	static 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 - 6C)t + 2B + 10C : (12 - 9B - 6C)t - 6 + 4B + 2C ))/(SS*S);
1079	}
1080	}
1081
1082	static NrrdKernel
1083	_nrrdKernelDDBC = {
1084	"BCcubicDD",
1085	3, _nrrdDDBCSup, returnZero,
1086	_nrrdDDBC1_f, _nrrdDDBCN_f, _nrrdDDBC1_d, _nrrdDDBCN_d
1087	};
1088	NrrdKernel *const
1089	nrrdKernelBCCubicDD = &_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 -20.0 - 80.5)x + 40.0/3 + 40.5 : ((2 - 30.0/2 - 0.5)x - 3 + 20.0 + 0.5)xx + 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 -20.0 - 80.5)x + 40.0/3 + 40.5 : ((2 - 30.0/2 - 0.5)x - 3 + 20.0 + 0.5)xx + 1 - 0.0/3))
1095
1096	static 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 -20.0 - 80.5)x + 40.0/3 + 40.5 : ((2 - 30.0/2 - 0.5)x - 3 + 20.0 + 0.5)x*x + 1 - 0.0/3));
1101	}
1102
1103	static 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 + 50.5)x -20.0 - 80.5)x + 40.0/3 + 40.5 : ((2 - 30.0/2 - 0.5)x - 3 + 20.0 + 0.5)x*x + 1 - 0.0/3))));
1108	}
1109
1110	static 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 -20.0 - 80.5)t + 40.0/3 + 40.5 : ((2 - 30.0/2 - 0.5)t - 3 + 20.0 + 0.5)t*t + 1 - 0.0/3));
1119	}
1120	}
1121
1122	static 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 + 50.5)t -20.0 - 80.5)t + 40.0/3 + 40.5 : ((2 - 30.0/2 - 0.5)t - 3 + 20.0 + 0.5)t*t + 1 - 0.0/3))));
1131	}
1132	}
1133
1134	static NrrdKernel
1135	_nrrdKernelCatmullRom = {
1136	"catmull-rom",
1137	0, returnTwo, returnOne,
1138	_nrrdCTMR1_f, _nrrdCTMRN_f, _nrrdCTMR1_d, _nrrdCTMRN_d
1139	};
1140	NrrdKernel *const
1141	nrrdKernelCatmullRom = &_nrrdKernelCatmullRom;
1142
1143
1144	static 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
1150	static NrrdKernel
1151	_nrrdKernelCatmullRomSupportDebug = {
1152	"ctmrsup",
1153	1, _nrrdCtmrSDSup, returnOne,
1154	_nrrdCTMR1_f, _nrrdCTMRN_f, _nrrdCTMR1_d, _nrrdCTMRN_d
1155	};
1156	NrrdKernel *const
1157	nrrdKernelCatmullRomSupportDebug = &_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 - 30.5)x + 20.0 + 100.5)x -20.0 - 80.5 : ((6 - 90.0/2 - 30.5)x - 6 + 4 0.0 + 20.5)x)) _DBCCUBIC(x, 0.0, 0.5)(x >= 2.0 ? 0.0 : (x >= 1.0 ? ((-0.0/2 - 30.5)x + 20.0 + 100.5)x -20.0 - 80.5 : ((6 - 90.0/2 - 30.5)x - 6 + 4 0.0 + 20.5)x))
1163
1164	static 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 - 30.5)x + 20.0 + 100.5)x -20.0 - 80.5 : ((6 - 90.0/2 - 30.5)x - 6 + 4 0.0 + 20.5)x));
1170	}
1171
1172	static 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 + 20.0 + 100.5)x -20.0 - 80.5 : ((6 - 90.0/2 - 3 0.5)x - 6 + 40.0 + 20.5)*x))));
1178	}
1179
1180	static 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 - 30.5)t + 20.0 + 100.5)t -20.0 - 80.5 : ((6 - 90.0/2 - 30.5)t - 6 + 4 0.0 + 20.5)t));
1190	}
1191	}
1192
1193	static 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 + 20.0 + 100.5)t -20.0 - 80.5 : ((6 - 90.0/2 - 3 0.5)t - 6 + 40.0 + 20.5)*t))));
1203	}
1204	}
1205
1206	static NrrdKernel
1207	_nrrdKernelCatmullRomD = {
1208	"catmull-romD",
1209	0, returnTwo, returnZero,
1210	_nrrdDCTMR1_f, _nrrdDCTMRN_f, _nrrdDCTMR1_d, _nrrdDCTMRN_d
1211	};
1212	NrrdKernel *const
1213	nrrdKernelCatmullRomD = &_nrrdKernelCatmullRomD;
1214
1215	static NrrdKernel
1216	_nrrdKernelCatmullRomSupportDebugD = {
1217	"ctmrsupD",
1218	1, _nrrdCtmrSDSup, returnZero,
1219	_nrrdDCTMR1_f, _nrrdDCTMRN_f, _nrrdDCTMR1_d, _nrrdDCTMRN_d
1220	};
1221	NrrdKernel *const
1222	nrrdKernelCatmullRomSupportDebugD = &_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 - 60.5)x + 20.0 + 10 0.5 : (12 - 90.0 - 60.5)x - 6 + 40.0 + 20.5 )) _DDBCCUBIC(x, 0.0, 0.5)(x >= 2.0 ? 0 : (x >= 1.0 ? (-0.0 - 60.5)x + 20.0 + 10 0.5 : (12 - 90.0 - 60.5)x - 6 + 40.0 + 20.5 ))
1228
1229	static 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 - 60.5)x + 20.0 + 10 0.5 : (12 - 90.0 - 60.5)x - 6 + 40.0 + 2*0.5 ));
1234	}
1235
1236	static 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 - 60.5)x + 20.0 + 100.5 : (12 - 90.0 - 60.5)x - 6 + 40.0 + 2*0.5 ) )));
1241	}
1242
1243	static 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 - 60.5)t + 20.0 + 10 0.5 : (12 - 90.0 - 60.5)t - 6 + 40.0 + 2*0.5 ));
1252	}
1253	}
1254
1255	static 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 - 60.5)t + 20.0 + 100.5 : (12 - 90.0 - 60.5)t - 6 + 40.0 + 2*0.5 ) )));
1264	}
1265	}
1266
1267	static NrrdKernel
1268	_nrrdKernelCatmullRomDD = {
1269	"catmull-romDD",
1270	0, returnTwo, returnZero,
1271	_nrrdDDCTMR1_f, _nrrdDDCTMRN_f, _nrrdDDCTMR1_d, _nrrdDDCTMRN_d
1272	};
1273	NrrdKernel *const
1274	nrrdKernelCatmullRomDD = &_nrrdKernelCatmullRomDD;
1275
1276	static NrrdKernel
1277	_nrrdKernelCatmullRomSupportDebugDD = {
1278	"ctmrsupDD",
1279	1, _nrrdCtmrSDSup, returnZero,
1280	_nrrdDDCTMR1_f, _nrrdDDCTMRN_f, _nrrdDDCTMR1_d, _nrrdDDCTMRN_d
1281	};
1282	NrrdKernel *const
1283	nrrdKernelCatmullRomSupportDebugDD = &_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 - 6A + x(-10 + 25A + x(9 - 33A + x(-3.5 + 17A + x(0.5 - 3A)))) : 1 + xx(-3 + 6A + x((2.5 - 10A) + x(-0.5 + 4A)))))) \
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 - 6A + x(-10 + 25A + x(9 - 33*A \
1293	+ x(-3.5 + 17A + x(0.5 - 3A)))) \
1294	: 1 + xx(-3 + 6A + x((2.5 - 10A) + x(-0.5 + 4*A))))))
1295
1296	static double
1297	_nrrdA4Sup(const double *parm) {
1298	double S;
1299
1300	S = parm[0];
1301	return 3*S;
1302	}
1303
1304	static 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 - 6A + x(-10 + 25A + x(9 - 33A + x(-3.5 + 17A + x(0.5 - 3A)))) : 1 + xx(-3 + 6A + x((2.5 - 10A) + x(-0.5 + 4A))))))/S;
1312	}
1313
1314	static 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 - 6A + x(-10 + 25* A + x(9 - 33A + x(-3.5 + 17A + x(0.5 - 3A)))) : 1 + xx (-3 + 6A + x((2.5 - 10A) + x(-0.5 + 4*A))))))/S));
1321	}
1322
1323	static 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 - 6A + t(-10 + 25A + t(9 - 33A + t(-3.5 + 17A + t(0.5 - 3A)))) : 1 + tt(-3 + 6A + t((2.5 - 10A) + t(-0.5 + 4A))))))/S;
1334	}
1335	}
1336
1337	static 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 - 6A + t(-10 + 25* A + t(9 - 33A + t(-3.5 + 17A + t(0.5 - 3A)))) : 1 + tt (-3 + 6A + t((2.5 - 10A) + t(-0.5 + 4*A))))))/S));
1347	}
1348	}
1349
1350	static NrrdKernel
1351	_nrrdKernelA4 = {
1352	"Aquartic",
1353	2, _nrrdA4Sup, returnOne,
1354	_nrrdA41_f, _nrrdA4N_f, _nrrdA41_d, _nrrdA4N_d
1355	};
1356	NrrdKernel *const
1357	nrrdKernelAQuartic = &_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 + 25A + x(18 - 66A + x(-10.5 + 51 A + x(2 - 12A))) : x(-6 + 12A + x(7.5 - 30A + x(-2 + 16 *A)))))) \
1362	(x >= 3.0 ? 0 : \
1363	(x >= 2.0 \
1364	? A(81 + x(-90 + x(33 - 4x))) \
1365	: (x >= 1.0 \
1366	? -10 + 25A + x(18 - 66A + x(-10.5 + 51A + x(2 - 12*A))) \
1367	: x(-6 + 12A + x(7.5 - 30A + x(-2 + 16A))))))
1368
1369	static double
1370	_nrrdDA4Sup(const double *parm) {
1371	double S;
1372
1373	S = parm[0];
1374	return 3*S;
1375	}
1376
1377	static 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 + 25A + x(18 - 66A + x(-10.5 + 51 A + x(2 - 12A))) : x(-6 + 12A + x(7.5 - 30A + x(-2 + 16 A))))))/(SS);
1387	}
1388
1389	static 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)/(SS))((float)(sgn(x >= 3.0 ? 0 : (x >= 2.0 ? A(81 + x(-90 + x(33 - 4x))) : (x >= 1.0 ? -10 + 25A + x(18 - 66A + x(-10.5 + 51A + x(2 - 12A))) : x(-6 + 12A + x(7.5 - 30 A + x(-2 + 16A))))))/(S*S)));
1398	}
1399
1400	static 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 + 25A + t(18 - 66A + t(-10.5 + 51 A + t(2 - 12A))) : t(-6 + 12A + t(7.5 - 30A + t(-2 + 16 A))))))/(SS);
1412	}
1413	}
1414
1415	static 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)/(SS))((float)(sgn(t >= 3.0 ? 0 : (t >= 2.0 ? A(81 + t(-90 + t(33 - 4t))) : (t >= 1.0 ? -10 + 25A + t(18 - 66A + t(-10.5 + 51A + t(2 - 12A))) : t(-6 + 12A + t(7.5 - 30 A + t(-2 + 16A))))))/(S*S)));
1426	}
1427	}
1428
1429	static NrrdKernel
1430	_nrrdKernelDA4 = {
1431	"AquarticD",
1432	2, _nrrdDA4Sup, returnZero,
1433	_nrrdDA41_f, _nrrdDA4N_f, _nrrdDA41_d, _nrrdDA4N_d
1434	};
1435	NrrdKernel *const
1436	nrrdKernelAQuarticD = &_nrrdKernelDA4;
1437
1438	/* ------------------------------------------------------------ */
1439
1440	#define _DDAQUARTIC(x, A)(x >= 3.0 ? 0 : (x >= 2.0 ? A(-90 + x(66 - x12)) : ( x >= 1.0 ? 18 - 66A + x(-21 + 102A + x(6 - 36A)) : -6 + 12A + x(15 - 60A + 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 - 66A + x(-21 + 102A + x(6 - 36*A)) \
1446	: -6 + 12A + x(15 - 60A + x(-6 + 48*A)))))
1447
1448	static double
1449	_nrrdDDA4Sup(const double *parm) {
1450	double S;
1451
1452	S = parm[0];
1453	return 3*S;
1454	}
1455
1456	static 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 - x12)) : ( x >= 1.0 ? 18 - 66A + x(-21 + 102A + x(6 - 36A)) : -6 + 12A + x(15 - 60A + x(-6 + 48A)))))/(SS*S);
1464	}
1465
1466	static 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 - x12)) : ( x >= 1.0 ? 18 - 66A + x(-21 + 102A + x(6 - 36A)) : -6 + 12A + x(15 - 60A + x(-6 + 48A)))))/(SS*S);
1473	}
1474
1475	static 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 - t12)) : ( t >= 1.0 ? 18 - 66A + t(-21 + 102A + t(6 - 36A)) : -6 + 12A + t(15 - 60A + t(-6 + 48A)))))/(SS*S);
1486	}
1487	}
1488
1489	static 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 - t12)) : ( t >= 1.0 ? 18 - 66A + t(-21 + 102A + t(6 - 36A)) : -6 + 12A + t(15 - 60A + t(-6 + 48A)))))/(SS*S);
1499	}
1500	}
1501
1502	static NrrdKernel
1503	_nrrdKernelDDA4 = {
1504	"AquarticDD",
1505	2, _nrrdDDA4Sup, returnZero,
1506	_nrrdDDA41_f, _nrrdDDA4N_f, _nrrdDDA41_d, _nrrdDDA4N_d
1507	};
1508	NrrdKernel *const
1509	nrrdKernelAQuarticDD = &_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 + xx(-2 + x(2 + x(- 0.75 + x0.1))) : 0.7 + xx(-1 + xx(0.75 - x0.3)))) \
1533	(x >= 2.0 ? 0 : \
1534	(x >= 1.0 \
1535	? 0.8 + xx(-2 + x(2 + x(-0.75 + x*0.1))) \
1536	: 0.7 + xx(-1 + xx(0.75 - x*0.3))))
1537
1538	static 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 + xx(-2 + x(2 + x(- 0.75 + x0.1))) : 0.7 + xx(-1 + xx(0.75 - x0.3))));
1543	}
1544
1545	static 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 + xx(-2 + x* (2 + x(-0.75 + x0.1))) : 0.7 + xx(-1 + xx(0.75 - x*0.3) )))));
1550	}
1551
1552	static 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 + tt(-2 + t(2 + t(- 0.75 + t0.1))) : 0.7 + tt(-1 + tt(0.75 - t0.3))));
1561	}
1562	}
1563
1564	static 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 + tt(-2 + t* (2 + t(-0.75 + t0.1))) : 0.7 + tt(-1 + tt(0.75 - t*0.3) )))));
1573	}
1574	}
1575
1576	static NrrdKernel
1577	_c3quint = {
1578	"C3Quintic",
1579	0, returnTwo, returnOne,
1580	_c3quint1_f, _c3quintN_f, _c3quint1_d, _c3quintN_d
1581	};
1582	NrrdKernel *const
1583	nrrdKernelC3Quintic = &_c3quint;
1584
1585	/* ------------------------------------------------------------ */
1586
1587	#define _DC3QUINTIC(x)(x >= 2.0 ? 0 : (x >= 1.0 ? x(-4 + x(6 + x(-3 + x0.5 ))) : x(-2 + xx(3 - x1.5)))) \
1588	(x >= 2.0 ? 0 : \
1589	(x >= 1.0 \
1590	? x(-4 + x(6 + x(-3 + x0.5))) \
1591	: x(-2 + xx(3 - x1.5))))
1592
1593	static 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 + x0.5 ))) : x(-2 + xx(3 - x*1.5))));
1599	}
1600
1601	static 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 + x0.5))) : x(-2 + xx(3 - x1.5))))));
1607	}
1608
1609	static 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 + t0.5 ))) : t(-2 + tt(3 - t*1.5))));
1619	}
1620	}
1621
1622	static 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 + t0.5))) : t(-2 + tt(3 - t1.5))))));
1632	}
1633	}
1634
1635	static NrrdKernel
1636	_nrrdKernelDC3Quintic = {
1637	"C3QuinticD",
1638	0, returnTwo, returnZero,
1639	_Dc3quint1_f, _Dc3quintN_f, _Dc3quint1_d, _Dc3quintN_d
1640	};
1641	NrrdKernel *const
1642	nrrdKernelC3QuinticD = &_nrrdKernelDC3Quintic;
1643
1644	/* ------------------------------------------------------------ */
1645
1646	#define _DDC3QUINTIC(x)(x >= 2.0 ? 0 : (x >= 1.0 ? -4 + x(12 + x(-9 + x2)) : -2 + xx(9 - x6))) \
1647	(x >= 2.0 ? 0 : \
1648	(x >= 1.0 \
1649	? -4 + x(12 + x(-9 + x*2)) \
1650	: -2 + xx(9 - x*6)))
1651
1652	static 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 + x2)) : -2 + xx(9 - x6)));
1657	}
1658
1659	static 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 + x2)) : -2 + xx(9 - x6)))));
1664	}
1665
1666	static 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 + t2)) : -2 + tt(9 - t6)));
1675	}
1676	}
1677
1678	static 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 + t2)) : -2 + tt(9 - t6)))));
1687	}
1688	}
1689
1690	static NrrdKernel
1691	_DDc3quint = {
1692	"C3QuinticDD",
1693	0, returnTwo, returnZero,
1694	_DDc3quint1_f, _DDc3quintN_f, _DDc3quint1_d, _DDc3quintN_d
1695	};
1696	NrrdKernel *const
1697	nrrdKernelC3QuinticDD = &_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 - x3.0/32.0) )))) : 69.0/80.0 + xx(-23.0/16.0 + xx(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 - x3.0/32.0))))) \
1719	: 69.0/80.0 + xx(-23.0/16.0 + xx(19.0/16.0 + x*(-7.0/12.0 + x/16.0))) )))
1720
1721	static 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 - x3.0/32.0) )))) : 69.0/80.0 + xx(-23.0/16.0 + xx(19.0/16.0 + x(-7.0 /12.0 + x/16.0))) )));
1726	}
1727
1728	static 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 - x3.0/32.0))))) : 69.0/80.0 + xx(-23.0/16.0 + xx(19.0/16.0 + x(-7.0/12.0 + x/16.0))) )))));
1733	}
1734
1735	static 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 - t3.0/32.0) )))) : 69.0/80.0 + tt(-23.0/16.0 + tt(19.0/16.0 + t(-7.0 /12.0 + t/16.0))) )));
1744	}
1745	}
1746
1747	static 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 - t3.0/32.0))))) : 69.0/80.0 + tt(-23.0/16.0 + tt(19.0/16.0 + t(-7.0/12.0 + t/16.0))) )))));
1756	}
1757	}
1758
1759	static NrrdKernel
1760	_c4hex = {
1761	"C4Hexic",
1762	0, returnThree, returnOne,
1763	_c4hex1_f, _c4hexN_f, _c4hex1_d, _c4hexN_d
1764	};
1765	NrrdKernel *const
1766	nrrdKernelC4Hexic = &_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 + xx (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 + xx(19.0/4.0 + x(-35.0/12.0 + x*(3.0/8.0)))) )))
1778
1779	static 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 + xx (19.0/4.0 + x(-35.0/12.0 + x(3.0/8.0)))) )));
1785	}
1786
1787	static 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 + xx(19.0/4.0 + x(-35.0/12.0 + x*(3.0/8.0)))) )))));
1793	}
1794
1795	static 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 + tt (19.0/4.0 + t(-35.0/12.0 + t(3.0/8.0)))) )));
1805	}
1806	}
1807
1808	static 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 + tt(19.0/4.0 + t(-35.0/12.0 + t*(3.0/8.0)))) )))));
1818	}
1819	}
1820
1821	static NrrdKernel
1822	_nrrdKernelDC4hexic = {
1823	"C4HexicD",
1824	0, returnThree, returnZero,
1825	_Dc4hex1_f, _Dc4hexN_f, _Dc4hex1_d, _Dc4hexN_d
1826	};
1827	NrrdKernel *const
1828	nrrdKernelC4HexicD = &_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 + x15.0/16.0))) : (x >= 1.0 ? -341.0 /16.0 + x(60 + x(-441.0/8.0 + x(125.0/6.0 - x45.0/16.0))) : -23.0/8.0 + xx(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 + x15.0/16.0))) \
1837	: (x >= 1.0 \
1838	? -341.0/16.0 + x(60 + x(-441.0/8.0 + x(125.0/6.0 - x45.0/16.0))) \
1839	: -23.0/8.0 + xx(57.0/4.0 + x(-35.0/3.0 + x(15.0/8.0))) )))
1840
1841	static 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 + x15.0/16.0))) : (x >= 1.0 ? -341.0 /16.0 + x(60 + x(-441.0/8.0 + x(125.0/6.0 - x45.0/16.0))) : -23.0/8.0 + xx(57.0/4.0 + x(-35.0/3.0 + x(15.0/8.0))) ) ));
1846	}
1847
1848	static 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 + x15.0/16.0))) : (x >= 1.0 ? -341.0/16.0 + x(60 + x(-441.0/8.0 + x(125.0/6.0 - x45.0 /16.0))) : -23.0/8.0 + xx(57.0/4.0 + x(-35.0/3.0 + x(15.0 /8.0))) )))));
1853	}
1854
1855	static 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 + t15.0/16.0))) : (t >= 1.0 ? -341.0 /16.0 + t(60 + t(-441.0/8.0 + t(125.0/6.0 - t45.0/16.0))) : -23.0/8.0 + tt(57.0/4.0 + t(-35.0/3.0 + t(15.0/8.0))) ) ));
1864	}
1865	}
1866
1867	static 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 + t15.0/16.0))) : (t >= 1.0 ? -341.0/16.0 + t(60 + t(-441.0/8.0 + t(125.0/6.0 - t45.0 /16.0))) : -23.0/8.0 + tt(57.0/4.0 + t(-35.0/3.0 + t(15.0 /8.0))) )))));
1876	}
1877	}
1878
1879	static NrrdKernel
1880	_DDc4hex = {
1881	"C4HexicDD",
1882	0, returnThree, returnZero,
1883	_DDc4hex1_f, _DDc4hexN_f, _DDc4hex1_d, _DDc4hexN_d
1884	};
1885	NrrdKernel *const
1886	nrrdKernelC4HexicDD = &_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 + 15x/4)) : (x >= 1.0 ? 60 + x(-441.0/4.0 + x( 125.0/2.0 - 45x/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
1900	static 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 + 15x/4)) : (x >= 1.0 ? 60 + x(-441.0/4.0 + x( 125.0/2.0 - 45x/4)) : x(57.0/2.0 + x(-35 + 15x/2)) )));
1906	}
1907
1908	static 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 + 15x/4)) : (x >= 1.0 ? 60 + x(-441.0 /4.0 + x(125.0/2.0 - 45x/4)) : x(57.0/2.0 + x(-35 + 15*x/ 2)) )))));
1914	}
1915
1916	static 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 + 15t/4)) : (t >= 1.0 ? 60 + t(-441.0/4.0 + t( 125.0/2.0 - 45t/4)) : t(57.0/2.0 + t(-35 + 15t/2)) )));
1926	}
1927	}
1928
1929	static 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 + 15t/4)) : (t >= 1.0 ? 60 + t(-441.0 /4.0 + t(125.0/2.0 - 45t/4)) : t(57.0/2.0 + t(-35 + 15*t/ 2)) )))));
1939	}
1940	}
1941
1942	static NrrdKernel
1943	_nrrdKernelDDDC4hexic = {
1944	"C4HexicDDD",
1945	0, returnThree, returnZero,
1946	_DDDc4hex1_f, _DDDc4hexN_f, _DDDc4hex1_d, _DDDc4hexN_d
1947	};
1948	NrrdKernel *const
1949	nrrdKernelC4HexicDDD = &_nrrdKernelDDDC4hexic;
1950
1951
1952	/* ------------- approximate inverse of c4h ------------------------- */
1953	/* see comments around "_bspl3_ANI" in bsplKernel.c */
1954
1955	static 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
1970	static 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
1984	static 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
1993	static 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
2002	static 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
2014	static 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
2026	static 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	};
2033	NrrdKernel *const
2034	nrrdKernelC4HexicApproxInverse = &_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 + xx(-1.654320987654321 + xx(1.073045267489712 + xx(-0.44997427983539096 + x0.13978909465020575)))) (0.9379776601998824 + xx(-1.654320987654321 + xx(1.073045267489712 + xx(-0.44997427983539096 + x0.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 - x0.054552469135802466))))))) (0.04651675485008818 + x(-0.7377829218106996 + x(0.9699074074074074 + x(0.18531378600823045 + x(-0.7839506172839507 + x(0.2357253086419753 + x(0.12021604938271604 - x0.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 - x0.0013631687242798354))))))) (-0.01860670194003527 + x(0.14022633744855967 + x(-0.16296296296296298 + x(-0.09825102880658436 + x(0.28858024691358025 + x(-0.18858024691358025 + x(0.04405864197530864 - x0.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 + x0.004383450911228689))))))) (0.003101116990005879 + x(-0.014223251028806585 + x(0.02021604938271605 + x(0.003729423868312757 + x(-0.0411522633744856 + x(0.04714506172839506 + x(-0.023199588477366254 + x0.004383450911228689)))))))
2049	#define _C5SEPT(i, x)(0 == i ? (0.9379776601998824 + xx(-1.654320987654321 + xx (1.073045267489712 + xx(-0.44997427983539096 + x0.13978909465020575 )))) : (1 == i ? (0.04651675485008818 + x(-0.7377829218106996 + x(0.9699074074074074 + x(0.18531378600823045 + x(-0.7839506172839507 + x(0.2357253086419753 + x(0.12021604938271604 - x0.054552469135802466 ))))))) : (2 == i ? (-0.01860670194003527 + x(0.14022633744855967 + x(-0.16296296296296298 + x(-0.09825102880658436 + x(0.28858024691358025 + x(-0.18858024691358025 + x(0.04405864197530864 - x0.0013631687242798354 ))))))) : (3 == i ? (0.003101116990005879 + x(-0.014223251028806585 + x(0.02021604938271605 + x(0.003729423868312757 + x(-0.0411522633744856 + x(0.04714506172839506 + x(-0.023199588477366254 + x0.004383450911228689 ))))))) : 0.0)))) \
2050	(0 == i ? _C5SEPT0(x)(0.9379776601998824 + xx(-1.654320987654321 + xx(1.073045267489712 + xx(-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
2056	static 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 + xx(-1.654320987654321 + x* x(1.073045267489712 + xx(-0.44997427983539096 + x0.13978909465020575 )))) : (1 == xi ? (0.04651675485008818 + x(-0.7377829218106996 + x(0.9699074074074074 + x(0.18531378600823045 + x(-0.7839506172839507 + x(0.2357253086419753 + x(0.12021604938271604 - x0.054552469135802466 ))))))) : (2 == xi ? (-0.01860670194003527 + x(0.14022633744855967 + x(-0.16296296296296298 + x(-0.09825102880658436 + x(0.28858024691358025 + x(-0.18858024691358025 + x(0.04405864197530864 - x0.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
2066	static 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 + xx(-1.654320987654321 + xx(1.073045267489712 + xx(-0.44997427983539096 + x0.13978909465020575 )))) : (1 == xi ? (0.04651675485008818 + x(-0.7377829218106996 + x(0.9699074074074074 + x(0.18531378600823045 + x(-0.7839506172839507 + x(0.2357253086419753 + x(0.12021604938271604 - x0.054552469135802466 ))))))) : (2 == xi ? (-0.01860670194003527 + x(0.14022633744855967 + x(-0.16296296296296298 + x(-0.09825102880658436 + x(0.28858024691358025 + x(-0.18858024691358025 + x(0.04405864197530864 - x0.0013631687242798354 ))))))) : (3 == xi ? (0.003101116990005879 + x(-0.014223251028806585 + x(0.02021604938271605 + x(0.003729423868312757 + x(-0.0411522633744856 + x(0.04714506172839506 + x(-0.023199588477366254 + x0.004383450911228689 ))))))) : 0.0))))));
2074	}
2075
2076	static 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 + tt(-1.654320987654321 + t* t(1.073045267489712 + tt(-0.44997427983539096 + t0.13978909465020575 )))) : (1 == ti ? (0.04651675485008818 + t(-0.7377829218106996 + t(0.9699074074074074 + t(0.18531378600823045 + t(-0.7839506172839507 + t(0.2357253086419753 + t(0.12021604938271604 - t0.054552469135802466 ))))))) : (2 == ti ? (-0.01860670194003527 + t(0.14022633744855967 + t(-0.16296296296296298 + t(-0.09825102880658436 + t(0.28858024691358025 + t(-0.18858024691358025 + t(0.04405864197530864 - t0.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
2091	static 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 + tt(-1.654320987654321 + tt(1.073045267489712 + tt(-0.44997427983539096 + t0.13978909465020575 )))) : (1 == ti ? (0.04651675485008818 + t(-0.7377829218106996 + t(0.9699074074074074 + t(0.18531378600823045 + t(-0.7839506172839507 + t(0.2357253086419753 + t(0.12021604938271604 - t0.054552469135802466 ))))))) : (2 == ti ? (-0.01860670194003527 + t(0.14022633744855967 + t(-0.16296296296296298 + t(-0.09825102880658436 + t(0.28858024691358025 + t(-0.18858024691358025 + t(0.04405864197530864 - t0.0013631687242798354 ))))))) : (3 == ti ? (0.003101116990005879 + t(-0.014223251028806585 + t(0.02021604938271605 + t(0.003729423868312757 + t(-0.0411522633744856 + t(0.04714506172839506 + t(-0.023199588477366254 + t0.004383450911228689 ))))))) : 0.0))))));
2103	}
2104	}
2105
2106	static NrrdKernel
2107	_c5sept = {
2108	"C5Septic",
2109	0, returnFour, returnOne,
2110	_c5sept1_f, _c5septN_f, _c5sept1_d, _c5septN_d
2111	};
2112	NrrdKernel *const
2113	nrrdKernelC5Septic = &_c5sept;
2114
2115	#define _DC5SEPT0(x)(x(-3.308641975308642 + xx(4.292181069958848 + xx(-2.6998456790123457 + x0.9785236625514403)))) (x(-3.308641975308642 + xx(4.292181069958848 + xx(-2.6998456790123457 + x0.9785236625514403))))
2116	#define _DC5SEPT1(x)(-0.7377829218106996 + x(1.9398148148148149 + x(0.5559413580246914 + x(-3.1358024691358026 + x(1.1786265432098766 + x(0.7212962962962963 - x0.3818672839506173)))))) (-0.7377829218106996 + x(1.9398148148148149 + x(0.5559413580246914 + x(-3.1358024691358026 + x(1.1786265432098766 + x(0.7212962962962963 - x0.3818672839506173))))))
2117	#define _DC5SEPT2(x)(0.14022633744855967 + x(-0.32592592592592595 + x(-0.29475308641975306 + x(1.154320987654321 + x(-0.9429012345679012 + x(0.26435185185185184 - x0.009542181069958848)))))) (0.14022633744855967 + x(-0.32592592592592595 + x(-0.29475308641975306 + x(1.154320987654321 + x(-0.9429012345679012 + x(0.26435185185185184 - x0.009542181069958848))))))
2118	#define _DC5SEPT3(x)(-0.014223251028806585 + x(0.0404320987654321 + x(0.011188271604938271 + x(-0.1646090534979424 + x(0.2357253086419753 + x(-0.13919753086419753 + x0.03068415637860082)))))) (-0.014223251028806585 + x(0.0404320987654321 + x(0.011188271604938271 + x(-0.1646090534979424 + x(0.2357253086419753 + x(-0.13919753086419753 + x0.03068415637860082))))))
2119	#define _DC5SEPT(i, x)(0 == i ? (x(-3.308641975308642 + xx(4.292181069958848 + x x(-2.6998456790123457 + x0.9785236625514403)))) : (1 == i ? (-0.7377829218106996 + x(1.9398148148148149 + x(0.5559413580246914 + x(-3.1358024691358026 + x(1.1786265432098766 + x(0.7212962962962963 - x0.3818672839506173)))))) : (2 == i ? (0.14022633744855967 + x(-0.32592592592592595 + x(-0.29475308641975306 + x(1.154320987654321 + x(-0.9429012345679012 + x(0.26435185185185184 - x0.009542181069958848 )))))) : (3 == i ? (-0.014223251028806585 + x(0.0404320987654321 + x(0.011188271604938271 + x(-0.1646090534979424 + x(0.2357253086419753 + x(-0.13919753086419753 + x0.03068415637860082)))))) : 0.0 )))) \
2120	(0 == i ? _DC5SEPT0(x)(x(-3.308641975308642 + xx(4.292181069958848 + xx(-2.6998456790123457 + x0.9785236625514403)))) \
2121	: (1 == i ? _DC5SEPT1(x)(-0.7377829218106996 + x(1.9398148148148149 + x(0.5559413580246914 + x(-3.1358024691358026 + x(1.1786265432098766 + x(0.7212962962962963 - x0.3818672839506173)))))) \
2122	: (2 == i ? _DC5SEPT2(x)(0.14022633744855967 + x(-0.32592592592592595 + x(-0.29475308641975306 + x(1.154320987654321 + x(-0.9429012345679012 + x(0.26435185185185184 - x0.009542181069958848)))))) \
2123	: (3 == i ? _DC5SEPT3(x)(-0.014223251028806585 + x(0.0404320987654321 + x(0.011188271604938271 + x(-0.1646090534979424 + x(0.2357253086419753 + x(-0.13919753086419753 + x0.03068415637860082)))))) \
2124	: 0.0))))
2125
2126	static 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 + xx(4.292181069958848 + x x(-2.6998456790123457 + x0.9785236625514403)))) : (1 == xi ? (-0.7377829218106996 + x(1.9398148148148149 + x(0.5559413580246914 + x(-3.1358024691358026 + x(1.1786265432098766 + x(0.7212962962962963 - x0.3818672839506173)))))) : (2 == xi ? (0.14022633744855967 + x(-0.32592592592592595 + x(-0.29475308641975306 + x(1.154320987654321 + x(-0.9429012345679012 + x(0.26435185185185184 - x0.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
2137	static 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 + xx(4.292181069958848 + xx(-2.6998456790123457 + x0.9785236625514403)))) : (1 == xi ? (-0.7377829218106996 + x(1.9398148148148149 + x(0.5559413580246914 + x(-3.1358024691358026 + x(1.1786265432098766 + x(0.7212962962962963 - x0.3818672839506173)))))) : (2 == xi ? (0.14022633744855967 + x(-0.32592592592592595 + x(-0.29475308641975306 + x(1.154320987654321 + x(-0.9429012345679012 + x(0.26435185185185184 - x0.009542181069958848 )))))) : (3 == xi ? (-0.014223251028806585 + x(0.0404320987654321 + x(0.011188271604938271 + x(-0.1646090534979424 + x(0.2357253086419753 + x(-0.13919753086419753 + x0.03068415637860082)))))) : 0.0 ))))));
2146	}
2147
2148	static 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 + tt(4.292181069958848 + t t(-2.6998456790123457 + t0.9785236625514403)))) : (1 == ti ? (-0.7377829218106996 + t(1.9398148148148149 + t(0.5559413580246914 + t(-3.1358024691358026 + t(1.1786265432098766 + t(0.7212962962962963 - t0.3818672839506173)))))) : (2 == ti ? (0.14022633744855967 + t(-0.32592592592592595 + t(-0.29475308641975306 + t(1.154320987654321 + t(-0.9429012345679012 + t(0.26435185185185184 - t0.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
2164	static 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 + tt(4.292181069958848 + tt(-2.6998456790123457 + t0.9785236625514403)))) : (1 == ti ? (-0.7377829218106996 + t(1.9398148148148149 + t(0.5559413580246914 + t(-3.1358024691358026 + t(1.1786265432098766 + t(0.7212962962962963 - t0.3818672839506173)))))) : (2 == ti ? (0.14022633744855967 + t(-0.32592592592592595 + t(-0.29475308641975306 + t(1.154320987654321 + t(-0.9429012345679012 + t(0.26435185185185184 - t0.009542181069958848 )))))) : (3 == ti ? (-0.014223251028806585 + t(0.0404320987654321 + t(0.011188271604938271 + t(-0.1646090534979424 + t(0.2357253086419753 + t(-0.13919753086419753 + t0.03068415637860082)))))) : 0.0 ))))));
2177	}
2178	}
2179
2180	static NrrdKernel
2181	_dc5sept = {
2182	"C5SepticD",
2183	0, returnFour, returnZero,
2184	_dc5sept1_f, _dc5septN_f, _dc5sept1_d, _dc5septN_d
2185	};
2186	NrrdKernel *const
2187	nrrdKernelC5SepticD = &_dc5sept;
2188
2189	#define _DDC5SEPT0(x)(-3.308641975308642 + xx(12.876543209876543 + xx(-13.499228395061728 + x5.871141975308642))) (-3.308641975308642 + xx(12.876543209876543 + xx(-13.499228395061728 + x5.871141975308642)))
2190	#define _DDC5SEPT1(x)(1.9398148148148149 + x(1.1118827160493827 + x(-9.407407407407407 + x(4.714506172839506 + x(3.6064814814814814 - x2.2912037037037036 ))))) (1.9398148148148149 + x(1.1118827160493827 + x(-9.407407407407407 + x(4.714506172839506 + x(3.6064814814814814 - x2.2912037037037036)))))
2191	#define _DDC5SEPT2(x)(-0.32592592592592595 + x(-0.5895061728395061 + x(3.462962962962963 + x(-3.771604938271605 + x(1.3217592592592593 - x0.05725308641975309 ))))) (-0.32592592592592595 + x(-0.5895061728395061 + x(3.462962962962963 + x(-3.771604938271605 + x(1.3217592592592593 - x0.05725308641975309)))))
2192	#define _DDC5SEPT3(x)(0.0404320987654321 + x(0.022376543209876542 + x(-0.49382716049382713 + x(0.9429012345679012 + x(-0.6959876543209876 + x0.18410493827160493 ))))) (0.0404320987654321 + x(0.022376543209876542 + x(-0.49382716049382713 + x(0.9429012345679012 + x(-0.6959876543209876 + x0.18410493827160493)))))
2193	#define _DDC5SEPT(i, x)(0 == i ? (-3.308641975308642 + xx(12.876543209876543 + xx (-13.499228395061728 + x5.871141975308642))) : (1 == i ? (1.9398148148148149 + x(1.1118827160493827 + x(-9.407407407407407 + x(4.714506172839506 + x(3.6064814814814814 - x2.2912037037037036))))) : (2 == i ? (-0.32592592592592595 + x(-0.5895061728395061 + x(3.462962962962963 + x(-3.771604938271605 + x(1.3217592592592593 - x0.05725308641975309 ))))) : (3 == i ? (0.0404320987654321 + x(0.022376543209876542 + x(-0.49382716049382713 + x(0.9429012345679012 + x(-0.6959876543209876 + x0.18410493827160493))))) : 0.0)))) \
2194	(0 == i ? _DDC5SEPT0(x)(-3.308641975308642 + xx(12.876543209876543 + xx(-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
2200	static 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 + xx(12.876543209876543 + x* x(-13.499228395061728 + x5.871141975308642))) : (1 == xi ? ( 1.9398148148148149 + x(1.1118827160493827 + x(-9.407407407407407 + x(4.714506172839506 + x(3.6064814814814814 - x2.2912037037037036 ))))) : (2 == xi ? (-0.32592592592592595 + x(-0.5895061728395061 + x(3.462962962962963 + x(-3.771604938271605 + x(1.3217592592592593 - x0.05725308641975309))))) : (3 == xi ? (0.0404320987654321 + x(0.022376543209876542 + x(-0.49382716049382713 + x(0.9429012345679012 + x(-0.6959876543209876 + x*0.18410493827160493))))) : 0.0) )));
2208	}
2209
2210	static 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 + xx(12.876543209876543 + xx(-13.499228395061728 + x5.871141975308642))) : (1 == xi ? (1.9398148148148149 + x(1.1118827160493827 + x(-9.407407407407407 + x(4.714506172839506 + x(3.6064814814814814 - x2.2912037037037036 ))))) : (2 == xi ? (-0.32592592592592595 + x(-0.5895061728395061 + x(3.462962962962963 + x(-3.771604938271605 + x(1.3217592592592593 - x0.05725308641975309))))) : (3 == xi ? (0.0404320987654321 + x(0.022376543209876542 + x(-0.49382716049382713 + x(0.9429012345679012 + x(-0.6959876543209876 + x0.18410493827160493))))) : 0.0) )))));
2218	}
2219
2220	static 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 + tt(12.876543209876543 + t* t(-13.499228395061728 + t5.871141975308642))) : (1 == ti ? ( 1.9398148148148149 + t(1.1118827160493827 + t(-9.407407407407407 + t(4.714506172839506 + t(3.6064814814814814 - t2.2912037037037036 ))))) : (2 == ti ? (-0.32592592592592595 + t(-0.5895061728395061 + t(3.462962962962963 + t(-3.771604938271605 + t(1.3217592592592593 - t0.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
2235	static 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 + tt(12.876543209876543 + tt(-13.499228395061728 + t5.871141975308642))) : (1 == ti ? (1.9398148148148149 + t(1.1118827160493827 + t(-9.407407407407407 + t(4.714506172839506 + t(3.6064814814814814 - t2.2912037037037036 ))))) : (2 == ti ? (-0.32592592592592595 + t(-0.5895061728395061 + t(3.462962962962963 + t(-3.771604938271605 + t(1.3217592592592593 - t0.05725308641975309))))) : (3 == ti ? (0.0404320987654321 + t(0.022376543209876542 + t(-0.49382716049382713 + t(0.9429012345679012 + t(-0.6959876543209876 + t0.18410493827160493))))) : 0.0) )))));
2247	}
2248	}
2249
2250	static NrrdKernel
2251	_ddc5sept = {
2252	"C5SepticDD",
2253	0, returnFour, returnZero,
2254	_ddc5sept1_f, _ddc5septN_f, _ddc5sept1_d, _ddc5septN_d
2255	};
2256	NrrdKernel *const
2257	nrrdKernelC5SepticDD = &_ddc5sept;
2258
2259	#define _DDDC5SEPT0(x)(x(25.75308641975309 + xx(-53.99691358024691 + x29.35570987654321 ))) (x(25.75308641975309 + xx(-53.99691358024691 + x29.35570987654321)))
2260	#define _DDDC5SEPT1(x)(1.111882716049383 + x(-18.81481481481481 + x(14.14351851851852 + x(14.42592592592593 - x11.45601851851852)))) (1.111882716049383 + x(-18.81481481481481 + x(14.14351851851852 + x(14.42592592592593 - x11.45601851851852))))
2261	#define _DDDC5SEPT2(x)(-0.5895061728395062 + x(6.925925925925926 + x(-11.31481481481481 + x(5.287037037037037 - x0.2862654320987654)))) (-0.5895061728395062 + x(6.925925925925926 + x(-11.31481481481481 + x(5.287037037037037 - x0.2862654320987654))))
2262	#define _DDDC5SEPT3(x)(0.02237654320987654 + x(-0.9876543209876543 + x(2.828703703703704 + x(-2.783950617283951 + x0.9205246913580247)))) (0.02237654320987654 + x(-0.9876543209876543 + x(2.828703703703704 + x(-2.783950617283951 + x0.9205246913580247))))
2263	#define _DDDC5SEPT(i, x)(0 == i ? (x(25.75308641975309 + xx(-53.99691358024691 + x 29.35570987654321))) : (1 == i ? (1.111882716049383 + x(-18.81481481481481 + x(14.14351851851852 + x(14.42592592592593 - x11.45601851851852 )))) : (2 == i ? (-0.5895061728395062 + x(6.925925925925926 + x(-11.31481481481481 + x(5.287037037037037 - x0.2862654320987654 )))) : (3 == i ? (0.02237654320987654 + x(-0.9876543209876543 + x(2.828703703703704 + x(-2.783950617283951 + x0.9205246913580247 )))) : 0.0)))) \
2264	(0 == i ? _DDDC5SEPT0(x)(x(25.75308641975309 + xx(-53.99691358024691 + x29.35570987654321 ))) \
2265	: (1 == i ? _DDDC5SEPT1(x)(1.111882716049383 + x(-18.81481481481481 + x(14.14351851851852 + x(14.42592592592593 - x11.45601851851852)))) \
2266	: (2 == i ? _DDDC5SEPT2(x)(-0.5895061728395062 + x(6.925925925925926 + x(-11.31481481481481 + x(5.287037037037037 - x0.2862654320987654)))) \
2267	: (3 == i ? _DDDC5SEPT3(x)(0.02237654320987654 + x(-0.9876543209876543 + x(2.828703703703704 + x(-2.783950617283951 + x0.9205246913580247)))) \
2268	: 0.0))))
2269
2270	static 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 + xx(-53.99691358024691 + x 29.35570987654321))) : (1 == xi ? (1.111882716049383 + x(-18.81481481481481 + x(14.14351851851852 + x(14.42592592592593 - x11.45601851851852 )))) : (2 == xi ? (-0.5895061728395062 + x(6.925925925925926 + x(-11.31481481481481 + x(5.287037037037037 - x0.2862654320987654 )))) : (3 == xi ? (0.02237654320987654 + x(-0.9876543209876543 + x(2.828703703703704 + x(-2.783950617283951 + x*0.9205246913580247 )))) : 0.0))));
2279	}
2280
2281	static 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 + xx(-53.99691358024691 + x29.35570987654321))) : (1 == xi ? (1.111882716049383 + x (-18.81481481481481 + x(14.14351851851852 + x(14.42592592592593 - x11.45601851851852)))) : (2 == xi ? (-0.5895061728395062 + x(6.925925925925926 + x(-11.31481481481481 + x(5.287037037037037 - x0.2862654320987654)))) : (3 == xi ? (0.02237654320987654 + x(-0.9876543209876543 + x(2.828703703703704 + x(-2.783950617283951 + x0.9205246913580247)))) : 0.0))))));
2290	}
2291
2292	static 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 + tt(-53.99691358024691 + t 29.35570987654321))) : (1 == ti ? (1.111882716049383 + t(-18.81481481481481 + t(14.14351851851852 + t(14.42592592592593 - t11.45601851851852 )))) : (2 == ti ? (-0.5895061728395062 + t(6.925925925925926 + t(-11.31481481481481 + t(5.287037037037037 - t0.2862654320987654 )))) : (3 == ti ? (0.02237654320987654 + t(-0.9876543209876543 + t(2.828703703703704 + t(-2.783950617283951 + t*0.9205246913580247 )))) : 0.0))));
2305	}
2306	}
2307
2308	static 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 + tt(-53.99691358024691 + t29.35570987654321))) : (1 == ti ? (1.111882716049383 + t (-18.81481481481481 + t(14.14351851851852 + t(14.42592592592593 - t11.45601851851852)))) : (2 == ti ? (-0.5895061728395062 + t(6.925925925925926 + t(-11.31481481481481 + t(5.287037037037037 - t0.2862654320987654)))) : (3 == ti ? (0.02237654320987654 + t(-0.9876543209876543 + t(2.828703703703704 + t(-2.783950617283951 + t0.9205246913580247)))) : 0.0))))));
2321	}
2322	}
2323
2324	static NrrdKernel
2325	_dddc5sept = {
2326	"C5SepticDDD",
2327	0, returnFour, returnZero,
2328	_dddc5sept1_f, _dddc5septN_f, _dddc5sept1_d, _dddc5septN_d
2329	};
2330	NrrdKernel *const
2331	nrrdKernelC5SepticDDD = &_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
2339	static 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
2368	static double
2369	_c5sept_ANI_sup(const double *parm) {
2370	AIR_UNUSED(parm)(void)(parm);
2371	return C5SEPT_AI_LEN26 + 0.5;
2372	}
2373
2374	static 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
2388	static 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
2398	static 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
2408	static 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
2421	static 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
2434	static 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	};
2441	NrrdKernel *const
2442	nrrdKernelC5SepticApproxInverse = &_nrrdKernelC5SepticApproxInverse;
2443
2444	/* ------------------------------------------------------------ */
2445
2446	#define _GAUSS(x, sig, cut)( x >= sigcut ? 0 : exp(-xx/(2.0sigsig))/(sig*2.50662827463100050241 )) ( \
2447	x >= sig*cut ? 0 \
2448	: exp(-xx/(2.0sigsig))/(sig2.50662827463100050241))
2449
2450	static double
2451	_nrrdGInt(const double *parm) {
2452	double cut;
2453
2454	cut = parm[1];
2455	return airErf(cut/sqrt(2.0));
2456	}
2457
2458	static 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
2467	static 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 >= sigcut ? 0 : exp(-xx/(2.0sigsig))/(sig*2.50662827463100050241 ));
2475	}
2476
2477	static 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 >= sigcut ? 0 : exp(-xx/(2.0sigsig))/(sig *2.50662827463100050241))));
2485	}
2486
2487	static 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 >= sigcut ? 0 : exp(-tt/(2.0sigsig))/(sig*2.50662827463100050241 ));
2498	}
2499	}
2500
2501	static 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 >= sigcut ? 0 : exp(-tt/(2.0sigsig))/(sig *2.50662827463100050241))));
2512	}
2513	}
2514
2515	static NrrdKernel
2516	_nrrdKernelG = {
2517	"gauss",
2518	2, _nrrdGSup, _nrrdGInt,
2519	_nrrdG1_f, _nrrdGN_f, _nrrdG1_d, _nrrdGN_d
2520	};
2521	NrrdKernel *const
2522	nrrdKernelGaussian = &_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
2544	static 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
2556	static 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
2566	static 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
2580	static 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
2590	static 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
2604	static 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
2617	static NrrdKernel
2618	_nrrdKernelDiscreteGaussian = {
2619	"discretegauss", 2,
2620	_nrrdDiscGaussianSup, _nrrdDiscGaussianInt,
2621	_nrrdDiscGaussian1_f, _nrrdDiscGaussianN_f,
2622	_nrrdDiscGaussian1_d, _nrrdDiscGaussianN_d
2623	};
2624	NrrdKernel *const
2625	nrrdKernelDiscreteGaussian = &_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) */
2636	const double nrrdKernelDiscreteGaussianGoodSigmaMax = 6.0;
2637
2638	/* ------------------------------------------------------------ */
2639
2640	#define _DGAUSS(x, sig, cut)( x >= sigcut ? 0 : -exp(-xx/(2.0sigsig))x/(sigsigsig 2.50662827463100050241)) ( \
2641	x >= sig*cut ? 0 \
2642	: -exp(-xx/(2.0sigsig))x/(sigsigsig*2.50662827463100050241))
2643
2644	static double
2645	_nrrdDGInt(const double *parm) {
2646	AIR_UNUSED(parm)(void)(parm);
2647	return 0;
2648	}
2649
2650	static 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
2659	static 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 >= sigcut ? 0 : -exp(-xx/(2.0sigsig))x/(sigsigsig *2.50662827463100050241));
2668	}
2669
2670	static 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 >= sigcut ? 0 : -exp(-xx/(2.0sigsig)) x/(sigsigsig2.50662827463100050241))));
2679	}
2680
2681	static 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 >= sigcut ? 0 : -exp(-tt/(2.0sigsig))t/(sigsigsig *2.50662827463100050241));
2693	}
2694	}
2695
2696	static 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 >= sigcut ? 0 : -exp(-tt/(2.0sigsig)) t/(sigsigsig2.50662827463100050241))));
2708	}
2709	}
2710
2711	static NrrdKernel
2712	_nrrdKernelDG = {
2713	"gaussD",
2714	2, _nrrdDGSup, _nrrdDGInt,
2715	_nrrdDG1_f, _nrrdDGN_f, _nrrdDG1_d, _nrrdDGN_d
2716	};
2717	NrrdKernel *const
2718	nrrdKernelGaussianD = &_nrrdKernelDG;
2719
2720	/* ------------------------------------------------------------ */
2721
2722	#define _DDGAUSS(x, sig, cut)( x >= sigcut ? 0 : exp(-xx/(2.0sigsig))(xx-sigsig) / (sigsigsigsigsig2.50662827463100050241)) ( \
2723	x >= sig*cut ? 0 \
2724	: exp(-xx/(2.0sigsig))(xx-sigsig) / \
2725	(sigsigsigsigsig*2.50662827463100050241))
2726
2727	static double
2728	_nrrdDDGInt(const double *parm) {
2729	double sig, cut;
2730
2731	sig = parm[0];
2732	cut = parm[1];
2733	return -0.79788456080286535587cutexp(-cutcut/2)/(sigsig);
2734	}
2735
2736	static 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
2745	static 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 >= sigcut ? 0 : exp(-xx/(2.0sigsig))(xx-sigsig) / (sigsigsigsigsig2.50662827463100050241));
2753	}
2754
2755	static 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 >= sigcut ? 0 : exp(-xx/(2.0sigsig))(xx -sigsig) / (sigsigsigsigsig2.50662827463100050241))));
2763	}
2764
2765	static 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 >= sigcut ? 0 : exp(-tt/(2.0sigsig))(tt-sigsig) / (sigsigsigsigsig2.50662827463100050241));
2776	}
2777	}
2778
2779	static 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 >= sigcut ? 0 : exp(-tt/(2.0sigsig))(tt -sigsig) / (sigsigsigsigsig2.50662827463100050241))));
2790	}
2791	}
2792
2793	static NrrdKernel
2794	_nrrdKernelDDG = {
2795	"gaussDD",
2796	2, _nrrdDDGSup, _nrrdDDGInt,
2797	_nrrdDDG1_f, _nrrdDDGN_f, _nrrdDDG1_d, _nrrdDDGN_d
2798	};
2799	NrrdKernel *const
2800	nrrdKernelGaussianDD = &_nrrdKernelDDG;
2801
2802
2803	/* ------------------------------------------------------------ */
2804
2805	NrrdKernel *
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 */
2956	int
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
2971	int
2972	nrrdKernelParse(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
3140	int
3141	nrrdKernelSpecParse(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	*/
3162	int
3163	nrrdKernelSpecSprint(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
3221	int
3222	nrrdKernelSprint(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	*/
3239	int
3240	nrrdKernelCompare(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)) {
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) {
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) {
3262	/* actually, we're done: no parms and kernels are equal */
3263	*differ = 0;
3264	return 0;
3265	}
3266	if (!(parmA && parmB)) {
3267	biffAddf(NRRDnrrdBiffKey, "%s: kernel %s needs %u parms but got NULL parm vectors",
3268	me, kernA->name ? kernA->name : "(unnamed)", pnum);
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	*/
3290	int
3291	nrrdKernelSpecCompare(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	*/
3359	int
3360	nrrdKernelCheck(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;
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) {
3382	biffAddf(NRRDnrrdBiffKey, "%s: got NULL kernel", me);
3383	airMopError(mop); return 1;
3384	}
3385	if (!(evalNum > 20)) {
3386	biffAddf(NRRDnrrdBiffKey, "%s: need evalNum > 20", me);
3387	airMopError(mop); return 1;
3388	}
3389	if (!(kern->support && kern->integral
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)
3404	\|\| nrrdKernelSpecSprint(kspstr, kspA)) {
3405	biffAddf(NRRDnrrdBiffKey, "%s: trouble", me);
3406	airMopError(mop); return 1;
3407	}
3408	if (strcmp(kstr, kspstr)) {
3409	biffAddf(NRRDnrrdBiffKey, "%s: sprinted kernel \|%s\| != kspec \|%s\|", me, kstr, kspstr);
3410	airMopError(mop); return 1;
3411	}
3412	if (nrrdKernelParse(&parsedkern, parsedparm, kstr)
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,
3419	&differ, explain)) {
3420	biffAddf(NRRDnrrdBiffKey, "%s: trouble comparing kern/parm pairs", me);
3421	airMopError(mop); return 1;
3422	}
3423	if (differ) {
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 = 2supp/(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
3592	int
3593	nrrdKernelParm0IsScale(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	}