GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: src/ten/tendEstim.c Lines: 33 196 16.8 %
Date: 2017-05-26 Branches: 1 125 0.8 %

Line Branch Exec Source
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
27
#define INFO "Estimate tensors from a set of DW images"
28
static const char *_tend_estimInfoL =
29
  (INFO
30
   ". The tensor coefficient weightings associated with "
31
   "each of the DWIs, the B-matrix, is given either as a separate array, "
32
   "(see \"tend bmat\" usage info for details), or by the key-value pairs "
33
   "in the DWI nrrd header.  A \"confidence\" value is computed with the "
34
   "tensor, based on a soft thresholding of the sum of all the DWIs, "
35
   "according to the threshold and softness parameters. ");
36
37
int
38
tend_estimMain(int argc, const char **argv, const char *me,
39
               hestParm *hparm) {
40
  int pret;
41
2
  hestOpt *hopt = NULL;
42
1
  char *perr, *err;
43
  airArray *mop;
44
45
1
  Nrrd **nin, *nin4d, *nbmat, *nterr, *nB0, *nout;
46
1
  char *outS, *terrS, *bmatS, *eb0S;
47
1
  float soft, scale, sigma;
48
1
  int dwiax, EE, knownB0, oldstuff, estmeth, verbose, fixneg;
49
1
  unsigned int ninLen, axmap[4], wlsi, *skip, skipNum, skipIdx;
50
1
  double valueMin, thresh;
51
52
1
  Nrrd *ngradKVP=NULL, *nbmatKVP=NULL;
53
1
  double bKVP, bval;
54
55
  tenEstimateContext *tec;
56
57
1
  hestOptAdd(&hopt, "old", NULL, airTypeInt, 0, 0, &oldstuff, NULL,
58
             "instead of the new tenEstimateContext code, use "
59
             "the old tenEstimateLinear code");
60
1
  hestOptAdd(&hopt, "sigma", "sigma", airTypeFloat, 1, 1, &sigma, "nan",
61
             "Rician noise parameter");
62
1
  hestOptAdd(&hopt, "v", "verbose", airTypeInt, 1, 1, &verbose, "0",
63
             "verbosity level");
64
1
  hestOptAdd(&hopt, "est", "estimate method", airTypeEnum, 1, 1, &estmeth,
65
             "lls",
66
             "estimation method to use. \"lls\": linear-least squares",
67
1
             NULL, tenEstimate1Method);
68
1
  hestOptAdd(&hopt, "wlsi", "WLS iters", airTypeUInt, 1, 1, &wlsi, "1",
69
             "when using weighted-least-squares (\"-est wls\"), how "
70
             "many iterations to do after the initial weighted fit.");
71
1
  hestOptAdd(&hopt, "fixneg", NULL, airTypeInt, 0, 0, &fixneg, NULL,
72
             "after estimating the tensor, ensure that there are no negative "
73
             "eigenvalues by adding (to all eigenvalues) the amount by which "
74
             "the smallest is negative (corresponding to increasing the "
75
             "non-DWI image value).");
76
1
  hestOptAdd(&hopt, "ee", "filename", airTypeString, 1, 1, &terrS, "",
77
             "Giving a filename here allows you to save out the tensor "
78
             "estimation error: a value which measures how much error there "
79
             "is between the tensor model and the given diffusion weighted "
80
             "measurements for each sample.  By default, no such error "
81
             "calculation is saved.");
82
1
  hestOptAdd(&hopt, "eb", "filename", airTypeString, 1, 1, &eb0S, "",
83
             "In those cases where there is no B=0 reference image given "
84
             "(\"-knownB0 false\"), "
85
             "giving a filename here allows you to save out the B=0 image "
86
             "which is estimated from the data.  By default, this image value "
87
             "is estimated but not saved.");
88
1
  hestOptAdd(&hopt, "t", "thresh", airTypeDouble, 1, 1, &thresh, "nan",
89
             "value at which to threshold the mean DWI value per pixel "
90
             "in order to generate the \"confidence\" mask.  By default, "
91
             "the threshold value is calculated automatically, based on "
92
             "histogram analysis.");
93
1
  hestOptAdd(&hopt, "soft", "soft", airTypeFloat, 1, 1, &soft, "0",
94
             "how fuzzy the confidence boundary should be.  By default, "
95
             "confidence boundary is perfectly sharp");
96
1
  hestOptAdd(&hopt, "scale", "scale", airTypeFloat, 1, 1, &scale, "1",
97
             "After estimating the tensor, scale all of its elements "
98
             "(but not the confidence value) by this amount.  Can help with "
99
             "downstream numerical precision if values are very large "
100
             "or small.");
101
1
  hestOptAdd(&hopt, "mv", "min val", airTypeDouble, 1, 1, &valueMin, "1.0",
102
             "minimum plausible value (especially important for linear "
103
             "least squares estimation)");
104
1
  hestOptAdd(&hopt, "B", "B-list", airTypeString, 1, 1, &bmatS, NULL,
105
             "6-by-N list of B-matrices characterizing "
106
             "the diffusion weighting for each "
107
             "image.  \"tend bmat\" is one source for such a matrix; see "
108
             "its usage info for specifics on how the coefficients of "
109
             "the B-matrix are ordered. "
110
             "An unadorned plain text file is a great way to "
111
             "specify the B-matrix.\n  **OR**\n "
112
             "Can say just \"-B kvp\" to try to learn B matrices from "
113
             "key/value pair information in input images.");
114
1
  hestOptAdd(&hopt, "b", "b", airTypeDouble, 1, 1, &bval, "nan",
115
             "\"b\" diffusion-weighting factor (units of sec/mm^2)");
116
1
  hestOptAdd(&hopt, "knownB0", "bool", airTypeBool, 1, 1, &knownB0, NULL,
117
             "Indicates if the B=0 non-diffusion-weighted reference image "
118
             "is known, or if it has to be estimated along with the tensor "
119
             "elements.\n "
120
             "\b\bo if \"true\": in the given list of diffusion gradients or "
121
             "B-matrices, there are one or more with zero norm, which are "
122
             "simply averaged to find the B=0 reference image value\n "
123
             "\b\bo if \"false\": there may or may not be diffusion-weighted "
124
             "images among the input; the B=0 image value is going to be "
125
             "estimated along with the diffusion model");
126
1
  hestOptAdd(&hopt, "i", "dwi0 dwi1", airTypeOther, 1, -1, &nin, "-",
127
             "all the diffusion-weighted images (DWIs), as separate 3D nrrds, "
128
             "**OR**: One 4D nrrd of all DWIs stacked along axis 0",
129
1
             &ninLen, NULL, nrrdHestNrrd);
130
1
  hestOptAdd(&hopt, "o", "nout", airTypeString, 1, 1, &outS, "-",
131
             "output tensor volume");
132
133
1
  mop = airMopNew();
134
1
  airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways);
