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 "teem/gage.h" |
25 |
|
|
#include <testDataPath.h> |
26 |
|
|
|
27 |
|
|
/* |
28 |
|
|
** Tests: |
29 |
|
|
** |
30 |
|
|
*/ |
31 |
|
|
|
32 |
|
|
#define INTERP_KERN_NUM 4 |
33 |
|
|
#define BLUR_KERN_NUM 5 |
34 |
|
|
|
35 |
|
|
int |
36 |
|
|
main(int argc, const char **argv) { |
37 |
|
|
const char *me; |
38 |
|
|
Nrrd *nscl; |
39 |
|
|
double *dscl; |
40 |
|
|
airArray *mop; |
41 |
|
|
char *fullname; |
42 |
|
2 |
gageContext *igctx[INTERP_KERN_NUM], *bgctx[BLUR_KERN_NUM]; |
43 |
|
1 |
const NrrdKernel *ikern[INTERP_KERN_NUM]; |
44 |
|
4 |
double ikparm[INTERP_KERN_NUM][NRRD_KERNEL_PARMS_NUM] = { |
45 |
|
1 |
{1.0}, |
46 |
|
1 |
{1.0}, |
47 |
|
1 |
{1.0, 0.0, 0.5}, |
48 |
|
1 |
{AIR_NAN}, |
49 |
|
|
}; |
50 |
|
1 |
const NrrdKernel *bkern[BLUR_KERN_NUM]; |
51 |
|
1 |
const NrrdKernel *bkernD[BLUR_KERN_NUM]; |
52 |
|
1 |
const NrrdKernel *bkernDD[BLUR_KERN_NUM]; |
53 |
|
5 |
double bkparm[BLUR_KERN_NUM][NRRD_KERNEL_PARMS_NUM] = { |
54 |
|
1 |
{1.0}, |
55 |
|
1 |
{AIR_NAN}, |
56 |
|
1 |
{AIR_NAN}, |
57 |
|
1 |
{2.0, 1.0, 0.0}, |
58 |
|
1 |
{1.2, 5.0}, |
59 |
|
|
}; |
60 |
|
1 |
const double *ivalAns[INTERP_KERN_NUM], *bvalAns[BLUR_KERN_NUM], |
61 |
|
|
*bgrdAns[BLUR_KERN_NUM], *bhesAns[BLUR_KERN_NUM]; |
62 |
|
|
int E; |
63 |
|
|
unsigned int sx, sy, sz, ki; |
64 |
|
|
|
65 |
|
|
/* C89 doesn't allow setting these in array declaration */ |
66 |
|
1 |
ikern[0] = nrrdKernelBox; |
67 |
|
1 |
ikern[1] = nrrdKernelTent; |
68 |
|
1 |
ikern[2] = nrrdKernelBCCubic; |
69 |
|
1 |
ikern[3] = nrrdKernelCatmullRom; |
70 |
|
1 |
bkern[0] = nrrdKernelTent; |
71 |
|
1 |
bkern[1] = nrrdKernelBSpline3; |
72 |
|
1 |
bkern[2] = nrrdKernelBSpline5; |
73 |
|
1 |
bkern[3] = nrrdKernelBCCubic; |
74 |
|
1 |
bkern[4] = nrrdKernelGaussian; |
75 |
|
1 |
bkernD[0] = nrrdKernelForwDiff; |
76 |
|
1 |
bkernD[1] = nrrdKernelBSpline3D; |
77 |
|
1 |
bkernD[2] = nrrdKernelBSpline5D; |
78 |
|
1 |
bkernD[3] = nrrdKernelBCCubicD; |
79 |
|
1 |
bkernD[4] = nrrdKernelGaussianD; |
80 |
|
1 |
bkernDD[0] = nrrdKernelZero; |
81 |
|
1 |
bkernDD[1] = nrrdKernelBSpline3DD; |
82 |
|
1 |
bkernDD[2] = nrrdKernelBSpline5DD; |
83 |
|
1 |
bkernDD[3] = nrrdKernelBCCubicDD; |
84 |
|
1 |
bkernDD[4] = nrrdKernelGaussianDD; |
85 |
|
|
|
86 |
|
|
AIR_UNUSED(argc); |
87 |
|
1 |
me = argv[0]; |
88 |
|
1 |
mop = airMopNew(); |
89 |
|
|
|
90 |
|
1 |
nscl = nrrdNew(); |
91 |
|
1 |
airMopAdd(mop, nscl, (airMopper)nrrdNuke, airMopAlways); |
92 |
|
1 |
fullname = testDataPathPrefix("fmob-c4h.nrrd"); |
93 |
|
1 |
airMopAdd(mop, fullname, airFree, airMopAlways); |
94 |
✗✓ |
1 |
if (nrrdLoad(nscl, fullname, NULL)) { |
95 |
|
|
char *err; |
96 |
|
|
airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways); |
97 |
|
|
fprintf(stderr, "%s: trouble reading data \"%s\":\n%s", |
98 |
|
|
me, fullname, err); |
99 |
|
|
airMopError(mop); return 1; |
100 |
|
|
} |
101 |
|
|
/* make sure its a double-type volume (assumed below) */ |
102 |
✗✓ |
1 |
if (nrrdTypeDouble != nscl->type) { |
103 |
|
|
fprintf(stderr, "%s: volume type %s != expected type %s\n", me, |
104 |
|
|
airEnumStr(nrrdType, nscl->type), |
105 |
|
|
airEnumStr(nrrdType, nrrdTypeDouble)); |
106 |
|
|
airMopError(mop); return 1; |
107 |
|
|
} |
108 |
|
1 |
dscl = AIR_CAST(double *, nscl->data); |
109 |
|
1 |
sx = AIR_CAST(unsigned int, nscl->axis[0].size); |
110 |
|
1 |
sy = AIR_CAST(unsigned int, nscl->axis[1].size); |
111 |
|
1 |
sz = AIR_CAST(unsigned int, nscl->axis[2].size); |
112 |
|
|
|
113 |
✓✓ |
10 |
for (ki=0; ki<INTERP_KERN_NUM; ki++) { |
114 |
|
|
gagePerVolume *gpvl; |
115 |
|
4 |
igctx[ki] = gageContextNew(); |
116 |
|
4 |
airMopAdd(mop, igctx[ki], (airMopper)gageContextNix, airMopAlways); |
117 |
|
4 |
gageParmSet(igctx[ki], gageParmRenormalize, AIR_FALSE); |
118 |
|
4 |
gageParmSet(igctx[ki], gageParmCheckIntegrals, AIR_TRUE); |
119 |
|
4 |
gageParmSet(igctx[ki], gageParmOrientationFromSpacing, AIR_FALSE); |
120 |
|
|
E = 0; |
121 |
✓✗ |
8 |
if (!E) E |= !(gpvl = gagePerVolumeNew(igctx[ki], nscl, gageKindScl)); |
122 |
✓✗ |
12 |
if (!E) E |= gageKernelSet(igctx[ki], gageKernel00, |
123 |
|
4 |
ikern[ki], ikparm[ki]); |
124 |
✓✗ |
8 |
if (!E) E |= gagePerVolumeAttach(igctx[ki], gpvl); |
125 |
✓✗ |
8 |
if (!E) E |= gageQueryItemOn(igctx[ki], gpvl, gageSclValue); |
126 |
✓✗ |
8 |
if (!E) E |= gageUpdate(igctx[ki]); |
127 |
✗✓ |
4 |
if (E) { |
128 |
|
|
char *err; |
129 |
|
|
airMopAdd(mop, err = biffGetDone(GAGE), airFree, airMopAlways); |
130 |
|
|
fprintf(stderr, "%s: trouble %s set-up:\n%s\n", me, |
131 |
|
|
ikern[ki]->name, err); |
132 |
|
|
airMopError(mop); return 1; |
133 |
|
|
} |
134 |
|
4 |
ivalAns[ki] = gageAnswerPointer(igctx[ki], gpvl, gageSclValue); |
135 |
|
4 |
} |
136 |
|
|
|
137 |
|
|
/* traverse all samples of volume, probing with the interpolating |
138 |
|
|
kernels, make sure we recover the original values */ |
139 |
|
|
{ |
140 |
|
|
unsigned int xi, yi, zi; |
141 |
|
1 |
double pval[INTERP_KERN_NUM], err, rval; |
142 |
|
|
int pret; |
143 |
✓✓ |
34 |
for (zi=0; zi<sz; zi++) { |
144 |
✓✓ |
1312 |
for (yi=0; yi<sy; yi++) { |
145 |
✓✓ |
42240 |
for (xi=0; xi<sx; xi++) { |
146 |
|
20480 |
rval = dscl[xi + sx*(yi + sy*zi)]; |
147 |
✓✓ |
204800 |
for (ki=0; ki<INTERP_KERN_NUM; ki++) { |
148 |
|
81920 |
pret = gageProbeSpace(igctx[ki], xi, yi, zi, |
149 |
|
|
AIR_TRUE /* indexSpace */, |
150 |
|
|
AIR_FALSE /* clamp */); |
151 |
✗✓ |
81920 |
if (pret) { |
152 |
|
|
fprintf(stderr, "%s: %s probe error(%d): %s\n", me, |
153 |
|
|
ikern[ki]->name, igctx[ki]->errNum, igctx[ki]->errStr); |
154 |
|
|
|
155 |
|
|
airMopError(mop); return 1; |
156 |
|
|
} |
157 |
|
81920 |
pval[ki] = *ivalAns[ki]; |
158 |
|
81920 |
err = AIR_ABS(rval - pval[ki]); |
159 |
✗✓ |
81920 |
if (err) { |
160 |
|
|
fprintf(stderr, "%s: interp's [%u,%u,%u] %s probe %f " |
161 |
|
|
"!= true %f (err %f)\n", me, xi, yi, zi, |
162 |
|
|
ikern[ki]->name, pval[ki], rval, err); |
163 |
|
|
airMopError(mop); return 1; |
164 |
|
|
} |
165 |
|
|
} |
166 |
|
|
} |
167 |
|
|
} |
168 |
|
|
} |
169 |
✓✗ |
2 |
} |
170 |
|
|
|
171 |
|
|
/* set up contexts for non-interpolating (blurring) kernels, |
172 |
|
|
and their first and second derivatives */ |
173 |
✓✓ |
12 |
for (ki=0; ki<BLUR_KERN_NUM; ki++) { |
174 |
|
|
gagePerVolume *gpvl; |
175 |
|
5 |
bgctx[ki] = gageContextNew(); |
176 |
|
5 |
airMopAdd(mop, bgctx[ki], (airMopper)gageContextNix, airMopAlways); |
177 |
|
5 |
gageParmSet(bgctx[ki], gageParmRenormalize, AIR_TRUE); |
178 |
|
5 |
gageParmSet(bgctx[ki], gageParmCheckIntegrals, AIR_TRUE); |
179 |
|
5 |
gageParmSet(bgctx[ki], gageParmOrientationFromSpacing, AIR_FALSE); |
180 |
|
|
E = 0; |
181 |
✓✗ |
10 |
if (!E) E |= !(gpvl = gagePerVolumeNew(bgctx[ki], nscl, gageKindScl)); |
182 |
✓✗ |
15 |
if (!E) E |= gageKernelSet(bgctx[ki], gageKernel00, |
183 |
|
5 |
bkern[ki], bkparm[ki]); |
184 |
✓✗ |
15 |
if (!E) E |= gageKernelSet(bgctx[ki], gageKernel11, |
185 |
|
5 |
bkernD[ki], bkparm[ki]); |
186 |
✓✗ |
15 |
if (!E) E |= gageKernelSet(bgctx[ki], gageKernel22, |
187 |
|
5 |
bkernDD[ki], bkparm[ki]); |
188 |
✓✗ |
10 |
if (!E) E |= gagePerVolumeAttach(bgctx[ki], gpvl); |
189 |
✓✗ |
10 |
if (!E) E |= gageQueryItemOn(bgctx[ki], gpvl, gageSclValue); |
190 |
✓✗ |
10 |
if (!E) E |= gageQueryItemOn(bgctx[ki], gpvl, gageSclGradVec); |
191 |
✓✗ |
10 |
if (!E) E |= gageQueryItemOn(bgctx[ki], gpvl, gageSclHessian); |
192 |
✓✗ |
10 |
if (!E) E |= gageUpdate(bgctx[ki]); |
193 |
✗✓ |
5 |
if (E) { |
194 |
|
|
char *err; |
195 |
|
|
airMopAdd(mop, err = biffGetDone(GAGE), airFree, airMopAlways); |
196 |
|
|
fprintf(stderr, "%s: trouble %s set-up:\n%s\n", me, |
197 |
|
|
bkern[ki]->name, err); |
198 |
|
|
airMopError(mop); return 1; |
199 |
|
|
} |
200 |
|
5 |
fprintf(stderr, "%s radius = %u\n", bkern[ki]->name, bgctx[ki]->radius); |
201 |
|
5 |
bvalAns[ki] = gageAnswerPointer(bgctx[ki], gpvl, gageSclValue); |
202 |
|
5 |
bgrdAns[ki] = gageAnswerPointer(bgctx[ki], gpvl, gageSclGradVec); |
203 |
|
5 |
bhesAns[ki] = gageAnswerPointer(bgctx[ki], gpvl, gageSclHessian); |
204 |
|
5 |
} |
205 |
|
|
|
206 |
|
|
{ |
207 |
|
|
#define POS_NUM 12 |
208 |
|
1 |
double xp[POS_NUM], yp[POS_NUM], zp[POS_NUM], |
209 |
|
|
pos[POS_NUM*POS_NUM*POS_NUM][3], *prbd, |
210 |
|
1 |
offs[POS_NUM/2] = {0, 1.22222, 2.444444, 3.777777, 5.88888, 7.55555}; |
211 |
|
|
Nrrd *nprbd, *ncorr; |
212 |
|
|
unsigned int ii, jj, kk, qlen = 1 + 3 + 9; |
213 |
|
1 |
char *corrfn, explain[AIR_STRLEN_LARGE]; |
214 |
|
1 |
int pret, differ; |
215 |
|
|
|
216 |
|
1 |
corrfn = testDataPathPrefix("test/probeSclAns.nrrd"); |
217 |
|
1 |
airMopAdd(mop, corrfn, airFree, airMopAlways); |
218 |
|
1 |
ncorr = nrrdNew(); |
219 |
|
1 |
airMopAdd(mop, ncorr, (airMopper)nrrdNuke, airMopAlways); |
220 |
✗✓ |
1 |
if (nrrdLoad(ncorr, corrfn, NULL)) { |
221 |
|
|
char *err; |
222 |
|
|
airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways); |
223 |
|
|
fprintf(stderr, "%s: trouble reading data \"%s\":\n%s", |
224 |
|
|
me, corrfn, err); |
225 |
|
|
airMopError(mop); return 1; |
226 |
|
|
} |
227 |
✓✓ |
14 |
for (ii=0; ii<POS_NUM/2; ii++) { |
228 |
|
6 |
xp[ii] = yp[ii] = zp[ii] = offs[ii]; |
229 |
|
6 |
xp[POS_NUM-1-ii] = AIR_CAST(double, sx)-1.0-offs[ii]; |
230 |
|
6 |
yp[POS_NUM-1-ii] = AIR_CAST(double, sy)-1.0-offs[ii]; |
231 |
|
6 |
zp[POS_NUM-1-ii] = AIR_CAST(double, sz)-1.0-offs[ii]; |
232 |
|
|
} |
233 |
✓✓ |
26 |
for (kk=0; kk<POS_NUM; kk++) { |
234 |
✓✓ |
312 |
for (jj=0; jj<POS_NUM; jj++) { |
235 |
✓✓ |
3744 |
for (ii=0; ii<POS_NUM; ii++) { |
236 |
|
1728 |
ELL_3V_SET(pos[ii + POS_NUM*(jj + POS_NUM*kk)], |
237 |
|
|
xp[ii], yp[jj], zp[kk]); |
238 |
|
|
} |
239 |
|
|
} |
240 |
|
|
} |
241 |
|
1 |
nprbd = nrrdNew(); |
242 |
|
1 |
airMopAdd(mop, nprbd, (airMopper)nrrdNuke, airMopAlways); |
243 |
✗✓ |
1 |
if (nrrdMaybeAlloc_va(nprbd, nrrdTypeDouble, 3, |
244 |
|
|
AIR_CAST(size_t, qlen), |
245 |
|
|
AIR_CAST(size_t, BLUR_KERN_NUM), |
246 |
|
|
AIR_CAST(size_t, POS_NUM*POS_NUM*POS_NUM))) { |
247 |
|
|
char *err; |
248 |
|
|
airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways); |
249 |
|
|
fprintf(stderr, "%s: trouble setting up prbd:\n%s", me, err); |
250 |
|
|
airMopError(mop); return 1; |
251 |
|
|
} |
252 |
|
1 |
prbd = AIR_CAST(double *, nprbd->data); |
253 |
✓✓ |
3458 |
for (ii=0; ii<POS_NUM*POS_NUM*POS_NUM; ii++) { |
254 |
✓✓ |
20736 |
for (ki=0; ki<BLUR_KERN_NUM; ki++) { |
255 |
|
8640 |
pret = gageProbeSpace(bgctx[ki], pos[ii][0], pos[ii][1], pos[ii][2], |
256 |
|
|
AIR_TRUE /* indexSpace */, |
257 |
|
|
AIR_FALSE /* clamp */); |
258 |
✗✓ |
8640 |
if (pret) { |
259 |
|
|
fprintf(stderr, "%s: %s probe error(%d): %s\n", me, |
260 |
|
|
bkern[ki]->name, bgctx[ki]->errNum, bgctx[ki]->errStr); |
261 |
|
|
airMopError(mop); return 1; |
262 |
|
|
} |
263 |
|
8640 |
prbd[0 + qlen*(ki + BLUR_KERN_NUM*(ii))] = bvalAns[ki][0]; |
264 |
|
8640 |
ELL_3V_COPY(prbd + 1 + qlen*(ki + BLUR_KERN_NUM*(ii)), bgrdAns[ki]); |
265 |
|
8640 |
ELL_9V_COPY(prbd + 4 + qlen*(ki + BLUR_KERN_NUM*(ii)), bhesAns[ki]); |
266 |
|
|
} |
267 |
|
|
} |
268 |
|
|
/* HEY: weirdly, so far its only on Windows (and more than 10 times worse |
269 |
|
|
on Cygwin) this epsilon needs to be larger than zero, and only for the |
270 |
|
|
radius 6 Gaussian? */ |
271 |
✗✓ |
1 |
if (nrrdCompare(ncorr, nprbd, AIR_FALSE /* onlyData */, |
272 |
|
1 |
8.0e-14 /* epsilon */, &differ, explain)) { |
273 |
|
|
char *err; |
274 |
|
|
airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways); |
275 |
|
|
fprintf(stderr, "%s: trouble comparing:\n%s", me, err); |
276 |
|
|
airMopError(mop); return 1; |
277 |
|
|
} |
278 |
✗✓✗✓
|
2 |
if (differ) { |
279 |
|
1 |
fprintf(stderr, "%s: probed values not correct: %s\n", me, explain); |
280 |
|
|
airMopError(mop); return 1; |
281 |
|
|
} else { |
282 |
|
1 |
fprintf(stderr, "all good\n"); |
283 |
|
|
} |
284 |
✓✗ |
2 |
} |
285 |
|
|
|
286 |
|
1 |
airMopOkay(mop); |
287 |
|
1 |
return 0; |
288 |
|
1 |
} |