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 const char *
28 tenModelPrefixStr = "DWMRI_model:";
29 
30 static const tenModel *
str2model(const char * str)31 str2model(const char *str) {
32   const tenModel *ret;
33 
34   if (!strcmp(str, TEN_MODEL_STR_ZERO)) {
35     ret = tenModelZero;
36   } else if (!strcmp(str, TEN_MODEL_STR_B0)) {
37     ret = tenModelB0;
38   } else if (!strcmp(str, TEN_MODEL_STR_BALL)) {
39     ret = tenModelBall;
40   } else if (!strcmp(str, TEN_MODEL_STR_1STICK)) {
41     ret = tenModel1Stick;
42   } else if (!strcmp(str, TEN_MODEL_STR_1VECTOR2D)) {
43     ret = tenModel1Vector2D;
44   } else if (!strcmp(str, TEN_MODEL_STR_1UNIT2D)) {
45     ret = tenModel1Unit2D;
46   } else if (!strcmp(str, TEN_MODEL_STR_2UNIT2D)) {
47     ret = tenModel2Unit2D;
48   } else if (!strcmp(str, TEN_MODEL_STR_BALL1STICKEMD)) {
49     ret = tenModelBall1StickEMD;
50   } else if (!strcmp(str, TEN_MODEL_STR_BALL1STICK)) {
51     ret = tenModelBall1Stick;
52   } else if (!strcmp(str, TEN_MODEL_STR_BALL1CYLINDER)) {
53     ret = tenModelBall1Cylinder;
54   } else if (!strcmp(str, TEN_MODEL_STR_1CYLINDER)) {
55     ret = tenModel1Cylinder;
56   } else if (!strcmp(str, TEN_MODEL_STR_1TENSOR2)) {
57     ret = tenModel1Tensor2;
58   } else {
59     /* we don't currently have a tenModelUnknown */
60     ret = NULL;
61   }
62   return ret;
63 }
64 
65 int
tenModelParse(const tenModel ** model,int * plusB0,int requirePrefix,const char * _str)66 tenModelParse(const tenModel **model, int *plusB0,
67               int requirePrefix, const char *_str) {
68   static const char me[]="tenModelParse";
69   char *str, *modstr, *pre;
70   airArray *mop;
71 
72   if (!( model && plusB0 && _str)) {
73     biffAddf(TEN, "%s: got NULL pointer", me);
74     return 1;
75   }
76   str = airStrdup(_str);
77   if (!str) {
78     biffAddf(TEN, "%s: couldn't strdup", me);
79     return 1;
80   }
81   mop = airMopNew();
82   airMopAdd(mop, str, airFree, airMopAlways);
83   pre = strstr(str, tenModelPrefixStr);
84   if (pre) {
85     str += strlen(tenModelPrefixStr);
86   } else {
87     if (requirePrefix) {
88       biffAddf(TEN, "%s: didn't see prefix \"%s\" in \"%s\"", me,
89                tenModelPrefixStr, _str);
90       airMopError(mop); return 1;
91     }
92   }
93   airToLower(str); /* for sake of "b0" and str2model below */
94 
95   if ((modstr = strchr(str, '+'))) {
96     *modstr = '\0';
97     ++modstr;
98     if (!strcmp(str, "b0")) {
99       *plusB0 = AIR_TRUE;
100     } else {
101       biffAddf(TEN, "%s: string (\"%s\") prior to \"+\" not \"b0\"", me, str);
102       airMopError(mop); return 1;
103     }
104   } else {
105     *plusB0 = AIR_FALSE;
106     modstr = str;
107   }
108   if (!(*model = str2model(modstr))) {
109     biffAddf(TEN, "%s: didn't recognize \"%s\" as model", me, modstr);
110     airMopError(mop); return 1;
111   }
112   airMopOkay(mop);
113   return 0;
114 }
115 
116 int
tenModelFromAxisLearnPossible(const NrrdAxisInfo * axinfo)117 tenModelFromAxisLearnPossible(const NrrdAxisInfo *axinfo) {
118 
119   /* HEY keep in synch with nrrdKind* code below */
120   return (nrrdKind3DSymMatrix == axinfo->kind
121           || nrrdKind3DMaskedSymMatrix == axinfo->kind
122           || airStrlen(axinfo->label));
123 }
124 
125 int
tenModelFromAxisLearn(const tenModel ** modelP,int * plusB0,const NrrdAxisInfo * axinfo)126 tenModelFromAxisLearn(const tenModel **modelP,
127                       int *plusB0,
128                       const NrrdAxisInfo *axinfo) {
129   static const char me[]="tenModelFromAxisLearn";
130 
131   if (!(modelP && plusB0 && axinfo)) {
132     biffAddf(TEN, "%s: got NULL pointer", me);
133     return 1;
134   }
135   *plusB0 = AIR_FALSE;
136   /* first try to learn model from axis "kind" */
137   /* HEY should probably also support 3 vector for stick? */
138   /* HEY keep in synch with nrrdKind* code above */
139   if (nrrdKind3DSymMatrix == axinfo->kind
140       || nrrdKind3DMaskedSymMatrix == axinfo->kind) {
141     *modelP = tenModel1Tensor2;
142   } else if (airStrlen(axinfo->label)) {
143     /* try to parse from label */
144     if (tenModelParse(modelP, plusB0, AIR_TRUE, axinfo->label)) {
145       biffAddf(TEN, "%s: couldn't parse label \"%s\"", me, axinfo->label);
146       *modelP = NULL;
147       return 1;
148     }
149   } else {
150     biffAddf(TEN, "%s: don't have kind or label info to learn model", me);
151     *modelP = NULL;
152     return 1;
153   }
154 
155   return 0;
156 }
157 
158 /*
159 ** If nB0 is given, then those B0 image values will be used.
160 ** In this case, either the parm vector can be short by one (seems to be
161 ** missing B0), or the parm vector includes B0, but these will be ignored
162 ** and over-written with the B0 values from nB0.
163 **
164 ** basic and axis info is derived from _nparm
165 */
166 int
tenModelSimulate(Nrrd * ndwi,int typeOut,tenExperSpec * espec,const tenModel * model,const Nrrd * _nB0,const Nrrd * _nparm,int keyValueSet)167 tenModelSimulate(Nrrd *ndwi, int typeOut,
168                  tenExperSpec *espec,
169                  const tenModel *model,
170                  const Nrrd *_nB0,
171                  const Nrrd *_nparm,
172                  int keyValueSet) {
173   static const char me[]="tenModelSimulate";
174   double *ddwi, *parm, (*ins)(void *v, size_t I, double d);
175   char *dwi;
176   size_t szOut[NRRD_DIM_MAX], II, numSamp;
177   const Nrrd *nB0,    /* B0 as doubles */
178     *ndparm,          /* parm as doubles */
179     *ndpparm,         /* parm as doubles, padded */
180     *nparm;           /* final parm as doubles, padded, w/ correct B0 values */
181   Nrrd *ntmp;         /* non-const pointer for working */
182   airArray *mop;
183   unsigned int gpsze, /* given parm size */
184     ii;
185   int useB0Img, needPad, axmap[NRRD_DIM_MAX];
186 
187   if (!(ndwi && espec && model /* _nB0 can be NULL */ && _nparm)) {
188     biffAddf(TEN, "%s: got NULL pointer", me);
189     return 1;
190   }
191   if (!espec->imgNum) {
192     biffAddf(TEN, "%s: given espec wants 0 images, unset?", me);
193     return 1;
194   }
195 
196   gpsze = _nparm->axis[0].size;
197   if (model->parmNum - 1 == gpsze) {
198     /* got one less than needed parm #, see if we got B0 */
199     if (!_nB0) {
200       biffAddf(TEN, "%s: got %u parms, need %u (for %s), "
201                "but didn't get B0 vol",
202                me, gpsze, model->parmNum, model->name);
203       return 1;
204     }
205     useB0Img = AIR_TRUE;
206     needPad = AIR_TRUE;
207   } else if (model->parmNum != gpsze) {
208     biffAddf(TEN, "%s: mismatch between getting %u parms, "
209              "needing %u (for %s)\n",
210              me, gpsze, model->parmNum, model->name);
211     return 1;
212   } else {
213     /* model->parmNum == gpsze */
214     needPad = AIR_FALSE;
215     useB0Img = !!_nB0;
216   }
217 
218   mop = airMopNew();
219   /* get parms as doubles */
220   if (nrrdTypeDouble == _nparm->type) {
221     ndparm = _nparm;
222   } else {
223     ntmp = nrrdNew();
224     airMopAdd(mop, ntmp, (airMopper)nrrdNuke, airMopAlways);
225     if (nrrdConvert(ntmp, _nparm, nrrdTypeDouble)) {
226       biffMovef(TEN, NRRD, "%s: couldn't convert parm to %s", me,
227                 airEnumStr(nrrdType, nrrdTypeDouble));
228       airMopError(mop); return 1;
229     }
230     ndparm = ntmp;
231   }
232   /* get parms the right length */
233   if (!needPad) {
234     ndpparm = ndparm;
235   } else {
236     ptrdiff_t min[NRRD_DIM_MAX], max[NRRD_DIM_MAX];
237     unsigned int ax;
238 
239     ntmp = nrrdNew();
240     airMopAdd(mop, ntmp, (airMopper)nrrdNuke, airMopAlways);
241     for (ax=0; ax<ndparm->dim; ax++) {
242       min[ax] = (!ax ? -1 : 0);
243       max[ax] = ndparm->axis[ax].size-1;
244     }
245     if (nrrdPad_nva(ntmp, ndparm, min, max, nrrdBoundaryBleed, 0.0)) {
246       biffMovef(TEN, NRRD, "%s: couldn't pad", me);
247       airMopError(mop); return 1;
248     }
249     ndpparm = ntmp;
250   }
251   /* put in B0 values if needed */
252   if (!useB0Img) {
253     nparm = ndpparm;
254   } else {
255     if (nrrdTypeDouble == _nB0->type) {
256       nB0 = _nB0;
257     } else {
258       ntmp = nrrdNew();
259       airMopAdd(mop, ntmp, (airMopper)nrrdNuke, airMopAlways);
260       if (nrrdConvert(ntmp, _nB0, nrrdTypeDouble)) {
261         biffMovef(TEN, NRRD, "%s: couldn't convert B0 to %s", me,
262                   airEnumStr(nrrdType, nrrdTypeDouble));
263         airMopError(mop); return 1;
264       }
265       nB0 = ntmp;
266     }
267     /* HEY: this is mostly likely a waste of memory,
268        but its all complicated by const-correctness */
269     ntmp = nrrdNew();
270     airMopAdd(mop, ntmp, (airMopper)nrrdNuke, airMopAlways);
271     if (nrrdSplice(ntmp, ndpparm, nB0, 0, 0)) {
272       biffMovef(TEN, NRRD, "%s: couldn't splice in B0", me);
273       airMopError(mop); return 1;
274     }
275     nparm = ntmp;
276   }
277 
278   /* allocate output (and set axmap) */
279   for (ii=0; ii<nparm->dim; ii++) {
280     szOut[ii] = (!ii
281                  ? espec->imgNum
282                  : nparm->axis[ii].size);
283     axmap[ii] = (!ii
284                  ? -1
285                  : AIR_CAST(int, ii));
286   }
287   if (nrrdMaybeAlloc_nva(ndwi, typeOut, nparm->dim, szOut)) {
288     biffMovef(TEN, NRRD, "%s: couldn't allocate output", me);
289     airMopError(mop); return 1;
290   }
291   if (!( ddwi = AIR_CALLOC(espec->imgNum, double))) {
292     biffAddf(TEN, "%s: couldn't allocate dwi buffer", me);
293     airMopError(mop); return 1;
294   }
295   airMopAdd(mop, ddwi, airFree, airMopAlways);
296   numSamp = nrrdElementNumber(nparm)/nparm->axis[0].size;
297 
298   /* set output */
299   ins = nrrdDInsert[typeOut];
300   parm = AIR_CAST(double *, nparm->data);
301   dwi = AIR_CAST(char *, ndwi->data);
302   for (II=0; II<numSamp; II++) {
303     model->simulate(ddwi, parm, espec);
304     for (ii=0; ii<espec->imgNum; ii++) {
305       ins(dwi, ii, ddwi[ii]);
306     }
307     parm += model->parmNum;
308     dwi += espec->imgNum*nrrdTypeSize[typeOut];
309   }
310 
311   if (keyValueSet) {
312     if (tenDWMRIKeyValueFromExperSpecSet(ndwi, espec)) {
313       biffAddf(TEN, "%s: trouble", me);
314       airMopError(mop); return 1;
315     }
316   }
317 
318   if (nrrdAxisInfoCopy(ndwi, _nparm, axmap, NRRD_AXIS_INFO_SIZE_BIT)
319       || nrrdBasicInfoCopy(ndwi, _nparm,
320                            NRRD_BASIC_INFO_DATA_BIT
321                            | NRRD_BASIC_INFO_TYPE_BIT
322                            | NRRD_BASIC_INFO_BLOCKSIZE_BIT
323                            | NRRD_BASIC_INFO_DIMENSION_BIT
324                            | NRRD_BASIC_INFO_CONTENT_BIT
325                            | NRRD_BASIC_INFO_COMMENTS_BIT
326                            | (nrrdStateKeyValuePairsPropagate
327                               ? 0
328                               : NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT))) {
329     biffMovef(TEN, NRRD, "%s: couldn't copy axis or basic info", me);
330     airMopError(mop); return 1;
331   }
332   ndwi->axis[0].kind = nrrdKindList;
333 
334   airMopOkay(mop);
335   return 0;
336 }
337 
338 /*
339 ** _tenModelSqeFitSingle
340 **
341 ** callable function (as opposed to tenModel method) for doing
342 ** sqe fitting.  Returns the sqe at the converged fit location
343 ** Requires PARM_NUM length buffers testParm and grad
344 */
345 double
_tenModelSqeFitSingle(const tenModel * model,double * testParm,double * grad,double * parm,double * convFrac,unsigned int * itersTaken,const tenExperSpec * espec,double * dwiBuff,const double * dwiMeas,const double * parmInit,int knownB0,unsigned int minIter,unsigned int maxIter,double convEps,int verbose)346 _tenModelSqeFitSingle(const tenModel *model,
347                       double *testParm, double *grad,
348                       double *parm, double *convFrac, unsigned int *itersTaken,
349                       const tenExperSpec *espec,
350                       double *dwiBuff, const double *dwiMeas,
351                       const double *parmInit, int knownB0,
352                       unsigned int minIter, unsigned int maxIter,
353                       double convEps, int verbose) {
354   static const char me[]="_tenModelSqeFitSingle";
355   unsigned int iter, subIter;
356   double step, bak, opp, val, testval, dist, td;
357   int done;
358   char pstr[AIR_STRLEN_MED];
359 
360   step = 1;
361   model->copy(parm, parmInit);
362   val = model->sqe(parm, espec, dwiBuff, dwiMeas, knownB0);
363   model->sqeGrad(grad, parm, espec, dwiBuff, dwiMeas, knownB0);
364   if (verbose > 1) {
365     model->sprint(pstr, parm);
366     fprintf(stderr, "\n");
367     fprintf(stderr, "%s(%s): minIter = %u, maxIter = %u\n", me, model->name,
368             minIter, maxIter);
369     fprintf(stderr, "%s(%s): starting at %s -> %g (step %g)\n", me,
370             model->name, pstr, val, step);
371   }
372 
373   opp = 1.2;  /* opportunistic step size increase */
374   bak = 0.5;  /* scaling back because of bad step */
375   iter = 0;
376   dist = convEps*8;
377   do {
378     subIter = 0;
379     do {
380       model->step(testParm, -step, grad, parm);
381       testval = model->sqe(testParm, espec, dwiBuff, dwiMeas, knownB0);
382       if (verbose > 1) {
383         model->sprint(pstr, testParm);
384         fprintf(stderr, "%s(%s): (iter %u/%u) tried %s -> %g (step %g)\n",
385                 me, model->name, iter, subIter, pstr, testval, step);
386       }
387       if (testval > val) {
388         step *= bak;
389       }
390       subIter++;
391     } while (testval > val && subIter <= maxIter);
392     if (subIter > maxIter) {
393       /* something went wrong with merely trying to find a downhill step;
394          this has occurred previously when (because of a bug) the
395          per-parameter bounds put the test location inside the bounding
396          box while the initial location was outside => could never converge.
397          Not using biff, so this is one way of trying to signal the problem */
398       model->copy(parm, parmInit);
399       *convFrac = AIR_POS_INF;
400       *itersTaken = maxIter;
401       return AIR_POS_INF;
402     }
403     td = model->dist(testParm, parm);
404     dist = (td + dist)/2;
405     val = testval;
406     model->copy(parm, testParm);
407     model->sqeGrad(grad, parm, espec, dwiBuff, dwiMeas, knownB0);
408     step *= opp;
409     iter++;
410     done = (iter < minIter
411             ? AIR_FALSE
412             : (iter > maxIter) || dist < convEps);
413   } while (!done);
414   *convFrac = dist/convEps;
415   *itersTaken = iter;
416   return val;
417 }
418 
419 int
tenModelSqeFit(Nrrd * nparm,Nrrd ** nsqeP,Nrrd ** nconvP,Nrrd ** niterP,const tenModel * model,const tenExperSpec * espec,const Nrrd * ndwi,int knownB0,int saveB0,int typeOut,unsigned int minIter,unsigned int maxIter,unsigned int starts,double convEps,airRandMTState * _rng,int verbose)420 tenModelSqeFit(Nrrd *nparm,
421                Nrrd **nsqeP, Nrrd **nconvP, Nrrd **niterP,
422                const tenModel *model,
423                const tenExperSpec *espec, const Nrrd *ndwi,
424                int knownB0, int saveB0, int typeOut,
425                unsigned int minIter, unsigned int maxIter,
426                unsigned int starts, double convEps,
427                airRandMTState *_rng, int verbose) {
428   static const char me[]="tenModelSqeFit";
429   char doneStr[13];
430   double *ddwi, *dwibuff, sqe, sqeBest,
431     *dparm, *dparmBest,
432     (*ins)(void *v, size_t I, double d),
433     (*lup)(const void *v, size_t I);
434   airArray *mop;
435   unsigned int saveParmNum, dwiNum, ii, lablen, itersTaken;
436   size_t szOut[NRRD_DIM_MAX], II, numSamp;
437   int axmap[NRRD_DIM_MAX], erraxmap[NRRD_DIM_MAX], fitVerbose;
438   const char *dwi;
439   char *parm;
440   airRandMTState *rng;
441   Nrrd *nsqe, *nconv, *niter;
442 
443   /* nsqeP, nconvP, niterP can be NULL */
444   if (!( nparm && model && espec && ndwi )) {
445     biffAddf(TEN, "%s: got NULL pointer", me);
446     return 1;
447   }
448   if (!( starts > 0 )) {
449     biffAddf(TEN, "%s: need non-zero starts", me);
450     return 1;
451   }
452   if (!( nrrdTypeFloat == typeOut || nrrdTypeDouble == typeOut )) {
453     biffAddf(TEN, "%s: typeOut must be %s or %s, not %s", me,
454              airEnumStr(nrrdType, nrrdTypeFloat),
455              airEnumStr(nrrdType, nrrdTypeDouble),
456              airEnumStr(nrrdType, typeOut));
457     return 1;
458   }
459   dwiNum = ndwi->axis[0].size;
460   if (espec->imgNum != dwiNum) {
461     biffAddf(TEN, "%s: espec expects %u images but dwi has %u on axis 0",
462              me, espec->imgNum, AIR_CAST(unsigned int, dwiNum));
463     return 1;
464   }
465 
466   /* allocate output (and set axmap) */
467   dparm = model->alloc();
468   dparmBest = model->alloc();
469   if (!( dparm && dparmBest )) {
470     biffAddf(TEN, "%s: couldn't allocate parm vecs", me);
471     return 1;
472   }
473   mop = airMopNew();
474   airMopAdd(mop, dparm, airFree, airMopAlways);
475   airMopAdd(mop, dparmBest, airFree, airMopAlways);
476   saveParmNum = saveB0 ? model->parmNum : model->parmNum-1;
477   for (ii=0; ii<ndwi->dim; ii++) {
478     szOut[ii] = (!ii
479                  ? saveParmNum
480                  : ndwi->axis[ii].size);
481     axmap[ii] = (!ii
482                  ? -1
483                  : AIR_CAST(int, ii));
484     if (ii) {
485       erraxmap[ii-1] = AIR_CAST(int, ii);
486     }
487   }
488   if (nrrdMaybeAlloc_nva(nparm, typeOut, ndwi->dim, szOut)) {
489     biffMovef(TEN, NRRD, "%s: couldn't allocate output "
490               "(saveB0 %d, knownB0 %d)", me, saveB0, knownB0);
491     airMopError(mop); return 1;
492   }
493   if (nsqeP) {
494     nsqe = *nsqeP;
495     if (!nsqe) {
496       nsqe = nrrdNew();
497       *nsqeP = nsqe;
498     }
499     if (nrrdMaybeAlloc_nva(nsqe, typeOut, ndwi->dim-1, szOut+1)) {
500       biffMovef(TEN, NRRD, "%s: couldn't allocate error output", me);
501       airMopError(mop); return 1;
502     }
503   } else {
504     nsqe = NULL;
505   }
506   if (nconvP) {
507     nconv = *nconvP;
508     if (!nconv) {
509       nconv = nrrdNew();
510       *nconvP = nconv;
511     }
512     if (nrrdMaybeAlloc_nva(nconv, nrrdTypeDouble, ndwi->dim-1, szOut+1)) {
513       biffMovef(TEN, NRRD, "%s: couldn't allocate conv output", me);
514       airMopError(mop); return 1;
515     }
516   } else {
517     nconv = NULL;
518   }
519   if (niterP) {
520     niter = *niterP;
521     if (!niter) {
522       niter = nrrdNew();
523       *niterP = niter;
524     }
525     if (nrrdMaybeAlloc_nva(niter, nrrdTypeUInt, ndwi->dim-1, szOut+1)) {
526       biffMovef(TEN, NRRD, "%s: couldn't allocate iter output", me);
527       airMopError(mop); return 1;
528     }
529   } else {
530     niter = NULL;
531   }
532   ddwi = AIR_CALLOC(espec->imgNum, double);
533   dwibuff = AIR_CALLOC(espec->imgNum, double);
534   if (!(ddwi && dwibuff)) {
535     biffAddf(TEN, "%s: couldn't allocate dwi buffers", me);
536     airMopError(mop); return 1;
537   }
538   airMopAdd(mop, ddwi, airFree, airMopAlways);
539   airMopAdd(mop, dwibuff, airFree, airMopAlways);
540 
541   /* set output */
542   if (_rng) {
543     rng = _rng;
544   } else {
545     airRandMTStateGlobalInit();
546     rng = airRandMTStateGlobal;
547   }
548   numSamp = nrrdElementNumber(ndwi)/ndwi->axis[0].size;
549   lup = nrrdDLookup[ndwi->type];
550   ins = nrrdDInsert[typeOut];
551   parm = AIR_CAST(char *, nparm->data);
552   dwi = AIR_CAST(char *, ndwi->data);
553   itersTaken = 0;
554   if (verbose) {
555     fprintf(stderr, "%s: fitting ...       ", me);
556     fflush(stderr);
557   }
558   for (II=0; II<numSamp; II++) {
559     double cvf, convFrac=0;
560     unsigned int ss, itak;
561     if (verbose) {
562       fprintf(stderr, "%s", airDoneStr(0, II, numSamp, doneStr));
563       fflush(stderr);
564     }
565     for (ii=0; ii<dwiNum; ii++) {
566       ddwi[ii] = lup(dwi, ii);
567     }
568     sqeBest = DBL_MAX; /* forces at least one improvement */
569     for (ss=0; ss<starts; ss++) {
570       /* can add other debugging conditions here */
571       fitVerbose = verbose;
572       if (knownB0) {
573         dparm[0] = tenExperSpecKnownB0Get(espec, ddwi);
574       }
575       model->rand(dparm, rng, knownB0);
576       sqe = model->sqeFit(dparm, &cvf, &itak,
577                           espec, dwibuff, ddwi,
578                           dparm, knownB0, minIter, maxIter,
579                           convEps, fitVerbose);
580       if (sqe <= sqeBest) {
581         sqeBest = sqe;
582         model->copy(dparmBest, dparm);
583         itersTaken = itak;
584         convFrac = cvf;
585       }
586     }
587     for (ii=0; ii<saveParmNum; ii++) {
588       ins(parm, ii, saveB0 ? dparmBest[ii] : dparmBest[ii+1]);
589     }
590     /* save things about fitting into nrrds */
591     if (nsqeP) {
592       ins(nsqe->data, II, sqeBest);
593     }
594     if (nconvP) {
595       nrrdDInsert[nrrdTypeDouble](nconv->data, II, convFrac);
596     }
597     if (niterP) {
598       nrrdDInsert[nrrdTypeUInt](niter->data, II, itersTaken);
599     }
600     parm += saveParmNum*nrrdTypeSize[typeOut];
601     dwi += espec->imgNum*nrrdTypeSize[ndwi->type];
602   }
603   if (verbose) {
604     fprintf(stderr, "%s\n", airDoneStr(0, II, numSamp, doneStr));
605   }
606 
607   if (nrrdAxisInfoCopy(nparm, ndwi, axmap, NRRD_AXIS_INFO_SIZE_BIT)
608       || nrrdBasicInfoCopy(nparm, ndwi,
609                            NRRD_BASIC_INFO_DATA_BIT
610                            | NRRD_BASIC_INFO_TYPE_BIT
611                            | NRRD_BASIC_INFO_BLOCKSIZE_BIT
612                            | NRRD_BASIC_INFO_DIMENSION_BIT
613                            | NRRD_BASIC_INFO_CONTENT_BIT
614                            | NRRD_BASIC_INFO_COMMENTS_BIT
615                            | (nrrdStateKeyValuePairsPropagate
616                               ? 0
617                               : NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT))) {
618     biffMovef(TEN, NRRD, "%s: couldn't copy axis or basic info", me);
619     airMopError(mop); return 1;
620   }
621   if (nsqeP) {
622     if (nrrdAxisInfoCopy(nsqe, ndwi, erraxmap, NRRD_AXIS_INFO_SIZE_BIT)
623         || nrrdBasicInfoCopy(nsqe, ndwi,
624                              NRRD_BASIC_INFO_DATA_BIT
625                              | NRRD_BASIC_INFO_TYPE_BIT
626                              | NRRD_BASIC_INFO_BLOCKSIZE_BIT
627                              | NRRD_BASIC_INFO_DIMENSION_BIT
628                              | NRRD_BASIC_INFO_CONTENT_BIT
629                              | NRRD_BASIC_INFO_COMMENTS_BIT
630                              | (nrrdStateKeyValuePairsPropagate
631                                 ? 0
632                                 : NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT))) {
633       biffMovef(TEN, NRRD,
634                 "%s: couldn't copy axis or basic info to error out", me);
635       airMopError(mop); return 1;
636     }
637   }
638   if (nconvP) {
639     if (nrrdAxisInfoCopy(nconv, ndwi, erraxmap, NRRD_AXIS_INFO_SIZE_BIT)
640         || nrrdBasicInfoCopy(nconv, ndwi,
641                              NRRD_BASIC_INFO_DATA_BIT
642                              | NRRD_BASIC_INFO_TYPE_BIT
643                              | NRRD_BASIC_INFO_BLOCKSIZE_BIT
644                              | NRRD_BASIC_INFO_DIMENSION_BIT
645                              | NRRD_BASIC_INFO_CONTENT_BIT
646                              | NRRD_BASIC_INFO_COMMENTS_BIT
647                              | (nrrdStateKeyValuePairsPropagate
648                                 ? 0
649                                 : NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT))) {
650       biffMovef(TEN, NRRD,
651                 "%s: couldn't copy axis or basic info to conv out", me);
652       airMopError(mop); return 1;
653     }
654   }
655   if (niterP) {
656     if (nrrdAxisInfoCopy(niter, ndwi, erraxmap, NRRD_AXIS_INFO_SIZE_BIT)
657         || nrrdBasicInfoCopy(niter, ndwi,
658                              NRRD_BASIC_INFO_DATA_BIT
659                              | NRRD_BASIC_INFO_TYPE_BIT
660                              | NRRD_BASIC_INFO_BLOCKSIZE_BIT
661                              | NRRD_BASIC_INFO_DIMENSION_BIT
662                              | NRRD_BASIC_INFO_CONTENT_BIT
663                              | NRRD_BASIC_INFO_COMMENTS_BIT
664                              | (nrrdStateKeyValuePairsPropagate
665                                 ? 0
666                                 : NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT))) {
667       biffMovef(TEN, NRRD,
668                 "%s: couldn't copy axis or basic info to iter out", me);
669       airMopError(mop); return 1;
670     }
671   }
672   lablen = (strlen(tenModelPrefixStr)
673             + (saveB0 ? strlen("B0+") : 0)
674             + strlen(model->name)
675             + 1);
676   nparm->axis[0].label = AIR_CALLOC(lablen, char);
677   sprintf(nparm->axis[0].label, "%s%s%s",
678           tenModelPrefixStr,
679           saveB0 ? "B0+" : "",
680           model->name);
681 
682   airMopOkay(mop);
683   return 0;
684 }
685 
686 int
tenModelNllFit(Nrrd * nparm,Nrrd ** nnllP,const tenModel * model,const tenExperSpec * espec,const Nrrd * ndwi,int rician,double sigma,int knownB0)687 tenModelNllFit(Nrrd *nparm, Nrrd **nnllP,
688                const tenModel *model,
689                const tenExperSpec *espec, const Nrrd *ndwi,
690                int rician, double sigma, int knownB0) {
691 
692   AIR_UNUSED(nparm);
693   AIR_UNUSED(nnllP);
694   AIR_UNUSED(model);
695   AIR_UNUSED(espec);
696   AIR_UNUSED(ndwi);
697   AIR_UNUSED(rician);
698   AIR_UNUSED(sigma);
699   AIR_UNUSED(knownB0);
700 
701   return 0;
702 }
703 
704 /*
705 ** copy the B0 info if we have it
706 ** use the same type on the way out.
707 */
708 int
tenModelConvert(Nrrd * nparmDst,int * convRetP,const tenModel * modelDst,const Nrrd * nparmSrc,const tenModel * _modelSrc)709 tenModelConvert(Nrrd *nparmDst, int *convRetP, const tenModel *modelDst,
710                 const Nrrd *nparmSrc, const tenModel *_modelSrc) {
711   static char me[]="tenModelConvert";
712   const tenModel *modelSrc;
713   double *dpdst, *dpsrc, (*lup)(const void *v, size_t I),
714     (*ins)(void *v, size_t I, double d);
715   size_t szOut[NRRD_DIM_MAX], II, NN, tsize;
716   airArray *mop;
717   int withB0, axmap[NRRD_DIM_MAX], convRet=0;
718   unsigned int parmNumDst, parmNumSrc, ii, lablen;
719   const char *parmSrc;
720   char *parmDst;
721 
722   if (!( nparmDst && modelDst && nparmSrc )) {
723     biffAddf(TEN, "%s: got NULL pointer", me);
724     return 1;
725   }
726   if (!_modelSrc) {
727     /* we have to try to learn the source model from the nrrd */
728     if (tenModelFromAxisLearn(&modelSrc, &withB0, nparmSrc->axis + 0)) {
729       biffAddf(TEN, "%s: couldn't learn model from src nparm", me);
730       return 1;
731     }
732   } else {
733     modelSrc = _modelSrc;
734     if (modelSrc->parmNum == nparmSrc->axis[0].size) {
735       withB0 = AIR_TRUE;
736     } if (modelSrc->parmNum-1 == nparmSrc->axis[0].size) {
737       withB0 = AIR_FALSE;
738     } else {
739       biffAddf(TEN, "%s: axis[0].size %u is not \"%s\" parmnum %u or 1 less",
740                me, AIR_CAST(unsigned int, nparmSrc->axis[0].size),
741                modelSrc->name, modelSrc->parmNum);
742       return 1;
743     }
744   }
745 
746   mop = airMopNew();
747   dpdst = modelDst->alloc();
748   airMopAdd(mop, dpdst, airFree, airMopAlways);
749   dpsrc = modelSrc->alloc();
750   airMopAdd(mop, dpsrc, airFree, airMopAlways);
751   lup = nrrdDLookup[nparmSrc->type];
752   ins = nrrdDInsert[nparmSrc->type];
753   parmNumDst = withB0 ? modelDst->parmNum : modelDst->parmNum-1;
754   parmNumSrc = nparmSrc->axis[0].size;
755   for (ii=0; ii<nparmSrc->dim; ii++) {
756     szOut[ii] = (!ii
757                  ? parmNumDst
758                  : nparmSrc->axis[ii].size);
759     axmap[ii] = (!ii
760                  ? -1
761                  : AIR_CAST(int, ii));
762   }
763   if (nrrdMaybeAlloc_nva(nparmDst, nparmSrc->type, nparmSrc->dim, szOut)) {
764     biffMovef(TEN, NRRD, "%s: couldn't allocate output", me);
765     airMopError(mop); return 1;
766   }
767 
768   NN = nrrdElementNumber(nparmSrc)/nparmSrc->axis[0].size;
769   tsize = nrrdTypeSize[nparmSrc->type];
770   parmSrc = AIR_CAST(char *, nparmSrc->data);
771   parmDst = AIR_CAST(char *, nparmDst->data);
772   if (!withB0) {
773     dpsrc[0] = 0;
774   }
775   for (II=0; II<NN; II++) {
776     for (ii=0; ii<parmNumSrc; ii++) {
777       dpsrc[withB0 ? ii : ii+1] = lup(parmSrc, ii);
778     }
779     convRet = modelDst->convert(dpdst, dpsrc, modelSrc);
780     if (2 == convRet) {  /* HEY should be enum for this value */
781       biffAddf(TEN, "%s: error converting from \"%s\" to \"%s\"", me,
782                modelSrc->name, modelDst->name);
783       airMopError(mop); return 1;
784     }
785     for (ii=0; ii<parmNumDst; ii++) {
786       ins(parmDst, ii, dpdst[withB0 ? ii : ii+1]);
787     }
788     parmSrc += parmNumSrc*tsize;
789     parmDst += parmNumDst*tsize;
790   }
791   if (convRetP) {
792     *convRetP = convRet;
793   }
794 
795   if (nrrdAxisInfoCopy(nparmDst, nparmSrc, axmap, NRRD_AXIS_INFO_SIZE_BIT)
796       || nrrdBasicInfoCopy(nparmDst, nparmSrc,
797                            NRRD_BASIC_INFO_DATA_BIT
798                            | NRRD_BASIC_INFO_TYPE_BIT
799                            | NRRD_BASIC_INFO_BLOCKSIZE_BIT
800                            | NRRD_BASIC_INFO_DIMENSION_BIT
801                            | NRRD_BASIC_INFO_CONTENT_BIT
802                            | NRRD_BASIC_INFO_COMMENTS_BIT
803                            | (nrrdStateKeyValuePairsPropagate
804                               ? 0
805                               : NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT))) {
806     biffMovef(TEN, NRRD, "%s: couldn't copy axis or basic info", me);
807     airMopError(mop); return 1;
808   }
809   /* HEY: COPY AND PASTE! from above. perhaps make helper functions? */
810   lablen = (strlen(tenModelPrefixStr)
811             + (withB0 ? strlen("B0+") : 0)
812             + strlen(modelDst->name)
813             + 1);
814   nparmDst->axis[0].label = AIR_CALLOC(lablen, char);
815   sprintf(nparmDst->axis[0].label, "%s%s%s",
816           tenModelPrefixStr,
817           withB0 ? "B0+" : "",
818           modelDst->name);
819 
820   airMopOkay(mop);
821   return 0;
822 }
823