135
2
  USAGE(_tend_estimInfoL);
136
  JUSTPARSE();
137
  airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways);
138
139
  nout = nrrdNew();
140
  airMopAdd(mop, nout, (airMopper)nrrdNuke, airMopAlways);
141
  nbmat = nrrdNew();
142
  airMopAdd(mop, nbmat, (airMopper)nrrdNuke, airMopAlways);
143
144
  /* figure out B-matrix */
145
  if (strcmp("kvp", airToLower(bmatS))) {
146
    /* its NOT coming from key/value pairs */
147
    if (!AIR_EXISTS(bval)) {
148
      fprintf(stderr, "%s: need to specify scalar b-value\n", me);
149
      airMopError(mop); return 1;
150
    }
151
    if (nrrdLoad(nbmat, bmatS, NULL)) {
152
      airMopAdd(mop, err=biffGetDone(NRRD), airFree, airMopAlways);
153
      fprintf(stderr, "%s: trouble loading B-matrix:\n%s\n", me, err);
154
      airMopError(mop); return 1;
155
    }
156
    nin4d = nin[0];
157
    skip = NULL;
158
    skipNum = 0;
159
  } else {
160
    /* it IS coming from key/value pairs */
161
    if (1 != ninLen) {
162
      fprintf(stderr, "%s: require a single 4-D DWI volume for "
163
              "key/value pair based calculation of B-matrix\n", me);
164
      airMopError(mop); return 1;
165
    }
166
    if (oldstuff) {
167
      if (knownB0) {
168
        fprintf(stderr, "%s: sorry, key/value-based DWI info not compatible "
169
                "with older implementation of knownB0\n", me);
170
        airMopError(mop); return 1;
171
      }
172
    }
173
    if (tenDWMRIKeyValueParse(&ngradKVP, &nbmatKVP, &bKVP,
174
                              &skip, &skipNum, nin[0])) {
175
      airMopAdd(mop, err=biffGetDone(TEN), airFree, airMopAlways);
176
      fprintf(stderr, "%s: trouble parsing DWI info:\n%s\n", me, err);
177
      airMopError(mop); return 1;
178
    }
179
    if (AIR_EXISTS(bval)) {
180
      fprintf(stderr, "%s: WARNING: key/value pair derived b-value %g "
181
              "over-riding %g from command-line", me, bKVP, bval);
182
    }
183
    bval = bKVP;
184
    if (ngradKVP) {
185
      airMopAdd(mop, ngradKVP, (airMopper)nrrdNuke, airMopAlways);
186
      if (tenBMatrixCalc(nbmat, ngradKVP)) {
187
        airMopAdd(mop, err=biffGetDone(TEN), airFree, airMopAlways);
188
        fprintf(stderr, "%s: trouble finding B-matrix:\n%s\n", me, err);
189
        airMopError(mop); return 1;
190
      }
191
    } else {
192
      airMopAdd(mop, nbmatKVP, (airMopper)nrrdNuke, airMopAlways);
193
      if (nrrdConvert(nbmat, nbmatKVP, nrrdTypeDouble)) {
194
        airMopAdd(mop, err=biffGetDone(NRRD), airFree, airMopAlways);
195
        fprintf(stderr, "%s: trouble converting B-matrix:\n%s\n", me, err);
196
        airMopError(mop); return 1;
197
      }
198
    }
199
    /* this will work because of the impositions of tenDWMRIKeyValueParse */
200
    dwiax = ((nrrdKindList == nin[0]->axis[0].kind ||
201
              nrrdKindVector == nin[0]->axis[0].kind)
202
             ? 0
203
             : ((nrrdKindList == nin[0]->axis[1].kind ||
204
                 nrrdKindVector == nin[0]->axis[1].kind)
205
                ? 1
206
                : ((nrrdKindList == nin[0]->axis[2].kind ||
207
                    nrrdKindVector == nin[0]->axis[2].kind)
208
                   ? 2
209
                   : 3)));
210
    if (0 == dwiax) {
211
      nin4d = nin[0];
212
    } else {
213
      axmap[0] = dwiax;
214
      axmap[1] = 1 > dwiax ? 1 : 0;
215
      axmap[2] = 2 > dwiax ? 2 : 1;
216
      axmap[3] = 3 > dwiax ? 3 : 2;
217
      nin4d = nrrdNew();
218
      airMopAdd(mop, nin4d, (airMopper)nrrdNuke, airMopAlways);
219
      if (nrrdAxesPermute(nin4d, nin[0], axmap)) {
220
        airMopAdd(mop, err=biffGetDone(NRRD), airFree, airMopAlways);
221
        fprintf(stderr, "%s: trouble creating DWI volume:\n%s\n", me, err);
222
        airMopError(mop); return 1;
223
      }
224
    }
225
  }
