ANIMA  4.0
animaT1RelaxometryEstimation.cxx
Go to the documentation of this file.
1 #include <iostream>
2 #include <tclap/CmdLine.h>
3 
4 #include <cmath>
5 
6 #include <itkImageFileReader.h>
7 #include <itkImageFileWriter.h>
10 #include <itkTimeProbe.h>
11 
12 int main(int argc, char **argv)
13 {
14  TCLAP::CmdLine cmd("INRIA / IRISA - VisAGeS/Empenn Team", ' ',ANIMA_VERSION);
15 
16  TCLAP::ValueArg<std::string> t1Arg("l","t1","List of T1 relaxometry images",true,"","T1 relaxometry images",cmd);
17  TCLAP::ValueArg<std::string> maskArg("m","maskname","Computation mask",false,"","computation mask",cmd);
18  TCLAP::ValueArg<std::string> b1Arg("b","b1","B1 inhomogenity map",false,"","B1 inhomogenity map",cmd);
19 
20  TCLAP::ValueArg<std::string> resT1Arg("o","out-t1","Result T1 image",true,"","result T1 image",cmd);
21  TCLAP::ValueArg<std::string> resM0Arg("O","out-m0","Result M0 image",false,"","result M0 image",cmd);
22 
23  TCLAP::ValueArg<double> trArg("","tr","Repetition time for T1 relaxometry (default: 15)",false,15,"repetition time",cmd);
24  TCLAP::ValueArg<std::string> flipAnglesArg("f","flip","Text file listing flip angles (in degrees)",true,"","T1 flip angles list",cmd);
25 
26  TCLAP::ValueArg<double> upperBoundT1Arg("u","upper-bound-t1","T1 value upper bound (default: 5000)",false,5000,"T1 upper bound",cmd);
27  TCLAP::ValueArg<double> upperBoundM0Arg("","upper-bound-m0","M0 value upper bound (default: 5000)",false,5000,"M0 upper bound",cmd);
28  TCLAP::ValueArg<double> backgroundSignalThresholdArg("t","signal-thr","Background signal threshold (default: 10)",false,10,"Background signal threshold",cmd);
29 
30  TCLAP::ValueArg<unsigned int> nbpArg("T","numberofthreads","Number of threads to run on (default : all cores)",false,itk::MultiThreaderBase::GetGlobalDefaultNumberOfThreads(),"number of threads",cmd);
31 
32  try
33  {
34  cmd.parse(argc,argv);
35  }
36  catch (TCLAP::ArgException& e)
37  {
38  std::cerr << "Error: " << e.error() << "for argument " << e.argId() << std::endl;
39  return EXIT_FAILURE;
40  }
41 
42  typedef itk::Image <double,3> InputImageType;
43  typedef InputImageType OutputImageType;
45  typedef itk::ImageFileReader <InputImageType> InputImageReaderType;
46  typedef itk::ImageFileWriter <OutputImageType> OutputImageWriterType;
47 
48  FilterType::Pointer mainFilter = FilterType::New();
49 
50  // Read and set T1 relaxometry images
51  anima::setMultipleImageFilterInputsFromFileName<InputImageType,FilterType>(t1Arg.getValue(),mainFilter);
52 
53  std::ifstream inputFlip(flipAnglesArg.getValue().c_str());
54  if (!inputFlip.is_open())
55  {
56  std::cerr << "Text file for flip angles not readable " << flipAnglesArg.getValue() << std::endl;
57  return EXIT_FAILURE;
58  }
59 
60  std::vector <double> flipAngles;
61  while (!inputFlip.eof())
62  {
63  char tmpStr[8192];
64  inputFlip.getline(tmpStr,8192);
65  std::string workStr(tmpStr);
66  if (workStr == "")
67  continue;
68 
69  workStr.erase(workStr.find_last_not_of(" \n\r\t")+1);
70  std::istringstream iss(workStr);
71  std::string shortStr;
72  iss >> shortStr;
73  flipAngles.push_back(std::stod(shortStr) * M_PI / 180.0);
74  }
75 
76  mainFilter->SetFlipAngles(flipAngles);
77  mainFilter->SetTRValue(trArg.getValue());
78  mainFilter->SetM0UpperBoundValue(upperBoundM0Arg.getValue());
79  mainFilter->SetT1UpperBoundValue(upperBoundT1Arg.getValue());
80 
81  if (b1Arg.getValue() != "")
82  {
83  InputImageReaderType::Pointer b1Read = InputImageReaderType::New();
84  b1Read->SetFileName(b1Arg.getValue().c_str());
85  b1Read->Update();
86 
87  mainFilter->SetB1Map(b1Read->GetOutput());
88  }
89 
90  if (maskArg.getValue() != "")
91  {
92  typedef itk::ImageFileReader < itk::Image <unsigned char, 3> > itkMaskReader;
93  itkMaskReader::Pointer maskRead = itkMaskReader::New();
94  maskRead->SetFileName(maskArg.getValue().c_str());
95  maskRead->Update();
96 
97  mainFilter->SetComputationMask(maskRead->GetOutput());
98  }
99 
100  mainFilter->SetAverageSignalThreshold(backgroundSignalThresholdArg.getValue());
101  mainFilter->SetNumberOfWorkUnits(nbpArg.getValue());
102 
103  itk::TimeProbe tmpTime;
104  tmpTime.Start();
105 
106  mainFilter->Update();
107 
108  tmpTime.Stop();
109 
110  std::cout << "Total computation time: " << tmpTime.GetTotal() << std::endl;
111 
112  OutputImageWriterType::Pointer resultT1Writer = OutputImageWriterType::New();
113  resultT1Writer->SetFileName(resT1Arg.getValue().c_str());
114  resultT1Writer->SetUseCompression(true);
115  resultT1Writer->SetInput(mainFilter->GetOutput(0));
116 
117  resultT1Writer->Update();
118 
119  if (resM0Arg.getValue() != "")
120  {
121  OutputImageWriterType::Pointer resultM0Writer = OutputImageWriterType::New();
122  resultM0Writer->SetFileName(resM0Arg.getValue().c_str());
123  resultM0Writer->SetUseCompression(true);
124  resultM0Writer->SetInput(mainFilter->GetOutput(1));
125 
126  resultM0Writer->Update();
127  }
128 
129  return EXIT_SUCCESS;
130 }
int main(int argc, char **argv)