ANIMA  4.0
animaMultiT2RelaxometryEstimation.cxx
Go to the documentation of this file.
1 #include <iostream>
2 #include <tclap/CmdLine.h>
3 
4 #include <itkImageFileReader.h>
5 #include <itkImageFileWriter.h>
6 #include <animaMultiT2RelaxometryEstimationImageFilter.h>
7 #include <animaReadWriteFunctions.h>
8 #include <itkTimeProbe.h>
9 
10 #include <boost/lexical_cast.hpp>
11 #include <boost/algorithm/string.hpp>
12 
13 #include <itkCommand.h>
14 
15 //Update progression of the process
16 void eventCallback (itk::Object* caller, const itk::EventObject& event, void* clientData)
17 {
18  itk::ProcessObject * processObject = (itk::ProcessObject*) caller;
19  std::cout<<"\033[K\rProgression: "<<(int)(processObject->GetProgress() * 100)<<"%"<<std::flush;
20 }
21 
22 int main(int argc, char **argv)
23 {
24  TCLAP::CmdLine cmd("INRIA / IRISA - VisAGeS/Empenn Team", ' ',ANIMA_VERSION);
25 
26  TCLAP::ValueArg<std::string> t2Arg("l","t2","List of T2 relaxometry images",true,"","T2 relaxometry images",cmd);
27  TCLAP::ValueArg<std::string> maskArg("m","maskname","Computation mask",false,"","computation mask",cmd);
28 
29  TCLAP::ValueArg<std::string> t1MapArg("","t1","T1 map",false,"","T1 map",cmd);
30 
31  TCLAP::ValueArg<std::string> resT2Arg("O","out-t2","Result T2 image",false,"","result T2 image",cmd);
32  TCLAP::ValueArg<std::string> resM0Arg("","out-m0","Result M0 image",false,"","result M0 image",cmd);
33  TCLAP::ValueArg<std::string> resMWFArg("o","out-mwf","Result MWF image",true,"","result MWF image",cmd);
34  TCLAP::ValueArg<std::string> resB1Arg("","out-b1","Result B1 image",false,"","result B1 image",cmd);
35  TCLAP::ValueArg<std::string> resCostArg("c","out-cost","Result cost image",false,"","result cost image",cmd);
36 
37  TCLAP::ValueArg<double> echoSpacingArg("e","echo-spacing","Spacing between two successive echoes (default: 10)",false,10,"Spacing between echoes",cmd);
38  TCLAP::ValueArg<double> excitationT2FlipAngleArg("","t2-ex-flip","Excitation flip angle for T2 (in degrees, default: 90)",false,90,"T2 excitation flip angle",cmd);
39  TCLAP::ValueArg<double> t2FlipAngleArg("","t2-flip","All flip angles for T2 (in degrees, default: 180)",false,180,"T2 flip angle",cmd);
40  TCLAP::ValueArg<double> backgroundSignalThresholdArg("t","signal-thr","Background signal threshold (default: 10)",false,10,"Background signal threshold",cmd);
41  TCLAP::ValueArg<double> regulIntensityArg("r","reg-int","Tikhonov regularization intensity (default: 1.08)",false,1.08,"regularization intensity",cmd);
42  TCLAP::ValueArg<double> nlRegulIntensityArg("","nl-reg-int","Tikhonov regularization intensity in NL mode (default: 1.5)",false,1.5,"regularization intensity",cmd);
43 
44  TCLAP::ValueArg<double> myelinThrArg("","mye-thr","T2 myelin threshold for MWF computation (default: 50)",false,50,"T2 myelin threshold",cmd);
45  TCLAP::ValueArg<unsigned int> numT2CompartmentsArg("n","num-t2","Number of T2 compartments (default: 40)",false,40,"Number of T2 compartments",cmd);
46  TCLAP::ValueArg<double> t2LowerBoundArg("","low-t2","Lower T2 value (default: 15)",false,15,"T2 lower value",cmd);
47  TCLAP::ValueArg<double> t2UpperBoundArg("","up-t2","Upper T2 value (default: 2000)",false,2000,"T2 upper value",cmd);
48 
49  TCLAP::ValueArg<unsigned int> b1NumOptimizerIterArg("","b1-opt-iter","Maximal number of optimizer iterations (default: 200)",false,200,"Maximal number of optimizer iterations",cmd);
50  TCLAP::ValueArg<double> b1OptimizerStopConditionArg("","b1-opt-stop","Optimizer stopping threshold (default: 1.0e-4)",false,1.0e-4,"B1 Optimizer stopping threshold",cmd);
51  TCLAP::ValueArg<double> b1OptimizerInitialStepArg("i","b1-opt-init","Optimizer initial step (default: 0.1)",false,0.1,"B1 Optimizer initial step",cmd);
52  TCLAP::ValueArg<double> b1ToleranceArg("","b1-tol","B1 tolerance threshold of convergence (default: 1.0e-4)",false,1.0e-4,"B1 tolerance for optimization convergence",cmd);
53 
54  //NL params
55  TCLAP::SwitchArg nlEstimationArg("N","nl-est","NL estimation of MWF",cmd, false);
56  TCLAP::ValueArg<double> weightThrArg("w","weightThr","Weight threshold: patches around have to be similar enough -> default: 0.0",false,0.0,"Weight threshold",cmd);
57  TCLAP::ValueArg<double> betaArg("b","beta","Beta parameter for local noise estimation -> default: 1",false,1,"Beta for local noise estimation",cmd);
58  TCLAP::ValueArg<double> meanMinArg("","meanMin","Minimun mean threshold (default: 0.95)",false,0.95,"Minimun mean threshold",cmd);
59  TCLAP::ValueArg<double> varMinArg("v","varMin","Minimun variance threshold -> default: 0.5",false,0.5,"Minimun variance threshold",cmd);
60  TCLAP::ValueArg<unsigned int> patchHSArg("S","patchHalfSize","Patch half size in each direction -> default: 1",false,1,"patch half size",cmd);
61  TCLAP::ValueArg<unsigned int> patchSSArg("s","patchStepSize","Patch step size for searching -> default: 1",false,1,"Patch search step size",cmd);
62  TCLAP::ValueArg<unsigned int> patchNeighArg("","patchNeighborhood","Patch half neighborhood size -> default: 5",false,5,"Patch search neighborhood size",cmd);
63 
64  TCLAP::ValueArg<unsigned int> nbpArg("T","numberofthreads","Number of threads to run on (default : all cores)",false,itk::MultiThreaderBase::GetGlobalDefaultNumberOfThreads(),"number of threads",cmd);
65 
66  try
67  {
68  cmd.parse(argc,argv);
69  }
70  catch (TCLAP::ArgException& e)
71  {
72  std::cerr << "Error: " << e.error() << "for argument " << e.argId() << std::endl;
73  return EXIT_FAILURE;
74  }
75 
76  typedef itk::Image <double,3> InputImageType;
77  typedef InputImageType OutputImageType;
78  typedef anima::MultiT2RelaxometryEstimationImageFilter <double> FilterType;
79  typedef FilterType::OutputImageType OutputImageType;
80  typedef FilterType::VectorOutputImageType VectorOutputImageType;
81  typedef itk::ImageFileReader <InputImageType> InputImageReaderType;
82 
83  FilterType::Pointer mainFilter = FilterType::New();
84 
85  unsigned int numInputs = anima::setMultipleImageFilterInputsFromFileName<InputImageType,FilterType>(t2Arg.getValue(),mainFilter);
86 
87  mainFilter->SetEchoSpacing(echoSpacingArg.getValue());
88  mainFilter->SetT2FlipAngles(t2FlipAngleArg.getValue() * M_PI / 180.0,numInputs);
89  mainFilter->SetT2ExcitationFlipAngle(excitationT2FlipAngleArg.getValue() * M_PI / 180.0);
90  mainFilter->SetB1MaximumOptimizerIterations(b1NumOptimizerIterArg.getValue());
91  mainFilter->SetB1OptimizerStopCondition(b1OptimizerStopConditionArg.getValue());
92  mainFilter->SetB1OptimizerInitialStep(b1OptimizerInitialStepArg.getValue());
93  mainFilter->SetB1Tolerance(b1ToleranceArg.getValue());
94 
95  mainFilter->SetLowerT2Bound(t2LowerBoundArg.getValue());
96  mainFilter->SetUpperT2Bound(t2UpperBoundArg.getValue());
97  mainFilter->SetMyelinThreshold(myelinThrArg.getValue());
98  mainFilter->SetNumberOfT2Compartments(numT2CompartmentsArg.getValue());
99  mainFilter->SetRegularizationIntensity(regulIntensityArg.getValue());
100  mainFilter->SetNLEstimation(false);
101 
102  if (t1MapArg.getValue() != "")
103  {
104  InputImageReaderType::Pointer t1MapRead = InputImageReaderType::New();
105  t1MapRead->SetFileName(t1MapArg.getValue().c_str());
106  t1MapRead->Update();
107 
108  mainFilter->SetT1Map(t1MapRead->GetOutput());
109  }
110 
111  if (maskArg.getValue() != "")
112  {
113  typedef itk::ImageFileReader < itk::Image <unsigned char, 3> > itkMaskReader;
114  itkMaskReader::Pointer maskRead = itkMaskReader::New();
115  maskRead->SetFileName(maskArg.getValue().c_str());
116  maskRead->Update();
117 
118  mainFilter->SetComputationMask(maskRead->GetOutput());
119  }
120 
121  mainFilter->SetAverageSignalThreshold(backgroundSignalThresholdArg.getValue());
122  mainFilter->SetNumberOfWorkUnits(nbpArg.getValue());
123 
124  itk::CStyleCommand::Pointer callback = itk::CStyleCommand::New();
125  callback->SetCallback(eventCallback);
126  mainFilter->AddObserver(itk::ProgressEvent(), callback );
127 
128  itk::TimeProbe tmpTime;
129  tmpTime.Start();
130 
131  try
132  {
133  mainFilter->Update();
134  }
135  catch (itk::ExceptionObject &e)
136  {
137  std::cerr << e << std::endl;
138  }
139 
140  tmpTime.Stop();
141 
142  std::cout << "Regular estimation computation time: " << tmpTime.GetTotal() << std::endl;
143 
144  FilterType::OutputImagePointer outMWFImage = mainFilter->GetMWFOutputImage();
145  outMWFImage->DisconnectPipeline();
146  FilterType::OutputImagePointer outM0Image = mainFilter->GetM0OutputImage();
147  outM0Image->DisconnectPipeline();
148  FilterType::OutputImagePointer outCostImage = mainFilter->GetCostOutputImage();
149  outCostImage->DisconnectPipeline();
150  FilterType::OutputImagePointer outB1Image = mainFilter->GetB1OutputImage();
151  outB1Image->DisconnectPipeline();
152  FilterType::VectorOutputImagePointer outT2Image = mainFilter->GetT2OutputImage();
153  outT2Image->DisconnectPipeline();
154 
155  if (nlEstimationArg.isSet())
156  {
157  FilterType::Pointer secondaryFilter = FilterType::New();
158  for (unsigned int i = 0;i < mainFilter->GetNumberOfIndexedInputs();++i)
159  secondaryFilter->SetInput(i,mainFilter->GetInput(i));
160 
161  secondaryFilter->SetEchoSpacing(echoSpacingArg.getValue());
162  secondaryFilter->SetT2FlipAngles(t2FlipAngleArg.getValue() * M_PI / 180.0,numInputs);
163  secondaryFilter->SetT2ExcitationFlipAngle(excitationT2FlipAngleArg.getValue() * M_PI / 180.0);
164  secondaryFilter->SetB1MaximumOptimizerIterations(b1NumOptimizerIterArg.getValue());
165  secondaryFilter->SetB1OptimizerStopCondition(b1OptimizerStopConditionArg.getValue());
166  secondaryFilter->SetB1OptimizerInitialStep(b1OptimizerInitialStepArg.getValue());
167  secondaryFilter->SetB1Tolerance(b1ToleranceArg.getValue());
168 
169  secondaryFilter->SetLowerT2Bound(t2LowerBoundArg.getValue());
170  secondaryFilter->SetUpperT2Bound(t2UpperBoundArg.getValue());
171  secondaryFilter->SetMyelinThreshold(myelinThrArg.getValue());
172  secondaryFilter->SetNumberOfT2Compartments(numT2CompartmentsArg.getValue());
173  secondaryFilter->SetRegularizationIntensity(nlRegulIntensityArg.getValue());
174 
175  secondaryFilter->SetT1Map(mainFilter->GetT1Map());
176  secondaryFilter->SetComputationMask(mainFilter->GetComputationMask());
177 
178  secondaryFilter->SetAverageSignalThreshold(backgroundSignalThresholdArg.getValue());
179  secondaryFilter->SetNumberOfWorkUnits(nbpArg.getValue());
180 
181  itk::CStyleCommand::Pointer secondaryCallback = itk::CStyleCommand::New();
182  secondaryCallback->SetCallback(eventCallback);
183  secondaryFilter->AddObserver(itk::ProgressEvent(), secondaryCallback);
184 
185  secondaryFilter->SetNLEstimation(true);
186  secondaryFilter->SetWeightThreshold(weightThrArg.getValue());
187  secondaryFilter->SetBetaParameter(betaArg.getValue());
188  secondaryFilter->SetMeanMinThreshold(meanMinArg.getValue());
189  secondaryFilter->SetVarMinThreshold(varMinArg.getValue());
190  secondaryFilter->SetPatchHalfSize(patchHSArg.getValue());
191  secondaryFilter->SetSearchStepSize(patchSSArg.getValue());
192  secondaryFilter->SetSearchNeighborhood(patchNeighArg.getValue());
193 
194  secondaryFilter->SetInitialB1Map(outB1Image);
195  secondaryFilter->SetInitialT2Map(outT2Image);
196  secondaryFilter->SetInitialM0Map(outM0Image);
197 
198  itk::TimeProbe tmpTimeNL;
199  tmpTimeNL.Start();
200 
201  try
202  {
203  secondaryFilter->Update();
204  }
205  catch (itk::ExceptionObject &e)
206  {
207  std::cerr << e << std::endl;
208  }
209 
210  tmpTimeNL.Stop();
211 
212  std::cout << "NL estimation computation time: " << tmpTimeNL.GetTotal() << std::endl;
213 
214  outMWFImage = secondaryFilter->GetMWFOutputImage();
215  outMWFImage->DisconnectPipeline();
216  outM0Image = secondaryFilter->GetM0OutputImage();
217  outM0Image->DisconnectPipeline();
218  outCostImage = secondaryFilter->GetCostOutputImage();
219  outCostImage->DisconnectPipeline();
220  outB1Image = secondaryFilter->GetB1OutputImage();
221  outB1Image->DisconnectPipeline();
222  outT2Image = secondaryFilter->GetT2OutputImage();
223  outT2Image->DisconnectPipeline();
224  }
225 
226  anima::writeImage<OutputImageType> (resMWFArg.getValue(),outMWFImage);
227 
228  if (resT2Arg.getValue() != "")
229  anima::writeImage<VectorOutputImageType> (resT2Arg.getValue(),outT2Image);
230 
231  if (resM0Arg.getValue() != "")
232  anima::writeImage<OutputImageType> (resM0Arg.getValue(),outM0Image);
233 
234  if (resB1Arg.getValue() != "")
235  anima::writeImage<OutputImageType> (resB1Arg.getValue(),outB1Image);
236 
237  if (resCostArg.getValue() != "")
238  anima::writeImage<OutputImageType> (resCostArg.getValue(),outCostImage);
239 
240  return EXIT_SUCCESS;
241 }
main
int main(int argc, char **argv)
Definition: animaMultiT2RelaxometryEstimation.cxx:22
anima::MultiT2RelaxometryEstimationImageFilter
Implements multi-peak T2 relaxometry estimation (with or without regularization)
Definition: animaMultiT2RelaxometryEstimationImageFilter.h:23
eventCallback
void eventCallback(itk::Object *caller, const itk::EventObject &event, void *clientData)
Definition: animaMultiT2RelaxometryEstimation.cxx:16
Generated by   doxygen 1.8.13