Bug Summary

File:src/ten/epireg.c
Location:line 224, column 5
Description:Value stored to 'range' 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 "ten.h"
25#include "privateTen.h"
26
27int
28_tenEpiRegSave(const char *fname, Nrrd *nsingle, Nrrd **nmulti,
29 int len, const char *desc) {
30 static const char me[]="_tenEpiRegSave";
31 Nrrd *nout;
32 airArray *mop;
33
34 mop = airMopNew();
35 if (nsingle) {
36 nout = nsingle;
37 } else {
38 airMopAdd(mop, nout=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
39 if (nrrdJoin(nout, (const Nrrd*const*)nmulti, len, 0, AIR_TRUE1)) {
40 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: couldn't join %s for output", me, desc);
41 airMopError(mop); return 1;
42 }
43 }
44 if (nrrdSave(fname, nout, NULL((void*)0))) {
45 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: trouble saving %s to \"%s\"", me, desc, fname);
46 airMopError(mop); return 1;
47 }
48 fprintf(stderr__stderrp, "%s: saved %s to \"%s\"\n", me, desc, fname);
49 airMopOkay(mop);
50 return 0;
51}
52
53
54int
55_tenEpiRegCheck(Nrrd **nout, Nrrd **ndwi, unsigned int dwiLen, Nrrd *ngrad,
56 int reference,
57 double bwX, double bwY, double DWthr,
58 const NrrdKernel *kern, double *kparm) {
59 static const char me[]="_tenEpiRegCheck";
60 unsigned int ni;
61
62 AIR_UNUSED(DWthr)(void)(DWthr);
63 if (!( nout && ndwi && ngrad && kern && kparm )) {
64 biffAddf(TENtenBiffKey, "%s: got NULL pointer", me);
65 return 1;
66 }
67 if (tenGradientCheck(ngrad, nrrdTypeDefault, 6 /* <-- HEY: not sure */)) {
68 biffAddf(TENtenBiffKey, "%s: problem with given gradient list", me);
69 return 1;
70 }
71 if (dwiLen != ngrad->axis[1].size) {
72 char stmp[AIR_STRLEN_SMALL(128+1)];
73 biffAddf(TENtenBiffKey, "%s: got %u DWIs, but %s gradient directions", me,
74 dwiLen, airSprintSize_t(stmp, ngrad->axis[1].size));
75 return 1;
76 }
77 for (ni=0; ni<dwiLen; ni++) {
78 if (!nout[ni]) {
79 biffAddf(TENtenBiffKey, "%s: nout[%d] is NULL", me, ni);
80 return 1;
81 }
82 if (nrrdCheck(ndwi[ni])) {
83 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey,
84 "%s: basic nrrd validity failed on ndwi[%d]", me, ni);
85 return 1;
86 }
87 if (!nrrdSameSize(ndwi[0], ndwi[ni], AIR_TRUE1)) {
88 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: ndwi[%d] is different from ndwi[0]", me, ni);
89 return 1;
90 }
91 }
92 if (!( 3 == ndwi[0]->dim )) {
93 biffAddf(TENtenBiffKey, "%s: didn't get a set of 3-D arrays (got %d-D)", me,
94 ndwi[0]->dim);
95 return 1;
96 }
97 if (!( AIR_IN_CL(-1, reference, (int)dwiLen-1)((-1) <= (reference) && (reference) <= ((int)dwiLen
-1)) )) {
98 biffAddf(TENtenBiffKey, "%s: reference index %d not in valid range [-1,%d]",
99 me, reference, dwiLen-1);
100 return 1;
101 }
102 if (!( AIR_EXISTS(bwX)(((int)(!((bwX) - (bwX))))) && AIR_EXISTS(bwY)(((int)(!((bwY) - (bwY))))) )) {
103 biffAddf(TENtenBiffKey, "%s: bwX, bwY don't both exist", me);
104 return 1;
105 }
106 if (!( bwX >= 0 && bwY >= 0 )) {
107 biffAddf(TENtenBiffKey, "%s: bwX (%g) and bwY (%g) are not both non-negative",
108 me, bwX, bwY);
109 return 1;
110 }
111 return 0;
112}
113
114/*
115** this assumes that all nblur[i] are valid nrrds, and does nothing
116** to manage them
117*/
118int
119_tenEpiRegBlur(Nrrd **nblur, Nrrd **ndwi, unsigned int dwiLen,
120 double bwX, double bwY, int verb) {
121 static const char me[]="_tenEpiRegBlur";
122 NrrdResampleInfo *rinfo;
123 airArray *mop;
124 size_t ni, sx, sy, sz;
125 double savemin[2], savemax[2];
126
127 if (!( bwX || bwY )) {
128 if (verb) {
129 fprintf(stderr__stderrp, "%s:\n ", me); fflush(stderr__stderrp);
130 }
131 for (ni=0; ni<dwiLen; ni++) {
132 if (verb) {
133 fprintf(stderr__stderrp, "%2u ", (unsigned int)ni); fflush(stderr__stderrp);
134 }
135 if (nrrdCopy(nblur[ni], ndwi[ni])) {
136 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: trouble copying ndwi[%u]",
137 me, (unsigned int)ni);
138 return 1;
139 }
140 }
141 if (verb) {
142 fprintf(stderr__stderrp, "done\n");
143 }
144 return 0;
145 }
146 /* else we need to blur */
147 sx = ndwi[0]->axis[0].size;
148 sy = ndwi[0]->axis[1].size;
149 sz = ndwi[0]->axis[2].size;
150 mop = airMopNew();
151 rinfo = nrrdResampleInfoNew();
152 airMopAdd(mop, rinfo, (airMopper)nrrdResampleInfoNix, airMopAlways);
153 if (bwX) {
154 rinfo->kernel[0] = nrrdKernelGaussian;
155 rinfo->parm[0][0] = bwX;
156 rinfo->parm[0][1] = 3.0; /* how many stnd devs do we cut-off at */
157 } else {
158 rinfo->kernel[0] = NULL((void*)0);
159 }
160 if (bwY) {
161 rinfo->kernel[1] = nrrdKernelGaussian;
162 rinfo->parm[1][0] = bwY;
163 rinfo->parm[1][1] = 3.0; /* how many stnd devs do we cut-off at */
164 } else {
165 rinfo->kernel[1] = NULL((void*)0);
166 }
167 rinfo->kernel[2] = NULL((void*)0);
168 ELL_3V_SET(rinfo->samples, sx, sy, sz)((rinfo->samples)[0] = (sx), (rinfo->samples)[1] = (sy)
, (rinfo->samples)[2] = (sz));
169 ELL_3V_SET(rinfo->min, 0, 0, 0)((rinfo->min)[0] = (0), (rinfo->min)[1] = (0), (rinfo->
min)[2] = (0));
170 ELL_3V_SET(rinfo->max, sx-1, sy-1, sz-1)((rinfo->max)[0] = (sx-1), (rinfo->max)[1] = (sy-1), (rinfo
->max)[2] = (sz-1));
171 rinfo->boundary = nrrdBoundaryBleed;
172 rinfo->type = nrrdTypeDefault;
173 rinfo->renormalize = AIR_TRUE1;
174 rinfo->clamp = AIR_TRUE1;
175 if (verb) {
176 fprintf(stderr__stderrp, "%s:\n ", me); fflush(stderr__stderrp);
177 }
178 for (ni=0; ni<dwiLen; ni++) {
179 if (verb) {
180 fprintf(stderr__stderrp, "%2u ", (unsigned int)ni); fflush(stderr__stderrp);
181 }
182 savemin[0] = ndwi[ni]->axis[0].min; savemax[0] = ndwi[ni]->axis[0].max;
183 savemin[1] = ndwi[ni]->axis[1].min; savemax[1] = ndwi[ni]->axis[1].max;
184 ndwi[ni]->axis[0].min = 0; ndwi[ni]->axis[0].max = sx-1;
185 ndwi[ni]->axis[1].min = 0; ndwi[ni]->axis[1].max = sy-1;
186 if (nrrdSpatialResample(nblur[ni], ndwi[ni], rinfo)) {
187 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: trouble blurring ndwi[%u]",
188 me, (unsigned int)ni);
189 airMopError(mop); return 1;
190 }
191 ndwi[ni]->axis[0].min = savemin[0]; ndwi[ni]->axis[0].max = savemax[0];
192 ndwi[ni]->axis[1].min = savemin[1]; ndwi[ni]->axis[1].max = savemax[1];
193 }
194 if (verb) {
195 fprintf(stderr__stderrp, "done\n");
196 }
197 airMopOkay(mop);
198 return 0;
199}
200
201int
202_tenEpiRegThresholdFind(double *DWthrP, Nrrd **nin, int ninLen,
203 int save, double expo) {
204 static const char me[]="_tenEpiRegThresholdFind";
205 Nrrd *nhist, *ntmp;
206 airArray *mop;
207 int ni, bins, E;
208 double min=0, max=0;
209 NrrdRange *range;
210
211 mop = airMopNew();
212 airMopAdd(mop, nhist=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
213 airMopAdd(mop, ntmp=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
214
215 for (ni=0; ni<ninLen; ni++) {
216 range = nrrdRangeNewSet(nin[ni], nrrdBlind8BitRangeFalse);
217 if (!ni) {
218 min = range->min;
219 max = range->max;
220 } else {
221 min = AIR_MIN(min, range->min)((min) < (range->min) ? (min) : (range->min));
222 max = AIR_MAX(max, range->max)((max) > (range->max) ? (max) : (range->max));
223 }
224 range = nrrdRangeNix(range);
Value stored to 'range' is never read
225 }
226 bins = AIR_MIN(1024, (int)(max - min + 1))((1024) < ((int)(max - min + 1)) ? (1024) : ((int)(max - min
+ 1)));
227 ntmp->axis[0].min = min;
228 ntmp->axis[0].max = max;
229 for (ni=0; ni<ninLen; ni++) {
230 if (nrrdHisto(ntmp, nin[ni], NULL((void*)0), NULL((void*)0), bins, nrrdTypeFloat)) {
231 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey,
232 "%s: problem forming histogram of DWI %d", me, ni);
233 airMopError(mop); return 1;
234 }
235 if (!ni) {
236 E = nrrdCopy(nhist, ntmp);
237 } else {
238 E = nrrdArithBinaryOp(nhist, nrrdBinaryOpAdd, nhist, ntmp);
239 }
240 if (E) {
241 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey,
242 "%s: problem updating histogram sum on DWI %d",
243 me, ni);
244 airMopError(mop); return 1;
245 }
246 }
247 if (save) {
248 nrrdSave("regtmp-dwihist.nrrd", nhist, NULL((void*)0));
249 }
250 /*
251 if (_tenFindValley(DWthrP, nhist, 0.65, save)) {
252 biffAddf(TEN, "%s: problem finding DWI histogram valley", me);
253 airMopError(mop); return 1;
254 }
255 */
256 if (nrrdHistoThresholdOtsu(DWthrP, nhist, expo)) {
257 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: problem finding DWI threshold", me);
258 airMopError(mop); return 1;
259 }
260
261 airMopOkay(mop);
262 return 0;
263}
264
265int
266_tenEpiRegThreshold(Nrrd **nthresh, Nrrd **nblur, unsigned int ninLen,
267 double DWthr, int verb, int progress, double expo) {
268 static const char me[]="_tenEpiRegThreshold";
269 airArray *mop;
270 size_t I, sx, sy, sz, ni;
271 float val;
272 unsigned char *thr;
273
274 if (!( AIR_EXISTS(DWthr)(((int)(!((DWthr) - (DWthr))))) )) {
275 if (_tenEpiRegThresholdFind(&DWthr, nblur, ninLen, progress, expo)) {
276 biffAddf(TENtenBiffKey, "%s: trouble with automatic threshold determination", me);
277 return 1;
278 }
279 fprintf(stderr__stderrp, "%s: using %g for DWI threshold\n", me, DWthr);
280 }
281
282 mop = airMopNew();
283 if (verb) {
284 fprintf(stderr__stderrp, "%s:\n ", me); fflush(stderr__stderrp);
285 }
286 sx = nblur[0]->axis[0].size;
287 sy = nblur[0]->axis[1].size;
288 sz = nblur[0]->axis[2].size;
289 for (ni=0; ni<ninLen; ni++) {
290 if (verb) {
291 fprintf(stderr__stderrp, "%2u ", (unsigned int)ni); fflush(stderr__stderrp);
292 }
293 if (nrrdMaybeAlloc_va(nthresh[ni], nrrdTypeUChar, 3,
294 sx, sy, sz)) {
295 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: trouble allocating threshold %u",
296 me, (unsigned int)ni);
297 airMopError(mop); return 1;
298 }
299 thr = (unsigned char *)(nthresh[ni]->data);
300 for (I=0; I<sx*sy*sz; I++) {
301 val = nrrdFLookup[nblur[ni]->type](nblur[ni]->data, I);
302 val -= AIR_CAST(float, DWthr)((float)(DWthr));
303 thr[I] = (val >= 0 ? 1 : 0);
304 }
305 }
306 if (verb) {
307 fprintf(stderr__stderrp, "done\n");
308 }
309
310 airMopOkay(mop);
311 return 0;
312}
313
314/*
315** _tenEpiRegBB: find the biggest bright CC
316*/
317int
318_tenEpiRegBB(Nrrd *nval, Nrrd *nsize) {
319 unsigned char *val;
320 int *size, big;
321 unsigned int ci;
322
323 val = (unsigned char *)(nval->data);
324 size = (int *)(nsize->data);
325 big = 0;
326 for (ci=0; ci<nsize->axis[0].size; ci++) {
327 big = val[ci] ? AIR_MAX(big, size[ci])((big) > (size[ci]) ? (big) : (size[ci])) : big;
328 }
329 return big;
330}
331
332int
333_tenEpiRegCC(Nrrd **nthr, int ninLen, int conny, int verb) {
334 static const char me[]="_tenEpiRegCC";
335 Nrrd *nslc, *ncc, *nval, *nsize;
336 airArray *mop;
337 int ni, z, sz, big;
338
339 mop = airMopNew();
340 airMopAdd(mop, nslc=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
341 airMopAdd(mop, nval=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
342 airMopAdd(mop, ncc=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
343 airMopAdd(mop, nsize=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
344 sz = nthr[0]->axis[2].size;
345 if (verb) {
346 fprintf(stderr__stderrp, "%s:\n ", me); fflush(stderr__stderrp);
347 }
348 for (ni=0; ni<ninLen; ni++) {
349 if (verb) {
350 fprintf(stderr__stderrp, "%2d ", ni); fflush(stderr__stderrp);
351 }
352 /* for each volume, we find the biggest bright 3-D CC, and merge
353 down (to dark) all smaller bright pieces. Then, within each
354 slice, we do 2-D CCs, find the biggest bright CC (size == big),
355 and merge up (to bright) all small dark pieces, where
356 (currently) small is big/2 */
357 big = -1;
358 if (nrrdCCFind(ncc, &nval, nthr[ni], nrrdTypeDefault, conny)
359 || nrrdCCSize(nsize, ncc)
360 || !(big = _tenEpiRegBB(nval, nsize))
361 || nrrdCCMerge(ncc, ncc, nval, -1, big-1, 0, conny)
362 || nrrdCCRevalue(nthr[ni], ncc, nval)) {
363 if (big) {
364 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey,
365 "%s: trouble with 3-D processing nthr[%d]", me, ni);
366 return 1;
367 } else {
368 biffAddf(TENtenBiffKey, "%s: got size 0 for biggest bright CC of nthr[%d]",
369 me, ni);
370 return 1;
371 }
372 }
373 for (z=0; z<sz; z++) {
374 big = -1;
375 if ( nrrdSlice(nslc, nthr[ni], 2, z)
376 || nrrdCCFind(ncc, &nval, nslc, nrrdTypeDefault, conny)
377 || nrrdCCSize(nsize, ncc)
378 || !(big = _tenEpiRegBB(nval, nsize))
379 || nrrdCCMerge(ncc, ncc, nval, 1, big/2, 0, conny)
380 || nrrdCCRevalue(nslc, ncc, nval)
381 || nrrdSplice(nthr[ni], nthr[ni], nslc, 2, z) ) {
382 if (big) {
383 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: trouble processing slice %d of nthr[%d]",
384 me, z, ni);
385 return 1;
386 } else {
387 /* biggest bright CC on this slice had size 0
388 <==> there was no bright CC on this slice, move on */
389 }
390 }
391 }
392 }
393 if (verb) {
394 fprintf(stderr__stderrp, "done\n");
395 }
396
397 airMopOkay(mop);
398 return 0;
399}
400
401#define MEAN_X0 0
402#define MEAN_Y1 1
403#define M_022 2
404#define M_113 3
405#define M_204 4
406
407/*
408** _tenEpiRegMoments()
409**
410** the moments are stored in (of course) a nrrd, one scanline per slice,
411** with each scanline containing:
412**
413** 0 1 2 3 4
414** mean(x) mean(y) M_02 M_11 M_20
415*/
416int
417_tenEpiRegMoments(Nrrd **nmom, Nrrd **nthresh, unsigned int ninLen,
418 int verb) {
419 static const char me[]="_tenEpiRegMoments";
420 size_t sx, sy, sz, xi, yi, zi, ni;
421 double N, mx, my, cx, cy, x, y, M02, M11, M20, *mom;
422 float val;
423 unsigned char *thr;
424
425 sx = nthresh[0]->axis[0].size;
426 sy = nthresh[0]->axis[1].size;
427 sz = nthresh[0]->axis[2].size;
428 if (verb) {
429 fprintf(stderr__stderrp, "%s:\n ", me); fflush(stderr__stderrp);
430 }
431 for (ni=0; ni<ninLen; ni++) {
432 if (verb) {
433 fprintf(stderr__stderrp, "%2u ", (unsigned int)ni); fflush(stderr__stderrp);
434 }
435 if (nrrdMaybeAlloc_va(nmom[ni], nrrdTypeDouble, 2,
436 AIR_CAST(size_t, 5)((size_t)(5)), sz)) {
437 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey,
438 "%s: couldn't allocate nmom[%u]", me, (unsigned int)ni);
439 return 1;
440 }
441 nrrdAxisInfoSet_va(nmom[ni], nrrdAxisInfoLabel, "mx,my,h,s,t", "z");
442 thr = (unsigned char *)(nthresh[ni]->data);
443 mom = (double *)(nmom[ni]->data);
444 for (zi=0; zi<sz; zi++) {
445 /* ------ find mx, my */
446 N = 0;
447 mx = my = 0.0;
448 for (yi=0; yi<sy; yi++) {
449 for (xi=0; xi<sx; xi++) {
450 val = thr[xi + sx*yi];
451 N += val;
452 mx += xi*val;
453 my += yi*val;
454 }
455 }
456 if (N == sx*sy) {
457 biffAddf(TENtenBiffKey, "%s: saw only non-zero pixels in nthresh[%u]; "
458 "DWI hreshold too low?", me, (unsigned int)ni);
459 return 1;
460 }
461 if (N) {
462 /* there were non-zero pixels */
463 mx /= N;
464 my /= N;
465 cx = sx/2.0;
466 cy = sy/2.0;
467 /* ------ find M02, M11, M20 */
468 M02 = M11 = M20 = 0.0;
469 for (yi=0; yi<sy; yi++) {
470 for (xi=0; xi<sx; xi++) {
471 val = thr[xi + sx*yi];
472 x = xi - cx;
473 y = yi - cy;
474 M02 += y*y*val;
475 M11 += x*y*val;
476 M20 += x*x*val;
477 }
478 }
479 M02 /= N;
480 M11 /= N;
481 M20 /= N;
482 /* ------ set output */
483 mom[MEAN_X0] = mx;
484 mom[MEAN_Y1] = my;
485 mom[M_022] = M02;
486 mom[M_113] = M11;
487 mom[M_204] = M20;
488 } else {
489 /* there were no non-zero pixels */
490 mom[MEAN_X0] = 0;
491 mom[MEAN_Y1] = 0;
492 mom[M_022] = 0;
493 mom[M_113] = 0;
494 mom[M_204] = 0;
495 }
496 thr += sx*sy;
497 mom += 5;
498 }
499 }
500 if (verb) {
501 fprintf(stderr__stderrp, "done\n");
502 }
503
504 return 0;
505}
506
507/*
508** _tenEpiRegPairXforms
509**
510** uses moment information to compute all pair-wise transforms, which are
511** stored in the 3 x ninLen x ninLen x sizeZ output. If xfr = npxfr->data,
512** xfr[0 + 3*(zi + sz*(A + ninLen*B))] is shear,
513** xfr[1 + " ] is scale, and
514** xfr[2 + " ] is translate in the transform
515** that maps slice zi from volume A to volume B.
516*/
517int
518_tenEpiRegPairXforms(Nrrd *npxfr, Nrrd **nmom, int ninLen) {
519 static const char me[]="_tenEpiRegPairXforms";
520 double *xfr, *A, *B, hh, ss, tt;
521 int ai, bi, zi, sz;
522
523 sz = nmom[0]->axis[1].size;
524 if (nrrdMaybeAlloc_va(npxfr, nrrdTypeDouble, 4,
525 AIR_CAST(size_t, 5)((size_t)(5)),
526 AIR_CAST(size_t, sz)((size_t)(sz)),
527 AIR_CAST(size_t, ninLen)((size_t)(ninLen)),
528 AIR_CAST(size_t, ninLen)((size_t)(ninLen)))) {
529 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: couldn't allocate transform nrrd", me);
530 return 1;
531 }
532 nrrdAxisInfoSet_va(npxfr, nrrdAxisInfoLabel,
533 "mx,my,h,s,t", "zi", "orig", "target");
534 xfr = (double *)(npxfr->data);
535 for (bi=0; bi<ninLen; bi++) {
536 for (ai=0; ai<ninLen; ai++) {
537 for (zi=0; zi<sz; zi++) {
538 A = (double*)(nmom[ai]->data) + 5*zi;
539 B = (double*)(nmom[bi]->data) + 5*zi;
540 ss = sqrt((A[M_204]*B[M_022] - B[M_113]*B[M_113]) /
541 (A[M_204]*A[M_022] - A[M_113]*A[M_113]));
542 hh = (B[M_113] - ss*A[M_113])/A[M_204];
543 tt = B[MEAN_Y1] - A[MEAN_Y1];
544 ELL_5V_SET(xfr + 5*(zi + sz*(ai + ninLen*bi)),((xfr + 5*(zi + sz*(ai + ninLen*bi)))[0]=(A[0]), (xfr + 5*(zi
+ sz*(ai + ninLen*bi)))[1]=(A[1]), (xfr + 5*(zi + sz*(ai + ninLen
*bi)))[2]=(hh), (xfr + 5*(zi + sz*(ai + ninLen*bi)))[3]=(ss),
(xfr + 5*(zi + sz*(ai + ninLen*bi)))[4]=(tt))
545 A[MEAN_X], A[MEAN_Y], hh, ss, tt)((xfr + 5*(zi + sz*(ai + ninLen*bi)))[0]=(A[0]), (xfr + 5*(zi
+ sz*(ai + ninLen*bi)))[1]=(A[1]), (xfr + 5*(zi + sz*(ai + ninLen
*bi)))[2]=(hh), (xfr + 5*(zi + sz*(ai + ninLen*bi)))[3]=(ss),
(xfr + 5*(zi + sz*(ai + ninLen*bi)))[4]=(tt));
546 }
547 }
548 }
549 return 0;
550}
551
552#define SHEAR2 2
553#define SCALE3 3
554#define TRAN4 4
555
556int
557_tenEpiRegEstimHST(Nrrd *nhst, Nrrd *npxfr, int ninLen, Nrrd *ngrad) {
558 static const char me[]="_tenEpiRegEstimHST";
559 double *hst, *grad, *mat, *vec, *ans, *pxfr, *gA, *gB;
560 int z, sz, A, B, npairs, ri;
561 Nrrd **nmat, *nvec, **ninv, *nans;
562 airArray *mop;
563 int order;
564
565 order = 1;
566
567 sz = npxfr->axis[1].size;
568 npairs = ninLen*(ninLen-1);
569
570 mop = airMopNew();
571 nmat = (Nrrd**)calloc(sz, sizeof(Nrrd*));
572 ninv = (Nrrd**)calloc(sz, sizeof(Nrrd*));
573 airMopAdd(mop, nmat, airFree, airMopAlways);
574 airMopAdd(mop, ninv, airFree, airMopAlways);
575 for (z=0; z<sz; z++) {
576 airMopAdd(mop, nmat[z]=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
577 if (nrrdMaybeAlloc_va(nmat[z], nrrdTypeDouble, 2,
578 AIR_CAST(size_t, (1 == order ? 3 : 9))((size_t)((1 == order ? 3 : 9))),
579 AIR_CAST(size_t, npairs)((size_t)(npairs)))) {
580 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: couldn't allocate fitting matrices", me);
581 airMopError(mop); return 1;
582 }
583 airMopAdd(mop, ninv[z]=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
584 }
585 airMopAdd(mop, nvec=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
586 airMopAdd(mop, nans=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
587 if (nrrdMaybeAlloc_va(nhst, nrrdTypeDouble, 2,
588 AIR_CAST(size_t, (1 == order ? 9 : 27))((size_t)((1 == order ? 9 : 27))),
589 AIR_CAST(size_t, sz)((size_t)(sz)))
590 || nrrdMaybeAlloc_va(nvec, nrrdTypeDouble, 2,
591 AIR_CAST(size_t, 1)((size_t)(1)),
592 AIR_CAST(size_t, npairs)((size_t)(npairs)))) {
593 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: couldn't allocate HST nrrd", me);
594 airMopError(mop); return 1;
595 }
596 nrrdAxisInfoSet_va(nhst, nrrdAxisInfoLabel,
597 (1 == order
598 ? "Hx,Hy,Hz,Sx,Sy,Sz,Tx,Ty,Tz"
599 : "HST parms"), "z");
600
601 /* ------ per-slice: compute model fitting matrix and its pseudo-inverse */
602 grad = (double *)(ngrad->data);
603 for (z=0; z<sz; z++) {
604 hst = (double *)(nhst->data) + (1 == order ? 9 : 27)*z;
605 mat = (double *)(nmat[z]->data);
606 ri = 0;
607 for (A=0; A<ninLen; A++) {
608 for (B=0; B<ninLen; B++) {
609 if (A == B) continue;
610 pxfr = (double *)(npxfr->data) + 0 + 5*(z + sz*(A + ninLen*B));
611 gA = grad + 0 + 3*A;
612 gB = grad + 0 + 3*B;
613 if (1 == order) {
614 ELL_3V_SET(mat + 3*ri,((mat + 3*ri)[0] = (gB[0] - pxfr[3]*gA[0]), (mat + 3*ri)[1] =
(gB[1] - pxfr[3]*gA[1]), (mat + 3*ri)[2] = (gB[2] - pxfr[3]*
gA[2]))
615 gB[0] - pxfr[SCALE]*gA[0],((mat + 3*ri)[0] = (gB[0] - pxfr[3]*gA[0]), (mat + 3*ri)[1] =
(gB[1] - pxfr[3]*gA[1]), (mat + 3*ri)[2] = (gB[2] - pxfr[3]*
gA[2]))
616 gB[1] - pxfr[SCALE]*gA[1],((mat + 3*ri)[0] = (gB[0] - pxfr[3]*gA[0]), (mat + 3*ri)[1] =
(gB[1] - pxfr[3]*gA[1]), (mat + 3*ri)[2] = (gB[2] - pxfr[3]*
gA[2]))
617 gB[2] - pxfr[SCALE]*gA[2])((mat + 3*ri)[0] = (gB[0] - pxfr[3]*gA[0]), (mat + 3*ri)[1] =
(gB[1] - pxfr[3]*gA[1]), (mat + 3*ri)[2] = (gB[2] - pxfr[3]*
gA[2]));
618 } else {
619 ELL_9V_SET(mat + 9*ri,((mat + 9*ri)[0]=(gB[0] - pxfr[3]*gA[0]), (mat + 9*ri)[1]=(gB
[1] - pxfr[3]*gA[1]), (mat + 9*ri)[2]=(gB[2] - pxfr[3]*gA[2])
, (mat + 9*ri)[3]=(gB[0]*gB[0]*gB[0] - pxfr[3]*gA[0]*gA[0]*gA
[0]), (mat + 9*ri)[4]=(gB[1]*gB[1]*gB[1] - pxfr[3]*gA[1]*gA[1
]*gA[1]), (mat + 9*ri)[5]=(gB[2]*gB[2]*gB[2] - pxfr[3]*gA[2]*
gA[2]*gA[2]), (mat + 9*ri)[6]=(gB[0]*gB[1]*gB[2] - pxfr[3]*gA
[0]*gA[1]*gA[2]), (mat + 9*ri)[7]=(1), (mat + 9*ri)[8]=(1))
620 gB[0] - pxfr[SCALE]*gA[0],((mat + 9*ri)[0]=(gB[0] - pxfr[3]*gA[0]), (mat + 9*ri)[1]=(gB
[1] - pxfr[3]*gA[1]), (mat + 9*ri)[2]=(gB[2] - pxfr[3]*gA[2])
, (mat + 9*ri)[3]=(gB[0]*gB[0]*gB[0] - pxfr[3]*gA[0]*gA[0]*gA
[0]), (mat + 9*ri)[4]=(gB[1]*gB[1]*gB[1] - pxfr[3]*gA[1]*gA[1
]*gA[1]), (mat + 9*ri)[5]=(gB[2]*gB[2]*gB[2] - pxfr[3]*gA[2]*
gA[2]*gA[2]), (mat + 9*ri)[6]=(gB[0]*gB[1]*gB[2] - pxfr[3]*gA
[0]*gA[1]*gA[2]), (mat + 9*ri)[7]=(1), (mat + 9*ri)[8]=(1))
621 gB[1] - pxfr[SCALE]*gA[1],((mat + 9*ri)[0]=(gB[0] - pxfr[3]*gA[0]), (mat + 9*ri)[1]=(gB
[1] - pxfr[3]*gA[1]), (mat + 9*ri)[2]=(gB[2] - pxfr[3]*gA[2])
, (mat + 9*ri)[3]=(gB[0]*gB[0]*gB[0] - pxfr[3]*gA[0]*gA[0]*gA
[0]), (mat + 9*ri)[4]=(gB[1]*gB[1]*gB[1] - pxfr[3]*gA[1]*gA[1
]*gA[1]), (mat + 9*ri)[5]=(gB[2]*gB[2]*gB[2] - pxfr[3]*gA[2]*
gA[2]*gA[2]), (mat + 9*ri)[6]=(gB[0]*gB[1]*gB[2] - pxfr[3]*gA
[0]*gA[1]*gA[2]), (mat + 9*ri)[7]=(1), (mat + 9*ri)[8]=(1))
622 gB[2] - pxfr[SCALE]*gA[2],((mat + 9*ri)[0]=(gB[0] - pxfr[3]*gA[0]), (mat + 9*ri)[1]=(gB
[1] - pxfr[3]*gA[1]), (mat + 9*ri)[2]=(gB[2] - pxfr[3]*gA[2])
, (mat + 9*ri)[3]=(gB[0]*gB[0]*gB[0] - pxfr[3]*gA[0]*gA[0]*gA
[0]), (mat + 9*ri)[4]=(gB[1]*gB[1]*gB[1] - pxfr[3]*gA[1]*gA[1
]*gA[1]), (mat + 9*ri)[5]=(gB[2]*gB[2]*gB[2] - pxfr[3]*gA[2]*
gA[2]*gA[2]), (mat + 9*ri)[6]=(gB[0]*gB[1]*gB[2] - pxfr[3]*gA
[0]*gA[1]*gA[2]), (mat + 9*ri)[7]=(1), (mat + 9*ri)[8]=(1))
623 gB[0]*gB[0]*gB[0] - pxfr[SCALE]*gA[0]*gA[0]*gA[0],((mat + 9*ri)[0]=(gB[0] - pxfr[3]*gA[0]), (mat + 9*ri)[1]=(gB
[1] - pxfr[3]*gA[1]), (mat + 9*ri)[2]=(gB[2] - pxfr[3]*gA[2])
, (mat + 9*ri)[3]=(gB[0]*gB[0]*gB[0] - pxfr[3]*gA[0]*gA[0]*gA
[0]), (mat + 9*ri)[4]=(gB[1]*gB[1]*gB[1] - pxfr[3]*gA[1]*gA[1
]*gA[1]), (mat + 9*ri)[5]=(gB[2]*gB[2]*gB[2] - pxfr[3]*gA[2]*
gA[2]*gA[2]), (mat + 9*ri)[6]=(gB[0]*gB[1]*gB[2] - pxfr[3]*gA
[0]*gA[1]*gA[2]), (mat + 9*ri)[7]=(1), (mat + 9*ri)[8]=(1))
624 gB[1]*gB[1]*gB[1] - pxfr[SCALE]*gA[1]*gA[1]*gA[1],((mat + 9*ri)[0]=(gB[0] - pxfr[3]*gA[0]), (mat + 9*ri)[1]=(gB
[1] - pxfr[3]*gA[1]), (mat + 9*ri)[2]=(gB[2] - pxfr[3]*gA[2])
, (mat + 9*ri)[3]=(gB[0]*gB[0]*gB[0] - pxfr[3]*gA[0]*gA[0]*gA
[0]), (mat + 9*ri)[4]=(gB[1]*gB[1]*gB[1] - pxfr[3]*gA[1]*gA[1
]*gA[1]), (mat + 9*ri)[5]=(gB[2]*gB[2]*gB[2] - pxfr[3]*gA[2]*
gA[2]*gA[2]), (mat + 9*ri)[6]=(gB[0]*gB[1]*gB[2] - pxfr[3]*gA
[0]*gA[1]*gA[2]), (mat + 9*ri)[7]=(1), (mat + 9*ri)[8]=(1))
625 gB[2]*gB[2]*gB[2] - pxfr[SCALE]*gA[2]*gA[2]*gA[2],((mat + 9*ri)[0]=(gB[0] - pxfr[3]*gA[0]), (mat + 9*ri)[1]=(gB
[1] - pxfr[3]*gA[1]), (mat + 9*ri)[2]=(gB[2] - pxfr[3]*gA[2])
, (mat + 9*ri)[3]=(gB[0]*gB[0]*gB[0] - pxfr[3]*gA[0]*gA[0]*gA
[0]), (mat + 9*ri)[4]=(gB[1]*gB[1]*gB[1] - pxfr[3]*gA[1]*gA[1
]*gA[1]), (mat + 9*ri)[5]=(gB[2]*gB[2]*gB[2] - pxfr[3]*gA[2]*
gA[2]*gA[2]), (mat + 9*ri)[6]=(gB[0]*gB[1]*gB[2] - pxfr[3]*gA
[0]*gA[1]*gA[2]), (mat + 9*ri)[7]=(1), (mat + 9*ri)[8]=(1))
626 gB[0]*gB[1]*gB[2] - pxfr[SCALE]*gA[0]*gA[1]*gA[2],((mat + 9*ri)[0]=(gB[0] - pxfr[3]*gA[0]), (mat + 9*ri)[1]=(gB
[1] - pxfr[3]*gA[1]), (mat + 9*ri)[2]=(gB[2] - pxfr[3]*gA[2])
, (mat + 9*ri)[3]=(gB[0]*gB[0]*gB[0] - pxfr[3]*gA[0]*gA[0]*gA
[0]), (mat + 9*ri)[4]=(gB[1]*gB[1]*gB[1] - pxfr[3]*gA[1]*gA[1
]*gA[1]), (mat + 9*ri)[5]=(gB[2]*gB[2]*gB[2] - pxfr[3]*gA[2]*
gA[2]*gA[2]), (mat + 9*ri)[6]=(gB[0]*gB[1]*gB[2] - pxfr[3]*gA
[0]*gA[1]*gA[2]), (mat + 9*ri)[7]=(1), (mat + 9*ri)[8]=(1))
627 1, 1)((mat + 9*ri)[0]=(gB[0] - pxfr[3]*gA[0]), (mat + 9*ri)[1]=(gB
[1] - pxfr[3]*gA[1]), (mat + 9*ri)[2]=(gB[2] - pxfr[3]*gA[2])
, (mat + 9*ri)[3]=(gB[0]*gB[0]*gB[0] - pxfr[3]*gA[0]*gA[0]*gA
[0]), (mat + 9*ri)[4]=(gB[1]*gB[1]*gB[1] - pxfr[3]*gA[1]*gA[1
]*gA[1]), (mat + 9*ri)[5]=(gB[2]*gB[2]*gB[2] - pxfr[3]*gA[2]*
gA[2]*gA[2]), (mat + 9*ri)[6]=(gB[0]*gB[1]*gB[2] - pxfr[3]*gA
[0]*gA[1]*gA[2]), (mat + 9*ri)[7]=(1), (mat + 9*ri)[8]=(1));
628 /*
629 gB[0]*gB[1] - pxfr[SCALE]*gA[0]*gA[1],
630 gB[0]*gB[2] - pxfr[SCALE]*gA[0]*gA[2],
631 gB[1]*gB[2] - pxfr[SCALE]*gA[1]*gA[2]);
632 */
633 }
634 ri += 1;
635 }
636 }
637 if (nrrdHasNonExist(nmat[z])) {
638 /* as happens if there were zero slices in the segmentation output */
639 if (1 == order) {
640 ELL_3V_SET(hst + 0*3, 0, 0, 0)((hst + 0*3)[0] = (0), (hst + 0*3)[1] = (0), (hst + 0*3)[2] =
(0));
641 ELL_3V_SET(hst + 1*3, 0, 0, 0)((hst + 1*3)[0] = (0), (hst + 1*3)[1] = (0), (hst + 1*3)[2] =
(0));
642 ELL_3V_SET(hst + 2*3, 0, 0, 0)((hst + 2*3)[0] = (0), (hst + 2*3)[1] = (0), (hst + 2*3)[2] =
(0));
643 } else {
644 ELL_9V_SET(hst + 0*9, 0, 0, 0, 0, 0, 0, 0, 0, 0)((hst + 0*9)[0]=(0), (hst + 0*9)[1]=(0), (hst + 0*9)[2]=(0), (
hst + 0*9)[3]=(0), (hst + 0*9)[4]=(0), (hst + 0*9)[5]=(0), (hst
+ 0*9)[6]=(0), (hst + 0*9)[7]=(0), (hst + 0*9)[8]=(0));
645 ELL_9V_SET(hst + 1*9, 0, 0, 0, 0, 0, 0, 0, 0, 0)((hst + 1*9)[0]=(0), (hst + 1*9)[1]=(0), (hst + 1*9)[2]=(0), (
hst + 1*9)[3]=(0), (hst + 1*9)[4]=(0), (hst + 1*9)[5]=(0), (hst
+ 1*9)[6]=(0), (hst + 1*9)[7]=(0), (hst + 1*9)[8]=(0));
646 ELL_9V_SET(hst + 2*9, 0, 0, 0, 0, 0, 0, 0, 0, 0)((hst + 2*9)[0]=(0), (hst + 2*9)[1]=(0), (hst + 2*9)[2]=(0), (
hst + 2*9)[3]=(0), (hst + 2*9)[4]=(0), (hst + 2*9)[5]=(0), (hst
+ 2*9)[6]=(0), (hst + 2*9)[7]=(0), (hst + 2*9)[8]=(0));
647 }
648 } else {
649 if (ell_Nm_pseudo_inv(ninv[z], nmat[z])) {
650 biffMovef(TENtenBiffKey, ELLell_biff_key, "%s: trouble estimating model (slice %d)", me, z);
651 airMopError(mop); return 1;
652 }
653 }
654 }
655 vec = (double *)(nvec->data);
656
657 /* ------ find Sx, Sy, Sz per slice */
658 for (z=0; z<sz; z++) {
659 if (nrrdHasNonExist(nmat[z])) {
660 /* we've already zero-ed out this row in the HST nrrd */
661 continue;
662 }
663 hst = (double *)(nhst->data) + (1 == order ? 9 : 27)*z;
664 ri = 0;
665 for (A=0; A<ninLen; A++) {
666 for (B=0; B<ninLen; B++) {
667 if (A == B) continue;
668 pxfr = (double *)(npxfr->data) + 0 + 5*(z + sz*(A + ninLen*B));
669 vec[ri] = pxfr[SCALE3] - 1;
670 ri += 1;
671 }
672 }
673 if (ell_Nm_mul(nans, ninv[z], nvec)) {
674 biffMovef(TENtenBiffKey, ELLell_biff_key,
675 "%s: trouble estimating model (slice %d): Sx, Sy, Sz",
676 me, z);
677 airMopError(mop); return 1;
678 }
679 ans = (double *)(nans->data);
680 if (1 == order) {
681 ELL_3V_COPY(hst + 1*3, ans)((hst + 1*3)[0] = (ans)[0], (hst + 1*3)[1] = (ans)[1], (hst +
1*3)[2] = (ans)[2]);
682 } else {
683 ELL_9V_COPY(hst + 1*9, ans)((hst + 1*9)[0] = (ans)[0], (hst + 1*9)[1] = (ans)[1], (hst +
1*9)[2] = (ans)[2], (hst + 1*9)[3] = (ans)[3], (hst + 1*9)[4
] = (ans)[4], (hst + 1*9)[5] = (ans)[5], (hst + 1*9)[6] = (ans
)[6], (hst + 1*9)[7] = (ans)[7], (hst + 1*9)[8] = (ans)[8]);
684 }
685 }
686
687 /* ------ find Hx, Hy, Hz per slice */
688 for (z=0; z<sz; z++) {
689 if (nrrdHasNonExist(nmat[z])) {
690 /* we've already zero-ed out this row in the HST nrrd */
691 continue;
692 }
693 hst = (double *)(nhst->data) + (1 == order ? 9 : 27)*z;
694 ri = 0;
695 for (A=0; A<ninLen; A++) {
696 for (B=0; B<ninLen; B++) {
697 if (A == B) continue;
698 pxfr = (double *)(npxfr->data) + 0 + 5*(z + sz*(A + ninLen*B));
699 vec[ri] = pxfr[SHEAR2];
700 ri += 1;
701 }
702 }
703 if (ell_Nm_mul(nans, ninv[z], nvec)) {
704 biffAddf(TENtenBiffKey, ELLell_biff_key, "%s: trouble estimating model (slice %d): Hx, Hy, Hz",
705 me, z);
706 airMopError(mop); return 1;
707 }
708 ans = (double *)(nans->data);
709 if (1 == order) {
710 ELL_3V_COPY(hst + 0*3, ans)((hst + 0*3)[0] = (ans)[0], (hst + 0*3)[1] = (ans)[1], (hst +
0*3)[2] = (ans)[2]);
711 } else {
712 ELL_9V_COPY(hst + 0*9, ans)((hst + 0*9)[0] = (ans)[0], (hst + 0*9)[1] = (ans)[1], (hst +
0*9)[2] = (ans)[2], (hst + 0*9)[3] = (ans)[3], (hst + 0*9)[4
] = (ans)[4], (hst + 0*9)[5] = (ans)[5], (hst + 0*9)[6] = (ans
)[6], (hst + 0*9)[7] = (ans)[7], (hst + 0*9)[8] = (ans)[8]);
713 }
714 }
715
716 /* ------ find Tx, Ty, Tz per slice */
717 for (z=0; z<sz; z++) {
718 if (nrrdHasNonExist(nmat[z])) {
719 /* we've already zero-ed out this row in the HST nrrd */
720 continue;
721 }
722 hst = (double *)(nhst->data) + (1 == order ? 9 : 27)*z;
723 ri = 0;
724 for (A=0; A<ninLen; A++) {
725 for (B=0; B<ninLen; B++) {
726 if (A == B) continue;
727 pxfr = (double *)(npxfr->data) + 0 + 5*(z + sz*(A + ninLen*B));
728 vec[ri] = pxfr[TRAN4];
729 ri += 1;
730 }
731 }
732 if (ell_Nm_mul(nans, ninv[z], nvec)) {
733 biffMovef(TENtenBiffKey, ELLell_biff_key,
734 "%s: trouble estimating model (slice %d): Tx, Ty, Tz",
735 me, z);
736 airMopError(mop); return 1;
737 }
738 ans = (double *)(nans->data);
739 if (1 == order) {
740 ELL_3V_COPY(hst + 2*3, ans)((hst + 2*3)[0] = (ans)[0], (hst + 2*3)[1] = (ans)[1], (hst +
2*3)[2] = (ans)[2]);
741 } else {
742 ELL_9V_COPY(hst + 2*9, ans)((hst + 2*9)[0] = (ans)[0], (hst + 2*9)[1] = (ans)[1], (hst +
2*9)[2] = (ans)[2], (hst + 2*9)[3] = (ans)[3], (hst + 2*9)[4
] = (ans)[4], (hst + 2*9)[5] = (ans)[5], (hst + 2*9)[6] = (ans
)[6], (hst + 2*9)[7] = (ans)[7], (hst + 2*9)[8] = (ans)[8]);
743 }
744 }
745
746 airMopOkay(mop);
747 return 0;
748}
749
750int
751_tenEpiRegFitHST(Nrrd *nhst, Nrrd **_ncc, int ninLen,
752 double goodFrac, int prog, int verb) {
753 static const char me[]="_tenEpiRegFitHST";
754 airArray *mop;
755 Nrrd *ncc, *ntA, *ntB, *nsd, *nl2;
756 unsigned int cc, sz, zi, sh, hi;
757 float *mess, *two, tmp;
758 double *hst, x, y, xx, xy, mm, bb;
759
760 mop = airMopNew();
761 airMopAdd(mop, ncc=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
762 airMopAdd(mop, ntA=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
763 airMopAdd(mop, ntB=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
764 airMopAdd(mop, nsd=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
765 airMopAdd(mop, nl2=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
766 /* do SD and L2 projections of the CCs along the DWI axis,
767 integrate these over the X and Y axes of the slices,
768 and define per-slice "messiness" as the quotient of the
769 SD integral with the L2 integral */
770 if (verb) {
771 fprintf(stderr__stderrp, "%s: measuring segmentation uncertainty ... ", me);
772 fflush(stderr__stderrp);
773 }
774 if (nrrdJoin(ncc, (const Nrrd*const*)_ncc, ninLen, 0, AIR_TRUE1)
775 || nrrdProject(ntA, ncc, 0, nrrdMeasureSD, nrrdTypeFloat)
776 || nrrdProject(ntB, ntA, 0, nrrdMeasureSum, nrrdTypeFloat)
777 || nrrdProject(nsd, ntB, 0, nrrdMeasureSum, nrrdTypeFloat)
778 || nrrdProject(ntA, ncc, 0, nrrdMeasureL2, nrrdTypeFloat)
779 || nrrdProject(ntB, ntA, 0, nrrdMeasureSum, nrrdTypeFloat)
780 || nrrdProject(nl2, ntB, 0, nrrdMeasureSum, nrrdTypeFloat)
781 || nrrdArithBinaryOp(ntA, nrrdBinaryOpDivide, nsd, nl2)) {
782 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: trouble doing CC projections", me);
783 airMopError(mop); return 1;
784 }
785 if (verb) {
786 fprintf(stderr__stderrp, "done\n");
787 }
788 if (prog && _tenEpiRegSave("regtmp-messy.txt", ntA,
789 NULL((void*)0), 0, "segmentation uncertainty")) {
790 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: EpiRegSave failed", me);
791 airMopError(mop); return 1;
792 }
793 /* now ntA stores the per-slice messiness */
794 mess = AIR_CAST(float*, ntA->data)((float*)(ntA->data));
795
796 /* allocate an array of 2 floats per slice */
797 sz = ntA->axis[0].size;
798 two = AIR_CAST(float*, calloc(2*sz, sizeof(float)))((float*)(calloc(2*sz, sizeof(float))));
799 if (!two) {
800 biffAddf(TENtenBiffKey, "%s: couldn't allocate tmp buffer", me);
801 airMopError(mop); return 1;
802 }
803 airMopAdd(mop, two, airFree, airMopAlways);
804 /* initial ordering: messiness, slice index */
805 for (zi=0; zi<sz; zi++) {
806 two[0 + 2*zi] = (AIR_EXISTS(mess[zi])(((int)(!((mess[zi]) - (mess[zi])))))
807 ? mess[zi]
808 : 666); /* don't use empty slices */
809 two[1 + 2*zi] = AIR_CAST(float, zi)((float)(zi));
810 }
811 /* sort into ascending messiness */
812 qsort(two, zi, 2*sizeof(float), nrrdValCompare[nrrdTypeFloat]);
813 /* flip ordering while thresholding messiness into usability */
814 for (zi=0; zi<sz; zi++) {
815 tmp = two[1 + 2*zi];
816 two[1 + 2*zi] = AIR_AFFINE(0, zi, sz-1, 0, 1)( ((double)(1)-(0))*((double)(zi)-(0)) / ((double)(sz-1)-(0))
+ (0)) <= goodFrac ? 1.0f : 0.0f;
817 two[0 + 2*zi] = tmp;
818 }
819 /* sort again, now into ascending slice order */
820 qsort(two, zi, 2*sizeof(float), nrrdValCompare[nrrdTypeFloat]);
821 if (verb) {
822 fprintf(stderr__stderrp, "%s: using slices", me);
823 for (zi=0; zi<sz; zi++) {
824 if (two[1 + 2*zi]) {
825 fprintf(stderr__stderrp, " %u", zi);
826 }
827 }
828 fprintf(stderr__stderrp, " for fitting\n");
829 }
830
831 /* perform fitting for each column in hst (regardless of
832 whether we're using a 1st or 2nd order model */
833 hst = (double*)(nhst->data);
834 sh = nhst->axis[0].size;
835 for (hi=0; hi<sh; hi++) {
836 x = y = xy = xx = 0;
837 cc = 0;
838 for (zi=0; zi<sz; zi++) {
839 if (!two[1 + 2*zi])
840 continue;
841 cc += 1;
842 x += zi;
843 xx += zi*zi;
844 y += hst[hi + sh*zi];
845 xy += zi*hst[hi + sh*zi];
846 }
847 x /= cc; xx /= cc; y /= cc; xy /= cc;
848 mm = (xy - x*y)/(xx - x*x);
849 bb = y - mm*x;
850 for (zi=0; zi<sz; zi++) {
851 hst[hi + sh*zi] = mm*zi + bb;
852 }
853 }
854
855 airMopOkay(mop);
856 return 0;
857}
858
859int
860_tenEpiRegGetHST(double *hhP, double *ssP, double *ttP,
861 int reference, int ni, int zi,
862 Nrrd *npxfr, Nrrd *nhst, Nrrd *ngrad) {
863 double *xfr, *hst, *_g, grad[9]; /* big enough for 2nd order */
864 int sz, ninLen;
865
866 int order;
867
868 order = 1;
869
870 /* these could also have been passed to us, but we can also discover them */
871 sz = npxfr->axis[1].size;
872 ninLen = npxfr->axis[2].size;
873
874 if (-1 == reference) {
875 /* we use the estimated H,S,T vectors to determine distortion
876 as a function of gradient direction, and then invert this */
877 _g = (double*)(ngrad->data) + 0 + 3*ni;
878 if (1 == order) {
879 hst = (double*)(nhst->data) + 0 + 9*zi;
880 *hhP = ELL_3V_DOT(_g, hst + 0*3)((_g)[0]*(hst + 0*3)[0] + (_g)[1]*(hst + 0*3)[1] + (_g)[2]*(hst
+ 0*3)[2]);
881 *ssP = 1 + ELL_3V_DOT(_g, hst + 1*3)((_g)[0]*(hst + 1*3)[0] + (_g)[1]*(hst + 1*3)[1] + (_g)[2]*(hst
+ 1*3)[2]);
882 *ttP = ELL_3V_DOT(_g, hst + 2*3)((_g)[0]*(hst + 2*3)[0] + (_g)[1]*(hst + 2*3)[1] + (_g)[2]*(hst
+ 2*3)[2]);
883 } else {
884 hst = (double*)(nhst->data) + 0 + 27*zi;
885 ELL_9V_SET(grad, _g[0], _g[1], _g[2],((grad)[0]=(_g[0]), (grad)[1]=(_g[1]), (grad)[2]=(_g[2]), (grad
)[3]=(_g[0]*_g[0]*_g[0]), (grad)[4]=(_g[1]*_g[1]*_g[1]), (grad
)[5]=(_g[2]*_g[2]*_g[2]), (grad)[6]=(_g[0]*_g[1]*_g[2]), (grad
)[7]=(0), (grad)[8]=(0))
886 _g[0]*_g[0]*_g[0], _g[1]*_g[1]*_g[1], _g[2]*_g[2]*_g[2],((grad)[0]=(_g[0]), (grad)[1]=(_g[1]), (grad)[2]=(_g[2]), (grad
)[3]=(_g[0]*_g[0]*_g[0]), (grad)[4]=(_g[1]*_g[1]*_g[1]), (grad
)[5]=(_g[2]*_g[2]*_g[2]), (grad)[6]=(_g[0]*_g[1]*_g[2]), (grad
)[7]=(0), (grad)[8]=(0))
887 _g[0]*_g[1]*_g[2],((grad)[0]=(_g[0]), (grad)[1]=(_g[1]), (grad)[2]=(_g[2]), (grad
)[3]=(_g[0]*_g[0]*_g[0]), (grad)[4]=(_g[1]*_g[1]*_g[1]), (grad
)[5]=(_g[2]*_g[2]*_g[2]), (grad)[6]=(_g[0]*_g[1]*_g[2]), (grad
)[7]=(0), (grad)[8]=(0))
888 0, 0)((grad)[0]=(_g[0]), (grad)[1]=(_g[1]), (grad)[2]=(_g[2]), (grad
)[3]=(_g[0]*_g[0]*_g[0]), (grad)[4]=(_g[1]*_g[1]*_g[1]), (grad
)[5]=(_g[2]*_g[2]*_g[2]), (grad)[6]=(_g[0]*_g[1]*_g[2]), (grad
)[7]=(0), (grad)[8]=(0));
889 /*
890 _g[0]*_g[1], _g[0]*_g[2], _g[1]*_g[2]);
891 */
892 *hhP = ELL_9V_DOT(grad, hst + 0*9)((grad)[0]*(hst + 0*9)[0] + (grad)[1]*(hst + 0*9)[1] + (grad)
[2]*(hst + 0*9)[2] + (grad)[3]*(hst + 0*9)[3] + (grad)[4]*(hst
+ 0*9)[4] + (grad)[5]*(hst + 0*9)[5] + (grad)[6]*(hst + 0*9)
[6] + (grad)[7]*(hst + 0*9)[7] + (grad)[8]*(hst + 0*9)[8]);
893 *ssP = 1 + ELL_9V_DOT(grad, hst + 1*9)((grad)[0]*(hst + 1*9)[0] + (grad)[1]*(hst + 1*9)[1] + (grad)
[2]*(hst + 1*9)[2] + (grad)[3]*(hst + 1*9)[3] + (grad)[4]*(hst
+ 1*9)[4] + (grad)[5]*(hst + 1*9)[5] + (grad)[6]*(hst + 1*9)
[6] + (grad)[7]*(hst + 1*9)[7] + (grad)[8]*(hst + 1*9)[8]);
894 *ttP = ELL_9V_DOT(grad, hst + 2*9)((grad)[0]*(hst + 2*9)[0] + (grad)[1]*(hst + 2*9)[1] + (grad)
[2]*(hst + 2*9)[2] + (grad)[3]*(hst + 2*9)[3] + (grad)[4]*(hst
+ 2*9)[4] + (grad)[5]*(hst + 2*9)[5] + (grad)[6]*(hst + 2*9)
[6] + (grad)[7]*(hst + 2*9)[7] + (grad)[8]*(hst + 2*9)[8]);
895 }
896 } else {
897 /* we register against a specific DWI */
898 xfr = (double*)(npxfr->data) + 0 + 5*(zi + sz*(reference + ninLen*ni));
899 *hhP = xfr[2];
900 *ssP = xfr[3];
901 *ttP = xfr[4];
902 }
903 return 0;
904}
905
906/*
907** _tenEpiRegSliceWarp
908**
909** Apply [hh,ss,tt] transform to nin, putting results in nout, but with
910** some trickiness:
911** - nwght and nidx are already allocated to the the weights (type float)
912** and indices for resampling nin with "kern" and "kparm"
913** - nout is already allocated to the correct size and type
914** - nin is type float, but output must be type nout->type
915** - nin is been transposed to have the resampled axis fastest in memory,
916** but nout output will not be transposed
917*/
918int
919_tenEpiRegSliceWarp(Nrrd *nout, Nrrd *nin, Nrrd *nwght, Nrrd *nidx,
920 const NrrdKernel *kern, double *kparm,
921 double hh, double ss, double tt, double cx, double cy) {
922 float *wght, *in, pp, pf, tmp;
923 int *idx;
924 unsigned int supp;
925 size_t sx, sy, xi, yi, pb, pi;
926 double (*ins)(void *, size_t, double), (*clamp)(double);
927
928 sy = nin->axis[0].size;
929 sx = nin->axis[1].size;
930 supp = AIR_CAST(unsigned int, kern->support(kparm))((unsigned int)(kern->support(kparm)));
931 ins = nrrdDInsert[nout->type];
932 clamp = nrrdDClamp[nout->type];
933
934 in = AIR_CAST(float*, nin->data)((float*)(nin->data));
935 for (xi=0; xi<sx; xi++) {
936 idx = AIR_CAST(int*, nidx->data)((int*)(nidx->data));
937 wght = AIR_CAST(float*, nwght->data)((float*)(nwght->data));
938 for (yi=0; yi<sy; yi++) {
939 pp = AIR_CAST(float, hh*(xi - cx) + ss*(yi - cy) + tt + cy)((float)(hh*(xi - cx) + ss*(yi - cy) + tt + cy));
940 pb = AIR_CAST(size_t, floor(pp))((size_t)(floor(pp)));
941 pf = pp - pb;
942 for (pi=0; pi<2*supp; pi++) {
943 idx[pi] = AIR_MIN(pb + pi - (supp-1), sy-1)((pb + pi - (supp-1)) < (sy-1) ? (pb + pi - (supp-1)) : (sy
-1));
944 wght[pi] = pi - (supp-1) - pf;
945 }
946 idx += 2*supp;
947 wght += 2*supp;
948 }
949 idx = AIR_CAST(int*, nidx->data)((int*)(nidx->data));
950 wght = AIR_CAST(float*, nwght->data)((float*)(nwght->data));
951 kern->evalN_f(wght, wght, 2*supp*sy, kparm);
952 for (yi=0; yi<sy; yi++) {
953 tmp = 0;
954 for (pi=0; pi<2*supp; pi++) {
955 tmp += in[idx[pi]]*wght[pi];
956 }
957 ins(nout->data, xi + sx*yi, clamp(ss*tmp));
958 idx += 2*supp;
959 wght += 2*supp;
960 }
961 in += sy;
962 }
963
964 return 0;
965}
966
967/*
968** _tenEpiRegWarp()
969**
970*/
971int
972_tenEpiRegWarp(Nrrd **ndone, Nrrd *npxfr, Nrrd *nhst, Nrrd *ngrad,
973 Nrrd **nin, int ninLen,
974 int reference, const NrrdKernel *kern, double *kparm,
975 int verb) {
976 static const char me[]="_tenEpiRegWarp";
977 Nrrd *ntmp, *nfin, *nslcA, *nslcB, *nwght, *nidx;
978 airArray *mop;
979 int sx, sy, sz, ni, zi, supp;
980 double hh, ss, tt, cx, cy;
981
982 mop = airMopNew();
983 airMopAdd(mop, ntmp=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
984 airMopAdd(mop, nfin=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
985 airMopAdd(mop, nslcA=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
986 airMopAdd(mop, nslcB=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
987 airMopAdd(mop, nwght=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
988 airMopAdd(mop, nidx=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
989
990 if (verb) {
991 fprintf(stderr__stderrp, "%s:\n ", me); fflush(stderr__stderrp);
992 }
993 sx = nin[0]->axis[0].size;
994 sy = nin[0]->axis[1].size;
995 sz = nin[0]->axis[2].size;
996 cx = sx/2.0;
997 cy = sy/2.0;
998 supp = (int)kern->support(kparm);
999 if (nrrdMaybeAlloc_va(nwght, nrrdTypeFloat, 2,
1000 AIR_CAST(size_t, 2*supp)((size_t)(2*supp)),
1001 AIR_CAST(size_t, sy)((size_t)(sy)))
1002 || nrrdMaybeAlloc_va(nidx, nrrdTypeInt, 2,
1003 AIR_CAST(size_t, 2*supp)((size_t)(2*supp)),
1004 AIR_CAST(size_t, sy)((size_t)(sy)))) {
1005 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: trouble allocating buffers", me);
1006 airMopError(mop); return 1;
1007 }
1008 for (ni=0; ni<ninLen; ni++) {
1009 if (verb) {
1010 fprintf(stderr__stderrp, "%2d ", ni); fflush(stderr__stderrp);
1011 }
1012 if (nrrdCopy(ndone[ni], nin[ni])
1013 || ((!ni) && nrrdSlice(nslcB, ndone[ni], 2, 0)) /* only when 0==ni */
1014 || nrrdAxesSwap(ntmp, nin[ni], 0, 1)
1015 || nrrdConvert(nfin, ntmp, nrrdTypeFloat)) {
1016 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: trouble prepping at ni=%d", me, ni);
1017 airMopError(mop); return 1;
1018 }
1019 for (zi=0; zi<sz; zi++) {
1020 if (_tenEpiRegGetHST(&hh, &ss, &tt, reference,
1021 ni, zi, npxfr, nhst, ngrad)
1022 || nrrdSlice(nslcA, nfin, 2, zi)
1023 || _tenEpiRegSliceWarp(nslcB, nslcA, nwght, nidx, kern, kparm,
1024 hh, ss, tt, cx, cy)
1025 || nrrdSplice(ndone[ni], ndone[ni], nslcB, 2, zi)) {
1026 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: trouble on slice %d if ni=%d", me, zi, ni);
1027 /* because the _tenEpiReg calls above don't use biff */
1028 airMopError(mop); return 1;
1029 }
1030 }
1031 }
1032 if (verb) {
1033 fprintf(stderr__stderrp, "done\n");
1034 }
1035
1036 airMopOkay(mop);
1037 return 0;
1038}
1039
1040int
1041tenEpiRegister3D(Nrrd **nout, Nrrd **nin, unsigned int ninLen, Nrrd *_ngrad,
1042 int reference,
1043 double bwX, double bwY, double fitFrac,
1044 double DWthr, int doCC,
1045 const NrrdKernel *kern, double *kparm,
1046 int progress, int verbose) {
1047 static const char me[]="tenEpiRegister3D";
1048 airArray *mop;
1049 Nrrd **nbuffA, **nbuffB, *npxfr, *nhst, *ngrad;
1050 int hack1, hack2;
1051 unsigned int i;
1052
1053 hack1 = nrrdStateAlwaysSetContent;
1054 hack2 = nrrdStateDisableContent;
1055 nrrdStateAlwaysSetContent = AIR_FALSE0;
1056 nrrdStateDisableContent = AIR_TRUE1;
1057
1058 mop = airMopNew();
1059 if (_tenEpiRegCheck(nout, nin, ninLen, _ngrad, reference,
1060 bwX, bwY, DWthr,
1061 kern, kparm)) {
1062 biffAddf(TENtenBiffKey, "%s: trouble with input", me);
1063 airMopError(mop); return 1;
1064 }
1065
1066 nbuffA = (Nrrd **)calloc(ninLen, sizeof(Nrrd*));
1067 nbuffB = (Nrrd **)calloc(ninLen, sizeof(Nrrd*));
1068 if (!( nbuffA && nbuffB )) {
1069 biffAddf(TENtenBiffKey, "%s: couldn't allocate tmp nrrd pointer arrays", me);
1070 airMopError(mop); return 1;
1071 }
1072 airMopAdd(mop, nbuffA, airFree, airMopAlways);
1073 airMopAdd(mop, nbuffB, airFree, airMopAlways);
1074 for (i=0; i<ninLen; i++) {
1075 airMopAdd(mop, nbuffA[i] = nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
1076 airMopAdd(mop, nbuffB[i] = nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
1077 nrrdAxisInfoCopy(nout[i], nin[i], NULL((void*)0), NRRD_AXIS_INFO_NONE0);
1078 }
1079 airMopAdd(mop, npxfr = nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
1080 airMopAdd(mop, nhst = nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
1081 airMopAdd(mop, ngrad = nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
1082 if (nrrdConvert(ngrad, _ngrad, nrrdTypeDouble)) {
1083 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: trouble converting gradients to doubles", me);
1084 airMopError(mop); return 1;
1085 }
1086
1087 /* ------ blur */
1088 if (_tenEpiRegBlur(nbuffA, nin, ninLen, bwX, bwY, verbose)) {
1089 biffAddf(TENtenBiffKey, "%s: trouble %s", me, (bwX || bwY) ? "blurring" : "copying");
1090 airMopError(mop); return 1;
1091 }
1092 if (progress && _tenEpiRegSave("regtmp-blur.nrrd", NULL((void*)0),
1093 nbuffA, ninLen, "blurred DWIs")) {
1094 airMopError(mop); return 1;
1095 }
1096
1097 /* ------ threshold */
1098 if (_tenEpiRegThreshold(nbuffB, nbuffA, ninLen,
1099 DWthr, verbose, progress, 1.5)) {
1100 biffAddf(TENtenBiffKey, "%s: trouble thresholding", me);
1101 airMopError(mop); return 1;
1102 }
1103 if (progress && _tenEpiRegSave("regtmp-thresh.nrrd", NULL((void*)0),
1104 nbuffB, ninLen, "thresholded DWIs")) {
1105 airMopError(mop); return 1;
1106 }
1107
1108 /* ------ connected components */
1109 if (doCC) {
1110 if (_tenEpiRegCC(nbuffB, ninLen, 1, verbose)) {
1111 biffAddf(TENtenBiffKey, "%s: trouble doing connected components", me);
1112 airMopError(mop); return 1;
1113 }
1114 if (progress && _tenEpiRegSave("regtmp-ccs.nrrd", NULL((void*)0),
1115 nbuffB, ninLen, "connected components")) {
1116 airMopError(mop); return 1;
1117 }
1118 }
1119
1120 /* ------ moments */
1121 if (_tenEpiRegMoments(nbuffA, nbuffB, ninLen, verbose)) {
1122 biffAddf(TENtenBiffKey, "%s: trouble finding moments", me);
1123 airMopError(mop); return 1;
1124 }
1125 if (progress && _tenEpiRegSave("regtmp-mom.nrrd", NULL((void*)0),
1126 nbuffA, ninLen, "moments")) {
1127 airMopError(mop); return 1;
1128 }
1129
1130 /* ------ transforms */
1131 if (_tenEpiRegPairXforms(npxfr, nbuffA, ninLen)) {
1132 biffAddf(TENtenBiffKey, "%s: trouble calculating transforms", me);
1133 airMopError(mop); return 1;
1134 }
1135 if (progress && _tenEpiRegSave("regtmp-pxfr.nrrd", npxfr,
1136 NULL((void*)0), 0, "pair-wise xforms")) {
1137 airMopError(mop); return 1;
1138 }
1139
1140 if (-1 == reference) {
1141 /* ------ HST estimation */
1142 if (_tenEpiRegEstimHST(nhst, npxfr, ninLen, ngrad)) {
1143 biffAddf(TENtenBiffKey, "%s: trouble estimating HST", me);
1144 airMopError(mop); return 1;
1145 }
1146 if (progress && _tenEpiRegSave("regtmp-hst.txt", nhst,
1147 NULL((void*)0), 0, "HST estimates")) {
1148 airMopError(mop); return 1;
1149 }
1150
1151 if (fitFrac) {
1152 /* ------ HST parameter fitting */
1153 if (_tenEpiRegFitHST(nhst, nbuffB, ninLen, fitFrac, progress, verbose)) {
1154 biffAddf(TENtenBiffKey, "%s: trouble fitting HST", me);
1155 airMopError(mop); return 1;
1156 }
1157 if (progress && _tenEpiRegSave("regtmp-fit-hst.txt", nhst,
1158 NULL((void*)0), 0, "fitted HST")) {
1159 airMopError(mop); return 1;
1160 }
1161 }
1162 }
1163
1164 /* ------ doit */
1165 if (_tenEpiRegWarp(nout, npxfr, nhst, ngrad, nin, ninLen,
1166 reference, kern, kparm, verbose)) {
1167 biffAddf(TENtenBiffKey, "%s: trouble performing final registration", me);
1168 airMopError(mop); return 1;
1169 }
1170
1171 airMopOkay(mop);
1172 nrrdStateAlwaysSetContent = hack1;
1173 nrrdStateDisableContent = hack2;
1174 return 0;
1175}
1176
1177int
1178tenEpiRegister4D(Nrrd *_nout, Nrrd *_nin, Nrrd *_ngrad,
1179 int reference,
1180 double bwX, double bwY, double fitFrac,
1181 double DWthr, int doCC,
1182 const NrrdKernel *kern, double *kparm,
1183 int progress, int verbose) {
1184 static const char me[]="tenEpiRegister4D";
1185 unsigned int ninIdx, ninLen,
1186 dwiAx, rangeAxisNum, rangeAxisIdx[NRRD_DIM_MAX16];
1187 int dwiIdx;
1188 Nrrd **nout, **nin, **ndwi, **ndwiOut, *ngrad, *ndwigrad;
1189 airArray *mop;
1190 double *grad, *dwigrad, glen;
1191
1192 if (!(_nout && _nin)) {
1193 biffAddf(TENtenBiffKey, "%s: got NULL pointer", me);
1194 return 1;
1195 }
1196 if (4 != _nin->dim) {
1197 biffAddf(TENtenBiffKey, "%s: need a 4-D input array, not %d-D", me, _nin->dim);
1198 return 1;
1199 }
1200 if (tenGradientCheck(_ngrad, nrrdTypeDefault, 6 /* <-- HEY: not sure */)) {
1201 biffAddf(TENtenBiffKey, "%s: problem with given gradient list", me);
1202 return 1;
1203 }
1204
1205 rangeAxisNum = nrrdRangeAxesGet(_nin, rangeAxisIdx);
1206 if (0 == rangeAxisNum) {
1207 /* we fall back on old behavior */
1208 dwiAx = 0;
1209 } else if (1 == rangeAxisNum) {
1210 /* thankfully there's exactly one range axis */
1211 dwiAx = rangeAxisIdx[0];
1212 } else {
1213 biffAddf(TENtenBiffKey, "%s: have %u range axes instead of 1, don't know which "
1214 "is DWI axis", me, rangeAxisNum);
1215 return 1;
1216 }
1217
1218 ninLen = _nin->axis[dwiAx].size;
1219 /* outdated
1220 if (!( AIR_IN_CL(6, ninLen, 120) )) {
1221 biffAddf(TEN, "%s: %u (size of axis %u, and # DWIs) is unreasonable",
1222 me, ninLen, dwiAx);
1223 return 1;
1224 }
1225 */
1226 if (ninLen != _ngrad->axis[1].size) {
1227 char stmp[AIR_STRLEN_SMALL(128+1)];
1228 biffAddf(TENtenBiffKey, "%s: ninLen %u != # grads %s", me, ninLen,
1229 airSprintSize_t(stmp, _ngrad->axis[1].size));
1230 return 1;
1231 }
1232
1233 mop = airMopNew();
1234 ngrad = nrrdNew();
1235 ndwigrad = nrrdNew();
1236 airMopAdd(mop, ngrad, (airMopper)nrrdNuke, airMopAlways);
1237 airMopAdd(mop, ndwigrad, (airMopper)nrrdNuke, airMopAlways);
1238 if (nrrdConvert(ngrad, _ngrad, nrrdTypeDouble)
1239 || nrrdConvert(ndwigrad, _ngrad, nrrdTypeDouble)) { /* HACK applies */
1240 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: trouble converting gradients to doubles", me);
1241 airMopError(mop); return 1;
1242 }
1243
1244 nin = (Nrrd **)calloc(ninLen, sizeof(Nrrd*));
1245 ndwi = (Nrrd **)calloc(ninLen, sizeof(Nrrd*));
1246 nout = (Nrrd **)calloc(ninLen, sizeof(Nrrd*));
1247 ndwiOut = (Nrrd **)calloc(ninLen, sizeof(Nrrd*));
1248 if (!(nin && ndwi && nout && ndwiOut)) {
1249 biffAddf(TENtenBiffKey, "%s: trouble allocating local arrays", me);
1250 airMopError(mop); return 1;
1251 }
1252 airMopAdd(mop, nin, airFree, airMopAlways);
1253 airMopAdd(mop, ndwi, airFree, airMopAlways);
1254 airMopAdd(mop, nout, airFree, airMopAlways);
1255 airMopAdd(mop, ndwiOut, airFree, airMopAlways);
1256 dwiIdx = -1;
1257 grad = AIR_CAST(double *, ngrad->data)((double *)(ngrad->data));
1258 dwigrad = AIR_CAST(double *, ndwigrad->data)((double *)(ndwigrad->data));
1259 for (ninIdx=0; ninIdx<ninLen; ninIdx++) {
1260 glen = ELL_3V_LEN(grad + 3*ninIdx)(sqrt((((grad + 3*ninIdx))[0]*((grad + 3*ninIdx))[0] + ((grad
+ 3*ninIdx))[1]*((grad + 3*ninIdx))[1] + ((grad + 3*ninIdx))
[2]*((grad + 3*ninIdx))[2])));
1261 if (-1 != reference || glen > 0.0) {
1262 /* we're not allowed to use only those images that are truly DWIs,
1263 or,
1264 (we are so allowed and) this is a true DWI */
1265 dwiIdx++;
1266 ndwi[dwiIdx] = nrrdNew();
1267 ndwiOut[dwiIdx] = nrrdNew();
1268 airMopAdd(mop, ndwi[dwiIdx], (airMopper)nrrdNuke, airMopAlways);
1269 airMopAdd(mop, ndwiOut[dwiIdx], (airMopper)nrrdNuke, airMopAlways);
1270 if (nrrdSlice(ndwi[dwiIdx], _nin, dwiAx,
1271 AIR_CAST(unsigned int, ninIdx)((unsigned int)(ninIdx)))) {
1272 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: trouble slicing at %d on axis %u",
1273 me, ninIdx, dwiAx);
1274 airMopError(mop); return 1;
1275 }
1276 /* NOTE: this works because dwiIdx <= ninIdx */
1277 ELL_3V_COPY(dwigrad + 3*dwiIdx, grad + 3*ninIdx)((dwigrad + 3*dwiIdx)[0] = (grad + 3*ninIdx)[0], (dwigrad + 3
*dwiIdx)[1] = (grad + 3*ninIdx)[1], (dwigrad + 3*dwiIdx)[2] =
(grad + 3*ninIdx)[2]);
1278 } else {
1279 /* we are only looking at true DWIs, and this isn't one of them */
1280 continue;
1281 }
1282 }
1283 if (-1 == dwiIdx) {
1284 biffAddf(TENtenBiffKey, "%s: somehow got no DWIs", me);
1285 airMopError(mop); return 1;
1286 }
1287 /* HEY: HACK! */
1288 ndwigrad->axis[1].size = 1 + AIR_CAST(unsigned int, dwiIdx)((unsigned int)(dwiIdx));
1289 if (tenEpiRegister3D(ndwiOut, ndwi, ndwigrad->axis[1].size, ndwigrad,
1290 reference,
1291 bwX, bwY, fitFrac, DWthr,
1292 doCC,
1293 kern, kparm,
1294 progress, verbose)) {
1295 biffAddf(TENtenBiffKey, "%s: trouble", me);
1296 airMopError(mop); return 1;
1297 }
1298 dwiIdx = -1;
1299 for (ninIdx=0; ninIdx<ninLen; ninIdx++) {
1300 glen = ELL_3V_LEN(grad + 3*ninIdx)(sqrt((((grad + 3*ninIdx))[0]*((grad + 3*ninIdx))[0] + ((grad
+ 3*ninIdx))[1]*((grad + 3*ninIdx))[1] + ((grad + 3*ninIdx))
[2]*((grad + 3*ninIdx))[2])));
1301 if (-1 != reference || glen > 0.0) {
1302 /* we're not allowed to use only those images that are truly DWIs,
1303 or,
1304 (we are so allowed and) this is a true DWI */
1305 dwiIdx++;
1306 nout[ninIdx] = ndwiOut[dwiIdx];
1307 /*
1308 fprintf(stderr, "!%s: (A) nout[%u] = %p\n", me, ninIdx, nout[ninIdx]);
1309 */
1310 } else {
1311 /* we are only looking at true DWIs, and this isn't one of them */
1312 nout[ninIdx] = nrrdNew();
1313 airMopAdd(mop, nout[ninIdx], (airMopper)nrrdNuke, airMopAlways);
1314 if (nrrdSlice(nout[ninIdx], _nin, dwiAx, ninIdx)) {
1315 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: trouble slicing at %d on axis %u",
1316 me, dwiIdx, dwiAx);
1317 airMopError(mop); return 1;
1318 }
1319 /*
1320 fprintf(stderr, "!%s: (B) nout[%u] = %p\n", me, ninIdx, nout[ninIdx]);
1321 */
1322 }
1323 }
1324 if (nrrdJoin(_nout, (const Nrrd*const*)nout, ninLen, dwiAx, AIR_TRUE1)) {
1325 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: trouble joining output", me);
1326 airMopError(mop); return 1;
1327 }
1328 nrrdAxisInfoCopy(_nout, _nin, NULL((void*)0), NRRD_AXIS_INFO_NONE0);
1329 if (nrrdBasicInfoCopy(_nout, _nin,
1330 NRRD_BASIC_INFO_DATA_BIT(1<< 1)
1331 | NRRD_BASIC_INFO_TYPE_BIT(1<< 2)
1332 | NRRD_BASIC_INFO_BLOCKSIZE_BIT(1<< 3)
1333 | NRRD_BASIC_INFO_DIMENSION_BIT(1<< 4)
1334 | NRRD_BASIC_INFO_CONTENT_BIT(1<< 5)
1335 | NRRD_BASIC_INFO_COMMENTS_BIT(1<<14))) {
1336 /* note that we're ALWAYS copying the key/value pairs- its just
1337 too annoying to have to set nrrdStateKeyValuePairsPropagate
1338 in order for the DWI-specific key/value pairs to be set */
1339 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s:", me);
1340 airMopError(mop); return 1;
1341 }
1342
1343 airMopOkay(mop);
1344 return 0;
1345}