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 "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
tend_msimMain(int argc,const char ** argv,const char * me,hestParm * hparm)40 tend_msimMain(int argc, const char **argv, const char *me,
41 hestParm *hparm) {
42 int pret;
43 hestOpt *hopt = NULL;
44 char *perr, *err;
45 airArray *mop;
46
47 tenExperSpec *espec;
48 const tenModel *model;
49 int E, seed, keyValueSet, outType, plusB0, insertB0;
50 Nrrd *nin, *nT2, *_ngrad, *ngrad, *nout;
51 char *outS, *modS;
52 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 hparm->elideSingleOtherDefault = AIR_TRUE;
57
58 hestOptAdd(&hopt, "sigma", "sigma", airTypeDouble, 1, 1, &sigma, "0.0",
59 "Gaussian/Rician noise parameter");
60 hestOptAdd(&hopt, "seed", "seed", airTypeInt, 1, 1, &seed, "42",
61 "seed value for RNG which creates noise");
62 hestOptAdd(&hopt, "g", "grad list", airTypeOther, 1, 1, &_ngrad, NULL,
63 "gradient list, one row per diffusion-weighted image",
64 NULL, NULL, nrrdHestNrrd);
65 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 NULL, NULL, nrrdHestNrrd);
69 hestOptAdd(&hopt, "i", "model image", airTypeOther, 1, 1, &nin, "-",
70 "input model image", NULL, NULL, nrrdHestNrrd);
71 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 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 hestOptAdd(&hopt, "b", "b", airTypeDouble, 1, 1, &bval, "1000",
80 "b value for simulated scan");
81 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 hestOptAdd(&hopt, "t", "type", airTypeEnum, 1, 1, &outType, "float",
85 "output type of DWIs", NULL, nrrdType);
86 hestOptAdd(&hopt, "o", "nout", airTypeString, 1, 1, &outS, "-",
87 "output dwis");
88
89 mop = airMopNew();
90 airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways);
91 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 }
162 /* TEND_CMD(msim, INFO); */
163 unrrduCmd tend_msimCmd = { "msim", INFO, tend_msimMain };
164