1 /*
2 Teem: Tools to process and visualize scientific data and images .
3 Copyright (C) 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
27 /*
28 ** learned: when it looks like good-old LLS estimation is producing
29 ** nothing but zero tensors, see if your tec->valueMin is larger
30 ** than (what are problably) floating-point DWIs
31 */
32
33 /*
34
35 http://www.mathworks.com/access/helpdesk/help/toolbox/curvefit/ch_fitt5.html#40515
36
37 */
38
39 /* ---------------------------------------------- */
40
41 int
_tenGaussian(double * retP,double m,double t,double s)42 _tenGaussian(double *retP, double m, double t, double s) {
43 static const char me[]="_tenGaussian";
44 double diff, earg, den;
45
46 if (!retP) {
47 biffAddf(TEN, "%s: got NULL pointer", me);
48 return 1;
49 }
50 diff = (m-t)/2;
51 earg = -diff*diff/2;
52 den = s*sqrt(2*AIR_PI);
53 *retP = exp(earg)/den;
54 if (!AIR_EXISTS(*retP)) {
55 biffAddf(TEN, "%s: m=%g, t=%g, s=%g", me, m, t, s);
56 biffAddf(TEN, "%s: diff=%g, earg=%g, den=%g", me, diff, earg, den);
57 biffAddf(TEN, "%s: failed with ret = exp(%g)/%g = %g/%g = %g",
58 me, earg, den, exp(earg), den, *retP);
59 *retP = AIR_NAN; return 1;
60 }
61 return 0;
62 }
63
64 int
_tenRicianTrue(double * retP,double m,double t,double s)65 _tenRicianTrue(double *retP,
66 double m /* measured */,
67 double t /* truth */,
68 double s /* sigma */) {
69 static const char me[]="_tenRicianTrue";
70 double mos, moss, mos2, tos2, tos, ss, earg, barg;
71
72 if (!retP) {
73 biffAddf(TEN, "%s: got NULL pointer", me);
74 return 1;
75 }
76
77 mos = m/s;
78 moss = mos/s;
79 tos = t/s;
80 ss = s*s;
81 mos2 = mos*mos;
82 tos2 = tos*tos;
83 earg = -(mos2 + tos2)/2;
84 barg = mos*tos;
85 *retP = exp(earg)*airBesselI0(barg)*moss;
86
87 if (!AIR_EXISTS(*retP)) {
88 biffAddf(TEN, "%s: m=%g, t=%g, s=%g", me, m, t, s);
89 biffAddf(TEN, "%s: mos=%g, moss=%g, tos=%g, ss=%g",
90 me, mos, moss, tos, ss);
91 biffAddf(TEN, "%s: mos2=%g, tos2=%g, earg=%g, barg=%g", me,
92 mos2, tos2, earg, barg);
93 biffAddf(TEN, "%s: failed: ret=exp(%g)*bessi0(%g)*%g = %g * %g * %g = %g",
94 me, earg, barg, moss, exp(earg), airBesselI0(barg), moss, *retP);
95 *retP = AIR_NAN; return 1;
96 }
97 return 0;
98 }
99
100 int
_tenRicianSafe(double * retP,double m,double t,double s)101 _tenRicianSafe(double *retP, double m, double t, double s) {
102 static const char me[]="_tenRicianSafe";
103 double diff, ric, gau, neer=10, faar=20;
104 int E;
105
106 if (!retP) {
107 biffAddf(TEN, "%s: got NULL pointer", me);
108 return 1;
109 }
110
111 diff = AIR_ABS(m-t)/s;
112 E = 0;
113 if (diff < neer) {
114 if (!E) E |= _tenRicianTrue(retP, m, t, s);
115 } else if (diff < faar) {
116 if (!E) E |= _tenRicianTrue(&ric, m, t, s);
117 if (!E) E |= _tenGaussian(&gau, m, t, s);
118 if (!E) *retP = AIR_AFFINE(neer, diff, faar, ric, gau);
119 } else {
120 if (!E) E |= _tenGaussian(retP, m, t, s);
121 }
122 if (E) {
123 biffAddf(TEN, "%s: failed with m=%g, t=%g, s=%g -> diff=%g",
124 me, m, t, s, diff);
125 *retP = AIR_NAN; return 1;
126 }
127 return 0;
128 }
129
130 int
_tenRician(double * retP,double m,double t,double s)131 _tenRician(double *retP,
132 double m /* measured */,
133 double t /* truth */,
134 double s /* sigma */) {
135 static const char me[]="_tenRician";
136 double tos, ric, gau, loSignal=4.0, hiSignal=8.0;
137 int E;
138
139 if (!retP) {
140 biffAddf(TEN, "%s: got NULL pointer", me);
141 return 1;
142 }
143 if (!( m >= 0 && t >= 0 && s > 0 )) {
144 biffAddf(TEN, "%s: got bad args: m=%g t=%g s=%g", me, m, t, s);
145 *retP = AIR_NAN; return 1;
146 }
147
148 tos = t/s;
149 E = 0;
150 if (tos < loSignal) {
151 if (!E) E |= _tenRicianSafe(retP, m, t, s);
152 } else if (tos < hiSignal) {
153 if (!E) E |= _tenRicianSafe(&ric, m, t, s);
154 if (!E) E |= _tenGaussian(&gau, m, t, s);
155 if (!E) *retP = AIR_AFFINE(loSignal, tos, hiSignal, ric, gau);
156 } else {
157 if (!E) E |= _tenGaussian(retP, m, t, s);
158 }
159 if (E) {
160 biffAddf(TEN, "%s: failed with m=%g, t=%g, s=%g -> tos=%g",
161 me, m, t, s, tos);
162 *retP = AIR_NAN; return 1;
163 }
164 return 0;
165 }
166
167 enum {
168 flagUnknown,
169 flagEstimateMethod,
170 flagBInfo,
171 flagAllNum,
172 flagDwiNum,
173 flagAllAlloc,
174 flagDwiAlloc,
175 flagAllSet,
176 flagDwiSet,
177 flagSkipSet,
178 flagWght,
179 flagEmat,
180 flagLast
181 };
182
183 void
_tenEstimateOutputInit(tenEstimateContext * tec)184 _tenEstimateOutputInit(tenEstimateContext *tec) {
185
186 tec->estimatedB0 = AIR_NAN;
187 TEN_T_SET(tec->ten, AIR_NAN,
188 AIR_NAN, AIR_NAN, AIR_NAN,
189 AIR_NAN, AIR_NAN,
190 AIR_NAN);
191 tec->conf = AIR_NAN;
192 tec->mdwi = AIR_NAN;
193 tec->time = AIR_NAN;
194 tec->errorDwi = AIR_NAN;
195 tec->errorLogDwi = AIR_NAN;
196 tec->likelihoodDwi = AIR_NAN;
197 }
198
199 tenEstimateContext *
tenEstimateContextNew()200 tenEstimateContextNew() {
201 tenEstimateContext *tec;
202 unsigned int fi;
203 airPtrPtrUnion appu;
204
205
206 tec = AIR_CAST(tenEstimateContext *, malloc(sizeof(tenEstimateContext)));
207 if (tec) {
208 tec->bValue = AIR_NAN;
209 tec->valueMin = AIR_NAN;
210 tec->sigma = AIR_NAN;
211 tec->dwiConfThresh = AIR_NAN;
212 tec->dwiConfSoft = AIR_NAN;
213 tec->_ngrad = NULL;
214 tec->_nbmat = NULL;
215 tec->skipList = NULL;
216 appu.ui = &(tec->skipList);
217 tec->skipListArr = airArrayNew(appu.v, NULL,
218 2*sizeof(unsigned int), 128);
219 tec->skipListArr->noReallocWhenSmaller = AIR_TRUE;
220 tec->all_f = NULL;
221 tec->all_d = NULL;
222 tec->simulate = AIR_FALSE;
223 tec->estimate1Method = tenEstimate1MethodUnknown;
224 tec->estimateB0 = AIR_TRUE;
225 tec->recordTime = AIR_FALSE;
226 tec->recordErrorDwi = AIR_FALSE;
227 tec->recordErrorLogDwi = AIR_FALSE;
228 tec->recordLikelihoodDwi = AIR_FALSE;
229 tec->verbose = 0;
230 tec->progress = AIR_FALSE;
231 tec->WLSIterNum = 3;
232 for (fi=flagUnknown+1; fi<flagLast; fi++) {
233 tec->flag[fi] = AIR_FALSE;
234 }
235 tec->allNum = 0;
236 tec->dwiNum = 0;
237 tec->nbmat = nrrdNew();
238 tec->nwght = nrrdNew();
239 tec->nemat = nrrdNew();
240 tec->knownB0 = AIR_NAN;
241 tec->all = NULL;
242 tec->bnorm = NULL;
243 tec->allTmp = NULL;
244 tec->dwiTmp = NULL;
245 tec->dwi = NULL;
246 tec->skipLut = NULL;
247 _tenEstimateOutputInit(tec);
248 }
249 return tec;
250 }
251
252 tenEstimateContext *
tenEstimateContextNix(tenEstimateContext * tec)253 tenEstimateContextNix(tenEstimateContext *tec) {
254
255 if (tec) {
256 nrrdNuke(tec->nbmat);
257 nrrdNuke(tec->nwght);
258 nrrdNuke(tec->nemat);
259 airArrayNuke(tec->skipListArr);
260 airFree(tec->all);
261 airFree(tec->bnorm);
262 airFree(tec->allTmp);
263 airFree(tec->dwiTmp);
264 airFree(tec->dwi);
265 airFree(tec->skipLut);
266 airFree(tec);
267 }
268 return NULL;
269 }
270
271 void
tenEstimateVerboseSet(tenEstimateContext * tec,int verbose)272 tenEstimateVerboseSet(tenEstimateContext *tec,
273 int verbose) {
274 if (tec) {
275 tec->verbose = verbose;
276 }
277 return;
278 }
279
280 void
tenEstimateNegEvalShiftSet(tenEstimateContext * tec,int doit)281 tenEstimateNegEvalShiftSet(tenEstimateContext *tec, int doit) {
282
283 if (tec) {
284 tec->negEvalShift = !!doit;
285 }
286 return;
287 }
288
289 int
tenEstimateMethodSet(tenEstimateContext * tec,int estimateMethod)290 tenEstimateMethodSet(tenEstimateContext *tec, int estimateMethod) {
291 static const char me[]="tenEstimateMethodSet";
292
293 if (!tec) {
294 biffAddf(TEN, "%s: got NULL pointer", me);
295 return 1;
296 }
297 if (airEnumValCheck(tenEstimate1Method, estimateMethod)) {
298 biffAddf(TEN, "%s: estimateMethod %d not a valid %s", me,
299 estimateMethod, tenEstimate1Method->name);
300 return 1;
301 }
302
303 if (tec->estimate1Method != estimateMethod) {
304 tec->estimate1Method = estimateMethod;
305 tec->flag[flagEstimateMethod] = AIR_TRUE;
306 }
307
308 return 0;
309 }
310
311 int
tenEstimateSigmaSet(tenEstimateContext * tec,double sigma)312 tenEstimateSigmaSet(tenEstimateContext *tec, double sigma) {
313 static const char me[]="tenEstimateSigmaSet";
314
315 if (!tec) {
316 biffAddf(TEN, "%s: got NULL pointer", me);
317 return 1;
318 }
319 if (!(AIR_EXISTS(sigma) && sigma >= 0.0)) {
320 biffAddf(TEN, "%s: given sigma (%g) not existent and >= 0.0", me, sigma);
321 return 1;
322 }
323
324 tec->sigma = sigma;
325
326 return 0;
327 }
328
329 int
tenEstimateValueMinSet(tenEstimateContext * tec,double valueMin)330 tenEstimateValueMinSet(tenEstimateContext *tec, double valueMin) {
331 static const char me[]="tenEstimateValueMinSet";
332
333 if (!tec) {
334 biffAddf(TEN, "%s: got NULL pointer", me);
335 return 1;
336 }
337 if (!(AIR_EXISTS(valueMin) && valueMin > 0.0)) {
338 biffAddf(TEN, "%s: given valueMin (%g) not existent and > 0.0",
339 me, valueMin);
340 return 1;
341 }
342
343 tec->valueMin = valueMin;
344
345 return 0;
346 }
347
348 int
tenEstimateGradientsSet(tenEstimateContext * tec,const Nrrd * ngrad,double bValue,int estimateB0)349 tenEstimateGradientsSet(tenEstimateContext *tec,
350 const Nrrd *ngrad, double bValue, int estimateB0) {
351 static const char me[]="tenEstimateGradientsSet";
352
353 if (!(tec && ngrad)) {
354 biffAddf(TEN, "%s: got NULL pointer", me);
355 return 1;
356 }
357 if (!AIR_EXISTS(bValue)) {
358 biffAddf(TEN, "%s: given b value doesn't exist", me);
359 return 1;
360 }
361 if (tenGradientCheck(ngrad, nrrdTypeDefault, 7)) {
362 biffAddf(TEN, "%s: problem with gradient list", me);
363 return 1;
364 }
365
366 tec->bValue = bValue;
367 tec->_ngrad = ngrad;
368 tec->_nbmat = NULL;
369 tec->estimateB0 = estimateB0;
370
371 tec->flag[flagBInfo] = AIR_TRUE;
372 return 0;
373 }
374
375 int
tenEstimateBMatricesSet(tenEstimateContext * tec,const Nrrd * nbmat,double bValue,int estimateB0)376 tenEstimateBMatricesSet(tenEstimateContext *tec,
377 const Nrrd *nbmat, double bValue, int estimateB0) {
378 static const char me[]="tenEstimateBMatricesSet";
379
380 if (!(tec && nbmat)) {
381 biffAddf(TEN, "%s: got NULL pointer", me);
382 return 1;
383 }
384 if (!AIR_EXISTS(bValue)) {
385 biffAddf(TEN, "%s: given b value doesn't exist", me);
386 return 1;
387 }
388 if (tenBMatrixCheck(nbmat, nrrdTypeDefault, 7)) {
389 biffAddf(TEN, "%s: problem with b-matrix list", me);
390 return 1;
391 }
392
393 tec->bValue = bValue;
394 tec->_ngrad = NULL;
395 tec->_nbmat = nbmat;
396 tec->estimateB0 = estimateB0;
397
398 tec->flag[flagBInfo] = AIR_TRUE;
399 return 0;
400 }
401
402 int
tenEstimateSkipSet(tenEstimateContext * tec,unsigned int valIdx,int doSkip)403 tenEstimateSkipSet(tenEstimateContext *tec,
404 unsigned int valIdx, int doSkip) {
405 static const char me[]="tenEstimateSkipSet";
406 unsigned int skipIdx;
407
408 if (!tec) {
409 biffAddf(TEN, "%s: got NULL pointer", me);
410 return 1;
411 }
412
413 skipIdx = airArrayLenIncr(tec->skipListArr, 1);
414 tec->skipList[0 + 2*skipIdx] = valIdx;
415 tec->skipList[1 + 2*skipIdx] = !!doSkip;
416
417 tec->flag[flagSkipSet] = AIR_TRUE;
418 return 0;
419 }
420
421 int
tenEstimateSkipReset(tenEstimateContext * tec)422 tenEstimateSkipReset(tenEstimateContext *tec) {
423 static const char me[]="tenEstimateSkipReset";
424
425 if (!tec) {
426 biffAddf(TEN, "%s: got NULL pointer", me);
427 return 1;
428 }
429
430 airArrayLenSet(tec->skipListArr, 0);
431
432 tec->flag[flagSkipSet] = AIR_TRUE;
433 return 0;
434 }
435
436 int
tenEstimateThresholdSet(tenEstimateContext * tec,double thresh,double soft)437 tenEstimateThresholdSet(tenEstimateContext *tec,
438 double thresh, double soft) {
439 static const char me[]="tenEstimateThresholdSet";
440
441 if (!tec) {
442 biffAddf(TEN, "%s: got NULL pointer", me);
443 return 1;
444 }
445 if (!(AIR_EXISTS(thresh) && AIR_EXISTS(soft))) {
446 biffAddf(TEN, "%s: not both threshold (%g) and softness (%g) exist", me,
447 thresh, soft);
448 return 1;
449 }
450
451 tec->dwiConfThresh = thresh;
452 tec->dwiConfSoft = soft;
453
454 return 0;
455 }
456
457 int
_tenEstimateCheck(tenEstimateContext * tec)458 _tenEstimateCheck(tenEstimateContext *tec) {
459 static const char me[]="_tenEstimateCheck";
460
461 if (!tec) {
462 biffAddf(TEN, "%s: got NULL pointer", me);
463 return 1;
464 }
465 if (!( AIR_EXISTS(tec->valueMin) && tec->valueMin > 0.0)) {
466 biffAddf(TEN, "%s: need a positive valueMin set (not %g)",
467 me, tec->valueMin);
468 return 1;
469 }
470 if (!tec->simulate) {
471 if (!AIR_EXISTS(tec->bValue)) {
472 biffAddf(TEN, "%s: b-value not set", me);
473 return 1;
474 }
475 if (airEnumValCheck(tenEstimate1Method, tec->estimate1Method)) {
476 biffAddf(TEN, "%s: estimation method not set", me);
477 return 1;
478 }
479 if (tenEstimate1MethodMLE == tec->estimate1Method
480 && !(AIR_EXISTS(tec->sigma) && (tec->sigma >= 0.0))
481 ) {
482 biffAddf(TEN, "%s: can't do %s estim w/out non-negative sigma set", me,
483 airEnumStr(tenEstimate1Method, tenEstimate1MethodMLE));
484 return 1;
485 }
486 if (!(AIR_EXISTS(tec->dwiConfThresh) && AIR_EXISTS(tec->dwiConfSoft))) {
487 biffAddf(TEN, "%s: not both threshold (%g) and softness (%g) exist", me,
488 tec->dwiConfThresh, tec->dwiConfSoft);
489 return 1;
490 }
491 }
492 if (!( tec->_ngrad || tec->_nbmat )) {
493 biffAddf(TEN, "%s: need to set either gradients or B-matrices", me);
494 return 1;
495 }
496
497 return 0;
498 }
499
500 /*
501 ** allNum includes the skipped images
502 ** dwiNum does not include the skipped images
503 */
504 int
_tenEstimateNumUpdate(tenEstimateContext * tec)505 _tenEstimateNumUpdate(tenEstimateContext *tec) {
506 static const char me[]="_tenEstimateNumUpdate";
507 unsigned int newAllNum, newDwiNum, allIdx,
508 skipListIdx, skipIdx, skipDo, skipNotNum;
509 double (*lup)(const void *, size_t), gg[3], bb[6];
510
511 if (tec->flag[flagBInfo]
512 || tec->flag[flagSkipSet]) {
513 if (tec->_ngrad) {
514 newAllNum = AIR_CAST(unsigned int, tec->_ngrad->axis[1].size);
515 lup = nrrdDLookup[tec->_ngrad->type];
516 } else {
517 newAllNum = AIR_CAST(unsigned int, tec->_nbmat->axis[1].size);
518 lup = nrrdDLookup[tec->_nbmat->type];
519 }
520 if (tec->allNum != newAllNum) {
521 tec->allNum = newAllNum;
522 tec->flag[flagAllNum] = AIR_TRUE;
523 }
524
525 /* HEY: this should probably be its own update function, but its very
526 convenient to allocate these allNum-length arrays here, immediately */
527 airFree(tec->skipLut);
528 tec->skipLut = AIR_CAST(unsigned char *, calloc(tec->allNum,
529 sizeof(unsigned char)));
530 airFree(tec->bnorm);
531 tec->bnorm = AIR_CAST(double *, calloc(tec->allNum, sizeof(double)));
532 if (!(tec->skipLut && tec->bnorm)) {
533 biffAddf(TEN, "%s: couldn't allocate skipLut, bnorm vectors length %u\n",
534 me, tec->allNum);
535 return 1;
536 }
537
538 for (skipListIdx=0; skipListIdx<tec->skipListArr->len; skipListIdx++) {
539 skipIdx = tec->skipList[0 + 2*skipListIdx];
540 skipDo = tec->skipList[1 + 2*skipListIdx];
541 if (!(skipIdx < tec->allNum)) {
542 biffAddf(TEN, "%s: skipList entry %u value index %u not < # vals %u",
543 me, skipListIdx, skipIdx, tec->allNum);
544 return 1;
545 }
546 tec->skipLut[skipIdx] = skipDo;
547 }
548 skipNotNum = 0;
549 for (skipIdx=0; skipIdx<tec->allNum; skipIdx++) {
550 skipNotNum += !tec->skipLut[skipIdx];
551 }
552 if (!(skipNotNum >= 7 )) {
553 biffAddf(TEN, "%s: number of not-skipped (used) values %u < minimum 7",
554 me, skipNotNum);
555 return 1;
556 }
557
558 newDwiNum = 0;
559 for (allIdx=0; allIdx<tec->allNum; allIdx++) {
560 if (tec->skipLut[allIdx]) {
561 tec->bnorm[allIdx] = AIR_NAN;
562 } else {
563 if (tec->_ngrad) {
564 gg[0] = lup(tec->_ngrad->data, 0 + 3*allIdx);
565 gg[1] = lup(tec->_ngrad->data, 1 + 3*allIdx);
566 gg[2] = lup(tec->_ngrad->data, 2 + 3*allIdx);
567 bb[0] = gg[0]*gg[0];
568 bb[1] = gg[1]*gg[0];
569 bb[2] = gg[2]*gg[0];
570 bb[3] = gg[1]*gg[1];
571 bb[4] = gg[2]*gg[1];
572 bb[5] = gg[2]*gg[2];
573 } else {
574 bb[0] = lup(tec->_nbmat->data, 0 + 6*allIdx);
575 bb[1] = lup(tec->_nbmat->data, 1 + 6*allIdx);
576 bb[2] = lup(tec->_nbmat->data, 2 + 6*allIdx);
577 bb[3] = lup(tec->_nbmat->data, 3 + 6*allIdx);
578 bb[4] = lup(tec->_nbmat->data, 4 + 6*allIdx);
579 bb[5] = lup(tec->_nbmat->data, 5 + 6*allIdx);
580 }
581 tec->bnorm[allIdx] = sqrt(bb[0]*bb[0] + 2*bb[1]*bb[1] + 2*bb[2]*bb[2]
582 + bb[3]*bb[3] + 2*bb[4]*bb[4]
583 + bb[5]*bb[5]);
584 if (tec->estimateB0) {
585 ++newDwiNum;
586 } else {
587 newDwiNum += (0.0 != tec->bnorm[allIdx]);
588 }
589 }
590 }
591 if (tec->dwiNum != newDwiNum) {
592 tec->dwiNum = newDwiNum;
593 tec->flag[flagDwiNum] = AIR_TRUE;
594 }
595 if (!tec->estimateB0 && (tec->allNum == tec->dwiNum)) {
596 biffAddf(TEN, "%s: don't want to estimate B0, but all values are DW", me);
597 return 1;
598 }
599 }
600 return 0;
601 }
602
603 int
_tenEstimateAllAllocUpdate(tenEstimateContext * tec)604 _tenEstimateAllAllocUpdate(tenEstimateContext *tec) {
605 static const char me[]="_tenEstimateAllAllocUpdate";
606
607 if (tec->flag[flagAllNum]) {
608 airFree(tec->all);
609 airFree(tec->allTmp);
610 tec->all = AIR_CAST(double *, calloc(tec->allNum, sizeof(double)));
611 tec->allTmp = AIR_CAST(double *, calloc(tec->allNum, sizeof(double)));
612 if (!( tec->all && tec->allTmp )) {
613 biffAddf(TEN, "%s: couldn't allocate \"all\" arrays (length %u)", me,
614 tec->allNum);
615 return 1;
616 }
617 tec->flag[flagAllAlloc] = AIR_TRUE;
618 }
619 return 0;
620 }
621
622 int
_tenEstimateDwiAllocUpdate(tenEstimateContext * tec)623 _tenEstimateDwiAllocUpdate(tenEstimateContext *tec) {
624 static const char me[]="_tenEstimateDwiAllocUpdate";
625 size_t size[2];
626 int E;
627
628 if (tec->flag[flagDwiNum]) {
629 airFree(tec->dwi);
630 airFree(tec->dwiTmp);
631 tec->dwi = AIR_CAST(double *, calloc(tec->dwiNum, sizeof(double)));
632 tec->dwiTmp = AIR_CAST(double *, calloc(tec->dwiNum, sizeof(double)));
633 if (!(tec->dwi && tec->dwiTmp)) {
634 biffAddf(TEN, "%s: couldn't allocate DWI arrays (length %u)", me,
635 tec->dwiNum);
636 return 1;
637 }
638 E = 0;
639 if (!E) size[0] = (tec->estimateB0 ? 7 : 6);
640 if (!E) size[1] = tec->dwiNum;
641 if (!E) E |= nrrdMaybeAlloc_nva(tec->nbmat, nrrdTypeDouble, 2, size);
642 if (!E) size[0] = tec->dwiNum;
643 if (!E) size[1] = tec->dwiNum;
644 if (!E) E |= nrrdMaybeAlloc_nva(tec->nwght, nrrdTypeDouble, 2, size);
645 if (E) {
646 biffMovef(TEN, NRRD, "%s: couldn't allocate dwi nrrds", me);
647 return 1;
648 }
649 /* nrrdSave("0-nbmat.txt", tec->nbmat, NULL); */
650 tec->flag[flagDwiAlloc] = AIR_TRUE;
651 }
652 return 0;
653 }
654
655 int
_tenEstimateAllSetUpdate(tenEstimateContext * tec)656 _tenEstimateAllSetUpdate(tenEstimateContext *tec) {
657 /* static const char me[]="_tenEstimateAllSetUpdate"; */
658 /* unsigned int skipListIdx, skipIdx, skip, dwiIdx */;
659
660 if (tec->flag[flagAllAlloc]
661 || tec->flag[flagDwiNum]) {
662
663 }
664 return 0;
665 }
666
667 int
_tenEstimateDwiSetUpdate(tenEstimateContext * tec)668 _tenEstimateDwiSetUpdate(tenEstimateContext *tec) {
669 /* static const char me[]="_tenEstimateDwiSetUpdate"; */
670 double (*lup)(const void *, size_t I), gg[3], *bmat;
671 unsigned int allIdx, dwiIdx, bmIdx;
672
673 if (tec->flag[flagAllNum]
674 || tec->flag[flagDwiAlloc]) {
675 if (tec->_ngrad) {
676 lup = nrrdDLookup[tec->_ngrad->type];
677 } else {
678 lup = nrrdDLookup[tec->_nbmat->type];
679 }
680 dwiIdx = 0;
681 bmat = AIR_CAST(double*, tec->nbmat->data);
682 for (allIdx=0; allIdx<tec->allNum; allIdx++) {
683 if (!tec->skipLut[allIdx]
684 && (tec->estimateB0 || tec->bnorm[allIdx])) {
685 if (tec->_ngrad) {
686 gg[0] = lup(tec->_ngrad->data, 0 + 3*allIdx);
687 gg[1] = lup(tec->_ngrad->data, 1 + 3*allIdx);
688 gg[2] = lup(tec->_ngrad->data, 2 + 3*allIdx);
689 bmat[0] = gg[0]*gg[0];
690 bmat[1] = gg[1]*gg[0];
691 bmat[2] = gg[2]*gg[0];
692 bmat[3] = gg[1]*gg[1];
693 bmat[4] = gg[2]*gg[1];
694 bmat[5] = gg[2]*gg[2];
695 } else {
696 for (bmIdx=0; bmIdx<6; bmIdx++) {
697 bmat[bmIdx] = lup(tec->_nbmat->data, bmIdx + 6*allIdx);
698 }
699 }
700 bmat[1] *= 2.0;
701 bmat[2] *= 2.0;
702 bmat[4] *= 2.0;
703 if (tec->estimateB0) {
704 bmat[6] = -1;
705 }
706 bmat += tec->nbmat->axis[0].size;
707 dwiIdx++;
708 }
709 }
710 }
711 return 0;
712 }
713
714 int
_tenEstimateWghtUpdate(tenEstimateContext * tec)715 _tenEstimateWghtUpdate(tenEstimateContext *tec) {
716 /* static const char me[]="_tenEstimateWghtUpdate"; */
717 unsigned int dwiIdx;
718 double *wght;
719
720 wght = AIR_CAST(double *, tec->nwght->data);
721 if (tec->flag[flagDwiAlloc]
722 || tec->flag[flagEstimateMethod]) {
723
724 /* HEY: this is only useful for linear LS, no? */
725 for (dwiIdx=0; dwiIdx<tec->dwiNum; dwiIdx++) {
726 wght[dwiIdx + tec->dwiNum*dwiIdx] = 1.0;
727 }
728
729 tec->flag[flagEstimateMethod] = AIR_FALSE;
730 tec->flag[flagWght] = AIR_TRUE;
731 }
732 return 0;
733 }
734
735 int
_tenEstimateEmatUpdate(tenEstimateContext * tec)736 _tenEstimateEmatUpdate(tenEstimateContext *tec) {
737 static const char me[]="tenEstimateEmatUpdate";
738
739 if (tec->flag[flagDwiSet]
740 || tec->flag[flagWght]) {
741
742 if (!tec->simulate) {
743 /* HEY: ignores weights! */
744 if (ell_Nm_pseudo_inv(tec->nemat, tec->nbmat)) {
745 biffMovef(TEN, ELL, "%s: trouble pseudo-inverting %ux%u B-matrix", me,
746 AIR_CAST(unsigned int, tec->nbmat->axis[1].size),
747 AIR_CAST(unsigned int, tec->nbmat->axis[0].size));
748 return 1;
749 }
750 }
751
752 tec->flag[flagDwiSet] = AIR_FALSE;
753 tec->flag[flagWght] = AIR_FALSE;
754 }
755 return 0;
756 }
757
758 int
tenEstimateUpdate(tenEstimateContext * tec)759 tenEstimateUpdate(tenEstimateContext *tec) {
760 static const char me[]="tenEstimateUpdate";
761 int EE;
762
763 EE = 0;
764 if (!EE) EE |= _tenEstimateCheck(tec);
765 if (!EE) EE |= _tenEstimateNumUpdate(tec);
766 if (!EE) EE |= _tenEstimateAllAllocUpdate(tec);
767 if (!EE) EE |= _tenEstimateDwiAllocUpdate(tec);
768 if (!EE) EE |= _tenEstimateAllSetUpdate(tec);
769 if (!EE) EE |= _tenEstimateDwiSetUpdate(tec);
770 if (!EE) EE |= _tenEstimateWghtUpdate(tec);
771 if (!EE) EE |= _tenEstimateEmatUpdate(tec);
772 if (EE) {
773 biffAddf(TEN, "%s: problem updating", me);
774 return 1;
775 }
776 return 0;
777 }
778
779 /*
780 ** from given tec->all_f or tec->all_d (whichever is non-NULL), sets:
781 ** tec->all[],
782 ** tec->dwi[]
783 ** tec->knownB0, if !tec->estimateB0,
784 ** tec->mdwi,
785 ** tec->conf (from tec->mdwi)
786 */
787 void
_tenEstimateValuesSet(tenEstimateContext * tec)788 _tenEstimateValuesSet(tenEstimateContext *tec) {
789 unsigned int allIdx, dwiIdx, B0Num;
790 double normSum;
791
792 if (!tec->estimateB0) {
793 tec->knownB0 = 0;
794 } else {
795 tec->knownB0 = AIR_NAN;
796 }
797 normSum = 0;
798 tec->mdwi = 0;
799 B0Num = 0;
800 dwiIdx = 0;
801 for (allIdx=0; allIdx<tec->allNum; allIdx++) {
802 if (!tec->skipLut[allIdx]) {
803 tec->all[allIdx] = (tec->all_f
804 ? tec->all_f[allIdx]
805 : tec->all_d[allIdx]);
806 tec->mdwi += tec->bnorm[allIdx]*tec->all[allIdx];
807 normSum += tec->bnorm[allIdx];
808 if (tec->estimateB0 || tec->bnorm[allIdx]) {
809 tec->dwi[dwiIdx++] = tec->all[allIdx];
810 } else {
811 tec->knownB0 += tec->all[allIdx];
812 B0Num += 1;
813 }
814 }
815 }
816 if (!tec->estimateB0) {
817 tec->knownB0 /= B0Num;
818 }
819 tec->mdwi /= normSum;
820 if (tec->dwiConfSoft > 0) {
821 tec->conf = AIR_AFFINE(-1, airErf((tec->mdwi - tec->dwiConfThresh)
822 /tec->dwiConfSoft), 1,
823 0, 1);
824 } else {
825 tec->conf = tec->mdwi > tec->dwiConfThresh;
826 }
827 return;
828 }
829
830 /*
831 ** ASSUMES THAT dwiTmp[] has been stuff with all values simulated from model
832 */
833 double
_tenEstimateErrorDwi(tenEstimateContext * tec)834 _tenEstimateErrorDwi(tenEstimateContext *tec) {
835 unsigned int dwiIdx;
836 double err, diff;
837
838 err = 0;
839 for (dwiIdx=0; dwiIdx<tec->dwiNum; dwiIdx++) {
840 diff = tec->dwi[dwiIdx] - tec->dwiTmp[dwiIdx];
841 /*
842 avg = (tec->dwi[dwiIdx] + tec->dwiTmp[dwiIdx])/2;
843 avg = AIR_ABS(avg);
844 if (avg) {
845 err += diff*diff/(avg*avg);
846 }
847 */
848 err += diff*diff;
849 }
850 err /= tec->dwiNum;
851 return sqrt(err);
852 }
853 double
_tenEstimateErrorLogDwi(tenEstimateContext * tec)854 _tenEstimateErrorLogDwi(tenEstimateContext *tec) {
855 unsigned int dwiIdx;
856 double err, diff;
857
858 err = 0;
859 for (dwiIdx=0; dwiIdx<tec->dwiNum; dwiIdx++) {
860 diff = (log(AIR_MAX(tec->valueMin, tec->dwi[dwiIdx]))
861 - log(AIR_MAX(tec->valueMin, tec->dwiTmp[dwiIdx])));
862 err += diff*diff;
863 }
864 err /= tec->dwiNum;
865 return sqrt(err);
866 }
867
868 /*
869 ** sets:
870 ** tec->dwiTmp[]
871 ** and sets of all of them, regardless of estimateB0
872 */
873 int
_tenEstimate1TensorSimulateSingle(tenEstimateContext * tec,double sigma,double bValue,double B0,const double ten[7])874 _tenEstimate1TensorSimulateSingle(tenEstimateContext *tec,
875 double sigma, double bValue, double B0,
876 const double ten[7]) {
877 static const char me[]="_tenEstimate1TensorSimulateSingle";
878 unsigned int dwiIdx, jj;
879 double nr, ni, vv;
880 const double *bmat;
881
882 if (!( ten && ten )) {
883 biffAddf(TEN, "%s: got NULL pointer", me);
884 return 1;
885 }
886 if (!( AIR_EXISTS(sigma) && sigma >= 0
887 && AIR_EXISTS(bValue) && AIR_EXISTS(B0) )) {
888 biffAddf(TEN, "%s: got bad args: sigma %g, bValue %g, B0 %g\n", me,
889 sigma, bValue, B0);
890 return 1;
891 }
892
893 bmat = AIR_CAST(const double *, tec->nbmat->data);
894 for (dwiIdx=0; dwiIdx<tec->dwiNum; dwiIdx++) {
895 vv = 0;
896 for (jj=0; jj<6; jj++) {
897 vv += bmat[jj]*ten[1+jj];
898 }
899 /*
900 fprintf(stderr, "!%s: sigma = %g, bValue = %g, B0 = %g\n", me,
901 sigma, bValue, B0);
902 fprintf(stderr, "!%s[%u]: bmat=(%g %g %g %g %g %g)."
903 "ten=(%g %g %g %g %g %g)\n",
904 me, dwiIdx,
905 bmat[0], bmat[1], bmat[2], bmat[3], bmat[4], bmat[5],
906 ten[1], ten[2], ten[3], ten[4], ten[5], ten[6]);
907 fprintf(stderr, "!%s: %g * exp(- %g * %g) = %g * exp(%g) = "
908 "%g * %g = ... \n", me,
909 B0, bValue, vv, B0, -bValue*vv, B0, exp(-bValue*vv));
910 */
911 /* need AIR_MAX(0, vv) because B:D might be negative */
912 vv = B0*exp(-bValue*AIR_MAX(0, vv));
913 /*
914 fprintf(stderr, "!%s: vv = %g\n", me, vv);
915 */
916 if (sigma > 0) {
917 airNormalRand(&nr, &ni);
918 nr *= sigma;
919 ni *= sigma;
920 vv = sqrt((vv+nr)*(vv+nr) + ni*ni);
921 }
922 tec->dwiTmp[dwiIdx] = vv;
923 if (!AIR_EXISTS(tec->dwiTmp[dwiIdx])) {
924 fprintf(stderr, "**********************************\n");
925
926 }
927 /*
928 if (tec->verbose) {
929 fprintf(stderr, "%s: dwi[%u] = %g\n", me, dwiIdx, tec->dwiTmp[dwiIdx]);
930 }
931 */
932 bmat += tec->nbmat->axis[0].size;
933 }
934 return 0;
935 }
936
937 int
tenEstimate1TensorSimulateSingle_f(tenEstimateContext * tec,float * simval,float sigma,float bValue,float B0,const float _ten[7])938 tenEstimate1TensorSimulateSingle_f(tenEstimateContext *tec,
939 float *simval,
940 float sigma, float bValue, float B0,
941 const float _ten[7]) {
942 static const char me[]="tenEstimate1TensorSimulateSingle_f";
943 unsigned int allIdx, dwiIdx;
944 double ten[7];
945
946 if (!(tec && simval && _ten)) {
947 biffAddf(TEN, "%s: got NULL pointer", me);
948 return 1;
949 }
950
951 TEN_T_COPY(ten, _ten);
952 if (_tenEstimate1TensorSimulateSingle(tec, sigma, bValue, B0, ten)) {
953 biffAddf(TEN, "%s: ", me);
954 return 1;
955 }
956 dwiIdx = 0;
957 for (allIdx=0; allIdx<tec->allNum; allIdx++) {
958 if (tec->estimateB0 || tec->bnorm[allIdx]) {
959 simval[allIdx] = AIR_CAST(float, tec->dwiTmp[dwiIdx++]);
960 } else {
961 simval[allIdx] = B0;
962 }
963 }
964 return 0;
965 }
966
967 int
tenEstimate1TensorSimulateSingle_d(tenEstimateContext * tec,double * simval,double sigma,double bValue,double B0,const double ten[7])968 tenEstimate1TensorSimulateSingle_d(tenEstimateContext *tec,
969 double *simval,
970 double sigma, double bValue, double B0,
971 const double ten[7]) {
972 static const char me[]="tenEstimate1TensorSimulateSingle_d";
973 unsigned int allIdx, dwiIdx;
974
975 if (!(tec && simval && ten)) {
976 biffAddf(TEN, "%s: got NULL pointer", me);
977 return 1;
978 }
979 if (!( AIR_EXISTS(sigma) && sigma >= 0
980 && AIR_EXISTS(bValue) && AIR_EXISTS(B0) )) {
981 biffAddf(TEN, "%s: got bad bargs sigma %g, bValue %g, B0 %g\n", me,
982 sigma, bValue, B0);
983 return 1;
984 }
985
986 if (_tenEstimate1TensorSimulateSingle(tec, sigma, bValue, B0, ten)) {
987 biffAddf(TEN, "%s: ", me);
988 return 1;
989 }
990 dwiIdx = 0;
991 for (allIdx=0; allIdx<tec->allNum; allIdx++) {
992 if (tec->estimateB0 || tec->bnorm[allIdx]) {
993 simval[allIdx] = tec->dwiTmp[dwiIdx++];
994 } else {
995 simval[allIdx] = B0;
996 }
997 }
998 return 0;
999 }
1000
1001 int
tenEstimate1TensorSimulateVolume(tenEstimateContext * tec,Nrrd * ndwi,double sigma,double bValue,const Nrrd * nB0,const Nrrd * nten,int outType,int keyValueSet)1002 tenEstimate1TensorSimulateVolume(tenEstimateContext *tec,
1003 Nrrd *ndwi,
1004 double sigma, double bValue,
1005 const Nrrd *nB0, const Nrrd *nten,
1006 int outType, int keyValueSet) {
1007 static const char me[]="tenEstimate1TensorSimulateVolume";
1008 size_t sizeTen, sizeX, sizeY, sizeZ, NN, II;
1009 double (*tlup)(const void *, size_t), (*blup)(const void *, size_t),
1010 (*lup)(const void *, size_t), ten_d[7], *dwi_d, B0;
1011 float *dwi_f, ten_f[7];
1012 unsigned int tt, allIdx;
1013 int axmap[4], E;
1014 airArray *mop;
1015 char stmp[3][AIR_STRLEN_SMALL];
1016
1017 if (!(tec && ndwi && nB0 && nten)) {
1018 biffAddf(TEN, "%s: got NULL pointer", me);
1019 return 1;
1020 }
1021 /* this should have been done by update(), but why not */
1022 if (_tenEstimateCheck(tec)) {
1023 biffAddf(TEN, "%s: problem in given context", me);
1024 return 1;
1025 }
1026 if (!(AIR_EXISTS(sigma) && sigma >= 0.0
1027 && AIR_EXISTS(bValue) && bValue >= 0.0)) {
1028 biffAddf(TEN, "%s: got invalid sigma (%g) or bValue (%g)\n", me,
1029 sigma, bValue);
1030 return 1;
1031 }
1032 if (airEnumValCheck(nrrdType, outType)) {
1033 biffAddf(TEN, "%s: requested output type %d not valid", me, outType);
1034 return 1;
1035 }
1036 if (!( nrrdTypeFloat == outType || nrrdTypeDouble == outType )) {
1037 biffAddf(TEN, "%s: requested output type (%s) not %s or %s", me,
1038 airEnumStr(nrrdType, outType),
1039 airEnumStr(nrrdType, nrrdTypeFloat),
1040 airEnumStr(nrrdType, nrrdTypeDouble));
1041 return 1;
1042 }
1043
1044 mop = airMopNew();
1045
1046 sizeTen = nrrdKindSize(nrrdKind3DMaskedSymMatrix);
1047 sizeX = nten->axis[1].size;
1048 sizeY = nten->axis[2].size;
1049 sizeZ = nten->axis[3].size;
1050 if (!(3 == nB0->dim &&
1051 sizeX == nB0->axis[0].size &&
1052 sizeY == nB0->axis[1].size &&
1053 sizeZ == nB0->axis[2].size)) {
1054 biffAddf(TEN, "%s: given B0 (%u-D) volume not 3-D %sx%sx%s", me, nB0->dim,
1055 airSprintSize_t(stmp[0], sizeX),
1056 airSprintSize_t(stmp[1], sizeY),
1057 airSprintSize_t(stmp[2], sizeZ));
1058 return 1;
1059 }
1060 if (nrrdMaybeAlloc_va(ndwi, outType, 4,
1061 AIR_CAST(size_t, tec->allNum), sizeX, sizeY, sizeZ)) {
1062 biffMovef(TEN, NRRD, "%s: couldn't allocate DWI output", me);
1063 airMopError(mop); return 1;
1064 }
1065 NN = sizeX * sizeY * sizeZ;
1066 tlup = nrrdDLookup[nten->type];
1067 blup = nrrdDLookup[nB0->type];
1068 dwi_d = AIR_CAST(double *, ndwi->data);
1069 dwi_f = AIR_CAST(float *, ndwi->data);
1070 E = 0;
1071 for (II=0; !E && II<NN; II++) {
1072 B0 = blup(nB0->data, II);
1073 if (nrrdTypeDouble == outType) {
1074 for (tt=0; tt<7; tt++) {
1075 ten_d[tt] = tlup(nten->data, tt + sizeTen*II);
1076 }
1077 E = tenEstimate1TensorSimulateSingle_d(tec, dwi_d, sigma,
1078 bValue, B0, ten_d);
1079 dwi_d += tec->allNum;
1080 } else {
1081 for (tt=0; tt<7; tt++) {
1082 ten_f[tt] = AIR_CAST(float, tlup(nten->data, tt + sizeTen*II));
1083 }
1084 E = tenEstimate1TensorSimulateSingle_f(tec, dwi_f,
1085 AIR_CAST(float, sigma),
1086 AIR_CAST(float, bValue),
1087 AIR_CAST(float, B0),
1088 ten_f);
1089 dwi_f += tec->allNum;
1090 }
1091 if (E) {
1092 biffAddf(TEN, "%s: failed at sample %s", me,
1093 airSprintSize_t(stmp[0], II));
1094 airMopError(mop); return 1;
1095 }
1096 }
1097
1098 ELL_4V_SET(axmap, -1, 1, 2, 3);
1099 nrrdAxisInfoCopy(ndwi, nten, axmap, NRRD_AXIS_INFO_NONE);
1100 ndwi->axis[0].kind = nrrdKindList;
1101 if (nrrdBasicInfoCopy(ndwi, nten,
1102 NRRD_BASIC_INFO_ALL ^ NRRD_BASIC_INFO_SPACE)) {
1103 biffMovef(TEN, NRRD, "%s:", me);
1104 airMopError(mop); return 1;
1105 }
1106 if (keyValueSet) {
1107 char keystr[AIR_STRLEN_MED], valstr[AIR_STRLEN_MED];
1108
1109 nrrdKeyValueAdd(ndwi, tenDWMRIModalityKey, tenDWMRIModalityVal);
1110 sprintf(valstr, "%g", bValue);
1111 nrrdKeyValueAdd(ndwi, tenDWMRIBValueKey, valstr);
1112 if (tec->_ngrad) {
1113 lup = nrrdDLookup[tec->_ngrad->type];
1114 for (allIdx=0; allIdx<tec->allNum; allIdx++) {
1115 sprintf(keystr, tenDWMRIGradKeyFmt, allIdx);
1116 sprintf(valstr, "%g %g %g",
1117 lup(tec->_ngrad->data, 0 + 3*allIdx),
1118 lup(tec->_ngrad->data, 1 + 3*allIdx),
1119 lup(tec->_ngrad->data, 2 + 3*allIdx));
1120 nrrdKeyValueAdd(ndwi, keystr, valstr);
1121 }
1122 } else {
1123 lup = nrrdDLookup[tec->_nbmat->type];
1124 for (allIdx=0; allIdx<tec->allNum; allIdx++) {
1125 sprintf(keystr, tenDWMRIBmatKeyFmt, allIdx);
1126 sprintf(valstr, "%g %g %g %g %g %g",
1127 lup(tec->_nbmat->data, 0 + 6*allIdx),
1128 lup(tec->_nbmat->data, 1 + 6*allIdx),
1129 lup(tec->_nbmat->data, 2 + 6*allIdx),
1130 lup(tec->_nbmat->data, 3 + 6*allIdx),
1131 lup(tec->_nbmat->data, 4 + 6*allIdx),
1132 lup(tec->_nbmat->data, 5 + 6*allIdx));
1133 nrrdKeyValueAdd(ndwi, keystr, valstr);
1134 }
1135 }
1136 }
1137 airMopOkay(mop);
1138 return 0;
1139 }
1140
1141 /*
1142 ** sets:
1143 ** tec->ten[1..6]
1144 ** tec->B0, if tec->estimateB0
1145 */
1146 int
_tenEstimate1Tensor_LLS(tenEstimateContext * tec)1147 _tenEstimate1Tensor_LLS(tenEstimateContext *tec) {
1148 static const char me[]="_tenEstimate1Tensor_LLS";
1149 double *emat, tmp, logB0;
1150 unsigned int ii, jj;
1151
1152 emat = AIR_CAST(double *, tec->nemat->data);
1153 if (tec->verbose) {
1154 fprintf(stderr, "!%s: estimateB0 = %d\n", me, tec->estimateB0);
1155 }
1156 if (tec->estimateB0) {
1157 for (ii=0; ii<tec->allNum; ii++) {
1158 tmp = AIR_MAX(tec->valueMin, tec->all[ii]);
1159 tec->allTmp[ii] = -log(tmp)/(tec->bValue);
1160 }
1161 for (jj=0; jj<7; jj++) {
1162 tmp = 0;
1163 for (ii=0; ii<tec->allNum; ii++) {
1164 tmp += emat[ii + tec->allNum*jj]*tec->allTmp[ii];
1165 }
1166 if (jj < 6) {
1167 tec->ten[1+jj] = tmp;
1168 if (!AIR_EXISTS(tmp)) {
1169 biffAddf(TEN, "%s: estimated non-existent tensor coef (%u) %g",
1170 me, jj, tmp);
1171 return 1;
1172 }
1173 } else {
1174 /* we're on seventh row, for finding B0 */
1175 tec->estimatedB0 = exp(tec->bValue*tmp);
1176 tec->estimatedB0 = AIR_MIN(FLT_MAX, tec->estimatedB0);
1177 if (!AIR_EXISTS(tec->estimatedB0)) {
1178 biffAddf(TEN, "%s: estimated non-existent B0 %g (b=%g, tmp=%g)",
1179 me, tec->estimatedB0, tec->bValue, tmp);
1180 return 1;
1181 }
1182 }
1183 }
1184 } else {
1185 logB0 = log(AIR_MAX(tec->valueMin, tec->knownB0));
1186 for (ii=0; ii<tec->dwiNum; ii++) {
1187 tmp = AIR_MAX(tec->valueMin, tec->dwi[ii]);
1188 tec->dwiTmp[ii] = (logB0 - log(tmp))/(tec->bValue);
1189 }
1190 for (jj=0; jj<6; jj++) {
1191 tmp = 0;
1192 for (ii=0; ii<tec->dwiNum; ii++) {
1193 tmp += emat[ii + tec->dwiNum*jj]*tec->dwiTmp[ii];
1194 if (tec->verbose > 5) {
1195 fprintf(stderr, "%s: emat[(%u,%u)=%u]*dwi[%u] = %g*%g --> %g\n", me,
1196 ii, jj, ii + tec->dwiNum*jj, ii,
1197 emat[ii + tec->dwiNum*jj], tec->dwiTmp[ii],
1198 tmp);
1199 }
1200 }
1201 tec->ten[1+jj] = tmp;
1202 }
1203 }
1204 return 0;
1205 }
1206
1207 int
_tenEstimate1Tensor_WLS(tenEstimateContext * tec)1208 _tenEstimate1Tensor_WLS(tenEstimateContext *tec) {
1209 static const char me[]="_tenEstimate1Tensor_WLS";
1210 unsigned int dwiIdx, iter;
1211 double *wght, dwi, sum;
1212
1213 if (!tec) {
1214 biffAddf(TEN, "%s: got NULL pointer", me);
1215 return 1;
1216 }
1217
1218 wght = AIR_CAST(double *, tec->nwght->data);
1219 sum = 0;
1220 for (dwiIdx=0; dwiIdx<tec->dwiNum; dwiIdx++) {
1221 dwi = tec->dwi[dwiIdx];
1222 dwi = AIR_MAX(tec->valueMin, dwi);
1223 sum += dwi*dwi;
1224 }
1225 for (dwiIdx=0; dwiIdx<tec->dwiNum; dwiIdx++) {
1226 dwi = tec->dwi[dwiIdx];
1227 dwi = AIR_MAX(tec->valueMin, dwi);
1228 wght[dwiIdx + tec->dwiNum*dwiIdx] = dwi*dwi/sum;
1229 }
1230 if (ell_Nm_wght_pseudo_inv(tec->nemat, tec->nbmat, tec->nwght)) {
1231 biffMovef(TEN, ELL,
1232 "%s(1): trouble wght-pseudo-inverting %ux%u B-matrix", me,
1233 AIR_CAST(unsigned int, tec->nbmat->axis[1].size),
1234 AIR_CAST(unsigned int, tec->nbmat->axis[0].size));
1235 return 1;
1236 }
1237 /*
1238 nrrdSave("wght.txt", tec->nwght, NULL);
1239 nrrdSave("bmat.txt", tec->nbmat, NULL);
1240 nrrdSave("emat.txt", tec->nemat, NULL);
1241 */
1242 if (_tenEstimate1Tensor_LLS(tec)) {
1243 biffAddf(TEN, "%s: initial weighted LLS failed", me);
1244 return 1;
1245 }
1246
1247 for (iter=0; iter<tec->WLSIterNum; iter++) {
1248 /*
1249 fprintf(stderr, "!%s: bValue = %g, B0 = %g, ten = %g %g %g %g %g %g\n",
1250 me,
1251 tec->bValue, (tec->estimateB0 ? tec->estimatedB0 : tec->knownB0),
1252 tec->ten[1], tec->ten[2], tec->ten[3],
1253 tec->ten[4], tec->ten[5], tec->ten[6]);
1254 */
1255 if (_tenEstimate1TensorSimulateSingle(tec, 0.0, tec->bValue,
1256 (tec->estimateB0 ?
1257 tec->estimatedB0
1258 : tec->knownB0), tec->ten)) {
1259 biffAddf(TEN, "%s: iter %u", me, iter);
1260 return 1;
1261 }
1262 for (dwiIdx=0; dwiIdx<tec->dwiNum; dwiIdx++) {
1263 dwi = tec->dwiTmp[dwiIdx];
1264 if (!AIR_EXISTS(dwi)) {
1265 biffAddf(TEN, "%s: bad simulated dwi[%u] == %g (iter %u)",
1266 me, dwiIdx, dwi, iter);
1267 return 1;
1268 }
1269 wght[dwiIdx + tec->dwiNum*dwiIdx] = AIR_MAX(FLT_MIN, dwi*dwi);
1270 }
1271 if (ell_Nm_wght_pseudo_inv(tec->nemat, tec->nbmat, tec->nwght)) {
1272 biffMovef(TEN, ELL, "%s(2): trouble w/ %ux%u B-matrix (iter %u)", me,
1273 AIR_CAST(unsigned int, tec->nbmat->axis[1].size),
1274 AIR_CAST(unsigned int, tec->nbmat->axis[0].size), iter);
1275 return 1;
1276 }
1277 _tenEstimate1Tensor_LLS(tec);
1278 }
1279
1280 return 0;
1281 }
1282
1283 int
_tenEstimate1TensorGradient(tenEstimateContext * tec,double * gradB0P,double gradTen[7],double B0,double ten[7],double epsilon,int (* gradientCB)(tenEstimateContext * tec,double * gradB0P,double gTen[7],double B0,double ten[7]),int (* badnessCB)(tenEstimateContext * tec,double * badP,double B0,double ten[7]))1284 _tenEstimate1TensorGradient(tenEstimateContext *tec,
1285 double *gradB0P, double gradTen[7],
1286 double B0, double ten[7],
1287 double epsilon,
1288 int (*gradientCB)(tenEstimateContext *tec,
1289 double *gradB0P, double gTen[7],
1290 double B0, double ten[7]),
1291 int (*badnessCB)(tenEstimateContext *tec,
1292 double *badP,
1293 double B0, double ten[7])) {
1294 static const char me[]="_tenEstimate1TensorGradper";
1295 double forwTen[7], backTen[7], forwBad, backBad;
1296 unsigned int ti;
1297
1298 if (!( tec && gradB0P && gradTen && badnessCB && ten)) {
1299 biffAddf(TEN, "%s: got NULL pointer", me);
1300 return 1;
1301 }
1302
1303 if (gradientCB) {
1304 if (gradientCB(tec, gradB0P, gradTen, B0, ten)) {
1305 biffAddf(TEN, "%s: problem with grad callback", me);
1306 return 1;
1307 }
1308 } else {
1309 /* we find gradient manually */
1310 gradTen[0] = 0;
1311 for (ti=0; ti<6; ti++) {
1312 TEN_T_COPY(forwTen, ten);
1313 TEN_T_COPY(backTen, ten);
1314 forwTen[ti+1] += epsilon;
1315 backTen[ti+1] -= epsilon;
1316 if (badnessCB(tec, &forwBad, B0, forwTen)
1317 || badnessCB(tec, &backBad, B0, backTen)) {
1318 biffAddf(TEN, "%s: trouble at ti=%u", me, ti);
1319 return 1;
1320 }
1321 gradTen[ti+1] = (forwBad - backBad)/(2*epsilon);
1322 }
1323 }
1324
1325 return 0;
1326 }
1327
1328 int
_tenEstimate1TensorDescent(tenEstimateContext * tec,int (* gradientCB)(tenEstimateContext * tec,double * gradB0,double gradTen[7],double B0,double ten[7]),int (* badnessCB)(tenEstimateContext * tec,double * badP,double B0,double ten[7]))1329 _tenEstimate1TensorDescent(tenEstimateContext *tec,
1330 int (*gradientCB)(tenEstimateContext *tec,
1331 double *gradB0,
1332 double gradTen[7],
1333 double B0,
1334 double ten[7]),
1335 int (*badnessCB)(tenEstimateContext *tec,
1336 double *badP,
1337 double B0,
1338 double ten[7])) {
1339 static const char me[]="_tenEstimate1TensorDescent";
1340 double currB0, lastB0, currTen[7], lastTen[7], gradB0, gradTen[7],
1341 epsilon,
1342 stepSize, badInit, bad, badDelta, stepSizeMin = 0.00000000001, badLast;
1343 unsigned int iter, iterMax = 100000;
1344
1345 /* start with WLS fit since its probably close */
1346 _tenEstimate1Tensor_WLS(tec);
1347 if (tec->verbose) {
1348 fprintf(stderr, "%s: WLS gave %g %g %g %g %g %g\n", me,
1349 tec->ten[1], tec->ten[2], tec->ten[3],
1350 tec->ten[4], tec->ten[5], tec->ten[6]);
1351 }
1352
1353 if (badnessCB(tec, &badInit,
1354 (tec->estimateB0 ? tec->estimatedB0 : tec->knownB0), tec->ten)
1355 || !AIR_EXISTS(badInit)) {
1356 biffAddf(TEN, "%s: problem getting initial bad", me);
1357 return 1;
1358 }
1359 if (tec->verbose) {
1360 fprintf(stderr, "\n%s: ________________________________________\n", me);
1361 fprintf(stderr, "%s: start: badInit = %g ---------------\n", me, badInit);
1362 }
1363
1364 epsilon = 0.0000001;
1365 newepsilon:
1366 if (_tenEstimate1TensorGradient(tec, &gradB0, gradTen,
1367 (tec->estimateB0
1368 ? tec->estimatedB0
1369 : tec->knownB0),
1370 tec->ten, epsilon,
1371 gradientCB, badnessCB)) {
1372 biffAddf(TEN, "%s: problem getting initial gradient", me);
1373 return 1;
1374 }
1375 if (!( AIR_EXISTS(gradB0) || 0 <= TEN_T_NORM(gradTen) )) {
1376 biffAddf(TEN, "%s: got bad gradB0 %g or zero-norm tensor grad",
1377 me, gradB0);
1378 return 1;
1379 }
1380 if (tec->verbose) {
1381 fprintf(stderr, "%s: gradTen (%s) = %g %g %g %g %g %g\n", me,
1382 gradientCB ? "analytic" : "cent-diff",
1383 gradTen[1], gradTen[2], gradTen[3],
1384 gradTen[4], gradTen[5], gradTen[6]);
1385 }
1386
1387 stepSize = 0.1;
1388 do {
1389 stepSize /= 10;
1390 TEN_T_SCALE_ADD2(currTen, 1.0, tec->ten, -stepSize, gradTen);
1391 if (tec->estimateB0) {
1392 currB0 = tec->estimatedB0 + -stepSize*gradB0;
1393 } else {
1394 currB0 = tec->knownB0;
1395 }
1396 if (badnessCB(tec, &bad, currB0, currTen)
1397 || !AIR_EXISTS(bad)) {
1398 biffAddf(TEN, "%s: problem getting badness for stepSize", me);
1399 return 1;
1400 }
1401 if (tec->verbose) {
1402 fprintf(stderr, "%s: ************ stepSize = %g --> bad = %g\n",
1403 me, stepSize, bad);
1404 }
1405 } while (bad > badInit && stepSize > stepSizeMin);
1406
1407 if (stepSize <= stepSizeMin) {
1408 if (epsilon > FLT_MIN) {
1409 epsilon /= 10;
1410 fprintf(stderr, "%s: re-trying initial step w/ eps %g\n", me, epsilon);
1411 goto newepsilon;
1412 } else {
1413 biffAddf(TEN, "%s: never found a usable step size", me);
1414 return 1;
1415 }
1416 } else if (tec->verbose) {
1417 biffAddf(TEN, "%s: using step size %g\n", me, stepSize);
1418 }
1419
1420 iter = 0;
1421 badLast = bad;
1422 do {
1423 iter++;
1424 TEN_T_COPY(lastTen, currTen);
1425 lastB0 = currB0;
1426 if (0 == (iter % 3)) {
1427 if (_tenEstimate1TensorGradient(tec, &gradB0, gradTen,
1428 currB0, currTen, stepSize/5,
1429 gradientCB, badnessCB)
1430 || !AIR_EXISTS(gradB0)) {
1431 biffAddf(TEN, "%s[%u]: problem getting iter grad", me, iter);
1432 return 1;
1433 }
1434 }
1435 TEN_T_SCALE_INCR(currTen, -stepSize, gradTen);
1436 if (tec->estimateB0) {
1437 currB0 -= stepSize*gradB0;
1438 }
1439 if (badnessCB(tec, &bad, currB0, currTen)
1440 || !AIR_EXISTS(bad)) {
1441 biffAddf(TEN, "%s[%u]: problem getting badness during grad", me, iter);
1442 return 1;
1443 }
1444 if (tec->verbose) {
1445 fprintf(stderr, "%s: %u bad = %g\n", me, iter, bad);
1446 }
1447 badDelta = bad - badLast;
1448 badLast = bad;
1449 if (badDelta > 0) {
1450 stepSize /= 10;
1451 if (tec->verbose) {
1452 fprintf(stderr, "%s: badDelta %g > 0 ---> stepSize = %g\n",
1453 me, badDelta, stepSize);
1454 }
1455 badDelta = -1; /* bogus improvement for loop continuation */
1456 TEN_T_COPY(currTen, lastTen);
1457 currB0 = lastB0;
1458 }
1459 } while (iter < iterMax && (iter < 2 || badDelta < -0.00005));
1460 if (iter >= iterMax) {
1461 biffAddf(TEN, "%s: didn't converge after %u iterations", me, iter);
1462 return 1;
1463 }
1464 if (tec->verbose) {
1465 fprintf(stderr, "%s: finished\n", me);
1466 }
1467
1468 ELL_6V_COPY(tec->ten+1, currTen+1);
1469 tec->estimatedB0 = currB0;
1470
1471 return 0;
1472 }
1473
1474 int
_tenEstimate1Tensor_GradientNLS(tenEstimateContext * tec,double * gradB0P,double gradTen[7],double currB0,double currTen[7])1475 _tenEstimate1Tensor_GradientNLS(tenEstimateContext *tec,
1476 double *gradB0P, double gradTen[7],
1477 double currB0, double currTen[7]) {
1478 static const char me[]="_tenEstimate1Tensor_GradientNLS";
1479 double *bmat, dot, tmp, diff, scl;
1480 unsigned int dwiIdx;
1481
1482 if (!(tec && gradB0P && gradTen && currTen)) {
1483 biffAddf(TEN, "%s: got NULL pointer", me);
1484 return 1;
1485 }
1486 *gradB0P = 0;
1487 TEN_T_SET(gradTen, 0, 0, 0, 0, 0, 0, 0);
1488 bmat = AIR_CAST(double *, tec->nbmat->data);
1489 for (dwiIdx=0; dwiIdx<tec->dwiNum; dwiIdx++) {
1490 dot = ELL_6V_DOT(bmat, currTen+1);
1491 tmp = currB0*exp(-(tec->bValue)*dot);
1492 diff = tec->dwi[dwiIdx] - tmp;
1493 scl = 2*diff*tmp*(tec->bValue);
1494 ELL_6V_SCALE_INCR(gradTen+1, scl, bmat);
1495 bmat += tec->nbmat->axis[0].size;
1496 /* HEY: increment *gradB0P */
1497 }
1498 ELL_6V_SCALE_INCR(gradTen+1, 1.0/tec->dwiNum, gradTen+1);
1499 return 0;
1500 }
1501
1502 int
_tenEstimate1Tensor_BadnessNLS(tenEstimateContext * tec,double * retP,double currB0,double currTen[7])1503 _tenEstimate1Tensor_BadnessNLS(tenEstimateContext *tec,
1504 double *retP,
1505 double currB0, double currTen[7]) {
1506 static const char me[]="_tenEstimate1Tensor_BadnessNLS";
1507
1508 if (!(retP && tec)) {
1509 biffAddf(TEN, "%s: got NULL pointer", me);
1510 return 1;
1511 }
1512 if (_tenEstimate1TensorSimulateSingle(tec, 0.0, tec->bValue,
1513 currB0, currTen)) {
1514 biffAddf(TEN, "%s: ", me);
1515 return 1;
1516 }
1517 if (tec->verbose > 2) {
1518 unsigned int di;
1519 fprintf(stderr, "%s: simdwi =", me);
1520 for (di=0; di<tec->dwiNum; di++) {
1521 fprintf(stderr, " %g", tec->dwiTmp[di]);
1522 }
1523 fprintf(stderr, "\n");
1524 }
1525 *retP = _tenEstimateErrorDwi(tec);
1526 if (tec->verbose > 2) {
1527 fprintf(stderr, "!%s: badness(%g, (%g) %g %g %g %g %g %g) = %g\n",
1528 me, currB0, currTen[0],
1529 currTen[1], currTen[2], currTen[3],
1530 currTen[4], currTen[5],
1531 currTen[6], *retP);
1532 }
1533 return 0;
1534 }
1535
1536 int
_tenEstimate1Tensor_NLS(tenEstimateContext * tec)1537 _tenEstimate1Tensor_NLS(tenEstimateContext *tec) {
1538 static const char me[]="_tenEstimate1Tensor_NLS";
1539
1540 if (_tenEstimate1TensorDescent(tec,
1541 NULL
1542 /* _tenEstimate1Tensor_GradientNLS */
1543 ,
1544 _tenEstimate1Tensor_BadnessNLS)) {
1545 biffAddf(TEN, "%s: ", me);
1546 return 1;
1547 }
1548 return 0;
1549 }
1550
1551 int
_tenEstimate1Tensor_GradientMLE(tenEstimateContext * tec,double * gradB0P,double gradTen[7],double currB0,double currTen[7])1552 _tenEstimate1Tensor_GradientMLE(tenEstimateContext *tec,
1553 double *gradB0P, double gradTen[7],
1554 double currB0, double currTen[7]) {
1555 static const char me[]="_tenEstimate1Tensor_GradientMLE";
1556 double *bmat, dot, barg, tmp, scl, dwi, sigma, bval;
1557 unsigned int dwiIdx;
1558
1559 if (!(tec && gradB0P && gradTen && currTen)) {
1560 biffAddf(TEN, "%s: got NULL pointer", me);
1561 return 1;
1562 }
1563 if (tec->verbose) {
1564 fprintf(stderr, "%s grad (currTen = %g %g %g %g %g %g)\n", me,
1565 currTen[1], currTen[2], currTen[3],
1566 currTen[4], currTen[5],
1567 currTen[6]);
1568 }
1569 TEN_T_SET(gradTen, 0, 0, 0, 0, 0, 0, 0);
1570 *gradB0P = 0;
1571 sigma = tec->sigma;
1572 bval = tec->bValue;
1573 bmat = AIR_CAST(double *, tec->nbmat->data);
1574 for (dwiIdx=0; dwiIdx<tec->dwiNum; dwiIdx++) {
1575 dwi = tec->dwi[dwiIdx];
1576 dot = ELL_6V_DOT(bmat, currTen+1);
1577 barg = exp(-bval*dot)*(dwi/sigma)*(currB0/sigma);
1578 tmp = (exp(bval*dot)/sigma)*dwi/airBesselI0(barg);
1579 if (tec->verbose) {
1580 fprintf(stderr, "%s[%u]: dot = %g, barg = %g, tmp = %g\n", me, dwiIdx,
1581 dot, barg, tmp);
1582 }
1583 if (tmp > DBL_MIN) {
1584 tmp = currB0/sigma - tmp*airBesselI1(barg);
1585 } else {
1586 tmp = currB0/sigma;
1587 }
1588 if (tec->verbose) {
1589 fprintf(stderr, " ---- tmp = %g\n", tmp);
1590 }
1591 scl = tmp*exp(-2*bval*dot)*bval*currB0/sigma;
1592 ELL_6V_SCALE_INCR(gradTen+1, scl, bmat);
1593 if (tec->verbose) {
1594 fprintf(stderr, "%s[%u]: bmat = %g %g %g %g %g %g\n",
1595 me, dwiIdx,
1596 bmat[0], bmat[1], bmat[2],
1597 bmat[3], bmat[4],
1598 bmat[5]);
1599 fprintf(stderr, "%s[%u]: scl = %g -> gradTen = %g %g %g %g %g %g\n",
1600 me, dwiIdx, scl,
1601 gradTen[1], gradTen[2], gradTen[3],
1602 gradTen[4], gradTen[5],
1603 gradTen[6]);
1604 }
1605 if (!AIR_EXISTS(scl)) {
1606 TEN_T_SET(gradTen, AIR_NAN,
1607 AIR_NAN, AIR_NAN, AIR_NAN,
1608 AIR_NAN, AIR_NAN, AIR_NAN);
1609 *gradB0P = AIR_NAN;
1610 biffAddf(TEN, "%s: scl = %g, very sorry", me, scl);
1611 return 1;
1612 }
1613 bmat += tec->nbmat->axis[0].size;
1614 /* HEY: increment gradB0 */
1615 }
1616 ELL_6V_SCALE_INCR(gradTen+1, 1.0/tec->dwiNum, gradTen+1);
1617 if (tec->verbose) {
1618 fprintf(stderr, "%s: final gradTen = %g %g %g %g %g %g\n", me,
1619 gradTen[1], gradTen[2], gradTen[3],
1620 gradTen[4], gradTen[5],
1621 gradTen[6]);
1622 }
1623 return 0;
1624 }
1625
1626 int
_tenEstimate1Tensor_BadnessMLE(tenEstimateContext * tec,double * retP,double currB0,double curt[7])1627 _tenEstimate1Tensor_BadnessMLE(tenEstimateContext *tec,
1628 double *retP,
1629 double currB0, double curt[7]) {
1630 static const char me[]="_tenEstimate1Tensor_BadnessMLE";
1631 unsigned int dwiIdx;
1632 double *bmat, sum, rice, logrice=0, mesdwi=0, simdwi=0, dot=0;
1633 int E;
1634
1635 E = 0;
1636 sum = 0;
1637 bmat = AIR_CAST(double *, tec->nbmat->data);
1638 for (dwiIdx=0; !E && dwiIdx<tec->dwiNum; dwiIdx++) {
1639 dot = ELL_6V_DOT(bmat, curt+1);
1640 simdwi = currB0*exp(-(tec->bValue)*dot);
1641 mesdwi = tec->dwi[dwiIdx];
1642 if (!E) E |= _tenRician(&rice, mesdwi, simdwi, tec->sigma);
1643 if (!E) E |= !AIR_EXISTS(rice);
1644 if (!E) logrice = log(rice + DBL_MIN);
1645 if (!E) sum += logrice;
1646 if (!E) E |= !AIR_EXISTS(sum);
1647 if (!E) bmat += tec->nbmat->axis[0].size;
1648 }
1649 if (E) {
1650 biffAddf(TEN, "%s[%u]: dot = (%g %g %g %g %g %g).(%g %g %g %g %g %g) = %g",
1651 me, dwiIdx,
1652 bmat[0], bmat[1], bmat[2], bmat[3], bmat[4], bmat[5],
1653 curt[1], curt[2], curt[3], curt[4], curt[5], curt[6], dot);
1654 biffAddf(TEN, "%s[%u]: simdwi = %g * exp(-%g * %g) = %g * exp(%g) "
1655 "= %g * %g = %g", me, dwiIdx,
1656 currB0, tec->bValue, dot,
1657 currB0, -(tec->bValue)*dot,
1658 currB0, exp(-(tec->bValue)*dot),
1659 currB0*exp(-(tec->bValue)*dot));
1660 biffAddf(TEN, "%s[%u]: mesdwi = %g, simdwi = %g, sigma = %g", me, dwiIdx,
1661 mesdwi, simdwi, tec->sigma);
1662 biffAddf(TEN, "%s[%u]: rice = %g, logrice = %g, sum = %g", me, dwiIdx,
1663 rice, logrice, sum);
1664 *retP = AIR_NAN;
1665 return 1;
1666 }
1667 *retP = -sum/tec->dwiNum;
1668 return 0;
1669 }
1670
1671 int
_tenEstimate1Tensor_MLE(tenEstimateContext * tec)1672 _tenEstimate1Tensor_MLE(tenEstimateContext *tec) {
1673 static const char me[]="_tenEstimate1Tensor_MLE";
1674
1675 if (_tenEstimate1TensorDescent(tec, NULL,
1676 _tenEstimate1Tensor_BadnessMLE)) {
1677 biffAddf(TEN, "%s: ", me);
1678 return 1;
1679 }
1680
1681 return 0;
1682 }
1683
1684 /*
1685 ** sets:
1686 ** tec->ten[0] (from tec->conf)
1687 ** tec->time, if tec->recordTime
1688 ** tec->errorDwi, if tec->recordErrorDwi
1689 ** tec->errorLogDwi, if tec->recordErrorLogDwi
1690 ** tec->likelihoodDwi, if tec->recordLikelihoodDwi
1691 */
1692 int
_tenEstimate1TensorSingle(tenEstimateContext * tec)1693 _tenEstimate1TensorSingle(tenEstimateContext *tec) {
1694 static const char me[]="_tenEstimate1TensorSingle";
1695 double time0, B0;
1696 int E;
1697
1698 _tenEstimateOutputInit(tec);
1699 time0 = tec->recordTime ? airTime() : 0;
1700 _tenEstimateValuesSet(tec);
1701 tec->ten[0] = tec->conf;
1702 switch(tec->estimate1Method) {
1703 case tenEstimate1MethodLLS:
1704 E = _tenEstimate1Tensor_LLS(tec);
1705 break;
1706 case tenEstimate1MethodWLS:
1707 E = _tenEstimate1Tensor_WLS(tec);
1708 break;
1709 case tenEstimate1MethodNLS:
1710 E = _tenEstimate1Tensor_NLS(tec);
1711 break;
1712 case tenEstimate1MethodMLE:
1713 E = _tenEstimate1Tensor_MLE(tec);
1714 break;
1715 default:
1716 biffAddf(TEN, "%s: estimation method %d unimplemented",
1717 me, tec->estimate1Method);
1718 return 1;
1719 }
1720 tec->time = tec->recordTime ? airTime() - time0 : 0;
1721 if (tec->negEvalShift) {
1722 double eval[3];
1723 tenEigensolve_d(eval, NULL, tec->ten);
1724 if (eval[2] < 0) {
1725 tec->ten[1] += -eval[2];
1726 tec->ten[4] += -eval[2];
1727 tec->ten[6] += -eval[2];
1728 }
1729 }
1730 if (E) {
1731 TEN_T_SET(tec->ten, AIR_NAN,
1732 AIR_NAN, AIR_NAN, AIR_NAN,
1733 AIR_NAN, AIR_NAN, AIR_NAN);
1734 if (tec->estimateB0) {
1735 tec->estimatedB0 = AIR_NAN;
1736 }
1737 biffAddf(TEN, "%s: estimation failed", me);
1738 return 1;
1739 }
1740
1741 if (tec->recordErrorDwi
1742 || tec->recordErrorLogDwi) {
1743 B0 = tec->estimateB0 ? tec->estimatedB0 : tec->knownB0;
1744 if (_tenEstimate1TensorSimulateSingle(tec, 0.0, tec->bValue,
1745 B0, tec->ten)) {
1746 biffAddf(TEN, "%s: simulation failed", me);
1747 return 1;
1748 }
1749 if (tec->recordErrorDwi) {
1750 tec->errorDwi = _tenEstimateErrorDwi(tec);
1751 }
1752 if (tec->recordErrorLogDwi) {
1753 tec->errorLogDwi = _tenEstimateErrorLogDwi(tec);
1754 }
1755 }
1756
1757 /* HEY: record likelihood! */
1758
1759 return 0;
1760 }
1761
1762 int
tenEstimate1TensorSingle_f(tenEstimateContext * tec,float ten[7],const float * all)1763 tenEstimate1TensorSingle_f(tenEstimateContext *tec,
1764 float ten[7], const float *all) {
1765 static const char me[]="tenEstimate1TensorSingle_f";
1766
1767 if (!(tec && ten && all)) {
1768 biffAddf(TEN, "%s: got NULL pointer", me);
1769 return 1;
1770 }
1771
1772 tec->all_f = all;
1773 tec->all_d = NULL;
1774 /*
1775 fprintf(stderr, "!%s(%u): B0 = %g,%g\n", me, __LINE__,
1776 tec->knownB0, tec->estimatedB0);
1777 */
1778 if (_tenEstimate1TensorSingle(tec)) {
1779 biffAddf(TEN, "%s: ", me);
1780 return 1;
1781 }
1782 /*
1783 fprintf(stderr, "!%s(%u): B0 = %g,%g\n", me, __LINE__,
1784 tec->knownB0, tec->estimatedB0);
1785 */
1786 TEN_T_COPY_TT(ten, float, tec->ten);
1787
1788 return 0;
1789 }
1790
1791 int
tenEstimate1TensorSingle_d(tenEstimateContext * tec,double ten[7],const double * all)1792 tenEstimate1TensorSingle_d(tenEstimateContext *tec,
1793 double ten[7], const double *all) {
1794 static const char me[]="tenEstimate1TensorSingle_d";
1795 unsigned int ii;
1796
1797 if (!(tec && ten && all)) {
1798 biffAddf(TEN, "%s: got NULL pointer", me);
1799 return 1;
1800 }
1801
1802 tec->all_f = NULL;
1803 tec->all_d = all;
1804 if (tec->verbose) {
1805 for (ii=0; ii<tec->allNum; ii++) {
1806 fprintf(stderr, "%s: dwi[%u] = %g\n", me, ii,
1807 tec->all_d ? tec->all_d[ii] : tec->all_f[ii]);
1808 }
1809 fprintf(stderr, "%s: will estimate by %d (%s) \n"
1810 " estimateB0 %d; valueMin %g\n", me,
1811 tec->estimate1Method,
1812 airEnumStr(tenEstimate1Method, tec->estimate1Method),
1813 tec->estimateB0, tec->valueMin);
1814 }
1815 if (_tenEstimate1TensorSingle(tec)) {
1816 biffAddf(TEN, "%s: ", me);
1817 return 1;
1818 }
1819 if (tec->verbose) {
1820 fprintf(stderr, "%s: ten = %g %g %g %g %g %g %g\n", me,
1821 tec->ten[0],
1822 tec->ten[1], tec->ten[2], tec->ten[3],
1823 tec->ten[4], tec->ten[5],
1824 tec->ten[6]);
1825 }
1826 TEN_T_COPY(ten, tec->ten);
1827 return 0;
1828 }
1829
1830 int
tenEstimate1TensorVolume4D(tenEstimateContext * tec,Nrrd * nten,Nrrd ** nB0P,Nrrd ** nterrP,const Nrrd * ndwi,int outType)1831 tenEstimate1TensorVolume4D(tenEstimateContext *tec,
1832 Nrrd *nten, Nrrd **nB0P, Nrrd **nterrP,
1833 const Nrrd *ndwi, int outType) {
1834 static const char me[]="tenEstimate1TensorVolume4D";
1835 char doneStr[20];
1836 size_t sizeTen, sizeX, sizeY, sizeZ, NN, II, tick;
1837 double *all, ten[7], (*lup)(const void *, size_t),
1838 (*ins)(void *v, size_t I, double d);
1839 unsigned int dd;
1840 airArray *mop;
1841 int axmap[4];
1842 char stmp[AIR_STRLEN_SMALL];
1843
1844 #if 0
1845 #define NUM 800
1846 double val[NUM], minVal=0, maxVal=10, arg;
1847 unsigned int valIdx;
1848 Nrrd *nval;
1849 for (valIdx=0; valIdx<NUM; valIdx++) {
1850 arg = AIR_AFFINE(0, valIdx, NUM-1, minVal, maxVal);
1851 if (_tenRician(val + valIdx, arg, 1, 1)) {
1852 biffAddf(TEN, "%s: you are out of luck", me);
1853 return 1;
1854 }
1855 }
1856 nval = nrrdNew();
1857 nrrdWrap(nval, val, nrrdTypeDouble, 1, AIR_CAST(size_t, NUM));
1858 nrrdSave("nval.nrrd", nval, NULL);
1859 nrrdNix(nval);
1860 #endif
1861
1862 if (!(tec && nten && ndwi)) {
1863 /* nerrP and _NB0P can be NULL */
1864 biffAddf(TEN, "%s: got NULL pointer", me);
1865 return 1;
1866 }
1867 if (nrrdCheck(ndwi)) {
1868 biffMovef(TEN, NRRD, "%s: DWI volume not valid", me);
1869 return 1;
1870 }
1871 if (!( 4 == ndwi->dim && 7 <= ndwi->axis[0].size )) {
1872 biffAddf(TEN, "%s: DWI volume should be 4-D with axis 0 size >= 7", me);
1873 return 1;
1874 }
1875 if (tec->allNum != ndwi->axis[0].size) {
1876 biffAddf(TEN, "%s: from %s info, expected %u values per sample, "
1877 "but have %s in volume", me,
1878 tec->_ngrad ? "gradient" : "B-matrix", tec->allNum,
1879 airSprintSize_t(stmp, ndwi->axis[0].size));
1880 return 1;
1881 }
1882 if (nrrdTypeBlock == ndwi->type) {
1883 biffAddf(TEN, "%s: DWI volume has non-scalar type %s", me,
1884 airEnumStr(nrrdType, ndwi->type));
1885 return 1;
1886 }
1887 if (airEnumValCheck(nrrdType, outType)) {
1888 biffAddf(TEN, "%s: requested output type %d not valid", me, outType);
1889 return 1;
1890 }
1891 if (!( nrrdTypeFloat == outType || nrrdTypeDouble == outType )) {
1892 biffAddf(TEN, "%s: requested output type (%s) not %s or %s", me,
1893 airEnumStr(nrrdType, outType),
1894 airEnumStr(nrrdType, nrrdTypeFloat),
1895 airEnumStr(nrrdType, nrrdTypeDouble));
1896 return 1;
1897 }
1898 if (nterrP) {
1899 int recE, recEL, recLK;
1900 recE = !!(tec->recordErrorDwi);
1901 recEL = !!(tec->recordErrorLogDwi);
1902 recLK = !!(tec->recordLikelihoodDwi);
1903 if (1 != recE + recEL + recLK) {
1904 biffAddf(TEN, "%s: requested error volume but need exactly one of "
1905 "recordErrorDwi, recordErrorLogDwi, recordLikelihoodDwi "
1906 "to be set", me);
1907 return 1;
1908 }
1909 }
1910
1911 mop = airMopNew();
1912
1913 sizeTen = nrrdKindSize(nrrdKind3DMaskedSymMatrix);
1914 sizeX = ndwi->axis[1].size;
1915 sizeY = ndwi->axis[2].size;
1916 sizeZ = ndwi->axis[3].size;
1917 all = AIR_CAST(double *, calloc(tec->allNum, sizeof(double)));
1918 if (!all) {
1919 biffAddf(TEN, "%s: couldn't allocate length %u array", me, tec->allNum);
1920 airMopError(mop); return 1;
1921 }
1922 airMopAdd(mop, all, airFree, airMopAlways);
1923
1924 if (nrrdMaybeAlloc_va(nten, outType, 4,
1925 sizeTen, sizeX, sizeY, sizeZ)) {
1926 biffMovef(TEN, NRRD, "%s: couldn't allocate tensor output", me);
1927 airMopError(mop); return 1;
1928 }
1929 if (nB0P) {
1930 *nB0P = nrrdNew();
1931 if (nrrdMaybeAlloc_va(*nB0P, outType, 3,
1932 sizeX, sizeY, sizeZ)) {
1933 biffMovef(TEN, NRRD, "%s: couldn't allocate B0 output", me);
1934 airMopError(mop); return 1;
1935 }
1936 airMopAdd(mop, *nB0P, (airMopper)nrrdNuke, airMopOnError);
1937 airMopAdd(mop, nB0P, (airMopper)airSetNull, airMopOnError);
1938 }
1939 if (nterrP) {
1940 *nterrP = nrrdNew();
1941 if (nrrdMaybeAlloc_va(*nterrP, outType, 3,
1942 sizeX, sizeY, sizeZ)
1943 || nrrdBasicInfoCopy(*nterrP, ndwi,
1944 NRRD_BASIC_INFO_DATA_BIT
1945 | NRRD_BASIC_INFO_TYPE_BIT
1946 | NRRD_BASIC_INFO_BLOCKSIZE_BIT
1947 | NRRD_BASIC_INFO_DIMENSION_BIT
1948 | NRRD_BASIC_INFO_CONTENT_BIT
1949 | NRRD_BASIC_INFO_MEASUREMENTFRAME_BIT
1950 | NRRD_BASIC_INFO_COMMENTS_BIT
1951 | (nrrdStateKeyValuePairsPropagate
1952 ? 0
1953 : NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT))) {
1954 biffMovef(TEN, NRRD, "%s: couldn't creatting fitting error output", me);
1955 airMopError(mop); return 1;
1956 }
1957 ELL_3V_SET(axmap, 1, 2, 3);
1958 nrrdAxisInfoCopy(*nterrP, ndwi, axmap, NRRD_AXIS_INFO_NONE);
1959 airMopAdd(mop, *nterrP, (airMopper)nrrdNuke, airMopOnError);
1960 airMopAdd(mop, nterrP, (airMopper)airSetNull, airMopOnError);
1961 }
1962 NN = sizeX * sizeY * sizeZ;
1963 lup = nrrdDLookup[ndwi->type];
1964 ins = nrrdDInsert[outType];
1965 if (tec->progress) {
1966 fprintf(stderr, "%s: ", me);
1967 }
1968 fflush(stderr);
1969 tick = NN / 200;
1970 tick = AIR_MAX(1, tick);
1971 for (II=0; II<NN; II++) {
1972 if (tec->progress && 0 == II%tick) {
1973 fprintf(stderr, "%s", airDoneStr(0, II, NN-1, doneStr));
1974 }
1975 for (dd=0; dd<tec->allNum; dd++) {
1976 all[dd] = lup(ndwi->data, dd + tec->allNum*II);
1977 }
1978 /*
1979 tec->verbose = 10*(II == 42509);
1980 */
1981 if (tec->verbose) {
1982 fprintf(stderr, "!%s: hello; II=%u\n", me, AIR_CAST(unsigned int, II));
1983 }
1984 if (tenEstimate1TensorSingle_d(tec, ten, all)) {
1985 biffAddf(TEN, "%s: failed at sample %s", me,
1986 airSprintSize_t(stmp, II));
1987 airMopError(mop); return 1;
1988 }
1989 ins(nten->data, 0 + sizeTen*II, ten[0]);
1990 ins(nten->data, 1 + sizeTen*II, ten[1]);
1991 ins(nten->data, 2 + sizeTen*II, ten[2]);
1992 ins(nten->data, 3 + sizeTen*II, ten[3]);
1993 ins(nten->data, 4 + sizeTen*II, ten[4]);
1994 ins(nten->data, 5 + sizeTen*II, ten[5]);
1995 ins(nten->data, 6 + sizeTen*II, ten[6]);
1996 if (nB0P) {
1997 ins((*nB0P)->data, II, (tec->estimateB0
1998 ? tec->estimatedB0
1999 : tec->knownB0));
2000 }
2001 if (nterrP) {
2002 /* this works because above we checked that only one of the
2003 tec->record* flags is set */
2004 if (tec->recordErrorDwi) {
2005 ins((*nterrP)->data, II, tec->errorDwi);
2006 } else if (tec->recordErrorLogDwi) {
2007 ins((*nterrP)->data, II, tec->errorLogDwi);
2008 } else if (tec->recordLikelihoodDwi) {
2009 ins((*nterrP)->data, II, tec->likelihoodDwi);
2010 }
2011 }
2012 }
2013 if (tec->progress) {
2014 fprintf(stderr, "%s\n", airDoneStr(0, II, NN-1, doneStr));
2015 }
2016
2017 ELL_4V_SET(axmap, -1, 1, 2, 3);
2018 nrrdAxisInfoCopy(nten, ndwi, axmap, NRRD_AXIS_INFO_NONE);
2019 nten->axis[0].kind = nrrdKind3DMaskedSymMatrix;
2020 if (nrrdBasicInfoCopy(nten, ndwi,
2021 NRRD_BASIC_INFO_ALL ^ NRRD_BASIC_INFO_SPACE)) {
2022 biffAddf(NRRD, "%s:", me);
2023 return 1;
2024 }
2025
2026 airMopOkay(mop);
2027 return 0;
2028 }
2029