226
227
  nterr = NULL;
228
  nB0 = NULL;
229
  if (!oldstuff) {
230
    if (1 != ninLen) {
231
      fprintf(stderr, "%s: sorry, currently need single 4D volume "
232
              "for new implementation\n", me);
233
      airMopError(mop); return 1;
234
    }
235
    if (!AIR_EXISTS(thresh)) {
236
      unsigned char *isB0 = NULL;
237
      double bten[7], bnorm, *bmat;
238
      unsigned int sl;
239
      /* from nbmat, create an array that indicates B0 images */
240
      if (tenBMatrixCheck(nbmat, nrrdTypeDouble, 6)) {
241
        biffAddf(TEN, "%s: problem within given b-matrix", me);
242
        airMopError(mop); return 1;
243
      }
244
      isB0 = AIR_CAST(unsigned char *, malloc(nbmat->axis[1].size));
245
      airMopAdd(mop, isB0, airFree, airMopAlways);
246
      bmat = (double*) nbmat->data;
247
      for (sl=0; sl<nbmat->axis[1].size; sl++) {
248
        TEN_T_SET(bten, 1.0,
249
                  bmat[0], bmat[1], bmat[2],
250
                  bmat[3], bmat[4],
251
                  bmat[5]);
252
        bnorm = TEN_T_NORM(bten);
253
        isB0[sl]=(bnorm==0.0);
254
        bmat+=6;
255
      }
256
      if (tenEstimateThresholdFind(&thresh, isB0, nin4d)) {
257
        airMopAdd(mop, err=biffGetDone(TEN), airFree, airMopAlways);
258
        fprintf(stderr, "%s: trouble finding threshold:\n%s\n", me, err);
259
        airMopError(mop); return 1;
260
      }
261
      /* HACK to lower threshold a titch */
262
      thresh *= 0.93;
263
      fprintf(stderr, "%s: using mean DWI threshold %g\n", me, thresh);
264
    }
265
    tec = tenEstimateContextNew();
266
    tec->progress = AIR_TRUE;
267
    airMopAdd(mop, tec, (airMopper)tenEstimateContextNix, airMopAlways);
268
    EE = 0;
269
    if (!EE) tenEstimateVerboseSet(tec, verbose);
270
    if (!EE) tenEstimateNegEvalShiftSet(tec, fixneg);
271
    if (!EE) EE |= tenEstimateMethodSet(tec, estmeth);
272
    if (!EE) EE |= tenEstimateBMatricesSet(tec, nbmat, bval, !knownB0);
273
    if (!EE) EE |= tenEstimateValueMinSet(tec, valueMin);
274
    for (skipIdx=0; skipIdx<skipNum; skipIdx++) {
275
      /* fprintf(stderr, "%s: skipping %u\n", me, skip[skipIdx]); */
276
      if (!EE) EE |= tenEstimateSkipSet(tec, skip[skipIdx], AIR_TRUE);
277
    }
278
    switch(estmeth) {
279
    case tenEstimate1MethodLLS:
280
      if (airStrlen(terrS)) {
281
        tec->recordErrorLogDwi = AIR_TRUE;
282
        /* tec->recordErrorDwi = AIR_TRUE; */
283
      }
284
      break;
285
    case tenEstimate1MethodNLS:
286
      if (airStrlen(terrS)) {
287
        tec->recordErrorDwi = AIR_TRUE;
288
      }
289
      break;
290
    case tenEstimate1MethodWLS:
291
      if (!EE) tec->WLSIterNum = wlsi;
292
      if (airStrlen(terrS)) {
293
        tec->recordErrorDwi = AIR_TRUE;
294
      }
295
      break;
296
    case tenEstimate1MethodMLE:
297
      if (!(AIR_EXISTS(sigma) && sigma > 0.0)) {
298
        fprintf(stderr, "%s: can't do %s w/out sigma > 0 (not %g)\n",
299
                me, airEnumStr(tenEstimate1Method, tenEstimate1MethodMLE),
300
                sigma);
301
        airMopError(mop); return 1;
302
      }
303
      if (!EE) EE |= tenEstimateSigmaSet(tec, sigma);
304
      if (airStrlen(terrS)) {
305
        tec->recordLikelihoodDwi = AIR_TRUE;
306
      }
307
      break;
308
    }
309
    if (!EE) EE |= tenEstimateThresholdSet(tec, thresh, soft);
310
    if (!EE) EE |= tenEstimateUpdate(tec);
311
    if (EE) {
312
      airMopAdd(mop, err=biffGetDone(TEN), airFree, airMopAlways);
313
      fprintf(stderr, "%s: trouble setting up estimation:\n%s\n", me, err);
314
      airMopError(mop); return 1;
315
    }
316
    if (tenEstimate1TensorVolume4D(tec, nout, &nB0,
317
                                   airStrlen(terrS)
318
                                   ? &nterr
319
                                   : NULL,
320
                                   nin4d, nrrdTypeFloat)) {
321
      airMopAdd(mop, err=biffGetDone(TEN), airFree, airMopAlways);
322
      fprintf(stderr, "%s: trouble doing estimation:\n%s\n", me, err);
323
      airMopError(mop); return 1;
324
    }
325
    if (airStrlen(terrS)) {
326
      airMopAdd(mop, nterr, (airMopper)nrrdNuke, airMopAlways);
327
    }
328
  } else {
329
    EE = 0;
330
    if (1 == ninLen) {
331
      EE = tenEstimateLinear4D(nout, airStrlen(terrS) ? &nterr : NULL, &nB0,
332
                               nin4d, nbmat, knownB0, thresh, soft, bval);
333
    } else {
334
      EE = tenEstimateLinear3D(nout, airStrlen(terrS) ? &nterr : NULL, &nB0,
335
                               (const Nrrd*const*)nin, ninLen, nbmat,
336
                               knownB0, thresh, soft, bval);
337
    }
338
    if (EE) {
339
      airMopAdd(mop, err=biffGetDone(TEN), airFree, airMopAlways);
340
      fprintf(stderr, "%s: trouble making tensor volume:\n%s\n", me, err);
341
      airMopError(mop); return 1;
342
    }
343
  }
