ANIMA  4.0
animaT1SERelaxometryEstimation.cxx
Go to the documentation of this file.
1 #include <iostream>
2 #include <tclap/CmdLine.h>
3 
6 #include <itkTimeProbe.h>
7 
8 int main(int argc, char **argv)
9 {
10  TCLAP::CmdLine cmd("INRIA / IRISA - VisAGeS/Empenn Team", ' ',ANIMA_VERSION);
11 
12  TCLAP::ValueArg<std::string> t1Arg("l","t1","List of T1 relaxometry images",true,"","T1 relaxometry images",cmd);
13  TCLAP::ValueArg<std::string> maskArg("m","maskname","Computation mask",false,"","computation mask",cmd);
14 
15  TCLAP::ValueArg<std::string> resT1Arg("o","out-t1","Result T1 image",true,"","result T1 image",cmd);
16  TCLAP::ValueArg<std::string> resM0Arg("O","out-m0","Result M0 image",false,"","result M0 image",cmd);
17 
18  TCLAP::ValueArg<double> upperBoundT1Arg("u","upper-bound-t1","T1 value upper bound (default: 5000)",false,5000,"T1 value upper bound",cmd);
19  TCLAP::ValueArg<double> upperBoundM0Arg("","upper-bound-m0","M0 value upper bound (default: 5000)",false,5000,"M0 value upper bound",cmd);
20 
21  TCLAP::ValueArg<std::string> trArg("","tr","Text file listing repetition times for input images",true,"","repetition times list",cmd);
22  TCLAP::ValueArg<double> backgroundSignalThresholdArg("t","signal-thr","Background signal threshold (default: 10)",false,10,"Background signal threshold",cmd);
23 
24  TCLAP::ValueArg<unsigned int> nbpArg("T","numberofthreads","Number of threads to run on (default : all cores)",false,itk::MultiThreaderBase::GetGlobalDefaultNumberOfThreads(),"number of threads",cmd);
25 
26  TCLAP::ValueArg<unsigned int> numOptimizerIterArg("","opt-iter","Maximal number of optimizer iterations (default: 200)",false,200,"Maximal number of optimizer iterations",cmd);
27  TCLAP::ValueArg<double> optimizerStopConditionArg("","opt-stop","Optimizer stopping threshold (default: 1.0e-4)",false,1.0e-4,"Optimizer stopping threshold",cmd);
28  TCLAP::ValueArg<double> optimizerInitialStepArg("i","opt-init","Optimizer initial step (default: 10)",false,10,"Optimizer initial step",cmd);
29 
30  try
31  {
32  cmd.parse(argc,argv);
33  }
34  catch (TCLAP::ArgException& e)
35  {
36  std::cerr << "Error: " << e.error() << "for argument " << e.argId() << std::endl;
37  return EXIT_FAILURE;
38  }
39 
40  typedef itk::Image <double,3> InputImageType;
41  typedef InputImageType OutputImageType;
43  typedef itk::ImageFileWriter <OutputImageType> OutputImageWriterType;
44 
45  FilterType::Pointer mainFilter = FilterType::New();
46 
47  mainFilter->SetM0UpperBound(upperBoundM0Arg.getValue());
48  mainFilter->SetT1UpperBound(upperBoundT1Arg.getValue());
49 
50  mainFilter->SetMaximumOptimizerIterations(numOptimizerIterArg.getValue());
51  mainFilter->SetOptimizerStopCondition(optimizerStopConditionArg.getValue());
52  mainFilter->SetOptimizerInitialStep(optimizerInitialStepArg.getValue());
53 
54  // Read and set T1 relaxometry images
55  anima::setMultipleImageFilterInputsFromFileName<InputImageType,FilterType>(t1Arg.getValue(),mainFilter);
56 
57  std::ifstream inputTRValues(trArg.getValue().c_str());
58  if (!inputTRValues.is_open())
59  {
60  std::cerr << "Text file for TR values not readable " << trArg.getValue() << std::endl;
61  return EXIT_FAILURE;
62  }
63 
64  std::vector <double> trValues;
65  while (!inputTRValues.eof())
66  {
67  char tmpStr[8192];
68  inputTRValues.getline(tmpStr,8192);
69  std::string workStr(tmpStr);
70  if (workStr == "")
71  continue;
72 
73  workStr.erase(workStr.find_last_not_of(" \n\r\t")+1);
74  std::istringstream iss(workStr);
75  std::string shortStr;
76  iss >> shortStr;
77  trValues.push_back(std::stod(shortStr));
78  }
79 
80  mainFilter->SetTRValues(trValues);
81 
82  if (maskArg.getValue() != "")
83  {
84  typedef itk::ImageFileReader < itk::Image <unsigned char, 3> > itkMaskReader;
85  itkMaskReader::Pointer maskRead = itkMaskReader::New();
86  maskRead->SetFileName(maskArg.getValue().c_str());
87  maskRead->Update();
88 
89  mainFilter->SetComputationMask(maskRead->GetOutput());
90  }
91 
92  mainFilter->SetAverageSignalThreshold(backgroundSignalThresholdArg.getValue());
93  mainFilter->SetNumberOfWorkUnits(nbpArg.getValue());
94 
95  itk::TimeProbe tmpTime;
96  tmpTime.Start();
97 
98  mainFilter->Update();
99 
100  tmpTime.Stop();
101 
102  std::cout << "Total computation time: " << tmpTime.GetTotal() << std::endl;
103 
104  anima::writeImage <OutputImageType> (resT1Arg.getValue(),mainFilter->GetOutput(1));
105 
106  if (resM0Arg.getValue() != "")
107  anima::writeImage <OutputImageType> (resM0Arg.getValue(),mainFilter->GetOutput(0));
108 
109  return EXIT_SUCCESS;
110 }
int main(int argc, char **argv)