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