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 "Simulate DW images from an image of models" |
28 |
|
|
static const char *_tend_msimInfoL = |
29 |
|
|
(INFO |
30 |
|
|
". The output will be in the same form as the input to \"tend estim\". " |
31 |
|
|
"The B-matrices (\"-B\") can be the output from \"tend bmat\", or the " |
32 |
|
|
"gradients can be given directly (\"-g\"); one of these is required. " |
33 |
|
|
"Note that the input tensor image (\"-i\") is the basis of the output " |
34 |
|
|
"per-axis fields and image orientation. NOTE: this includes the " |
35 |
|
|
"measurement frame used in the input tensor image, which implies that " |
36 |
|
|
"the given gradients or B-matrices are already expressed in that " |
37 |
|
|
"measurement frame. "); |
38 |
|
|
|
39 |
|
|
int |
40 |
|
|
tend_msimMain(int argc, const char **argv, const char *me, |
41 |
|
|
hestParm *hparm) { |
42 |
|
|
int pret; |
43 |
|
2 |
hestOpt *hopt = NULL; |
44 |
|
1 |
char *perr, *err; |
45 |
|
|
airArray *mop; |
46 |
|
|
|
47 |
|
|
tenExperSpec *espec; |
48 |
|
1 |
const tenModel *model; |
49 |
|
1 |
int E, seed, keyValueSet, outType, plusB0, insertB0; |
50 |
|
1 |
Nrrd *nin, *nT2, *_ngrad, *ngrad, *nout; |
51 |
|
1 |
char *outS, *modS; |
52 |
|
1 |
double bval, sigma; |
53 |
|
|
|
54 |
|
|
/* maybe this can go in tend.c, but for some reason its explicitly |
55 |
|
|
set to AIR_FALSE there */ |
56 |
|
1 |
hparm->elideSingleOtherDefault = AIR_TRUE; |
57 |
|
|
|
58 |
|
1 |
hestOptAdd(&hopt, "sigma", "sigma", airTypeDouble, 1, 1, &sigma, "0.0", |
59 |
|
|
"Gaussian/Rician noise parameter"); |
60 |
|
1 |
hestOptAdd(&hopt, "seed", "seed", airTypeInt, 1, 1, &seed, "42", |
61 |
|
|
"seed value for RNG which creates noise"); |
62 |
|
1 |
hestOptAdd(&hopt, "g", "grad list", airTypeOther, 1, 1, &_ngrad, NULL, |
63 |
|
|
"gradient list, one row per diffusion-weighted image", |
64 |
|
1 |
NULL, NULL, nrrdHestNrrd); |
65 |
|
1 |
hestOptAdd(&hopt, "b0", "b0 image", airTypeOther, 1, 1, &nT2, "", |
66 |
|
|
"reference non-diffusion-weighted (\"B0\") image, which " |
67 |
|
|
"may be needed if it isn't part of give model param image", |
68 |
|
1 |
NULL, NULL, nrrdHestNrrd); |
69 |
|
1 |
hestOptAdd(&hopt, "i", "model image", airTypeOther, 1, 1, &nin, "-", |
70 |
|
1 |
"input model image", NULL, NULL, nrrdHestNrrd); |
71 |
|
1 |
hestOptAdd(&hopt, "m", "model", airTypeString, 1, 1, &modS, NULL, |
72 |
|
|
"model with which to simulate DWIs, which must be specified if " |
73 |
|
|
"it is not indicated by the first axis in input model image."); |
74 |
|
1 |
hestOptAdd(&hopt, "ib0", "bool", airTypeBool, 1, 1, &insertB0, "false", |
75 |
|
|
"insert a non-DW B0 image at the beginning of the experiment " |
76 |
|
|
"specification (useful if the given gradient list doesn't " |
77 |
|
|
"already have one) and hence also insert a B0 image at the " |
78 |
|
|
"beginning of the output simulated DWIs"); |
79 |
|
1 |
hestOptAdd(&hopt, "b", "b", airTypeDouble, 1, 1, &bval, "1000", |
80 |
|
|
"b value for simulated scan"); |
81 |
|
1 |
hestOptAdd(&hopt, "kvp", "bool", airTypeBool, 1, 1, &keyValueSet, "true", |
82 |
|
|
"generate key/value pairs in the NRRD header corresponding " |
83 |
|
|
"to the input b-value and gradients."); |
84 |
|
1 |
hestOptAdd(&hopt, "t", "type", airTypeEnum, 1, 1, &outType, "float", |
85 |
|
1 |
"output type of DWIs", NULL, nrrdType); |
86 |
|
1 |
hestOptAdd(&hopt, "o", "nout", airTypeString, 1, 1, &outS, "-", |
87 |
|
|
"output dwis"); |
88 |
|
|
|
89 |
|
1 |
mop = airMopNew(); |
90 |
|
1 |
airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways); |
91 |
✓✗ |
2 |
USAGE(_tend_msimInfoL); |
92 |
|
|
PARSE(); |
93 |
|
|
airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways); |
94 |
|
|
|
95 |
|
|
nout = nrrdNew(); |
96 |
|
|
airMopAdd(mop, nout, (airMopper)nrrdNuke, airMopAlways); |
97 |
|
|
espec = tenExperSpecNew(); |
98 |
|
|
airMopAdd(mop, espec, (airMopper)tenExperSpecNix, airMopAlways); |
99 |
|
|
|
100 |
|
|
airSrandMT(seed); |
101 |
|
|
if (nrrdTypeDouble == _ngrad->type) { |
102 |
|
|
ngrad = _ngrad; |
103 |
|
|
} else { |
104 |
|
|
ngrad = nrrdNew(); |
105 |
|
|
airMopAdd(mop, ngrad, (airMopper)nrrdNuke, airMopAlways); |
106 |
|
|
if (nrrdConvert(ngrad, _ngrad, nrrdTypeDouble)) { |
107 |
|
|
airMopAdd(mop, err=biffGetDone(NRRD), airFree, airMopAlways); |
108 |
|
|
fprintf(stderr, "%s: trouble converting grads to %s:\n%s\n", me, |
109 |
|
|
airEnumStr(nrrdType, nrrdTypeDouble), err); |
110 |
|
|
airMopError(mop); return 1; |
111 |
|
|
} |
112 |
|
|
} |
113 |
|
|
plusB0 = AIR_FALSE; |
114 |
|
|
if (airStrlen(modS)) { |
115 |
|
|
if (tenModelParse(&model, &plusB0, AIR_FALSE, modS)) { |
116 |
|
|
airMopAdd(mop, err=biffGetDone(TEN), airFree, airMopAlways); |
117 |
|
|
fprintf(stderr, "%s: trouble parsing model \"%s\":\n%s\n", |
118 |
|
|
me, modS, err); |
119 |
|
|
airMopError(mop); return 1; |
120 |
|
|
} |
121 |
|
|
} else if (tenModelFromAxisLearnPossible(nin->axis + 0)) { |
122 |
|
|
if (tenModelFromAxisLearn(&model, &plusB0, nin->axis + 0)) { |
123 |
|
|
airMopAdd(mop, err=biffGetDone(TEN), airFree, airMopAlways); |
124 |
|
|
fprintf(stderr, "%s: trouble parsing model frmo axis 0 of nin:\n%s\n", |
125 |
|
|
me, err); |
126 |
|
|
airMopError(mop); return 1; |
127 |
|
|
} |
128 |
|
|
} else { |
129 |
|
|
fprintf(stderr, "%s: need model specified either via \"-m\" or input " |
130 |
|
|
"model image axis 0\n", me); |
131 |
|
|
airMopError(mop); return 1; |
132 |
|
|
} |
133 |
|
|
/* we have learned plusB0, but we don't actually need it; |
134 |
|
|
either: it describes the given model param image |
135 |
|
|
(which is courteous but not necessary since the logic inside |
136 |
|
|
tenModeSimulate will see this), |
137 |
|
|
or: it is trying to say something about including B0 amongst |
138 |
|
|
model parameters (which isn't actually meaningful in the |
139 |
|
|
context of simulated DWIs */ |
140 |
|
|
E = 0; |
141 |
|
|
if (!E) E |= tenGradientCheck(ngrad, nrrdTypeDouble, 1); |
142 |
|
|
if (!E) E |= tenExperSpecGradSingleBValSet(espec, insertB0, bval, |
143 |
|
|
AIR_CAST(const double *, |
144 |
|
|
ngrad->data), |
145 |
|
|
ngrad->axis[1].size); |
146 |
|
|
if (!E) E |= tenModelSimulate(nout, outType, espec, |
147 |
|
|
model, nT2, nin, keyValueSet); |
148 |
|
|
if (E) { |
149 |
|
|
airMopAdd(mop, err=biffGetDone(TEN), airFree, airMopAlways); |
150 |
|
|
fprintf(stderr, "%s: trouble:\n%s\n", me, err); |
151 |
|
|
airMopError(mop); return 1; |
152 |
|
|
} |
153 |
|
|
if (nrrdSave(outS, nout, NULL)) { |
154 |
|
|
airMopAdd(mop, err=biffGetDone(NRRD), airFree, airMopAlways); |
155 |
|
|
fprintf(stderr, "%s: trouble writing:\n%s\n", me, err); |
156 |
|
|
airMopError(mop); return 1; |
157 |
|
|
} |
158 |
|
|
|
159 |
|
|
airMopOkay(mop); |
160 |
|
|
return 0; |
161 |
|
1 |
} |
162 |
|
|
TEND_CMD(msim, INFO); |