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 |
|
|
|
26 |
|
|
/* |
27 |
|
|
** Tests: |
28 |
|
|
** nrrdQuantize (to 8-bits), nrrdUnquantize |
29 |
|
|
** nrrdArithBinaryOp (with nrrdBinaryOpSubtract) |
30 |
|
|
*/ |
31 |
|
|
|
32 |
|
|
#define KERN_SIZE_MAX 10 |
33 |
|
|
#define PROBE_NUM 300 |
34 |
|
|
|
35 |
|
|
static void |
36 |
|
|
errPrefix(char *dst, |
37 |
|
|
int typi, unsigned int supi, unsigned int prbi, |
38 |
|
|
unsigned int probePass, double dxi, double dyi, double dzi, |
39 |
|
|
unsigned int xi, unsigned int yi, unsigned int zi) { |
40 |
|
|
sprintf(dst, "%s[%s][%u] #%u pp %u: (%g,%g,%g)->(%u,%u,%u): ", |
41 |
|
|
"probeMulti", |
42 |
|
|
airEnumStr(nrrdType, typi), supi, prbi, probePass, |
43 |
|
|
dxi, dyi, dzi, xi, yi, zi); |
44 |
|
|
return; |
45 |
|
|
} |
46 |
|
|
|
47 |
|
|
|
48 |
|
|
int |
49 |
|
|
main() { |
50 |
|
|
airArray *mop, *submop; |
51 |
|
|
char *err; |
52 |
|
|
|
53 |
|
|
int typi; |
54 |
|
2 |
unsigned int supi, probePass, cti /* context copy index */, |
55 |
|
|
pvlIdx[NRRD_TYPE_MAX+1], sx, sy, sz, subnum; |
56 |
|
1 |
size_t sizes[3] = {42,61,50} /* one of these must be even */, |
57 |
|
|
ii, nn; |
58 |
|
1 |
Nrrd *norigScl, *nucharScl, *nunquant, *nqdiff, |
59 |
|
|
*nconvScl[NRRD_TYPE_MAX+1]; |
60 |
|
|
unsigned char *ucharScl; |
61 |
|
1 |
gageContext *gctx[2][KERN_SIZE_MAX+1]; |
62 |
|
1 |
gagePerVolume *gpvl[2][NRRD_TYPE_MAX+1][KERN_SIZE_MAX+1]; |
63 |
|
1 |
const double *vansScl[2][NRRD_TYPE_MAX+1][KERN_SIZE_MAX+1], |
64 |
|
|
*gansScl[2][NRRD_TYPE_MAX+1][KERN_SIZE_MAX+1], |
65 |
|
|
*hansScl[2][NRRD_TYPE_MAX+1][KERN_SIZE_MAX+1]; |
66 |
|
2 |
double *origScl, omin, omax, dsx, dsy, dsz, |
67 |
|
1 |
spcOrig[NRRD_SPACE_DIM_MAX] = {0.0, 0.0, 0.0}, |
68 |
|
1 |
spcVec[3][NRRD_SPACE_DIM_MAX] = { |
69 |
|
|
{1.1, 0.0, 0.0}, |
70 |
|
|
{0.0, 2.2, 0.0}, |
71 |
|
|
{0.0, 0.0, 3.3}}; |
72 |
|
|
|
73 |
|
1 |
mop = airMopNew(); |
74 |
|
|
|
75 |
|
|
#define NRRD_NEW(name, mop) \ |
76 |
|
|
(name) = nrrdNew(); \ |
77 |
|
|
airMopAdd((mop), (name), (airMopper)nrrdNuke, airMopAlways) |
78 |
|
|
|
79 |
|
|
/* --------------------------------------------------------------- */ |
80 |
|
|
/* Creating initial volume */ |
81 |
|
1 |
NRRD_NEW(norigScl, mop); |
82 |
✗✓ |
1 |
if (nrrdMaybeAlloc_nva(norigScl, nrrdTypeDouble, 3, sizes)) { |
83 |
|
|
airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways); |
84 |
|
|
fprintf(stderr, "trouble allocating:\n%s", err); |
85 |
|
|
airMopError(mop); return 1; |
86 |
|
|
} |
87 |
|
1 |
origScl = AIR_CAST(double *, norigScl->data); |
88 |
|
1 |
nn = nrrdElementNumber(norigScl); |
89 |
|
1 |
airSrandMT(42*42); |
90 |
✓✓ |
128102 |
for (ii=0; ii<nn/2; ii++) { |
91 |
|
64050 |
airNormalRand(origScl + 2*ii + 0, origScl + 2*ii + 1); |
92 |
|
|
} |
93 |
|
|
/* learn real range */ |
94 |
|
1 |
omin = omax = origScl[0]; |
95 |
✓✓ |
256200 |
for (ii=1; ii<nn; ii++) { |
96 |
✓✓ |
384297 |
omin = AIR_MIN(omin, origScl[ii]); |
97 |
✓✓ |
384297 |
omax = AIR_MAX(omax, origScl[ii]); |
98 |
|
|
} |
99 |
|
1 |
ELL_3V_SET(spcOrig, 0.0, 0.0, 0.0); |
100 |
✗✓ |
2 |
if (nrrdSpaceSet(norigScl, nrrdSpaceRightAnteriorSuperior) |
101 |
✓✗ |
2 |
|| nrrdSpaceOriginSet(norigScl, spcOrig)) { |
102 |
|
|
airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways); |
103 |
|
|
fprintf(stderr, "trouble setting space:\n%s", err); |
104 |
|
|
airMopError(mop); return 1; |
105 |
|
|
} |
106 |
|
1 |
nrrdAxisInfoSet_nva(norigScl, nrrdAxisInfoSpaceDirection, spcVec); |
107 |
|
1 |
dsx = AIR_CAST(double, sizes[0]); |
108 |
|
1 |
dsy = AIR_CAST(double, sizes[1]); |
109 |
|
1 |
dsz = AIR_CAST(double, sizes[2]); |
110 |
|
1 |
sx = AIR_CAST(unsigned int, sizes[0]); |
111 |
|
1 |
sy = AIR_CAST(unsigned int, sizes[1]); |
112 |
|
1 |
sz = AIR_CAST(unsigned int, sizes[2]); |
113 |
|
|
subnum = AIR_CAST(unsigned int, PROBE_NUM*0.9); |
114 |
|
|
|
115 |
|
|
|
116 |
|
|
/* --------------------------------------------------------------- */ |
117 |
|
|
/* Quantizing to 8-bits and checking */ |
118 |
|
1 |
submop = airMopNew(); |
119 |
|
1 |
NRRD_NEW(nucharScl, mop); |
120 |
|
1 |
NRRD_NEW(nunquant, submop); |
121 |
|
1 |
NRRD_NEW(nqdiff, submop); |
122 |
✗✓ |
2 |
if (nrrdQuantize(nucharScl, norigScl, NULL, 8) |
123 |
✓✗ |
2 |
|| nrrdUnquantize(nunquant, nucharScl, nrrdTypeDouble) |
124 |
✓✗ |
2 |
|| nrrdArithBinaryOp(nqdiff, nrrdBinaryOpSubtract, |
125 |
|
|
norigScl, nunquant)) { |
126 |
|
|
airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways); |
127 |
|
|
fprintf(stderr, "trouble quantizing and back:\n%s", err); |
128 |
|
|
airMopError(submop); airMopError(mop); return 1; |
129 |
|
|
} |
130 |
✗✓ |
2 |
if (!( nucharScl->oldMin == omin |
131 |
✓✗ |
2 |
&& nucharScl->oldMax == omax )) { |
132 |
|
|
fprintf(stderr, "quantization range [%g,%g] != real range [%g,%g]\n", |
133 |
|
|
nucharScl->oldMin, nucharScl->oldMax, omin, omax); |
134 |
|
|
airMopError(submop); airMopError(mop); return 1; |
135 |
|
|
} |
136 |
|
|
{ |
137 |
|
|
double *qdiff, *unquant; |
138 |
|
|
/* empirically determined tolerance, which had to be increased in |
139 |
|
|
order to work under valgrind (!)- perhaps because of a |
140 |
|
|
difference in the use of 80-bit registers */ |
141 |
|
|
double epsilon=0.50000000000004; |
142 |
|
1 |
qdiff = AIR_CAST(double *, nqdiff->data); |
143 |
|
1 |
unquant = AIR_CAST(double *, nunquant->data); |
144 |
✓✓ |
256202 |
for (ii=0; ii<nn; ii++) { |
145 |
|
|
double dd; |
146 |
|
|
/* with infinite precision, the max difference between original and |
147 |
|
|
quantized values should be exactly half the width (in value) |
148 |
|
|
of 1/256 of value range ==> dd = 0.5 */ |
149 |
|
128100 |
dd = qdiff[ii]*256/(omax - omin); |
150 |
✗✓ |
128100 |
if (AIR_ABS(dd) > epsilon) { |
151 |
|
|
unsigned int ui; |
152 |
|
|
ui = AIR_CAST(unsigned int, ii); |
153 |
|
|
fprintf(stderr, "|orig[%u]=%.17g - unquant=%.17g|*256/%.17g " |
154 |
|
|
"= %.17g > %.17g!\n", ui, origScl[ii], unquant[ii], |
155 |
|
|
omax - omin, AIR_ABS(dd), epsilon); |
156 |
|
|
airMopError(submop); airMopError(mop); return 1; |
157 |
|
|
} |
158 |
|
128100 |
} |
159 |
|
1 |
} |
160 |
|
1 |
airMopOkay(submop); |
161 |
|
1 |
ucharScl = AIR_CAST(unsigned char *, nucharScl->data); |
162 |
|
|
|
163 |
|
|
/* --------------------------------------------------------------- */ |
164 |
|
|
/* Converting to all other types */ |
165 |
✓✓ |
24 |
for (typi=nrrdTypeUnknown+1; typi<nrrdTypeLast; typi++) { |
166 |
✓✓ |
11 |
if (nrrdTypeBlock == typi) { |
167 |
|
1 |
nconvScl[typi] = NULL; |
168 |
|
1 |
continue; |
169 |
|
|
} |
170 |
|
10 |
NRRD_NEW(nconvScl[typi], mop); |
171 |
✗✓ |
10 |
if (nrrdConvert(nconvScl[typi], nucharScl, typi)) { |
172 |
|
|
airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways); |
173 |
|
|
fprintf(stderr, "trouble converting:\n%s", err); |
174 |
|
|
airMopError(mop); return 1; |
175 |
|
|
} |
176 |
|
|
} |
177 |
✓✓ |
22 |
for (supi=1; supi<=KERN_SIZE_MAX; supi++) { |
178 |
|
|
unsigned int pvii; |
179 |
|
|
int E; |
180 |
|
10 |
double kparm[1]; |
181 |
|
10 |
gctx[0][supi] = gageContextNew(); |
182 |
|
10 |
airMopAdd(mop, gctx[0][supi], (airMopper)gageContextNix, airMopAlways); |
183 |
|
10 |
gageParmSet(gctx[0][supi], gageParmRenormalize, AIR_FALSE); |
184 |
|
10 |
gageParmSet(gctx[0][supi], gageParmCheckIntegrals, AIR_TRUE); |
185 |
|
10 |
kparm[0] = supi; |
186 |
|
|
E = 0; |
187 |
✓✗ |
30 |
if (!E) E |= gageKernelSet(gctx[0][supi], gageKernel00, |
188 |
|
10 |
nrrdKernelBoxSupportDebug, kparm); |
189 |
|
|
pvii = 0; |
190 |
✓✓ |
240 |
for (typi=nrrdTypeUnknown+1; typi<nrrdTypeLast; typi++) { |
191 |
✓✓ |
110 |
if (nrrdTypeBlock == typi) { |
192 |
|
10 |
gpvl[0][typi][supi] = NULL; |
193 |
|
10 |
continue; |
194 |
|
|
} |
195 |
✓✗ |
300 |
if (!E) E |= !(gpvl[0][typi][supi] |
196 |
|
300 |
= gagePerVolumeNew(gctx[0][supi], nconvScl[typi], |
197 |
|
100 |
gageKindScl)); |
198 |
✓✗ |
200 |
if (!E) E |= gagePerVolumeAttach(gctx[0][supi], gpvl[0][typi][supi]); |
199 |
✓✓ |
100 |
if (1 == supi) { |
200 |
|
|
/* first time through this typi loop; its the occasion to |
201 |
|
|
set the pvlIdx array, which records the index into the |
202 |
|
|
gageContext->pvl[] array for the per-type perVolume. |
203 |
|
|
Having to do this is a symptom of bad API design in gage */ |
204 |
|
10 |
pvlIdx[typi] = pvii++; |
205 |
|
10 |
} |
206 |
✓✗ |
200 |
if (!E) E |= gageQueryItemOn(gctx[0][supi], gpvl[0][typi][supi], |
207 |
|
|
gageSclValue); |
208 |
✓✗ |
100 |
if (E) { |
209 |
|
|
break; |
210 |
|
|
} |
211 |
|
|
} |
212 |
✓✗ |
20 |
if (!E) E |= gageUpdate(gctx[0][supi]); |
213 |
✗✓ |
10 |
if (E) { |
214 |
|
|
airMopAdd(mop, err = biffGetDone(GAGE), airFree, airMopAlways); |
215 |
|
|
fprintf(stderr, "trouble (supi=%u, %d/%s) set-up:\n%s\n", |
216 |
|
|
supi, typi, airEnumStr(nrrdType, typi), err); |
217 |
|
|
airMopError(mop); return 1; |
218 |
|
|
} |
219 |
✗✓ |
10 |
if (gctx[0][supi]->radius != supi) { |
220 |
|
|
fprintf(stderr, "supi %u != gageContext->radius %u\n", |
221 |
|
|
supi, gctx[0][supi]->radius); |
222 |
|
|
airMopError(mop); return 1; |
223 |
|
|
} |
224 |
✓✓ |
240 |
for (typi=nrrdTypeUnknown+1; typi<nrrdTypeLast; typi++) { |
225 |
✓✓ |
110 |
if (nrrdTypeBlock == typi) { |
226 |
|
10 |
vansScl[0][typi][supi] = NULL; |
227 |
|
10 |
gansScl[0][typi][supi] = NULL; |
228 |
|
10 |
hansScl[0][typi][supi] = NULL; |
229 |
|
10 |
continue; |
230 |
|
|
} |
231 |
|
100 |
vansScl[0][typi][supi] = |
232 |
|
100 |
gageAnswerPointer(gctx[0][supi], gpvl[0][typi][supi], gageSclValue); |
233 |
|
100 |
gansScl[0][typi][supi] = |
234 |
|
100 |
gageAnswerPointer(gctx[0][supi], gpvl[0][typi][supi], gageSclGradVec); |
235 |
|
100 |
hansScl[0][typi][supi] = |
236 |
|
100 |
gageAnswerPointer(gctx[0][supi], gpvl[0][typi][supi], gageSclHessian); |
237 |
|
100 |
} |
238 |
✓✗ |
20 |
} |
239 |
|
|
|
240 |
|
|
/* --------------------------------------------------------------- */ |
241 |
|
|
/* Exercising gageContextCopy */ |
242 |
✓✓ |
22 |
for (supi=1; supi<=KERN_SIZE_MAX; supi++) { |
243 |
✗✓ |
10 |
if (!(gctx[1][supi] = gageContextCopy(gctx[0][supi]))) { |
244 |
|
|
airMopAdd(mop, err = biffGetDone(GAGE), airFree, airMopAlways); |
245 |
|
|
fprintf(stderr, "trouble copying gctx[%u]:\n%s\n", |
246 |
|
|
supi, err); |
247 |
|
|
airMopError(mop); return 1; |
248 |
|
|
} |
249 |
✓✓ |
240 |
for (typi=nrrdTypeUnknown+1; typi<nrrdTypeLast; typi++) { |
250 |
✓✓ |
110 |
if (nrrdTypeBlock == typi) { |
251 |
|
10 |
vansScl[1][typi][supi] = NULL; |
252 |
|
10 |
gansScl[1][typi][supi] = NULL; |
253 |
|
10 |
hansScl[1][typi][supi] = NULL; |
254 |
|
10 |
continue; |
255 |
|
|
} |
256 |
|
100 |
gpvl[1][typi][supi] = gctx[1][supi]->pvl[pvlIdx[typi]]; |
257 |
|
100 |
vansScl[1][typi][supi] = |
258 |
|
100 |
gageAnswerPointer(gctx[1][supi], gpvl[1][typi][supi], gageSclValue); |
259 |
|
100 |
gansScl[1][typi][supi] = |
260 |
|
100 |
gageAnswerPointer(gctx[1][supi], gpvl[1][typi][supi], gageSclGradVec); |
261 |
|
100 |
hansScl[1][typi][supi] = |
262 |
|
100 |
gageAnswerPointer(gctx[1][supi], gpvl[1][typi][supi], gageSclHessian); |
263 |
|
100 |
} |
264 |
|
|
} |
265 |
|
|
|
266 |
|
|
/* --------------------------------------------------------------- */ |
267 |
|
|
/* the two different probing passes are just to use two different |
268 |
|
|
sets of kernels (first nrrdKernelBoxSupportDebug on pass 0, then |
269 |
|
|
nrrdKernelCos4SupportDebug and its derivatives on pass 1). |
270 |
|
|
Because nrrdKernelBoxSupportDebug has already been set prior to |
271 |
|
|
the first pass, there is some some second-time only "if (1 == |
272 |
|
|
probePass)" logic that confuses the clarity of this */ |
273 |
✓✓ |
6 |
for (probePass=0; probePass<=1; probePass++) { |
274 |
|
|
unsigned int prbi, lastjj=0, xi, yi, zi; |
275 |
|
2 |
double thet, xu, yu, zu, dxi, dyi, dzi, |
276 |
|
|
elapsed[2][KERN_SIZE_MAX+1], time0; |
277 |
|
2 |
char errpre[AIR_STRLEN_LARGE]; |
278 |
|
|
|
279 |
✓✓ |
2 |
if (1 == probePass) { |
280 |
|
|
/* switch to cos^4 kernel, turn on gradient and hessian */ |
281 |
✓✓ |
6 |
for (cti=0; cti<2; cti++) { |
282 |
✓✓ |
44 |
for (supi=1; supi<=KERN_SIZE_MAX; supi++) { |
283 |
|
|
int E; |
284 |
|
20 |
double kparm[1]; |
285 |
|
20 |
gageParmSet(gctx[cti][supi], gageParmRenormalize, AIR_FALSE); |
286 |
|
20 |
gageParmSet(gctx[cti][supi], gageParmCheckIntegrals, AIR_TRUE); |
287 |
|
20 |
kparm[0] = supi; |
288 |
|
|
E = 0; |
289 |
✓✗ |
60 |
if (!E) E |= gageKernelSet(gctx[cti][supi], gageKernel00, |
290 |
|
20 |
nrrdKernelCos4SupportDebug, kparm); |
291 |
✓✗ |
60 |
if (!E) E |= gageKernelSet(gctx[cti][supi], gageKernel11, |
292 |
|
20 |
nrrdKernelCos4SupportDebugD, kparm); |
293 |
✓✗ |
60 |
if (!E) E |= gageKernelSet(gctx[cti][supi], gageKernel22, |
294 |
|
20 |
nrrdKernelCos4SupportDebugDD, kparm); |
295 |
✓✓ |
480 |
for (typi=nrrdTypeUnknown+1; typi<nrrdTypeLast; typi++) { |
296 |
✓✓ |
220 |
if (nrrdTypeBlock == typi) { |
297 |
|
|
continue; |
298 |
|
|
} |
299 |
✓✗ |
600 |
if (!E) E |= gageQueryItemOn(gctx[cti][supi], |
300 |
|
200 |
gpvl[cti][typi][supi], |
301 |
|
|
gageSclGradVec); |
302 |
✓✗ |
600 |
if (!E) E |= gageQueryItemOn(gctx[cti][supi], |
303 |
|
200 |
gpvl[cti][typi][supi], |
304 |
|
|
gageSclHessian); |
305 |
✓✗ |
200 |
if (E) { |
306 |
|
|
break; |
307 |
|
|
} |
308 |
|
|
} |
309 |
✓✗ |
40 |
if (!E) E |= gageUpdate(gctx[cti][supi]); |
310 |
✗✓ |
20 |
if (E) { |
311 |
|
|
airMopAdd(mop, err = biffGetDone(GAGE), airFree, airMopAlways); |
312 |
|
|
fprintf(stderr, "trouble (cti=%u, supi=%u, %d/%s) " |
313 |
|
|
"set-up:\n%s\n", cti, |
314 |
|
|
supi, typi, airEnumStr(nrrdType, typi), err); |
315 |
|
|
airMopError(mop); return 1; |
316 |
|
|
} |
317 |
✓✗ |
40 |
} |
318 |
|
|
} |
319 |
|
|
} |
320 |
✓✓ |
12 |
for (cti=0; cti<2; cti++) { |
321 |
✓✓ |
88 |
for (supi=1; supi<=KERN_SIZE_MAX; supi++) { |
322 |
|
40 |
elapsed[cti][supi] = 0.0; |
323 |
|
|
} |
324 |
|
|
} |
325 |
|
|
/* do the probes along a curvy path */ |
326 |
✓✓ |
1204 |
for (prbi=0; prbi<PROBE_NUM; prbi++) { |
327 |
|
|
unsigned int jj; |
328 |
|
600 |
jj = airIndex(0, prbi, PROBE_NUM-1, subnum); |
329 |
|
600 |
thet = AIR_AFFINE(0, jj, subnum-1, 0.0, AIR_PI); |
330 |
|
600 |
xu = -cos(5*thet); |
331 |
|
600 |
yu = -cos(3*thet); |
332 |
|
600 |
zu = -cos(thet); |
333 |
|
600 |
dxi = AIR_AFFINE(-1.0, xu, 1.0, -0.5, dsx-0.5); |
334 |
|
600 |
dyi = AIR_AFFINE(-1.0, yu, 1.0, -0.5, dsy-0.5); |
335 |
|
600 |
dzi = AIR_AFFINE(-1.0, zu, 1.0, -0.5, dsz-0.5); |
336 |
✓✓✓✓
|
1198 |
if (prbi && lastjj == jj) { |
337 |
|
|
/* this occasionally tests the logic in gage that seeks to |
338 |
|
|
re-compute convolution weights only when necessary */ |
339 |
|
60 |
dxi += airSgn(xu); |
340 |
|
60 |
dyi += airSgn(yu); |
341 |
|
60 |
dzi += airSgn(zu); |
342 |
|
60 |
} |
343 |
|
600 |
xi = airIndexClamp(-0.5, dxi, dsx-0.5, sx); |
344 |
|
600 |
yi = airIndexClamp(-0.5, dyi, dsy-0.5, sy); |
345 |
|
600 |
zi = airIndexClamp(-0.5, dzi, dsz-0.5, sz); |
346 |
|
|
lastjj = jj; |
347 |
✓✓ |
13200 |
for (supi=1; supi<=KERN_SIZE_MAX; supi++) { |
348 |
|
6000 |
double truevalOrig[NRRD_TYPE_MAX+1]; |
349 |
✓✓ |
36000 |
for (cti=0; cti<2; cti++) { |
350 |
|
12000 |
time0 = airTime(); |
351 |
✗✓ |
12000 |
if (gageProbeSpace(gctx[cti][supi], dxi, dyi, dzi, |
352 |
|
|
AIR_TRUE /* indexSpace */, |
353 |
|
|
AIR_TRUE /* clamp */)) { |
354 |
|
|
fprintf(stderr, "probe (cti %u support %u) error (%d): %s\n", |
355 |
|
|
cti, supi, |
356 |
|
|
gctx[cti][supi]->errNum, gctx[cti][supi]->errStr); |
357 |
|
|
airMopError(mop); return 1; |
358 |
|
|
} |
359 |
|
12000 |
elapsed[cti][supi] = airTime() - time0; |
360 |
✓✓ |
288000 |
for (typi=nrrdTypeUnknown+1; typi<nrrdTypeLast; typi++) { |
361 |
|
|
double arrayval, trueval, probeval; |
362 |
|
132000 |
if (nrrdTypeBlock == typi |
363 |
✓✓✓✓
|
252000 |
|| (1 == probePass && nrrdTypeChar == typi)) { |
364 |
|
|
/* can't easily correct interpolation on signed char |
365 |
|
|
values to make it match interpolation on unsigned char |
366 |
|
|
values, prior to wrap-around */ |
367 |
|
18000 |
continue; |
368 |
|
|
} |
369 |
|
|
/* probeval is what we learned by probing with gage */ |
370 |
|
114000 |
probeval = vansScl[cti][typi][supi][0]; |
371 |
✓✓ |
114000 |
if (0 == probePass) { |
372 |
|
|
/* arrayval is the value directly from array of same type |
373 |
|
|
(converted from original uchar) */ |
374 |
|
120000 |
arrayval = (nrrdDLookup[typi])(nconvScl[typi]->data, |
375 |
|
60000 |
xi + sx*(yi + sy*zi)); |
376 |
|
|
/* when using box, gage-reconstructed value should |
377 |
|
|
match value from probing */ |
378 |
✗✓ |
60000 |
if (arrayval != probeval) { |
379 |
|
|
#define SPRINT_ERR_PREFIX \ |
380 |
|
|
errPrefix(errpre, typi, supi, prbi, probePass, \ |
381 |
|
|
dxi, dyi, dzi, xi, yi, zi) |
382 |
|
|
SPRINT_ERR_PREFIX; |
383 |
|
|
fprintf(stderr, "%s: (cti %u) probed %g != conv %g\n", errpre, |
384 |
|
|
cti, probeval, arrayval); |
385 |
|
|
airMopError(mop); return 1; |
386 |
|
|
} |
387 |
|
|
/* trueval on pass 0 is the original uchar value */ |
388 |
|
60000 |
trueval = AIR_CAST(double, ucharScl[xi + sx*(yi + sy*zi)]); |
389 |
|
60000 |
} else { |
390 |
|
|
/* trueval on pass 1 is the value from probing uchar volume */ |
391 |
|
54000 |
trueval = vansScl[cti][nrrdTypeUChar][supi][0]; |
392 |
|
|
} |
393 |
✓✓ |
114000 |
if (nrrdTypeChar == typi && trueval > 127) { |
394 |
|
|
/* recreate value wrapping of signed char */ |
395 |
|
2740 |
trueval -= 256; |
396 |
|
2740 |
} |
397 |
✓✓ |
114000 |
if (0 == cti) { |
398 |
|
57000 |
truevalOrig[typi] = trueval; |
399 |
|
57000 |
} else { |
400 |
|
|
/* make sure that result from gageContextCopy'd context |
401 |
|
|
(trueval) is same result as original (truevalOrig[typi]) */ |
402 |
✗✓ |
57000 |
if (truevalOrig[typi] != trueval) { |
403 |
|
|
SPRINT_ERR_PREFIX; |
404 |
|
|
fprintf(stderr, "%s: original->%g, gageContextCopy->%g\n", |
405 |
|
|
errpre, truevalOrig[typi], trueval); |
406 |
|
|
airMopError(mop); return 1; |
407 |
|
|
} |
408 |
|
|
} |
409 |
|
|
/* regardless of the volume (excepting where we've continue'd, |
410 |
|
|
above) the reconstructed value probeval should match trueval */ |
411 |
✗✓ |
114000 |
if (trueval != probeval) { |
412 |
|
|
SPRINT_ERR_PREFIX; |
413 |
|
|
fprintf(stderr, "%s: (cti %u) probed %g != true %g\n", errpre, |
414 |
|
|
cti, probeval, trueval); |
415 |
|
|
airMopError(mop); return 1; |
416 |
|
|
} |
417 |
✓✓ |
114000 |
if (1 == probePass) { |
418 |
|
|
/* and when we use a differentiable kernel, the gradient |
419 |
|
|
and Hessian had better agree too */ |
420 |
|
|
double diff3[3], diff9[9]; |
421 |
|
54000 |
ELL_3V_SUB(diff3, |
422 |
|
|
gansScl[cti][nrrdTypeUChar][supi], |
423 |
|
|
gansScl[cti][typi][supi]); |
424 |
✗✓ |
54000 |
if (ELL_3V_LEN(diff3) > 0.0) { |
425 |
|
|
SPRINT_ERR_PREFIX; |
426 |
|
|
fprintf(stderr, "%s: (cti %u) probed gradient error len %f\n", |
427 |
|
|
errpre, cti, ELL_3V_LEN(diff3)); |
428 |
|
|
airMopError(mop); return 1; |
429 |
|
|
} |
430 |
|
54000 |
ELL_9V_SUB(diff9, |
431 |
|
|
hansScl[cti][nrrdTypeUChar][supi], |
432 |
|
|
hansScl[cti][typi][supi]); |
433 |
✗✓ |
54000 |
if (ELL_9V_LEN(diff9) > 0.0) { |
434 |
|
|
SPRINT_ERR_PREFIX; |
435 |
|
|
fprintf(stderr, "%s: (cti %u) probed hessian error len %f\n", |
436 |
|
|
errpre, cti, ELL_9V_LEN(diff9)); |
437 |
|
|
airMopError(mop); return 1; |
438 |
|
|
} |
439 |
|
54000 |
} |
440 |
|
114000 |
} |
441 |
|
|
} |
442 |
✓✗ |
12000 |
} |
443 |
|
600 |
} |
444 |
✓✓ |
12 |
for (cti=0; cti<2; cti++) { |
445 |
✓✓ |
88 |
for (supi=1; supi<=KERN_SIZE_MAX; supi++) { |
446 |
|
80 |
fprintf(stderr, "elapsed[%u][%u] = %g ms\n", |
447 |
|
40 |
cti, supi, 1000*elapsed[cti][supi]); |
448 |
|
|
} |
449 |
|
|
} |
450 |
✓✗ |
4 |
} |
451 |
|
|
|
452 |
|
|
#undef NRRD_NEW |
453 |
|
|
#undef SPRINT_ERR_PREFIX |
454 |
|
1 |
airMopOkay(mop); |
455 |
|
1 |
exit(0); |
456 |
|
|
} |