Bug Summary

File:Testing/meet/probeSS.c
Location:line 532, column 16
Description:Call to 'calloc' has an allocation size of 0 bytes

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 "teem/meet.h"
25
26/*
27** Tests:
28** ... lots of gage stuff ...
29**
30** The main point of all this is to make sure of two (or maybe 2.5) separate
31** things about gage, the testing of which requires so much in common that one
32** program might as well do them all. The tasks are:
33**
34** 1) ensure that values and their derivatives (where the gageKind supports
35** it) are correctly handled in the multi-value gageKinds (gageKindVec,
36** tenGageKind, tenDwiGageKind), relative to the gageKindScl ground-truth: the
37** per-component values and derivatives have to match those reconstructed from
38** sets of scalar volumes of the components.
39** ==> Also, that gageContextCopy works even on dynamic kinds (like DWIs)
40**
41** 1.5) that there is expected consistency between the scalar, vector, tensor,
42** and DWI properties of a related set of volumes. So the test starts with a
43** diffusion tensor field and generates DWIs, scalars (norm squared of the
44** tensor), and vectors (gradient of norm squared) from this.
45**
46** 2) that scale-space reconstruction works: that sets of pre-blurred volumes
47** can be generated and saved via the utilities in meet, that the smart
48** hermite-spline based scale-space interpolation is working (for all kinds),
49** and that gageContextCopy works on scale-space contexts
50*/
51
52#define BKEY"probeSS" "probeSS"
53#define KIND_NUM4 4
54#define KI_SCL0 0
55#define KI_VEC1 1
56#define KI_TEN2 2
57#define KI_DWI3 3
58
59/*
60static const char *kindStr[KIND_NUM] = {
61 "scalar",
62 "vector",
63 "tensor",
64 " DWI "
65};
66*/
67
68#define HALTON_BASE100 100
69
70static unsigned int volSize[3] = {45, 46, 47};
71
72#define NRRD_NEW(nn, mm) \
73 (nn) = nrrdNew(); \
74 airMopAdd((mm), (nn), (airMopper)nrrdNuke, airMopAlways)
75
76/* the weird function names of the various local functions here (created in a
77 rather adhoc and organic way) should be read in usual Teem order: from
78 right to left */
79
80static int
81engageGenTensor(gageContext *gctx, Nrrd *ncten,
82 /* out/in */
83 double noiseStdv, unsigned int seed,
84 unsigned int sx, unsigned int sy, unsigned int sz) {
85 static const char me[]="engageGenTensor";
86 hestParm *hparm;
87 airArray *smop;
88 char tmpStr[4][AIR_STRLEN_SMALL(128+1)];
89 Nrrd *nclean;
90 NrrdIter *narg0, *narg1;
91 const char *helixArgv[] =
92 /* 0 1 2 3 4 5 6 7 8 9 */
93 {"-s", NULL((void*)0), NULL((void*)0), NULL((void*)0), "-v", "0", "-r", "40", "-o", NULL((void*)0),
94 "-ev", "0.00086", "0.00043", "0.00021", "-bg", "0.003", "-b", "3",
95 NULL((void*)0)};
96 int helixArgc;
97 gagePerVolume *pvl;
98
99 smop = airMopNew();
100 /* NOTE: this is currently the only place where a unrrduCmd
101 is called from within C code; it was educational to get working.
102 Learned:
103 * hest does NOT tolerate having empty or NULL elements of
104 its argv[]! More error checking for this in hest is needed.
105 * the "const char **argv" type is not very convenient to
106 set up in a dynamic way; the per-element setting done below
107 is certainly awkward
108 */
109 hparm = hestParmNew();
110 airMopAdd(smop, hparm, (airMopper)hestParmFree, airMopAlways);
111 sprintf(tmpStr[0], "%u", sx)__builtin___sprintf_chk (tmpStr[0], 0, __builtin_object_size (
tmpStr[0], 2 > 1 ? 1 : 0), "%u", sx); helixArgv[1] = tmpStr[0];
112 sprintf(tmpStr[1], "%u", sy)__builtin___sprintf_chk (tmpStr[1], 0, __builtin_object_size (
tmpStr[1], 2 > 1 ? 1 : 0), "%u", sy); helixArgv[2] = tmpStr[1];
113 sprintf(tmpStr[2], "%u", sz)__builtin___sprintf_chk (tmpStr[2], 0, __builtin_object_size (
tmpStr[2], 2 > 1 ? 1 : 0), "%u", sz); helixArgv[3] = tmpStr[2];
114 sprintf(tmpStr[3], "tmp-ten.nrrd")__builtin___sprintf_chk (tmpStr[3], 0, __builtin_object_size (
tmpStr[3], 2 > 1 ? 1 : 0), "tmp-ten.nrrd");
115 helixArgv[9] = tmpStr[3];
116 helixArgc = AIR_CAST(int, sizeof(helixArgv)/sizeof(char *))((int)(sizeof(helixArgv)/sizeof(char *))) - 1;
117 if (tend_helixCmd.main(helixArgc, helixArgv, me, hparm)) {
118 /* error already went to stderr, not to any biff container */
119 biffAddf(BKEY"probeSS", "%s: problem running tend %s", me, tend_helixCmd.name);
120 airMopError(smop); return 1;
121 }
122 NRRD_NEW(nclean, smop);
123 if (nrrdLoad(nclean, tmpStr[3], NULL((void*)0))) {
124 biffAddf(BKEY"probeSS", "%s: trouble loading from new vol %s", me, tmpStr[3]);
125 airMopError(smop); return 1;
126 }
127
128 /* add some noise to tensor value; no, this isn't really physical;
129 since we're adding noise to the tensor and then simulating DWIs,
130 rather than adding noise to DWIs and then estimating tensor,
131 but for the purposes of gage testing its fine */
132 narg0 = nrrdIterNew();
133 narg1 = nrrdIterNew();
134 airMopAdd(smop, narg0, (airMopper)nrrdIterNix, airMopAlways);
135 airMopAdd(smop, narg1, (airMopper)nrrdIterNix, airMopAlways);
136 nrrdIterSetNrrd(narg0, nclean);
137 nrrdIterSetValue(narg1, noiseStdv);
138 airSrandMT(seed);
139 if (nrrdArithIterBinaryOp(ncten, nrrdBinaryOpNormalRandScaleAdd,
140 narg0, narg1)) {
141 biffMovef(BKEY"probeSS", NRRDnrrdBiffKey, "%s: trouble noisying output", me);
142 airMopError(smop); return 1;
143 }
144
145 /* wrap in gage context */
146 if ( !(pvl = gagePerVolumeNew(gctx, ncten, tenGageKind))
147 || gagePerVolumeAttach(gctx, pvl) ) {
148 biffMovef(BKEY"probeSS", GAGEgageBiffKey, "%s: trouble engaging tensor", me);
149 return 1;
150 }
151
152 airMopOkay(smop);
153 return 0;
154}
155
156/*
157** takes in ncten, measures "S" to nscl, copies that to nsclCopy,
158** wraps those in gctx and gctxComp
159*/
160static int
161engageGenScalar(gageContext *gctx, Nrrd *nscl,
162 gageContext *gctxComp, Nrrd *nsclCopy,
163 /* out/in */
164 const Nrrd *ncten) {
165 static const char me[]="engageGenScalar";
166 gagePerVolume *pvl;
167
168 if (tenAnisoVolume(nscl, ncten, tenAniso_S, 0)) {
169 biffMovef(BKEY"probeSS", TENtenBiffKey, "%s: trouble creating scalar volume", me);
170 return 1;
171 }
172 if (nrrdCopy(nsclCopy, nscl)) {
173 biffMovef(BKEY"probeSS", NRRDnrrdBiffKey, "%s: trouble copying scalar volume", me);
174 return 1;
175 }
176
177 /* wrap both in gage context */
178 if ( !(pvl = gagePerVolumeNew(gctx, nscl, gageKindScl))
179 || gagePerVolumeAttach(gctx, pvl)
180 || !(pvl = gagePerVolumeNew(gctxComp, nsclCopy, gageKindScl))
181 || gagePerVolumeAttach(gctxComp, pvl)) {
182 biffMovef(BKEY"probeSS", GAGEgageBiffKey, "%s: trouble engaging scalars", me);
183 return 1;
184 }
185
186 return 0;
187}
188
189/*
190** Makes a vector volume by measuring the gradient
191** Being the gradient of the given scalar volume is just to make
192** something vaguely interesting, and to test consistency with
193** the the same gradient as measured in the tensor field.
194*/
195static int
196engageGenVector(gageContext *gctx, Nrrd *nvec,
197 /* out/in */
198 const Nrrd *nscl) {
199 static const char me[]="engageGenVector";
200 ptrdiff_t padMin[4] = {0, 0, 0, 0}, padMax[4];
201 Nrrd *ntmp;
202 airArray *smop;
203 float *vec, *scl;
204 size_t sx, sy, sz, xi, yi, zi, Px, Mx, Py, My, Pz, Mz;
205 double dnmX, dnmY, dnmZ;
206 gagePerVolume *pvl;
207
208 smop = airMopNew();
209 if (nrrdTypeFloat != nscl->type) {
210 biffAddf(BKEY"probeSS", "%s: expected %s not %s type", me,
211 airEnumStr(nrrdType, nrrdTypeFloat),
212 airEnumStr(nrrdType, nscl->type));
213 airMopError(smop); return 1;
214 }
215 NRRD_NEW(ntmp, smop);
216 sx = nscl->axis[0].size;
217 sy = nscl->axis[1].size;
218 sz = nscl->axis[2].size;
219 ELL_4V_SET(padMax, 2,((padMax)[0] = (2), (padMax)[1] = (((ptrdiff_t)(sx-1))), (padMax
)[2] = (((ptrdiff_t)(sy-1))), (padMax)[3] = (((ptrdiff_t)(sz-
1))))
220 AIR_CAST(ptrdiff_t, sx-1),((padMax)[0] = (2), (padMax)[1] = (((ptrdiff_t)(sx-1))), (padMax
)[2] = (((ptrdiff_t)(sy-1))), (padMax)[3] = (((ptrdiff_t)(sz-
1))))
221 AIR_CAST(ptrdiff_t, sy-1),((padMax)[0] = (2), (padMax)[1] = (((ptrdiff_t)(sx-1))), (padMax
)[2] = (((ptrdiff_t)(sy-1))), (padMax)[3] = (((ptrdiff_t)(sz-
1))))
222 AIR_CAST(ptrdiff_t, sz-1))((padMax)[0] = (2), (padMax)[1] = (((ptrdiff_t)(sx-1))), (padMax
)[2] = (((ptrdiff_t)(sy-1))), (padMax)[3] = (((ptrdiff_t)(sz-
1))));
223 /* we do axinsert and pad in order to keep all the per-axis info */
224 if (nrrdAxesInsert(ntmp, nscl, 0)
225 || nrrdPad_nva(nvec, ntmp, padMin, padMax,
226 nrrdBoundaryPad, 0.0)) {
227 biffMovef(BKEY"probeSS", NRRDnrrdBiffKey, "%s: trouble", me);
228 airMopError(smop); return 1;
229 }
230 dnmX = 0.5/nrrdSpaceVecNorm(nscl->spaceDim, nscl->axis[0].spaceDirection);
231 dnmY = 0.5/nrrdSpaceVecNorm(nscl->spaceDim, nscl->axis[1].spaceDirection);
232 dnmZ = 0.5/nrrdSpaceVecNorm(nscl->spaceDim, nscl->axis[2].spaceDirection);
233 vec = AIR_CAST(float *, nvec->data)((float *)(nvec->data));
234 scl = AIR_CAST(float *, nscl->data)((float *)(nscl->data));
235#define INDEX(xj, yj, zj) (xj + sx*(yj + sy*zj))
236 for (zi=0; zi<sz; zi++) {
237 Mz = (zi == 0 ? 0 : zi-1);
238 Pz = AIR_MIN(zi+1, sz-1)((zi+1) < (sz-1) ? (zi+1) : (sz-1));
239 for (yi=0; yi<sy; yi++) {
240 My = (yi == 0 ? 0 : yi-1);
241 Py = AIR_MIN(yi+1, sy-1)((yi+1) < (sy-1) ? (yi+1) : (sy-1));
242 for (xi=0; xi<sx; xi++) {
243 Mx = (xi == 0 ? 0 : xi-1);
244 Px = AIR_MIN(xi+1, sx-1)((xi+1) < (sx-1) ? (xi+1) : (sx-1));
245 vec[0 + 3*INDEX(xi, yi, zi)] =
246 AIR_CAST(float, (scl[INDEX(Px, yi, zi)] - scl[INDEX(Mx, yi, zi)])*dnmX)((float)((scl[INDEX(Px, yi, zi)] - scl[INDEX(Mx, yi, zi)])*dnmX
));
247 vec[1 + 3*INDEX(xi, yi, zi)] =
248 AIR_CAST(float, (scl[INDEX(xi, Py, zi)] - scl[INDEX(xi, My, zi)])*dnmY)((float)((scl[INDEX(xi, Py, zi)] - scl[INDEX(xi, My, zi)])*dnmY
));
249 vec[2 + 3*INDEX(xi, yi, zi)] =
250 AIR_CAST(float, (scl[INDEX(xi, yi, Pz)] - scl[INDEX(xi, yi, Mz)])*dnmZ)((float)((scl[INDEX(xi, yi, Pz)] - scl[INDEX(xi, yi, Mz)])*dnmZ
));
251 }
252 }
253 }
254#undef INDEX
255
256 /* wrap in gage context */
257 if ( !(pvl = gagePerVolumeNew(gctx, nvec, gageKindVec))
258 || gagePerVolumeAttach(gctx, pvl) ) {
259 biffMovef(BKEY"probeSS", GAGEgageBiffKey, "%s: trouble engaging vectors", me);
260 return 1;
261 }
262 airMopError(smop);
263 return 0;
264}
265
266/*
267** make a DWI volume by simulating DWIs from given tensor
268** this includes generating a new gradient set
269*/
270static int
271genDwi(Nrrd *ndwi, Nrrd *ngrad,
272 /* out/in */
273 unsigned int gradNum, double bval, const Nrrd *ncten) {
274 static const char me[]="genDwi";
275 tenGradientParm *gparm;
276 tenExperSpec *espec;
277 Nrrd *nten, *nb0;
278 NrrdIter *narg0, *narg1;
279 size_t cropMin[4] = {1, 0, 0, 0}, cropMax[4];
280 airArray *smop;
281
282 smop = airMopNew();
283 gparm = tenGradientParmNew();
284 airMopAdd(smop, gparm, (airMopper)tenGradientParmNix, airMopAlways);
285 espec = tenExperSpecNew();
286 airMopAdd(smop, espec, (airMopper)tenExperSpecNix, airMopAlways);
287 gparm->verbose = 0;
288 gparm->minMean = 0.002;
289 gparm->seed = 4242;
290 gparm->insertZeroVec = AIR_TRUE1;
291 if (tenGradientGenerate(ngrad, gradNum, gparm)
292 || tenExperSpecGradSingleBValSet(espec, AIR_FALSE0 /* insertB0 */,
293 bval,
294 AIR_CAST(double *, ngrad->data)((double *)(ngrad->data)),
295 AIR_CAST(unsigned int,((unsigned int)(ngrad->axis[1].size))
296 ngrad->axis[1].size)((unsigned int)(ngrad->axis[1].size)))) {
297 biffMovef(BKEY"probeSS", TENtenBiffKey, "%s: trouble generating grads or espec", me);
298 airMopError(smop); return 1;
299 }
300 NRRD_NEW(nten, smop);
301 NRRD_NEW(nb0, smop);
302 ELL_4V_SET(cropMax,((cropMax)[0] = (ncten->axis[0].size-1), (cropMax)[1] = (ncten
->axis[1].size-1), (cropMax)[2] = (ncten->axis[2].size-
1), (cropMax)[3] = (ncten->axis[3].size-1))
303 ncten->axis[0].size-1,((cropMax)[0] = (ncten->axis[0].size-1), (cropMax)[1] = (ncten
->axis[1].size-1), (cropMax)[2] = (ncten->axis[2].size-
1), (cropMax)[3] = (ncten->axis[3].size-1))
304 ncten->axis[1].size-1,((cropMax)[0] = (ncten->axis[0].size-1), (cropMax)[1] = (ncten
->axis[1].size-1), (cropMax)[2] = (ncten->axis[2].size-
1), (cropMax)[3] = (ncten->axis[3].size-1))
305 ncten->axis[2].size-1,((cropMax)[0] = (ncten->axis[0].size-1), (cropMax)[1] = (ncten
->axis[1].size-1), (cropMax)[2] = (ncten->axis[2].size-
1), (cropMax)[3] = (ncten->axis[3].size-1))
306 ncten->axis[3].size-1)((cropMax)[0] = (ncten->axis[0].size-1), (cropMax)[1] = (ncten
->axis[1].size-1), (cropMax)[2] = (ncten->axis[2].size-
1), (cropMax)[3] = (ncten->axis[3].size-1));
307 if (nrrdSlice(nb0, ncten, 0, 0)
308 || nrrdCrop(nten, ncten, cropMin, cropMax)) {
309 biffMovef(BKEY"probeSS", NRRDnrrdBiffKey, "%s: trouble slicing or cropping ten vol", me);
310 airMopError(smop); return 1;
311 }
312 narg0 = nrrdIterNew();
313 narg1 = nrrdIterNew();
314 airMopAdd(smop, narg0, (airMopper)nrrdIterNix, airMopAlways);
315 airMopAdd(smop, narg1, (airMopper)nrrdIterNix, airMopAlways);
316 nrrdIterSetValue(narg1, 50000.0);
317 nrrdIterSetNrrd(narg0, nb0);
318 if (nrrdArithIterBinaryOp(nb0, nrrdBinaryOpMultiply,
319 narg0, narg1)) {
320 biffMovef(BKEY"probeSS", NRRDnrrdBiffKey, "%s: trouble generating b0 vol", me);
321 airMopError(smop); return 1;
322 }
323 if (tenModelSimulate(ndwi, nrrdTypeUShort, espec,
324 tenModel1Tensor2,
325 nb0, nten, AIR_TRUE1 /* keyValueSet */)) {
326 biffMovef(BKEY"probeSS", TENtenBiffKey, "%s: trouble simulating DWI vol", me);
327 airMopError(smop); return 1;
328 }
329
330 airMopOkay(smop);
331 return 0;
332}
333
334int
335engageMopDiceVector(gageContext *gctx, Nrrd *nvecComp[3],
336 /* out/in */
337 airArray *mop, const Nrrd* nvec) {
338 static const char me[]="engageMopDiceVector";
339 gagePerVolume *pvl;
340 unsigned int ci;
341 char stmp[AIR_STRLEN_SMALL(128+1)];
342
343 if (!( 4 == nvec->dim && 3 == nvec->axis[0].size )) {
344 biffAddf(BKEY"probeSS", "%s: expected 4-D 3-by-X nrrd (not %u-D %s-by-X)",
345 me, nvec->dim, airSprintSize_t(stmp, nvec->axis[0].size));
346 return 1;
347 }
348 for (ci=0; ci<3; ci++) {
349 NRRD_NEW(nvecComp[ci], mop);
350 if (nrrdSlice(nvecComp[ci], nvec, 0, ci)) {
351 biffMovef(BKEY"probeSS", NRRDnrrdBiffKey, "%s: trouble getting component %u", me, ci);
352 return 1;
353 }
354 if ( !(pvl = gagePerVolumeNew(gctx, nvecComp[ci], gageKindScl))
355 || gagePerVolumeAttach(gctx, pvl) ) {
356 biffMovef(BKEY"probeSS", GAGEgageBiffKey, "%s: trouble engaging component %u", me, ci);
357 return 1;
358 }
359 }
360
361 return 0;
362}
363
364int
365engageMopDiceTensor(gageContext *gctx, Nrrd *nctenComp[7],
366 /* out/in */
367 airArray *mop, const Nrrd* ncten) {
368 static const char me[]="engageMopDiceTensor";
369 gagePerVolume *pvl;
370 unsigned int ci;
371
372 if (tenTensorCheck(ncten, nrrdTypeFloat, AIR_TRUE1 /* want4F */,
373 AIR_TRUE1 /* useBiff */)) {
374 biffMovef(BKEY"probeSS", TENtenBiffKey, "%s: didn't get tensor volume", me);
375 return 1;
376 }
377 for (ci=0; ci<7; ci++) {
378 NRRD_NEW(nctenComp[ci], mop);
379 if (nrrdSlice(nctenComp[ci], ncten, 0, ci)) {
380 biffMovef(BKEY"probeSS", NRRDnrrdBiffKey, "%s: trouble getting component %u", me, ci);
381 return 1;
382 }
383 if ( !(pvl = gagePerVolumeNew(gctx, nctenComp[ci], gageKindScl))
384 || gagePerVolumeAttach(gctx, pvl) ) {
385 biffMovef(BKEY"probeSS", GAGEgageBiffKey, "%s: trouble engaging component %u", me, ci);
386 return 1;
387 }
388 }
389
390 return 0;
391}
392
393int
394engageMopDiceDwi(gageContext *gctx, Nrrd ***ndwiCompP,
395 /* out/in */
396 airArray *mop, const Nrrd* ndwi) {
397 static const char me[]="mopDiceDwi";
398 Nrrd **ndwiComp;
399 size_t dwiNum;
400 char stmp[AIR_STRLEN_SMALL(128+1)];
401 gagePerVolume *pvl;
402 unsigned int ci;
403
404 if (!( 4 == ndwi->dim )) {
405 biffAddf(BKEY"probeSS", "%s: wanted 4D volume (not %u)", me, ndwi->dim);
406 return 1;
407 }
408 if (!( nrrdKindList == ndwi->axis[0].kind &&
409 nrrdKindSpace == ndwi->axis[1].kind &&
410 nrrdKindSpace == ndwi->axis[2].kind &&
411 nrrdKindSpace == ndwi->axis[3].kind )) {
412 biffAddf(BKEY"probeSS", "%s: wanted kinds %s,3x%s, not %s,%s,%s,%s", me,
413 airEnumStr(nrrdKind, nrrdKindList),
414 airEnumStr(nrrdKind, nrrdKindSpace),
415 airEnumStr(nrrdKind, ndwi->axis[0].kind),
416 airEnumStr(nrrdKind, ndwi->axis[1].kind),
417 airEnumStr(nrrdKind, ndwi->axis[2].kind),
418 airEnumStr(nrrdKind, ndwi->axis[3].kind));
419 return 1;
420 }
421 dwiNum = ndwi->axis[0].size;
422 if (!(ndwiComp = AIR_CALLOC(dwiNum, Nrrd *)(Nrrd **)(calloc((dwiNum), sizeof(Nrrd *))))) {
423 biffAddf(BKEY"probeSS", "%s: couldn't alloc %s Nrrd*", me,
424 airSprintSize_t(stmp, dwiNum));
425 return 1;
426 }
427 airMopAdd(mop, ndwiComp, airFree, airMopAlways);
428 *ndwiCompP = ndwiComp;
429 for (ci=0; ci<dwiNum; ci++) {
430 NRRD_NEW(ndwiComp[ci], mop);
431 if (nrrdSlice(ndwiComp[ci], ndwi, 0, ci)) {
432 biffMovef(BKEY"probeSS", NRRDnrrdBiffKey, "%s: trouble getting component %u", me, ci);
433 return 1;
434 }
435 if ( !(pvl = gagePerVolumeNew(gctx, ndwiComp[ci], gageKindScl))
436 || gagePerVolumeAttach(gctx, pvl) ) {
437 biffMovef(BKEY"probeSS", GAGEgageBiffKey, "%s: trouble engaging component %u", me, ci);
438 return 1;
439 }
440 }
441 return 0;
442}
443
444typedef struct {
445 char name[AIR_STRLEN_SMALL(128+1)];
446 const double **aptr; /* array of pointers to (const) answers */
447 gageItemSpec *ispec; /* array of gageItemSpecs (not pointers to them) */
448 unsigned int *alen; /* array of answer lengths */
449 unsigned int *sidx; /* array of index (within san) of copied answer */
450 unsigned int anum; /* lenth of aptr, ispec, alen, sidx arrays */
451 double *san; /* single buffer for copy + concating answers */
452 unsigned int slen; /* length of single buffer (sum of all alen[i]) */
453} multiAnswer;
454
455static void
456multiAnswerInit(multiAnswer *man) {
457
458 /* HEY sloppy: not resetting name */
459 man->aptr = airFree(AIR_VOIDP(man->aptr)((void *)(man->aptr)));
460 man->ispec = airFree(man->ispec);
461 man->alen = airFree(man->alen);
462 man->sidx = airFree(man->sidx);
463 man->anum = 0;
464 man->san = airFree(man->san);
465 man->slen = 0;
466}
467
468static multiAnswer*
469multiAnswerNew(char *name) {
470 multiAnswer *man;
471
472 man = AIR_CALLOC(1, multiAnswer)(multiAnswer*)(calloc((1), sizeof(multiAnswer)));
473 airStrcpy(man->name, AIR_STRLEN_SMALL(128+1), name);
474 multiAnswerInit(man);
475 return man;
476}
477
478void
479multiAnswerLenSet(multiAnswer *man, unsigned int anum) {
480 /*
481 static const char me[]="multiAnswerLenSet";
482
483 fprintf(stderr, "!%s: %s hello -> answer number = %u\n", me,
484 man->name, anum);
485 */
486 man->aptr = AIR_CALLOC(anum, const double *)(const double **)(calloc((anum), sizeof(const double *)));
487 man->ispec = AIR_CALLOC(anum, gageItemSpec)(gageItemSpec*)(calloc((anum), sizeof(gageItemSpec)));
488 man->alen = AIR_CALLOC(anum, unsigned int)(unsigned int*)(calloc((anum), sizeof(unsigned int)));
489 man->sidx = AIR_CALLOC(anum, unsigned int)(unsigned int*)(calloc((anum), sizeof(unsigned int)));
490 man->anum = anum;
491 /* don't know answer lengths yet */
492 man->san = airFree(man->san);
493 man->slen = 0;
494 return;
495}
496
497static multiAnswer*
498multiAnswerNix(multiAnswer *man) {
499
500 airFree(AIR_VOIDP(man->aptr)((void *)(man->aptr)));
501 airFree(man->ispec);
502 airFree(man->alen);
503 airFree(man->sidx);
504 airFree(man->san);
505 airFree(man);
506 return NULL((void*)0);
507}
508
509static void
510multiAnswerAdd(multiAnswer *man, unsigned int ansIdx,
511 const gageContext *gctx, const gagePerVolume *pvl,
512 unsigned int item) {
513 /*
514 static const char me[]="multiAnswerAdd";
515
516 fprintf(stderr, "!%s: %s hello (aidx = %u)\n", me, man->name, ansIdx);
517 */
518 man->aptr[ansIdx] = gageAnswerPointer(gctx, pvl, item);
519 man->ispec[ansIdx].kind = pvl->kind;
520 man->ispec[ansIdx].item = item;
521 man->alen[ansIdx] = gageAnswerLength(gctx, pvl, item);
522 /* HEY hack: doing this tally and allocation only on the last time
523 we're called, presuming calls with increasing ansIdx */
524 man->slen = 0;
1
The value 0 is assigned to field 'slen'
→
525 if (ansIdx == man->anum-1) {
2
←
Taking true branch
→
526 unsigned int ai;
527 /* fprintf(stderr, "!%s: %s hello aidx reached anum-1 = %u\n", me, man->name, man->anum-1); */
528 for (ai=0; ai<man->anum; ai++) {
3
←
Loop condition is false. Execution continues on line 532
→
529 man->sidx[ai] = man->slen;
530 man->slen += man->alen[ai];
531 }
532 man->san = AIR_CALLOC(man->slen, double)(double*)(calloc((man->slen), sizeof(double)));
4
←
Within the expansion of the macro 'AIR_CALLOC':
a
Call to 'calloc' has an allocation size of 0 bytes
533 /*
534 fprintf(stderr, "!%s: (%s) slen = %u\n", "multiAnswerAdd",
535 pvl->kind->name, man->slen);
536 */
537 }
538 return;
539}
540
541static void
542multiAnswerCollect(multiAnswer *man) {
543 unsigned int ai;
544 /*
545 static const char me[]="multiAnswerCollect";
546
547 fprintf(stderr, "%s: %s hi\n", me, man->name);
548 */
549 for (ai=0; ai<man->anum; ai++) {
550 /*
551 fprintf(stderr, "!%s: (%s/%s) ai=%u/%u to %p+%u=%p for %u doubles\n", "multiAnswerCollect",
552 man->ispec[ai].kind->name,
553 airEnumStr(man->ispec[ai].kind->enm, man->ispec[ai].item),
554 ai, man->anum, man->san, man->sidx[ai], man->san + man->sidx[ai],
555 man->alen[ai]);
556 */
557 memcpy(man->san + man->sidx[ai],__builtin___memcpy_chk (man->san + man->sidx[ai], man->
aptr[ai], man->alen[ai]*sizeof(double), __builtin_object_size
(man->san + man->sidx[ai], 0))
558 man->aptr[ai],__builtin___memcpy_chk (man->san + man->sidx[ai], man->
aptr[ai], man->alen[ai]*sizeof(double), __builtin_object_size
(man->san + man->sidx[ai], 0))
559 man->alen[ai]*sizeof(double))__builtin___memcpy_chk (man->san + man->sidx[ai], man->
aptr[ai], man->alen[ai]*sizeof(double), __builtin_object_size
(man->san + man->sidx[ai], 0));
560 }
561 return;
562}
563
564static int
565multiAnswerCompare(multiAnswer *man1, multiAnswer *man2) {
566 static const char me[]="multiAnswerCompare";
567 unsigned int si, slen;
568
569
570#if 1
571 if (man1->slen != man2->slen) {
572 biffAddf(BKEY"probeSS", "%s: man1->slen %u != man2->slen %u\n", me,
573 man1->slen, man2->slen);
574 return 1;
575 }
576 slen = man1->slen;
577#else
578 slen = AIR_MIN(man1->slen, man2->slen)((man1->slen) < (man2->slen) ? (man1->slen) : (man2
->slen));
579#endif
580 for (si=0; si<slen; si++) {
581 if (man1->san[si] != man2->san[si]) {
582 /* HEY should track down which part of which answer,
583 in man1 and man2, is different, which was the
584 purpose of recording ispec and sidx */
585 biffAddf(BKEY"probeSS", "%s: man1->san[si] %.17g != man2->san[si] %.17g", me,
586 man1->san[si], man2->san[si]);
587 return 1;
588 }
589 }
590 return 0;
591}
592
593/*
594** setting up gageContexts for the first of the two tasks listed above: making
595** sure per-component information is handled correctly. NOTE: The combination
596** of function calls made here is very atypical for a Teem-using program
597*/
598static int
599updateQueryKernelSetTask1(gageContext *gctxComp[KIND_NUM4],
600 gageContext *gctx[KIND_NUM4], int gSetup,
601 multiAnswer *manComp[KIND_NUM4],
602 multiAnswer *man[KIND_NUM4],
603 NrrdKernel *kpack[3],
604 double support) {
605 static const char me[]="updateQueryKernelSetTask1";
606 double parm1[NRRD_KERNEL_PARMS_NUM8], parmV[NRRD_KERNEL_PARMS_NUM8];
607 unsigned int kindIdx;
608
609 if (4 != KIND_NUM4) {
610 biffAddf(BKEY"probeSS", "%s: sorry, confused: KIND_NUM %u != 4",
611 me, KIND_NUM4);
612 return 1;
613 }
614 parm1[0] = 1.0;
615 parmV[0] = support;
616 if (gSetup) {
617 for (kindIdx=0; kindIdx<KIND_NUM4; kindIdx++) {
618 if (0 /* (no-op for formatting) */
619 || gageKernelSet(gctxComp[kindIdx], gageKernel00, kpack[0], parm1)
620 || gageKernelSet(gctxComp[kindIdx], gageKernel11, kpack[1], parm1)
621 || gageKernelSet(gctxComp[kindIdx], gageKernel22, kpack[2], parm1)
622 || gageKernelSet(gctx[kindIdx], gageKernel00, kpack[0], parmV)
623 || gageKernelSet(gctx[kindIdx], gageKernel11, kpack[1], parmV)
624 || gageKernelSet(gctx[kindIdx], gageKernel22, kpack[2], parmV)) {
625 biffMovef(BKEY"probeSS", GAGEgageBiffKey, "%s: trouble setting kernel (kind %u)",
626 me, kindIdx);
627 return 1;
628 }
629 }
630 /* Note that the contexts for the diced-up volumes of components are always
631 of kind gageKindScl, and all the items are from the gageScl* enum */
632 for (kindIdx=0; kindIdx<KIND_NUM4; kindIdx++) {
633 unsigned int pvi;
634 /*
635 fprintf(stderr, "!%s: gctx[%u] has %u pvls\n", me, kindIdx, gctx[kindIdx]->pvlNum);
636 fprintf(stderr, "!%s: gctxComp[%u] has %u pvls\n", me, kindIdx, gctxComp[kindIdx]->pvlNum);
637 */
638 for (pvi=0; pvi<gctxComp[kindIdx]->pvlNum; pvi++) {
639 if (0 /* (no-op for formatting) */
640 || gageQueryItemOn(gctxComp[kindIdx], gctxComp[kindIdx]->pvl[pvi], gageSclValue)
641 || gageQueryItemOn(gctxComp[kindIdx], gctxComp[kindIdx]->pvl[pvi], gageSclGradVec)
642 || gageQueryItemOn(gctxComp[kindIdx], gctxComp[kindIdx]->pvl[pvi], gageSclHessian)) {
643 biffMovef(BKEY"probeSS", GAGEgageBiffKey, "%s: trouble setting query (kind %u) on pvi %u",
644 me, kindIdx, pvi);
645 return 1;
646 }
647 }
648 if (gageUpdate(gctxComp[kindIdx])) {
649 biffMovef(BKEY"probeSS", GAGEgageBiffKey, "%s: trouble updating comp gctx %u",
650 me, kindIdx);
651 return 1;
652 }
653 }
654 /* For the original contexts, we have to use the kind-specific items that
655 correspond to the values and derivatives */
656 if (0 /* (no-op for formatting) */
657 || gageQueryItemOn(gctx[KI_SCL0], gctx[KI_SCL0]->pvl[0], gageSclValue)
658 || gageQueryItemOn(gctx[KI_SCL0], gctx[KI_SCL0]->pvl[0], gageSclGradVec)
659 || gageQueryItemOn(gctx[KI_SCL0], gctx[KI_SCL0]->pvl[0], gageSclHessian)
660
661 || gageQueryItemOn(gctx[KI_VEC1], gctx[KI_VEC1]->pvl[0], gageVecVector)
662 || gageQueryItemOn(gctx[KI_VEC1], gctx[KI_VEC1]->pvl[0], gageVecJacobian)
663 || gageQueryItemOn(gctx[KI_VEC1], gctx[KI_VEC1]->pvl[0], gageVecHessian)
664
665 || gageQueryItemOn(gctx[KI_TEN2], gctx[KI_TEN2]->pvl[0], tenGageTensor)
666 || gageQueryItemOn(gctx[KI_TEN2], gctx[KI_TEN2]->pvl[0], tenGageTensorGrad)
667 || gageQueryItemOn(gctx[KI_TEN2], gctx[KI_TEN2]->pvl[0], tenGageS)
668 || gageQueryItemOn(gctx[KI_TEN2], gctx[KI_TEN2]->pvl[0], tenGageSGradVec)
669 || gageQueryItemOn(gctx[KI_TEN2], gctx[KI_TEN2]->pvl[0], tenGageHessian)
670 || gageQueryItemOn(gctx[KI_TEN2], gctx[KI_TEN2]->pvl[0], tenGageSHessian)
671
672 || gageQueryItemOn(gctx[KI_DWI3], gctx[KI_DWI3]->pvl[0], tenDwiGageAll)
673 || gageQueryItemOn(gctx[KI_DWI3], gctx[KI_DWI3]->pvl[0], tenDwiGageTensorLLS)) {
674 biffMovef(BKEY"probeSS", GAGEgageBiffKey, "%s: trouble setting item", me);
675 return 1;
676 }
677 for (kindIdx=0; kindIdx<KIND_NUM4; kindIdx++) {
678 if (gageUpdate(gctx[kindIdx])) {
679 biffMovef(BKEY"probeSS", GAGEgageBiffKey, "%s: trouble updating single gctx %u",
680 me, kindIdx);
681 return 1;
682 }
683 }
684 } /* if (gSetup) */
685
686 for (kindIdx=0; kindIdx<KIND_NUM4; kindIdx++) {
687 multiAnswerInit(man[kindIdx]);
688 multiAnswerInit(manComp[kindIdx]);
689 }
690 for (kindIdx=0; kindIdx<KIND_NUM4; kindIdx++) {
691 unsigned int pvi;
692 multiAnswerLenSet(manComp[kindIdx], (KI_DWI3 != kindIdx ? 3 : 1)*gctxComp[kindIdx]->pvlNum);
693 for (pvi=0; pvi<gctxComp[kindIdx]->pvlNum; pvi++) {
694 multiAnswerAdd(manComp[kindIdx], pvi,
695 gctxComp[kindIdx], gctxComp[kindIdx]->pvl[pvi],
696 gageSclValue);
697 }
698 if (KI_DWI3 != kindIdx) {
699 for (pvi=0; pvi<gctxComp[kindIdx]->pvlNum; pvi++) {
700 multiAnswerAdd(manComp[kindIdx], pvi + gctxComp[kindIdx]->pvlNum,
701 gctxComp[kindIdx], gctxComp[kindIdx]->pvl[pvi],
702 gageSclGradVec);
703 }
704 for (pvi=0; pvi<gctxComp[kindIdx]->pvlNum; pvi++) {
705 multiAnswerAdd(manComp[kindIdx], pvi + 2*gctxComp[kindIdx]->pvlNum,
706 gctxComp[kindIdx], gctxComp[kindIdx]->pvl[pvi],
707 gageSclHessian);
708 }
709 }
710 multiAnswerLenSet(man[kindIdx], KI_DWI3 != kindIdx ? 3 : 1);
711 switch(kindIdx) {
712 case KI_SCL0:
713 multiAnswerAdd(man[kindIdx], 0, gctx[kindIdx], gctx[kindIdx]->pvl[0], gageSclValue);
714 multiAnswerAdd(man[kindIdx], 1, gctx[kindIdx], gctx[kindIdx]->pvl[0], gageSclGradVec);
715 multiAnswerAdd(man[kindIdx], 2, gctx[kindIdx], gctx[kindIdx]->pvl[0], gageSclHessian);
716 break;
717 case KI_VEC1:
718 multiAnswerAdd(man[kindIdx], 0, gctx[kindIdx], gctx[kindIdx]->pvl[0], gageVecVector);
719 multiAnswerAdd(man[kindIdx], 1, gctx[kindIdx], gctx[kindIdx]->pvl[0], gageVecJacobian);
720 multiAnswerAdd(man[kindIdx], 2, gctx[kindIdx], gctx[kindIdx]->pvl[0], gageVecHessian);
721 break;
722 case KI_TEN2:
723 multiAnswerAdd(man[kindIdx], 0, gctx[kindIdx], gctx[kindIdx]->pvl[0], tenGageTensor);
724 multiAnswerAdd(man[kindIdx], 1, gctx[kindIdx], gctx[kindIdx]->pvl[0], tenGageTensorGrad);
725 multiAnswerAdd(man[kindIdx], 2, gctx[kindIdx], gctx[kindIdx]->pvl[0], tenGageHessian);
726 break;
727 case KI_DWI3:
728 multiAnswerAdd(man[kindIdx], 0, gctx[kindIdx], gctx[kindIdx]->pvl[0], tenDwiGageAll);
729 break;
730 }
731 }
732
733 return 0;
734}
735
736static int
737probeTask1(gageContext *gctxComp[KIND_NUM4],
738 gageContext *gctx[KIND_NUM4],
739 multiAnswer *manComp[KIND_NUM4],
740 multiAnswer *man[KIND_NUM4],
741 unsigned int probeNum, NrrdKernel *kpack[3],
742 double ksupport) {
743 static const char me[]="probeTask1";
744 unsigned int kindIdx, probeIdx, errNumMax, tenErrNum, sclErrNum;
745 double pos[3], upos[3], minp[3], maxp[3],
746 tenDiff[7], tenAvg[7], tenErr, tenErrMax,
747 vecDiff[3], vecAvg[3], vecErr, vecErrMax, vecErrNum,
748 sclDiff, sclAvg, sclErr, sclErrMax, errNumFrac;
749 const double *dwiTenEstP,
750 *tenTenP, *tenTenNormP, *tenTenNormGradP,
751 *sclSclP, *sclGradP, *vecVecP;
752
753 ELL_3V_SET(minp, ksupport, ksupport, ksupport)((minp)[0] = (ksupport), (minp)[1] = (ksupport), (minp)[2] = (
ksupport));
754 ELL_3V_SET(maxp, volSize[0]-1-ksupport, volSize[1]-1-ksupport, volSize[2]-1-ksupport)((maxp)[0] = (volSize[0]-1-ksupport), (maxp)[1] = (volSize[1]
-1-ksupport), (maxp)[2] = (volSize[2]-1-ksupport));
755
756 /* this is all for task 1.5 */
757 sclSclP = gageAnswerPointer(gctx[KI_SCL0], gctx[KI_SCL0]->pvl[0], gageSclValue);
758 sclGradP = gageAnswerPointer(gctx[KI_SCL0], gctx[KI_SCL0]->pvl[0], gageSclGradVec);
759 vecVecP = gageAnswerPointer(gctx[KI_VEC1], gctx[KI_VEC1]->pvl[0], gageVecVector);
760 tenTenP = gageAnswerPointer(gctx[KI_TEN2], gctx[KI_TEN2]->pvl[0], tenGageTensor);
761 tenTenNormP = gageAnswerPointer(gctx[KI_TEN2], gctx[KI_TEN2]->pvl[0], tenGageS);
762 tenTenNormGradP = gageAnswerPointer(gctx[KI_TEN2], gctx[KI_TEN2]->pvl[0], tenGageSGradVec);
763 dwiTenEstP = gageAnswerPointer(gctx[KI_DWI3], gctx[KI_DWI3]->pvl[0], tenDwiGageTensorLLS);
764 tenErrNum = sclErrNum = 0;
765 vecErrNum = 0.0;
766 errNumFrac = 0.02;
767
768 if (nrrdKernelBoxSupportDebug == kpack[0]) {
769 /* not actually here for any derivatives, mainly to check on
770 tensor esimation (in addition to usual check of multivariate as
771 set of scalars) */
772 tenErrMax = 0.0002;
773 vecErrMax = AIR_NAN(airFloatQNaN.f);
774 sclErrMax = AIR_NAN(airFloatQNaN.f);
775 } else if (nrrdKernelCos4SupportDebug == kpack[0]) {
776 tenErrMax = AIR_NAN(airFloatQNaN.f);
777 vecErrMax = AIR_NAN(airFloatQNaN.f);
778 sclErrMax = AIR_NAN(airFloatQNaN.f);
779 } else if (nrrdKernelCatmullRomSupportDebug == kpack[0]) {
780 tenErrMax = 0.14; /* honestly, not sure how meaningful this test
781 is, given how significant we're allowing the
782 error to be. might be more meaningful if
783 this was with comparison to a pre-computed
784 non-linear least-squares fit, instead of the
785 (log-based) linear least-squares. */
786 vecErrMax = 0.0016;
787 sclErrMax = 0.0008;
788 } else {
789 biffAddf(BKEY"probeSS", "%s: unexpected kpack[0] %s\n", me, kpack[0]->name);
790 return 1;
791 }
792 errNumMax = AIR_CAST(unsigned int, errNumFrac*probeNum)((unsigned int)(errNumFrac*probeNum));
793 for (probeIdx=0; probeIdx<probeNum; probeIdx++) {
794 airHalton(upos, probeIdx+HALTON_BASE100, airPrimeList + 0, 3);
795 pos[0] = AIR_AFFINE(0.0, upos[0], 1.0, minp[0], maxp[0])( ((double)(maxp[0])-(minp[0]))*((double)(upos[0])-(0.0)) / (
(double)(1.0)-(0.0)) + (minp[0]));
796 pos[1] = AIR_AFFINE(0.0, upos[1], 1.0, minp[1], maxp[1])( ((double)(maxp[1])-(minp[1]))*((double)(upos[1])-(0.0)) / (
(double)(1.0)-(0.0)) + (minp[1]));
797 pos[2] = AIR_AFFINE(0.0, upos[2], 1.0, minp[2], maxp[2])( ((double)(maxp[2])-(minp[2]))*((double)(upos[2])-(0.0)) / (
(double)(1.0)-(0.0)) + (minp[2]));
798 for (kindIdx=0; kindIdx<KIND_NUM4; kindIdx++) {
799#define PROBE(gctx, what, pos) \
800 if (gageProbeSpace((gctx), (pos)[0], (pos)[1], (pos)[2], \
801 AIR_TRUE1 /* indexSpace */, \
802 AIR_FALSE0 /* clamp */)) { \
803 biffAddf(BKEY"probeSS", "%s: probe (%s) error(%d): %s", me, \
804 what, (gctx)->errNum, (gctx)->errStr); \
805 return 1; \
806 }
807
808 PROBE(gctxComp[kindIdx], "comp", pos);
809 PROBE(gctx[kindIdx], "single", pos);
810
811#undef PROBE
812
813 multiAnswerCollect(man[kindIdx]);
814 multiAnswerCollect(manComp[kindIdx]);
815 if (multiAnswerCompare(manComp[kindIdx], man[kindIdx])) {
816 biffAddf(BKEY"probeSS", "%s: on point %u of kindIdx %u", me,
817 probeIdx, kindIdx);
818 return 1;
819 }
820 }
821
822 /* see if probed tensors mostly agree with
823 tensors estimated from probed DWIs */
824 if (AIR_EXISTS(tenErrMax)(((int)(!((tenErrMax) - (tenErrMax)))))) {
825 TEN_T_SUB(tenDiff, dwiTenEstP, tenTenP)( (tenDiff)[0] = ((dwiTenEstP)[0] + (tenTenP)[0])/2.0, (tenDiff
)[1] = (dwiTenEstP)[1] - (tenTenP)[1], (tenDiff)[2] = (dwiTenEstP
)[2] - (tenTenP)[2], (tenDiff)[3] = (dwiTenEstP)[3] - (tenTenP
)[3], (tenDiff)[4] = (dwiTenEstP)[4] - (tenTenP)[4], (tenDiff
)[5] = (dwiTenEstP)[5] - (tenTenP)[5], (tenDiff)[6] = (dwiTenEstP
)[6] - (tenTenP)[6]);
826 TEN_T_LERP(tenAvg, 0.5, dwiTenEstP, tenTenP)( (tenAvg)[0] = (((0.5))*(((tenTenP)[0]) - ((dwiTenEstP)[0]))
+ ((dwiTenEstP)[0])), (tenAvg)[1] = (((0.5))*(((tenTenP)[1])
- ((dwiTenEstP)[1])) + ((dwiTenEstP)[1])), (tenAvg)[2] = (((
0.5))*(((tenTenP)[2]) - ((dwiTenEstP)[2])) + ((dwiTenEstP)[2]
)), (tenAvg)[3] = (((0.5))*(((tenTenP)[3]) - ((dwiTenEstP)[3]
)) + ((dwiTenEstP)[3])), (tenAvg)[4] = (((0.5))*(((tenTenP)[4
]) - ((dwiTenEstP)[4])) + ((dwiTenEstP)[4])), (tenAvg)[5] = (
((0.5))*(((tenTenP)[5]) - ((dwiTenEstP)[5])) + ((dwiTenEstP)[
5])), (tenAvg)[6] = (((0.5))*(((tenTenP)[6]) - ((dwiTenEstP)[
6])) + ((dwiTenEstP)[6])));
827 tenErr = TEN_T_NORM(tenDiff)(sqrt(( (tenDiff)[1]*(tenDiff)[1] + 2*(tenDiff)[2]*(tenDiff)[
2] + 2*(tenDiff)[3]*(tenDiff)[3] + (tenDiff)[4]*(tenDiff)[4] +
2*(tenDiff)[5]*(tenDiff)[5] + (tenDiff)[6]*(tenDiff)[6] )))/TEN_T_NORM(tenAvg)(sqrt(( (tenAvg)[1]*(tenAvg)[1] + 2*(tenAvg)[2]*(tenAvg)[2] +
2*(tenAvg)[3]*(tenAvg)[3] + (tenAvg)[4]*(tenAvg)[4] + 2*(tenAvg
)[5]*(tenAvg)[5] + (tenAvg)[6]*(tenAvg)[6] )));
828 if (tenErr > tenErrMax) {
829 tenErrNum++;
830 if (tenErrNum > errNumMax) {
831 biffAddf(BKEY"probeSS", "%s: (probe %u) tenErr %g > %g too many times %u > %u",
832 me, probeIdx, tenErr, tenErrMax, tenErrNum, errNumMax);
833 return 1;
834 }
835 }
836 }
837
838 /* see if invariant gradient learned from tensor volume agrees with
839 volume of pre-computed invariant gradients, and with
840 gradient of pre-computed invariants */
841 if (AIR_EXISTS(vecErrMax)(((int)(!((vecErrMax) - (vecErrMax)))))) {
842 ELL_3V_SUB(vecDiff, sclGradP, vecVecP)((vecDiff)[0] = (sclGradP)[0] - (vecVecP)[0], (vecDiff)[1] = (
sclGradP)[1] - (vecVecP)[1], (vecDiff)[2] = (sclGradP)[2] - (
vecVecP)[2]);
843 ELL_3V_LERP(vecAvg, 0.5, sclGradP, vecVecP)((vecAvg)[0] = (((0.5))*(((vecVecP)[0]) - ((sclGradP)[0])) + (
(sclGradP)[0])), (vecAvg)[1] = (((0.5))*(((vecVecP)[1]) - ((sclGradP
)[1])) + ((sclGradP)[1])), (vecAvg)[2] = (((0.5))*(((vecVecP)
[2]) - ((sclGradP)[2])) + ((sclGradP)[2])));
844 vecErr = ELL_3V_LEN(vecDiff)(sqrt((((vecDiff))[0]*((vecDiff))[0] + ((vecDiff))[1]*((vecDiff
))[1] + ((vecDiff))[2]*((vecDiff))[2])))/sqrt(ELL_3V_LEN(vecAvg)(sqrt((((vecAvg))[0]*((vecAvg))[0] + ((vecAvg))[1]*((vecAvg))
[1] + ((vecAvg))[2]*((vecAvg))[2]))));
845 if (vecErr > vecErrMax) {
846 vecErrNum += 0.5;
847 if (vecErrNum > errNumMax) {
848 biffAddf(BKEY"probeSS", "%s: (probe %u) (A) vecErr %g > %g too many times %g > %u",
849 me, probeIdx, vecErr, vecErrMax, vecErrNum, errNumMax);
850 return 1;
851 }
852 }
853 ELL_3V_SUB(vecDiff, tenTenNormGradP, vecVecP)((vecDiff)[0] = (tenTenNormGradP)[0] - (vecVecP)[0], (vecDiff
)[1] = (tenTenNormGradP)[1] - (vecVecP)[1], (vecDiff)[2] = (tenTenNormGradP
)[2] - (vecVecP)[2]);
854 ELL_3V_LERP(vecAvg, 0.5, tenTenNormGradP, vecVecP)((vecAvg)[0] = (((0.5))*(((vecVecP)[0]) - ((tenTenNormGradP)[
0])) + ((tenTenNormGradP)[0])), (vecAvg)[1] = (((0.5))*(((vecVecP
)[1]) - ((tenTenNormGradP)[1])) + ((tenTenNormGradP)[1])), (vecAvg
)[2] = (((0.5))*(((vecVecP)[2]) - ((tenTenNormGradP)[2])) + (
(tenTenNormGradP)[2])));
855 vecErr = ELL_3V_LEN(vecDiff)(sqrt((((vecDiff))[0]*((vecDiff))[0] + ((vecDiff))[1]*((vecDiff
))[1] + ((vecDiff))[2]*((vecDiff))[2])))/sqrt(ELL_3V_LEN(vecAvg)(sqrt((((vecAvg))[0]*((vecAvg))[0] + ((vecAvg))[1]*((vecAvg))
[1] + ((vecAvg))[2]*((vecAvg))[2]))));
856 if (vecErr > vecErrMax) {
857 vecErrNum += 0.5;
858 if (vecErrNum > errNumMax) {
859 biffAddf(BKEY"probeSS", "%s: (probe %u) (B) vecErr %g > %g too many times %g > %u",
860 me, probeIdx, vecErr, vecErrMax, vecErrNum, errNumMax);
861 return 1;
862 }
863 }
864 }
865
866 /* see if invariant learned from tensor volume agrees with
867 volume of precomputed invariants */
868 if (AIR_EXISTS(sclErrMax)(((int)(!((sclErrMax) - (sclErrMax)))))) {
869 sclDiff = sclSclP[0] - tenTenNormP[0];
870 sclAvg = (sclSclP[0] + tenTenNormP[0])/2;
871 sclErr = AIR_ABS(sclDiff)((sclDiff) > 0.0f ? (sclDiff) : -(sclDiff))/sqrt(AIR_ABS(sclAvg)((sclAvg) > 0.0f ? (sclAvg) : -(sclAvg)));
872 if (sclErr > sclErrMax) {
873 sclErrNum++;
874 if (sclErrNum > errNumMax) {
875 biffAddf(BKEY"probeSS", "%s: (probe %u) (B) sclErr %g > %g too many times %u > %u",
876 me, probeIdx, sclErr, sclErrMax, sclErrNum, errNumMax);
877 return 1;
878 }
879 }
880 }
881 }
882
883 return 0;
884}
885
886static const char *testInfo = "for testing gage";
887
888int
889main(int argc, const char **argv) {
890 const char *me;
891 char *err = NULL((void*)0);
892 hestOpt *hopt=NULL((void*)0);
893 hestParm *hparm;
894 airArray *mop;
895
896 const gageKind *kind[KIND_NUM4];
897 char name[KIND_NUM4][AIR_STRLEN_SMALL(128+1)] = { "scl", "vec", "ten", "dwi" };
898 char nameComp[KIND_NUM4][AIR_STRLEN_SMALL(128+1)] = { "sclComp", "vecComp", "tenComp", "dwiComp" };
899 char *kernS;
900 gageKind *dwikind = NULL((void*)0);
901 gageContext *gctxComp[KIND_NUM4], *gctxCompCopy[KIND_NUM4],
902 *gctx[KIND_NUM4], *gctxCopy[KIND_NUM4];
903 multiAnswer *manComp[KIND_NUM4], *man[KIND_NUM4];
904 Nrrd *nin[KIND_NUM4],
905 /* these are the volumes that are used in gctxComp[] */
906 *nsclCopy, *nvecComp[3], *nctenComp[7], **ndwiComp,
907 *ngrad; /* need access to list of gradients used to make DWIs;
908 (this is not the gradient of a scalar field) */
909 double bval = 1000, noiseStdv=0.0001, ksupport;
910 unsigned int kindIdx, probeNum,
911 gradNum = 10; /* small number so that missing one will produce
912 a big reconstruction error */
913 NrrdKernel *kpack[3];
914
915 kind[0] = gageKindScl;
916 kind[1] = gageKindVec;
917 kind[2] = tenGageKind;
918 kind[3] = NULL((void*)0); /* dwi */
919
920 me = argv[0];
921 mop = airMopNew();
922 hparm = hestParmNew();
923 airMopAdd(mop, hparm, (airMopper)hestParmFree, airMopAlways);
924 /* learn things from hest */
925 hestOptAdd(&hopt, "supp", "r", airTypeDouble, 1, 1, &ksupport, "1.0",
926 "kernel support");
927 hestOptAdd(&hopt, "pnum", "N", airTypeUInt, 1, 1, &probeNum, "100",
928 "# of probes");
929 hestOptAdd(&hopt, "k", "kern", airTypeString, 1, 1, &kernS, NULL((void*)0),
930 "kernel to use; can be: box, cos, or ctmr");
931 hestParseOrDie(hopt, argc-1, argv+1, hparm, me, testInfo,
932 AIR_TRUE1, AIR_TRUE1, AIR_TRUE1);
933 airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways);
934 airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways);
935
936 if (!strcmp(kernS, "box")) {
937 ELL_3V_SET(kpack,((kpack)[0] = (nrrdKernelBoxSupportDebug), (kpack)[1] = (nrrdKernelZero
), (kpack)[2] = (nrrdKernelZero))
938 nrrdKernelBoxSupportDebug,((kpack)[0] = (nrrdKernelBoxSupportDebug), (kpack)[1] = (nrrdKernelZero
), (kpack)[2] = (nrrdKernelZero))
939 nrrdKernelZero,((kpack)[0] = (nrrdKernelBoxSupportDebug), (kpack)[1] = (nrrdKernelZero
), (kpack)[2] = (nrrdKernelZero))
940 nrrdKernelZero)((kpack)[0] = (nrrdKernelBoxSupportDebug), (kpack)[1] = (nrrdKernelZero
), (kpack)[2] = (nrrdKernelZero));
941 } else if (!strcmp(kernS, "cos")) {
942 ELL_3V_SET(kpack,((kpack)[0] = (nrrdKernelCos4SupportDebug), (kpack)[1] = (nrrdKernelCos4SupportDebugD
), (kpack)[2] = (nrrdKernelCos4SupportDebugDD))
943 nrrdKernelCos4SupportDebug,((kpack)[0] = (nrrdKernelCos4SupportDebug), (kpack)[1] = (nrrdKernelCos4SupportDebugD
), (kpack)[2] = (nrrdKernelCos4SupportDebugDD))
944 nrrdKernelCos4SupportDebugD,((kpack)[0] = (nrrdKernelCos4SupportDebug), (kpack)[1] = (nrrdKernelCos4SupportDebugD
), (kpack)[2] = (nrrdKernelCos4SupportDebugDD))
945 nrrdKernelCos4SupportDebugDD)((kpack)[0] = (nrrdKernelCos4SupportDebug), (kpack)[1] = (nrrdKernelCos4SupportDebugD
), (kpack)[2] = (nrrdKernelCos4SupportDebugDD));
946 } else if (!strcmp(kernS, "ctmr")) {
947 ELL_3V_SET(kpack,((kpack)[0] = (nrrdKernelCatmullRomSupportDebug), (kpack)[1] =
(nrrdKernelCatmullRomSupportDebugD), (kpack)[2] = (nrrdKernelCatmullRomSupportDebugDD
))
948 nrrdKernelCatmullRomSupportDebug,((kpack)[0] = (nrrdKernelCatmullRomSupportDebug), (kpack)[1] =
(nrrdKernelCatmullRomSupportDebugD), (kpack)[2] = (nrrdKernelCatmullRomSupportDebugDD
))
949 nrrdKernelCatmullRomSupportDebugD,((kpack)[0] = (nrrdKernelCatmullRomSupportDebug), (kpack)[1] =
(nrrdKernelCatmullRomSupportDebugD), (kpack)[2] = (nrrdKernelCatmullRomSupportDebugDD
))
950 nrrdKernelCatmullRomSupportDebugDD)((kpack)[0] = (nrrdKernelCatmullRomSupportDebug), (kpack)[1] =
(nrrdKernelCatmullRomSupportDebugD), (kpack)[2] = (nrrdKernelCatmullRomSupportDebugDD
));
951 } else {
952 fprintf(stderr__stderrp, "%s: kernel \"%s\" not recognized", me, kernS);
953 airMopError(mop); return 1;
954 }
955 fprintf(stderr__stderrp, "kpack = %s, %s, %s\n", kpack[0]->name, kpack[1]->name, kpack[2]->name);
956
957 /* This was a tricky bug: Adding gageContextNix(gctx) to the mop
958 as soon as a gctx is created makes sense. But, in the first
959 version of this code, gageContextNix was added to the mop
960 BEFORE tenDwiGageKindNix was added, which meant that gageContextNix
961 was being called AFTER tenDwiGageKindNix. However, that meant
962 that existing pvls had their pvl->kind being freed from under them,
963 but gagePerVolumeNix certainly needs to look at and possibly call
964 pvl->kind->pvlDataNix in order clean up pvl->data. *Therefore*,
965 we have to add tenDwiGageKindNix to the mop prior to gageContextNix,
966 even though nothing about the (dyanamic) kind has been set yet.
967 If we had something like smart pointers, then tenDwiGageKindNix
968 would not have freed something that an extant pvl was using */
969 dwikind = tenDwiGageKindNew();
970 airMopAdd(mop, dwikind, (airMopper)tenDwiGageKindNix, airMopAlways);
971
972#define GAGE_CTX_NEW(gg, mm) \
973 (gg) = gageContextNew(); \
974 airMopAdd((mm), (gg), (airMopper)gageContextNix, airMopAlways); \
975 gageParmSet((gg), gageParmRenormalize, AIR_FALSE0); \
976 gageParmSet((gg), gageParmCheckIntegrals, AIR_TRUE1)
977
978 for (kindIdx=0; kindIdx<KIND_NUM4; kindIdx++) {
979 GAGE_CTX_NEW(gctx[kindIdx], mop);
980 GAGE_CTX_NEW(gctxComp[kindIdx], mop);
981 /*
982 fprintf(stderr, "%s: kind %s (%u): gctx=%p, gctxComp=%p\n", me,
983 kindStr[kindIdx], kindIdx,
984 AIR_VOIDP(gctx[kindIdx]),
985 AIR_VOIDP(gctxComp[kindIdx]));
986 */
987 man[kindIdx] = multiAnswerNew(name[kindIdx]);
988 manComp[kindIdx] = multiAnswerNew(nameComp[kindIdx]);
989 airMopAdd(mop, man[kindIdx], (airMopper)multiAnswerNix, airMopAlways);
990 airMopAdd(mop, manComp[kindIdx], (airMopper)multiAnswerNix, airMopAlways);
991 }
992#undef GAGE_CTX_NEW
993
994 for (kindIdx=0; kindIdx<KIND_NUM4; kindIdx++) {
995 NRRD_NEW(nin[kindIdx], mop);
996 }
997 NRRD_NEW(ngrad, mop);
998 NRRD_NEW(nsclCopy, mop);
999 if (engageGenTensor(gctx[KI_TEN2], nin[KI_TEN2],
1000 /* out/in */
1001 noiseStdv, 42 /* RNG seed, could be a command-line option */,
1002 volSize[0], volSize[1], volSize[2])
1003 || engageGenScalar(gctx[KI_SCL0], nin[KI_SCL0],
1004 gctxComp[KI_SCL0], nsclCopy,
1005 /* out/in */
1006 nin[KI_TEN2])
1007 || engageGenVector(gctx[KI_VEC1], nin[KI_VEC1],
1008 /* out/in */
1009 nin[KI_SCL0])
1010 /* engage'ing of nin[KI_DWI] done below */
1011 || genDwi(nin[KI_DWI3], ngrad,
1012 /* out/in */
1013 gradNum /* for B0 */, bval, nin[KI_TEN2])
1014 || engageMopDiceVector(gctxComp[KI_VEC1], nvecComp,
1015 /* out/in */
1016 mop, nin[KI_VEC1])
1017 || engageMopDiceTensor(gctxComp[KI_TEN2], nctenComp,
1018 /* out/in */
1019 mop, nin[KI_TEN2])
1020 || engageMopDiceDwi(gctxComp[KI_DWI3], &ndwiComp,
1021 /* out/in */
1022 mop, nin[KI_DWI3])) {
1023 airMopAdd(mop, err = biffGetDone(BKEY"probeSS"), airFree, airMopAlways);
1024 fprintf(stderr__stderrp, "%s: trouble creating volumes:\n%s", me, err);
1025 airMopError(mop); return 1;
1026 }
1027 if (tenDwiGageKindSet(dwikind, -1 /* thresh */, 0 /* soft */,
1028 bval, 1 /* valueMin */,
1029 ngrad, NULL((void*)0),
1030 tenEstimate1MethodLLS,
1031 tenEstimate2MethodPeled,
1032 424242)) {
1033 airMopAdd(mop, err = biffGetDone(TENtenBiffKey), airFree, airMopAlways);
1034 fprintf(stderr__stderrp, "%s: trouble creating DWI kind:\n%s", me, err);
1035 airMopError(mop); return 1;
1036 }
1037 /* access through kind[] is const, but not through dwikind */
1038 kind[KI_DWI3] = dwikind;
1039 /* engage dwi vol */
1040 {
1041 gagePerVolume *pvl;
1042 if ( !(pvl = gagePerVolumeNew(gctx[KI_DWI3], nin[KI_DWI3], kind[KI_DWI3]))
1043 || gagePerVolumeAttach(gctx[KI_DWI3], pvl) ) {
1044 airMopAdd(mop, err = biffGetDone(GAGEgageBiffKey), airFree, airMopAlways);
1045 fprintf(stderr__stderrp, "%s: trouble creating DWI context:\n%s", me, err);
1046 airMopError(mop); return 1;
1047 }
1048 /*
1049 fprintf(stderr, "%s: %s gctx=%p gets pvl %p\n", me,
1050 kindStr[KI_DWI], AIR_VOIDP(gctx[KI_DWI]), AIR_VOIDP(pvl));
1051 */
1052 }
1053
1054 /* make sure kinds can parse back to themselves */
1055 /* the messiness here is in part because of problems with the gage
1056 API, and the weirdness of how meetGageKindParse is either allocating
1057 something, or not, depending on its input. There is no way to
1058 refer to the "name" of a dwi kind without having allocated something. */
1059 for (kindIdx=0; kindIdx<KIND_NUM4; kindIdx++) {
1060 if (!(kind[kindIdx]->dynamicAlloc)) {
1061 if (kind[kindIdx] != meetGageKindParse(kind[kindIdx]->name)) {
1062 fprintf(stderr__stderrp, "%s: static kind[%u]->name %s wasn't parsed\n", me,
1063 kindIdx, kind[kindIdx]->name);
1064 airMopError(mop); return 1;
1065 }
1066 } else {
1067 gageKind *ktmp;
1068 ktmp = meetGageKindParse(kind[kindIdx]->name);
1069 if (!ktmp) {
1070 fprintf(stderr__stderrp, "%s: dynamic kind[%u]->name %s wasn't parsed\n", me,
1071 kindIdx, kind[kindIdx]->name);
1072 airMopError(mop); return 1;
1073 }
1074 if (airStrcmp(ktmp->name, kind[kindIdx]->name)) {
1075 fprintf(stderr__stderrp, "%s: parsed dynamic kind[%u]->name %s didn't match %s\n", me,
1076 kindIdx, ktmp->name, kind[kindIdx]->name);
1077 airMopError(mop); return 1;
1078 }
1079 if (!airStrcmp(TEN_DWI_GAGE_KIND_NAME"dwi", ktmp->name)) {
1080 /* HEY: there needs to be a generic way of freeing
1081 a dynamic kind, perhaps a new nixer field in gageKind */
1082 ktmp = tenDwiGageKindNix(ktmp);
1083 }
1084 }
1085 }
1086 /*
1087 for (kindIdx=0; kindIdx<KIND_NUM; kindIdx++) {
1088 printf("%s: %s (%s) -> (%u) %u\n", me, kind[kindIdx]->name,
1089 kind[kindIdx]->dynamicAlloc ? "dynamic" : "static",
1090 kind[kindIdx]->baseDim, kind[kindIdx]->valLen);
1091 }
1092 */
1093
1094 /*
1095 nrrdSave("tmp-cten.nrrd", nin[KI_TEN], NULL);
1096 nrrdSave("tmp-scl.nrrd", nin[KI_SCL], NULL);
1097 nrrdSave("tmp-vec.nrrd", nin[KI_VEC], NULL);
1098 nrrdSave("tmp-dwi.nrrd", nin[KI_DWI], NULL);
1099 */
1100
1101 /* ========================== TASK 1 */
1102 if (updateQueryKernelSetTask1(gctxComp, gctx, AIR_TRUE1, manComp, man, kpack, ksupport)
1103 || probeTask1(gctxComp, gctx, manComp, man, probeNum, kpack, ksupport)) {
1104 airMopAdd(mop, err = biffGetDone(BKEY"probeSS"), airFree, airMopAlways);
1105 fprintf(stderr__stderrp, "%s: trouble (orig, orig):\n%s", me, err);
1106 airMopError(mop); return 1;
1107 }
1108 /* testing gageContextCopy on multi-variate volumes */
1109 for (kindIdx=0; kindIdx<KIND_NUM4; kindIdx++) {
1110 if (!( gctxCopy[kindIdx] = gageContextCopy(gctx[kindIdx]) )) {
1111 airMopAdd(mop, err = biffGetDone(GAGEgageBiffKey), airFree, airMopAlways);
1112 fprintf(stderr__stderrp, "%s: trouble w/ multi-variate contexts:\n%s", me, err);
1113 airMopError(mop); return 1;
1114 }
1115 airMopAdd(mop, gctxCopy[kindIdx], (airMopper)gageContextNix, airMopAlways);
1116 }
1117 if (updateQueryKernelSetTask1(gctxComp, gctxCopy, AIR_FALSE0, manComp, man, kpack, ksupport)
1118 || probeTask1(gctxComp, gctxCopy, manComp, man, probeNum, kpack, ksupport)) {
1119 airMopAdd(mop, err = biffGetDone(BKEY"probeSS"), airFree, airMopAlways);
1120 fprintf(stderr__stderrp, "%s: trouble (orig, copy):\n%s", me, err);
1121 airMopError(mop); return 1;
1122 }
1123 for (kindIdx=0; kindIdx<KIND_NUM4; kindIdx++) {
1124 if (!( gctxCompCopy[kindIdx] = gageContextCopy(gctxComp[kindIdx]) )) {
1125 airMopAdd(mop, err = biffGetDone(GAGEgageBiffKey), airFree, airMopAlways);
1126 fprintf(stderr__stderrp, "%s: trouble w/ component contexts:\n%s", me, err);
1127 airMopError(mop); return 1;
1128 }
1129 airMopAdd(mop, gctxCompCopy[kindIdx], (airMopper)gageContextNix, airMopAlways);
1130 }
1131 if (updateQueryKernelSetTask1(gctxCompCopy, gctxCopy, AIR_FALSE0, manComp, man, kpack, ksupport)
1132 || probeTask1(gctxCompCopy, gctxCopy, manComp, man, probeNum, kpack, ksupport)) {
1133 airMopAdd(mop, err = biffGetDone(BKEY"probeSS"), airFree, airMopAlways);
1134 fprintf(stderr__stderrp, "%s: trouble (copy, copy):\n%s", me, err);
1135 airMopError(mop); return 1;
1136 }
1137
1138 /* ========================== TASK 2 */
1139 /* create scale-space stacks with tent, ctmr, and hermite */
1140 /* save them, free, and read them back in */
1141 /* gageContextCopy on stack */
1142 /* pick a scale in-between tau samples */
1143 /* for all the tau's half-way between tau samples in scale:
1144 blur at that tau to get correct values
1145 check that error with hermite is lower than ctmr is lower than tent */
1146 /* for all tau samples:
1147 blur at that tau to (re-)get correct values
1148 check that everything agrees there */
1149
1150
1151 /* single probe with high verbosity */
1152
1153 airMopOkay(mop);
1154 return 0;
1155}
1156#undef NRRD_NEW