Bug Summary

File:src/ten/tenModel.c
Location:line 535, column 5
Description:Potential leak of memory pointed to by 'dwibuff'

Annotated Source Code

1/*
2 Teem: Tools to process and visualize scientific data and images .
3 Copyright (C) 2013, 2012, 2011, 2010, 2009 University of Chicago
4 Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann
5 Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah
6
7 This library is free software; you can redistribute it and/or
8 modify it under the terms of the GNU Lesser General Public License
9 (LGPL) as published by the Free Software Foundation; either
10 version 2.1 of the License, or (at your option) any later version.
11 The terms of redistributing and/or modifying this software also
12 include exceptions to the LGPL that facilitate static linking.
13
14 This library is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 Lesser General Public License for more details.
18
19 You should have received a copy of the GNU Lesser General Public License
20 along with this library; if not, write to Free Software Foundation, Inc.,
21 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
22*/
23
24#include "ten.h"
25#include "privateTen.h"
26
27const char *
28tenModelPrefixStr = "DWMRI_model:";
29
30static const tenModel *
31str2model(const char *str) {
32 const tenModel *ret;
33
34 if (!strcmp(str, TEN_MODEL_STR_ZERO"zero")) {
35 ret = tenModelZero;
36 } else if (!strcmp(str, TEN_MODEL_STR_B0"b0")) {
37 ret = tenModelB0;
38 } else if (!strcmp(str, TEN_MODEL_STR_BALL"ball")) {
39 ret = tenModelBall;
40 } else if (!strcmp(str, TEN_MODEL_STR_1STICK"1stick")) {
41 ret = tenModel1Stick;
42 } else if (!strcmp(str, TEN_MODEL_STR_1VECTOR2D"1vector2d")) {
43 ret = tenModel1Vector2D;
44 } else if (!strcmp(str, TEN_MODEL_STR_1UNIT2D"1unit2d")) {
45 ret = tenModel1Unit2D;
46 } else if (!strcmp(str, TEN_MODEL_STR_2UNIT2D"2unit2d")) {
47 ret = tenModel2Unit2D;
48 } else if (!strcmp(str, TEN_MODEL_STR_BALL1STICKEMD"ball1stickemd")) {
49 ret = tenModelBall1StickEMD;
50 } else if (!strcmp(str, TEN_MODEL_STR_BALL1STICK"ball1stick")) {
51 ret = tenModelBall1Stick;
52 } else if (!strcmp(str, TEN_MODEL_STR_BALL1CYLINDER"ball1cylinder")) {
53 ret = tenModelBall1Cylinder;
54 } else if (!strcmp(str, TEN_MODEL_STR_1CYLINDER"1cylinder")) {
55 ret = tenModel1Cylinder;
56 } else if (!strcmp(str, TEN_MODEL_STR_1TENSOR2"1tensor2")) {
57 ret = tenModel1Tensor2;
58 } else {
59 /* we don't currently have a tenModelUnknown */
60 ret = NULL((void*)0);
61 }
62 return ret;
63}
64
65int
66tenModelParse(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(TENtenBiffKey, "%s: got NULL pointer", me);
74 return 1;
75 }
76 str = airStrdup(_str);
77 if (!str) {
78 biffAddf(TENtenBiffKey, "%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(TENtenBiffKey, "%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_TRUE1;
100 } else {
101 biffAddf(TENtenBiffKey, "%s: string (\"%s\") prior to \"+\" not \"b0\"", me, str);
102 airMopError(mop); return 1;
103 }
104 } else {
105 *plusB0 = AIR_FALSE0;
106 modstr = str;
107 }
108 if (!(*model = str2model(modstr))) {
109 biffAddf(TENtenBiffKey, "%s: didn't recognize \"%s\" as model", me, modstr);
110 airMopError(mop); return 1;
111 }
112 airMopOkay(mop);
113 return 0;
114}
115
116int
117tenModelFromAxisLearnPossible(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
125int
126tenModelFromAxisLearn(const tenModel **modelP,
127 int *plusB0,
128 const NrrdAxisInfo *axinfo) {
129 static const char me[]="tenModelFromAxisLearn";
130
131 if (!(modelP && plusB0 && axinfo)) {
132 biffAddf(TENtenBiffKey, "%s: got NULL pointer", me);
133 return 1;
134 }
135 *plusB0 = AIR_FALSE0;
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_TRUE1, axinfo->label)) {
145 biffAddf(TENtenBiffKey, "%s: couldn't parse label \"%s\"", me, axinfo->label);
146 *modelP = NULL((void*)0);
147 return 1;
148 }
149 } else {
150 biffAddf(TENtenBiffKey, "%s: don't have kind or label info to learn model", me);
151 *modelP = NULL((void*)0);
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*/
166int
167tenModelSimulate(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_MAX16], 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_MAX16];
186
187 if (!(ndwi && espec && model /* _nB0 can be NULL */ && _nparm)) {
188 biffAddf(TENtenBiffKey, "%s: got NULL pointer", me);
189 return 1;
190 }
191 if (!espec->imgNum) {
192 biffAddf(TENtenBiffKey, "%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(TENtenBiffKey, "%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_TRUE1;
206 needPad = AIR_TRUE1;
207 } else if (model->parmNum != gpsze) {
208 biffAddf(TENtenBiffKey, "%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_FALSE0;
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(TENtenBiffKey, NRRDnrrdBiffKey, "%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_MAX16], max[NRRD_DIM_MAX16];
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(TENtenBiffKey, NRRDnrrdBiffKey, "%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(TENtenBiffKey, NRRDnrrdBiffKey, "%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(TENtenBiffKey, NRRDnrrdBiffKey, "%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)((int)(ii)));
286 }
287 if (nrrdMaybeAlloc_nva(ndwi, typeOut, nparm->dim, szOut)) {
288 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: couldn't allocate output", me);
289 airMopError(mop); return 1;
290 }
291 if (!( ddwi = AIR_CALLOC(espec->imgNum, double)(double*)(calloc((espec->imgNum), sizeof(double))))) {
292 biffAddf(TENtenBiffKey, "%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)((double *)(nparm->data));
301 dwi = AIR_CAST(char *, ndwi->data)((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(TENtenBiffKey, "%s: trouble", me);
314 airMopError(mop); return 1;
315 }
316 }
317
318 if (nrrdAxisInfoCopy(ndwi, _nparm, axmap, NRRD_AXIS_INFO_SIZE_BIT(1<< 1))
319 || nrrdBasicInfoCopy(ndwi, _nparm,
320 NRRD_BASIC_INFO_DATA_BIT(1<< 1)
321 | NRRD_BASIC_INFO_TYPE_BIT(1<< 2)
322 | NRRD_BASIC_INFO_BLOCKSIZE_BIT(1<< 3)
323 | NRRD_BASIC_INFO_DIMENSION_BIT(1<< 4)
324 | NRRD_BASIC_INFO_CONTENT_BIT(1<< 5)
325 | NRRD_BASIC_INFO_COMMENTS_BIT(1<<14)
326 | (nrrdStateKeyValuePairsPropagate
327 ? 0
328 : NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT(1<<15)))) {
329 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%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*/
345double
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(256+1)];
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__stderrp, "\n");
367 fprintf(stderr__stderrp, "%s(%s): minIter = %u, maxIter = %u\n", me, model->name,
368 minIter, maxIter);
369 fprintf(stderr__stderrp, "%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__stderrp, "%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(airFloatPosInf.f);
400 *itersTaken = maxIter;
401 return AIR_POS_INF(airFloatPosInf.f);
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_FALSE0
412 : (iter > maxIter) || dist < convEps);
413 } while (!done);
414 *convFrac = dist/convEps;
415 *itersTaken = iter;
416 return val;
417}
418
419int
420tenModelSqeFit(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_MAX16], II, numSamp;
437 int axmap[NRRD_DIM_MAX16], erraxmap[NRRD_DIM_MAX16], 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 )) {
1
Taking false branch
→
445 biffAddf(TENtenBiffKey, "%s: got NULL pointer", me);
446 return 1;
447 }
448 if (!( starts > 0 )) {
2
←
Assuming 'starts' is > 0
→
3
←
Taking false branch
→
449 biffAddf(TENtenBiffKey, "%s: need non-zero starts", me);
450 return 1;
451 }
452 if (!( nrrdTypeFloat == typeOut || nrrdTypeDouble == typeOut )) {
4
←
Assuming 'typeOut' is not equal to nrrdTypeFloat
→
5
←
Assuming 'typeOut' is equal to nrrdTypeDouble
→
6
←
Taking false branch
→
453 biffAddf(TENtenBiffKey, "%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) {
7
←
Taking false branch
→
461 biffAddf(TENtenBiffKey, "%s: espec expects %u images but dwi has %u on axis 0",
462 me, espec->imgNum, AIR_CAST(unsigned int, dwiNum)((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 )) {
8
←
Taking false branch
→
470 biffAddf(TENtenBiffKey, "%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;
9
←
Assuming 'saveB0' is 0
→
10
←
'?' condition is false
→
477 for (ii=0; ii<ndwi->dim; ii++) {
11
←
Loop condition is false. Execution continues on line 488
→
478 szOut[ii] = (!ii
479 ? saveParmNum
480 : ndwi->axis[ii].size);
481 axmap[ii] = (!ii
482 ? -1
483 : AIR_CAST(int, ii)((int)(ii)));
484 if (ii) {
485 erraxmap[ii-1] = AIR_CAST(int, ii)((int)(ii));
486 }
487 }
488 if (nrrdMaybeAlloc_nva(nparm, typeOut, ndwi->dim, szOut)) {
12
←
Taking false branch
→
489 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: couldn't allocate output "
490 "(saveB0 %d, knownB0 %d)", me, saveB0, knownB0);
491 airMopError(mop); return 1;
492 }
493 if (nsqeP) {
13
←
Assuming 'nsqeP' is null
→
14
←
Taking false branch
→
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(TENtenBiffKey, NRRDnrrdBiffKey, "%s: couldn't allocate error output", me);
501 airMopError(mop); return 1;
502 }
503 } else {
504 nsqe = NULL((void*)0);
505 }
506 if (nconvP) {
15
←
Assuming 'nconvP' is null
→
16
←
Taking false branch
→
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(TENtenBiffKey, NRRDnrrdBiffKey, "%s: couldn't allocate conv output", me);
514 airMopError(mop); return 1;
515 }
516 } else {
517 nconv = NULL((void*)0);
518 }
519 if (niterP) {
17
←
Assuming 'niterP' is null
→
18
←
Taking false branch
→
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(TENtenBiffKey, NRRDnrrdBiffKey, "%s: couldn't allocate iter output", me);
527 airMopError(mop); return 1;
528 }
529 } else {
530 niter = NULL((void*)0);
531 }
532 ddwi = AIR_CALLOC(espec->imgNum, double)(double*)(calloc((espec->imgNum), sizeof(double)));
533 dwibuff = AIR_CALLOC(espec->imgNum, double)(double*)(calloc((espec->imgNum), sizeof(double)));
19
←
Within the expansion of the macro 'AIR_CALLOC':
→
a
Memory is allocated
534 if (!(ddwi && dwibuff)) {
20
←
Taking true branch
→
535 biffAddf(TENtenBiffKey, "%s: couldn't allocate dwi buffers", me);
21
←
Potential leak of memory pointed to by 'dwibuff'
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)((char *)(nparm->data));
552 dwi = AIR_CAST(char *, ndwi->data)((char *)(ndwi->data));
553 itersTaken = 0;
554 if (verbose) {
555 fprintf(stderr__stderrp, "%s: fitting ... ", me);
556 fflush(stderr__stderrp);
557 }
558 for (II=0; II<numSamp; II++) {
559 double cvf, convFrac=0;
560 unsigned int ss, itak;
561 if (verbose) {
562 fprintf(stderr__stderrp, "%s", airDoneStr(0, II, numSamp, doneStr));
563 fflush(stderr__stderrp);
564 }
565 for (ii=0; ii<dwiNum; ii++) {
566 ddwi[ii] = lup(dwi, ii);
567 }
568 sqeBest = DBL_MAX1.7976931348623157e+308; /* 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__stderrp, "%s\n", airDoneStr(0, II, numSamp, doneStr));
605 }
606
607 if (nrrdAxisInfoCopy(nparm, ndwi, axmap, NRRD_AXIS_INFO_SIZE_BIT(1<< 1))
608 || nrrdBasicInfoCopy(nparm, ndwi,
609 NRRD_BASIC_INFO_DATA_BIT(1<< 1)
610 | NRRD_BASIC_INFO_TYPE_BIT(1<< 2)
611 | NRRD_BASIC_INFO_BLOCKSIZE_BIT(1<< 3)
612 | NRRD_BASIC_INFO_DIMENSION_BIT(1<< 4)
613 | NRRD_BASIC_INFO_CONTENT_BIT(1<< 5)
614 | NRRD_BASIC_INFO_COMMENTS_BIT(1<<14)
615 | (nrrdStateKeyValuePairsPropagate
616 ? 0
617 : NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT(1<<15)))) {
618 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%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(1<< 1))
623 || nrrdBasicInfoCopy(nsqe, ndwi,
624 NRRD_BASIC_INFO_DATA_BIT(1<< 1)
625 | NRRD_BASIC_INFO_TYPE_BIT(1<< 2)
626 | NRRD_BASIC_INFO_BLOCKSIZE_BIT(1<< 3)
627 | NRRD_BASIC_INFO_DIMENSION_BIT(1<< 4)
628 | NRRD_BASIC_INFO_CONTENT_BIT(1<< 5)
629 | NRRD_BASIC_INFO_COMMENTS_BIT(1<<14)
630 | (nrrdStateKeyValuePairsPropagate
631 ? 0
632 : NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT(1<<15)))) {
633 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey,
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(1<< 1))
640 || nrrdBasicInfoCopy(nconv, ndwi,
641 NRRD_BASIC_INFO_DATA_BIT(1<< 1)
642 | NRRD_BASIC_INFO_TYPE_BIT(1<< 2)
643 | NRRD_BASIC_INFO_BLOCKSIZE_BIT(1<< 3)
644 | NRRD_BASIC_INFO_DIMENSION_BIT(1<< 4)
645 | NRRD_BASIC_INFO_CONTENT_BIT(1<< 5)
646 | NRRD_BASIC_INFO_COMMENTS_BIT(1<<14)
647 | (nrrdStateKeyValuePairsPropagate
648 ? 0
649 : NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT(1<<15)))) {
650 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey,
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(1<< 1))
657 || nrrdBasicInfoCopy(niter, ndwi,
658 NRRD_BASIC_INFO_DATA_BIT(1<< 1)
659 | NRRD_BASIC_INFO_TYPE_BIT(1<< 2)
660 | NRRD_BASIC_INFO_BLOCKSIZE_BIT(1<< 3)
661 | NRRD_BASIC_INFO_DIMENSION_BIT(1<< 4)
662 | NRRD_BASIC_INFO_CONTENT_BIT(1<< 5)
663 | NRRD_BASIC_INFO_COMMENTS_BIT(1<<14)
664 | (nrrdStateKeyValuePairsPropagate
665 ? 0
666 : NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT(1<<15)))) {
667 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey,
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)(char*)(calloc((lablen), sizeof(char)));
677 sprintf(nparm->axis[0].label, "%s%s%s",__builtin___sprintf_chk (nparm->axis[0].label, 0, __builtin_object_size
(nparm->axis[0].label, 2 > 1 ? 1 : 0), "%s%s%s", tenModelPrefixStr
, saveB0 ? "B0+" : "", model->name)
678 tenModelPrefixStr,__builtin___sprintf_chk (nparm->axis[0].label, 0, __builtin_object_size
(nparm->axis[0].label, 2 > 1 ? 1 : 0), "%s%s%s", tenModelPrefixStr
, saveB0 ? "B0+" : "", model->name)
679 saveB0 ? "B0+" : "",__builtin___sprintf_chk (nparm->axis[0].label, 0, __builtin_object_size
(nparm->axis[0].label, 2 > 1 ? 1 : 0), "%s%s%s", tenModelPrefixStr
, saveB0 ? "B0+" : "", model->name)
680 model->name)__builtin___sprintf_chk (nparm->axis[0].label, 0, __builtin_object_size
(nparm->axis[0].label, 2 > 1 ? 1 : 0), "%s%s%s", tenModelPrefixStr
, saveB0 ? "B0+" : "", model->name);
681
682 airMopOkay(mop);
683 return 0;
684}
685
686int
687tenModelNllFit(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)(void)(nparm);
693 AIR_UNUSED(nnllP)(void)(nnllP);
694 AIR_UNUSED(model)(void)(model);
695 AIR_UNUSED(espec)(void)(espec);
696 AIR_UNUSED(ndwi)(void)(ndwi);
697 AIR_UNUSED(rician)(void)(rician);
698 AIR_UNUSED(sigma)(void)(sigma);
699 AIR_UNUSED(knownB0)(void)(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*/
708int
709tenModelConvert(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_MAX16], II, NN, tsize;
716 airArray *mop;
717 int withB0, axmap[NRRD_DIM_MAX16], convRet=0;
718 unsigned int parmNumDst, parmNumSrc, ii, lablen;
719 const char *parmSrc;
720 char *parmDst;
721
722 if (!( nparmDst && modelDst && nparmSrc )) {
723 biffAddf(TENtenBiffKey, "%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(TENtenBiffKey, "%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_TRUE1;
736 } if (modelSrc->parmNum-1 == nparmSrc->axis[0].size) {
737 withB0 = AIR_FALSE0;
738 } else {
739 biffAddf(TENtenBiffKey, "%s: axis[0].size %u is not \"%s\" parmnum %u or 1 less",
740 me, AIR_CAST(unsigned int, nparmSrc->axis[0].size)((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)((int)(ii)));
762 }
763 if (nrrdMaybeAlloc_nva(nparmDst, nparmSrc->type, nparmSrc->dim, szOut)) {
764 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%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)((char *)(nparmSrc->data));
771 parmDst = AIR_CAST(char *, nparmDst->data)((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(TENtenBiffKey, "%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(1<< 1))
796 || nrrdBasicInfoCopy(nparmDst, nparmSrc,
797 NRRD_BASIC_INFO_DATA_BIT(1<< 1)
798 | NRRD_BASIC_INFO_TYPE_BIT(1<< 2)
799 | NRRD_BASIC_INFO_BLOCKSIZE_BIT(1<< 3)
800 | NRRD_BASIC_INFO_DIMENSION_BIT(1<< 4)
801 | NRRD_BASIC_INFO_CONTENT_BIT(1<< 5)
802 | NRRD_BASIC_INFO_COMMENTS_BIT(1<<14)
803 | (nrrdStateKeyValuePairsPropagate
804 ? 0
805 : NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT(1<<15)))) {
806 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%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)(char*)(calloc((lablen), sizeof(char)));
815 sprintf(nparmDst->axis[0].label, "%s%s%s",__builtin___sprintf_chk (nparmDst->axis[0].label, 0, __builtin_object_size
(nparmDst->axis[0].label, 2 > 1 ? 1 : 0), "%s%s%s", tenModelPrefixStr
, withB0 ? "B0+" : "", modelDst->name)
816 tenModelPrefixStr,__builtin___sprintf_chk (nparmDst->axis[0].label, 0, __builtin_object_size
(nparmDst->axis[0].label, 2 > 1 ? 1 : 0), "%s%s%s", tenModelPrefixStr
, withB0 ? "B0+" : "", modelDst->name)
817 withB0 ? "B0+" : "",__builtin___sprintf_chk (nparmDst->axis[0].label, 0, __builtin_object_size
(nparmDst->axis[0].label, 2 > 1 ? 1 : 0), "%s%s%s", tenModelPrefixStr
, withB0 ? "B0+" : "", modelDst->name)
818 modelDst->name)__builtin___sprintf_chk (nparmDst->axis[0].label, 0, __builtin_object_size
(nparmDst->axis[0].label, 2 > 1 ? 1 : 0), "%s%s%s", tenModelPrefixStr
, withB0 ? "B0+" : "", modelDst->name);
819
820 airMopOkay(mop);
821 return 0;
822}