ANIMA  4.0
animaMCMEstimator.cxx
Go to the documentation of this file.
1 #include <fstream>
2 #include <itkTimeProbe.h>
3 
4 #include <tclap/CmdLine.h>
5 
10 
11 //Update progression of the process
12 void eventCallback (itk::Object* caller, const itk::EventObject& event, void* clientData)
13 {
14  itk::ProcessObject * processObject = (itk::ProcessObject*) caller;
15  std::cout<<"\033[K\rProgression: "<<(int)(processObject->GetProgress() * 100)<<"%"<<std::flush;
16 }
17 
18 int main(int argc, char **argv)
19 {
20  TCLAP::CmdLine cmd("INRIA / IRISA - VisAGeS/Empenn Team", ' ',ANIMA_VERSION);
21 
22  // Required arguments
23  TCLAP::ValueArg<std::string> dwiArg("i", "dwi", "DWI 4D Volume", true, "", "DWI images", cmd);
24  TCLAP::ValueArg<std::string> gradsArg("g", "grads", "Gradient table", true, "", "gradients", cmd);
25  TCLAP::ValueArg<std::string> bvalsArg("b", "bvals", "B-value list", true, "", "b-values", cmd);
26  TCLAP::SwitchArg bvalueScaleArg("B","b-no-scale","Do not scale b-values according to gradient norm",cmd);
27  TCLAP::ValueArg<double> smallDeltaArg("", "small-delta", "Diffusion small delta (in seconds)", false, anima::DiffusionSmallDelta, "small delta", cmd);
28  TCLAP::ValueArg<double> bigDeltaArg("", "big-delta", "Diffusion big delta (in seconds)", false, anima::DiffusionBigDelta, "big delta", cmd);
29 
30  // Outputs
31  TCLAP::ValueArg<std::string> outArg("o", "out-mcm", "MCM output volume", true, "", "MCM output", cmd);
32  TCLAP::ValueArg<std::string> aicArg("a", "out-aic", "Output estimated AICu image", false, "", "AICu output image", cmd);
33  TCLAP::ValueArg<std::string> outB0Arg("", "out-b0", "Output estimated B0 image", false, "", "B0 output image", cmd);
34  TCLAP::ValueArg<std::string> outSigmaArg("", "out-sig", "Output estimated noise sigma square image", false, "", "noise sigma square output image", cmd);
35  TCLAP::ValueArg<std::string> outMoseArg("", "out-mose", "Output model selection map", false, "", "model selection output image", cmd);
36  TCLAP::ValueArg<std::string> inMoseArg("I", "in-mose", "Input model selection map (overrides model selection step)", false, "", "model selection input image", cmd);
37 
38  // Optional arguments
39  TCLAP::ValueArg<std::string> computationMaskArg("m", "mask", "Computation mask", false, "", "computation mask", cmd);
40  TCLAP::ValueArg<double> b0thrArg("", "b0-thr", "Background threshold on B0 value (default: 10)", false, 10.0, "B0 theshold", cmd);
41 
42  TCLAP::ValueArg<unsigned int> nbFasciclesArg("n", "nb-fascicles", "Number of computed fascicles (default: 2)", false, 2, "number of fascicles", cmd);
43  TCLAP::ValueArg<unsigned int> compartmentTypeArg("c", "comp-type", "Compartment type for fascicles: 1: stick, 2: zeppelin, 3: tensor, 4: noddi (default: 3)", false, 3, "fascicles type", cmd);
44  TCLAP::SwitchArg aicSelectNbCompartmentsArg("M", "opt-nb-comp", "Activate AICC-based number of compartments selection", cmd, false);
45 
46  TCLAP::SwitchArg freeWaterCompartmentArg("F", "free-water", "Model with free water", cmd, false);
47  TCLAP::SwitchArg stationaryWaterCompartmentArg("S", "stationary-water", "Model with stationary water", cmd, false);
48  TCLAP::SwitchArg restrictedWaterCompartmentArg("R", "restricted-water", "Model with restricted water", cmd, false);
49  TCLAP::SwitchArg staniszCompartmentArg("Z", "stanisz", "Model with stanisz isotropic compartment", cmd, false);
50 
51  TCLAP::SwitchArg optFWDiffArg("", "opt-free-water-diff", "Optimize free water diffusivity value", cmd, false);
52  TCLAP::SwitchArg optIRWDiffArg("", "opt-ir-water-diff", "Optimize isotropic restricted water diffusivity value", cmd, false);
53  TCLAP::SwitchArg optStaniszRadiusArg("", "opt-stanisz-radius", "Optimize isotropic Stanisz radius value", cmd, false);
54  TCLAP::SwitchArg optStaniszDiffArg("", "opt-stanisz-diff", "Optimize isotropic Stanisz diffusivity value", cmd, false);
55 
56  TCLAP::SwitchArg fixDiffArg("", "fix-diff", "Fix diffusivity value", cmd, false);
57  TCLAP::SwitchArg fixKappaArg("", "fix-kappa", "Fix orientation concentration values", cmd, false);
58  TCLAP::SwitchArg fixEAFArg("", "fix-eaf", "Fix extra axonal fraction values", cmd, false);
59 
60  TCLAP::SwitchArg commonDiffusivitiesArg("", "common-diffusivities", "Share diffusivity values among compartments", cmd, false);
61  TCLAP::SwitchArg commonKappaArg("", "common-kappa", "Share orientation concentration values among compartments", cmd, false);
62  TCLAP::SwitchArg commonEAFArg("", "common-eaf", "Share extra axonal fraction values among compartments", cmd, false);
63 
64  //Initial values for diffusivities
65  TCLAP::ValueArg<double> initAxialDiffArg("", "init-axial-diff", "Initial axial diffusivity (default: 1.71e-3)", false, 1.71e-3, "initial axial diffusivity", cmd);
66  TCLAP::ValueArg<double> initRadialDiff1Arg("", "init-radial-diff1", "Initial first radial diffusivity (default: 1.9e-4)", false, 1.9e-4, "initial first radial diffusivity", cmd);
67  TCLAP::ValueArg<double> initRadialDiff2Arg("", "init-radial-diff2", "Initial second radial diffusivity (default: 1.5e-4)", false, 1.5e-4, "initial second radial diffusivity", cmd);
68  TCLAP::ValueArg<double> initIRWDiffArg("", "init-irw-diff", "Initial isotropic restricted diffusivity (default: 7.5e-4)", false, 7.5e-4, "initial IRW diffusivity", cmd);
69  TCLAP::ValueArg<double> initStaniszDiffArg("", "init-stanisz-diff", "Initial Stanisz diffusivity (default: 1.71e-3)", false, 1.71e-3, "initial Stanisz diffusivity", cmd);
70 
71  // Optimization parameters
72  TCLAP::ValueArg<std::string> optimizerArg("", "optimizer", "Optimizer for estimation: bobyqa (default), ccsaq, bfgs or levenberg", false, "bobyqa", "optimizer", cmd);
73  TCLAP::ValueArg<double> absCostChangeArg("", "abs-cost-change", "Cost function change to stop estimation (default: 0.01)", false, 0.01, "cost change threshold", cmd);
74  TCLAP::ValueArg <unsigned int> mlModeArg("", "ml-mode", "ML estimation strategy: marginal likelihood (0), profile likelihood (1, default), Variable projection (2)", false, 1, "ML mode", cmd);
75  TCLAP::ValueArg<double> xTolArg("x", "x-tol", "Tolerance for relative position in optimization (default: 0 -> 1.0e-4 or 1.0e-7 for bobyqa)", false, 0, "position relative tolerance", cmd);
76  TCLAP::ValueArg<double> fTolArg("", "f-tol", "Tolerance for relative cost in optimization (default: 0 -> function of position tolerance)", false, 0, "cost relative tolerance", cmd);
77  TCLAP::ValueArg<unsigned int> maxEvalArg("e", "max-eval", "Maximum evaluations (default: 0 -> function of number of unknowns)", false, 0, "max evaluations", cmd);
78 
79  TCLAP::ValueArg<unsigned int> nbThreadsArg("T", "nb-threads", "Number of threads to run on (default: all cores)", false, itk::MultiThreaderBase::GetGlobalDefaultNumberOfThreads(), "number of threads", cmd);
80 
81  try
82  {
83  cmd.parse(argc,argv);
84  }
85  catch (TCLAP::ArgException& e)
86  {
87  std::cerr << "Error: " << e.error() << "for argument " << e.argId() << std::endl;
88  return EXIT_FAILURE;
89  }
90 
92  typedef FilterType::InputImageType InputImageType;
93  typedef FilterType::MaskImageType MaskImageType;
94  typedef FilterType::Pointer FilterPointer;
95 
96  itk::CStyleCommand::Pointer callback = itk::CStyleCommand::New();
97  callback->SetCallback(eventCallback);
98 
99  FilterPointer filter = FilterType::New();
100 
101  filter->SetUseConstrainedOrientationConcentration(fixKappaArg.isSet());
102  if (!fixKappaArg.isSet())
103  filter->SetUseCommonConcentrations(commonKappaArg.isSet());
104  else
105  filter->SetUseCommonConcentrations(false);
106 
107  filter->SetUseConstrainedExtraAxonalFraction(fixEAFArg.isSet());
108  if (!fixEAFArg.isSet())
109  filter->SetUseCommonExtraAxonalFractions(commonEAFArg.isSet());
110  else
111  filter->SetUseCommonExtraAxonalFractions(false);
112 
113  std::cout << "Loading input DWI image..." << std::endl;
114 
115  anima::setMultipleImageFilterInputsFromFileName<InputImageType,FilterType>(dwiArg.getValue(), filter);
116 
117  // Load gradient table and b-value list
118  std::cout << "Importing gradient table and b-values..." << std::endl;
119 
120  typedef anima::GradientFileReader < vnl_vector_fixed<double,3>, double > GFReaderType;
121  GFReaderType gfReader;
122  gfReader.SetGradientFileName(gradsArg.getValue());
123  gfReader.SetBValueBaseString(bvalsArg.getValue());
124  gfReader.SetGradientIndependentNormalization(bvalueScaleArg.isSet());
125  gfReader.SetSmallDelta(smallDeltaArg.getValue());
126  gfReader.SetBigDelta(bigDeltaArg.getValue());
127  gfReader.Update();
128 
129  GFReaderType::GradientVectorType directions = gfReader.GetGradients();
130  for(unsigned int i = 0;i < directions.size();++i)
131  filter->AddGradientDirection(i,directions[i]);
132 
133  GFReaderType::BValueVectorType mb = gfReader.GetGradientStrengths();
134 
135  filter->SetGradientStrengths(mb);
136  filter->SetSmallDelta(smallDeltaArg.getValue());
137  filter->SetBigDelta(bigDeltaArg.getValue());
138 
139  if (computationMaskArg.getValue() != "")
140  filter->SetComputationMask(anima::readImage<MaskImageType>(computationMaskArg.getValue()));
141 
142  if (inMoseArg.getValue() != "")
143  filter->SetMoseVolume(anima::readImage<FilterType::MoseImageType>(inMoseArg.getValue()));
144 
145  filter->SetB0Threshold(b0thrArg.getValue());
146 
147  filter->SetModelWithFreeWaterComponent(freeWaterCompartmentArg.isSet());
148  filter->SetModelWithStationaryWaterComponent(stationaryWaterCompartmentArg.isSet());
149  filter->SetModelWithRestrictedWaterComponent(restrictedWaterCompartmentArg.isSet());
150  filter->SetModelWithStaniszComponent(staniszCompartmentArg.isSet());
151 
152  switch (compartmentTypeArg.getValue())
153  {
154  case 1:
155  filter->SetCompartmentType(anima::Stick);
156  break;
157 
158  case 2:
159  filter->SetCompartmentType(anima::Zeppelin);
160  break;
161 
162  case 3:
163  filter->SetCompartmentType(anima::Tensor);
164  break;
165 
166  case 4:
167  filter->SetCompartmentType(anima::NODDI);
168  break;
169 
170  default:
171  std::cerr << "Unsupported compartment type" << std::endl;
172  return EXIT_FAILURE;
173  }
174 
175  filter->SetAxialDiffusivityValue(initAxialDiffArg.getValue());
176  filter->SetRadialDiffusivity1Value(initRadialDiff1Arg.getValue());
177  filter->SetRadialDiffusivity2Value(initRadialDiff2Arg.getValue());
178  filter->SetIRWDiffusivityValue(initIRWDiffArg.getValue());
179  filter->SetStaniszDiffusivityValue(initStaniszDiffArg.getValue());
180 
181  filter->SetNumberOfCompartments(nbFasciclesArg.getValue());
182  filter->SetFindOptimalNumberOfCompartments(aicSelectNbCompartmentsArg.isSet());
183 
184  filter->SetOptimizer(optimizerArg.getValue());
185  filter->SetAbsoluteCostChange(absCostChangeArg.getValue());
186 
187  filter->SetNoiseType(FilterType::Gaussian);
188  filter->SetMLEstimationStrategy((FilterType::MaximumLikelihoodEstimationMode)mlModeArg.getValue());
189 
190  filter->SetXTolerance(xTolArg.getValue());
191  filter->SetFTolerance(fTolArg.getValue());
192  filter->SetMaxEval(maxEvalArg.getValue());
193 
194  filter->SetUseConstrainedDiffusivity(fixDiffArg.isSet());
195  filter->SetUseConstrainedFreeWaterDiffusivity(!optFWDiffArg.isSet());
196  filter->SetUseConstrainedIRWDiffusivity(!optIRWDiffArg.isSet());
197  filter->SetUseConstrainedStaniszDiffusivity(!optStaniszDiffArg.isSet());
198  filter->SetUseConstrainedStaniszRadius(!optStaniszRadiusArg.isSet());
199 
200  if (!fixDiffArg.isSet())
201  filter->SetUseCommonDiffusivities(commonDiffusivitiesArg.isSet());
202  else
203  filter->SetUseCommonDiffusivities(false);
204 
205  filter->SetNumberOfWorkUnits(nbThreadsArg.getValue());
206  filter->AddObserver(itk::ProgressEvent(), callback);
207 
208  itk::TimeProbe tmpTimer;
209  tmpTimer.Start();
210 
211  try
212  {
213  filter->Update();
214  }
215  catch (itk::ExceptionObject &e)
216  {
217  std::cerr << e << std::endl;
218  return EXIT_FAILURE;
219  }
220 
221  tmpTimer.Stop();
222 
223  std::cout << "\nEstimation done in " << tmpTimer.GetTotal() << " s" << std::endl;
224  std::cout << "Writing MCM to: " << outArg.getValue() << std::endl;
225 
226  try
227  {
228  filter->WriteMCMOutput(outArg.getValue());
229  }
230  catch (std::exception &e)
231  {
232  std::cerr << "Error writing output " << e.what() << std::endl;
233  }
234 
235  if (aicArg.getValue() != "")
236  {
237  std::cout << "Writing AICu image to: " << aicArg.getValue() << std::endl;
238  anima::writeImage(aicArg.getValue(),filter->GetAICcVolume());
239  }
240 
241  if (outB0Arg.getValue() != "")
242  {
243  std::cout << "Writing B0 image to: " << outB0Arg.getValue() << std::endl;
244  anima::writeImage(outB0Arg.getValue(),filter->GetB0Volume());
245  }
246 
247  if (outSigmaArg.getValue() != "")
248  {
249  std::cout << "Writing noise sigma square image to: " << outSigmaArg.getValue() << std::endl;
250  anima::writeImage(outSigmaArg.getValue(),filter->GetSigmaSquareVolume());
251  }
252 
253  if (outMoseArg.getValue() != "")
254  {
255  std::cout << "Writing model selection image to: " << outMoseArg.getValue() << std::endl;
256  anima::writeImage(outMoseArg.getValue(),filter->GetMoseVolume());
257  }
258 
259  return EXIT_SUCCESS;
260 }
int main(int argc, char **argv)
void eventCallback(itk::Object *caller, const itk::EventObject &event, void *clientData)
const double DiffusionBigDelta
Default big delta value (classical values)
const double DiffusionSmallDelta
Default small delta value (classical values)
void writeImage(std::string filename, OutputImageType *img)
void SetGradientFileName(std::string fName)