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