344
  if (nterr) {
345
    /* it was allocated by tenEstimate*, we have to clean it up */
346
    airMopAdd(mop, nterr, (airMopper)nrrdNuke, airMopAlways);
347
  }
348
  if (nB0) {
349
    /* it was allocated by tenEstimate*, we have to clean it up */
350
    airMopAdd(mop, nB0, (airMopper)nrrdNuke, airMopAlways);
351
  }
352
  if (1 != scale) {
353
    if (tenSizeScale(nout, nout, scale)) {
354
      airMopAdd(mop, err=biffGetDone(TEN), airFree, airMopAlways);
355
      fprintf(stderr, "%s: trouble doing scaling:\n%s\n", me, err);
356
      airMopError(mop); return 1;
357
    }
358
  }
359
  if (nterr) {
360
    if (nrrdSave(terrS, nterr, NULL)) {
361
      airMopAdd(mop, err=biffGetDone(NRRD), airFree, airMopAlways);
362
      fprintf(stderr, "%s: trouble writing error image:\n%s\n", me, err);
363
      airMopError(mop); return 1;
364
    }
365
  }
366
  if (!knownB0 && airStrlen(eb0S)) {
367
    if (nrrdSave(eb0S, nB0, NULL)) {
368
      airMopAdd(mop, err=biffGetDone(NRRD), airFree, airMopAlways);
369
      fprintf(stderr, "%s: trouble writing estimated B=0 image:\n%s\n",
370
              me, err);
371
      airMopError(mop); return 1;
372
    }
373
  }
374
375
  if (nrrdSave(outS, nout, NULL)) {
376
    airMopAdd(mop, err=biffGetDone(NRRD), airFree, airMopAlways);
377
    fprintf(stderr, "%s: trouble writing:\n%s\n", me, err);
378
    airMopError(mop); return 1;
379
  }
380
381
  airMopOkay(mop);
382
  return 0;
383
1
}
384
TEND_CMD(estim, INFO);