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