ANIMA  4.0
animaDTIEstimator.cxx
Go to the documentation of this file.
3 
4 #include <tclap/CmdLine.h>
5 #include <itkTimeProbe.h>
6 #include <fstream>
7 
9 #include <animaReorientation.h>
10 
11 //Update progression of the process
12 void eventCallback (itk::Object* caller, const itk::EventObject& event, void* clientData)
13 {
14  itk::ProcessObject * processObject = static_cast<itk::ProcessObject *> (caller);
15  std::cout<<"\033[K\rProgression: " << static_cast<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  TCLAP::ValueArg<std::string> inArg("i","inputdwi","dwi_volume",true,"","DWI volume",cmd);
23  TCLAP::ValueArg<std::string> resArg("o","output","dit_volume",true,"","result DTI image",cmd);
24  TCLAP::ValueArg<std::string> b0OutArg("O","output-b0","output_b0",false,"","result B0 image",cmd);
25  TCLAP::ValueArg<std::string> varOutArg("N","output-variance","output_variance",false,"","result noise variance image",cmd);
26 
27  TCLAP::ValueArg<std::string> gradsArg("g","grad","input_gradients",true,"","Input gradients",cmd);
28  TCLAP::ValueArg<std::string> bvalArg("b","bval","input_b-values",true,"","Input b-values",cmd);
29  TCLAP::SwitchArg bvalueScaleArg("B","b-no-scale","Do not scale b-values according to gradient norm",cmd);
30  TCLAP::ValueArg<std::string> computationMaskArg("m","mask","Computation mask", false,"","computation mask",cmd);
31 
32  TCLAP::ValueArg<unsigned int> b0ThrArg("t","b0thr","bot_treshold",false,0,"B0 threshold (default : 0)",cmd);
33  TCLAP::ValueArg<unsigned int> nbpArg("p","numberofthreads","nb_thread",false,itk::MultiThreaderBase::GetGlobalDefaultNumberOfThreads(),"Number of threads to run on (default: all cores)",cmd);
34  TCLAP::ValueArg<std::string> reorientArg("r","reorient","dwi_reoriented",false,"","Reorient DWI given as input",cmd);
35  TCLAP::ValueArg<std::string> reorientGradArg("R","reorient-G","gradient reoriented output",false,"","Reorient gradients so that they are in MrTrix format (in image coordinates)",cmd);
36 
37  try
38  {
39  cmd.parse(argc,argv);
40  }
41  catch (TCLAP::ArgException& e)
42  {
43  std::cerr << "Error: " << e.error() << "for argument " << e.argId() << std::endl;
44  return EXIT_FAILURE;
45  }
46 
48  typedef FilterType::MaskImageType MaskImageType;
49  typedef FilterType::InputImageType InputImageType;
50  typedef FilterType::Image4DType Image4DType;
51  typedef FilterType::OutputImageType VectorImageType;
52 
53  itk::CStyleCommand::Pointer callback = itk::CStyleCommand::New();
54  callback->SetCallback(eventCallback);
55 
56  FilterType::Pointer mainFilter = FilterType::New();
57 
58  typedef anima::GradientFileReader < vnl_vector_fixed<double,3>, double > GFReaderType;
59  GFReaderType gfReader;
60  gfReader.SetGradientFileName(gradsArg.getValue());
61  gfReader.SetBValueBaseString(bvalArg.getValue());
62  gfReader.SetGradientIndependentNormalization(bvalueScaleArg.isSet());
63 
64  gfReader.Update();
65 
66  GFReaderType::GradientVectorType directions = gfReader.GetGradients();
67  GFReaderType::BValueVectorType mb = gfReader.GetBValues();
68 
69  std::string inputFile = inArg.getValue();
70 
71  Image4DType::Pointer input = anima::readImage<Image4DType>(inArg.getValue());
72 
73  if (reorientArg.getValue() != "")
74  {
75  input = anima::reorientImage<Image4DType>(input, itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RAI);
76  anima::writeImage<Image4DType>(reorientArg.getValue(), input);
77 
78  inputFile = reorientArg.getValue();
79  }
80 
81  if(reorientGradArg.getValue() != "")
82  {
83  anima::reorientGradients <Image4DType, vnl_vector_fixed<double,3> > (input, directions);
84 
85  std::ofstream outGrads(reorientGradArg.getValue().c_str());
86  outGrads.precision(15);
87  for (unsigned int i = 0;i < 3;++i)
88  {
89  for (unsigned int j = 0;j < directions.size();++j)
90  outGrads << directions[j][i] << " ";
91  outGrads << std::endl;
92  }
93 
94  outGrads.close();
95  }
96 
97  mainFilter->SetBValuesList(mb);
98  for(unsigned int i = 0;i < directions.size();++i)
99  mainFilter->AddGradientDirection(i, directions[i]);
100 
101  std::vector <InputImageType::Pointer> inputData;
102  inputData = anima::getImagesFromHigherDimensionImage<Image4DType,InputImageType>(input);
103 
104  for (unsigned int i = 0;i < inputData.size();++i)
105  mainFilter->SetInput(i,inputData[i]);
106 
107  if (computationMaskArg.getValue() != "")
108  mainFilter->SetComputationMask(anima::readImage<MaskImageType>(computationMaskArg.getValue()));
109 
110  mainFilter->SetB0Threshold(b0ThrArg.getValue());
111  mainFilter->SetNumberOfWorkUnits(nbpArg.getValue());
112  mainFilter->AddObserver(itk::ProgressEvent(), callback);
113 
114  itk::TimeProbe tmpTimer;
115 
116  tmpTimer.Start();
117 
118  try
119  {
120  mainFilter->Update();
121  } catch (itk::ExceptionObject &e)
122  {
123  std::cerr << e << std::endl;
124  return EXIT_FAILURE;
125  }
126 
127  tmpTimer.Stop();
128 
129  std::cout << "\nEstimation done in " << tmpTimer.GetTotal() << " s" << std::endl;
130  std::cout << "Writing result to : " << resArg.getValue() << std::endl;
131 
132  VectorImageType::Pointer output = mainFilter->GetOutput();
133  output->DisconnectPipeline();
134 
135  anima::writeImage <VectorImageType> (resArg.getValue(),output);
136 
137  if (b0OutArg.getValue() != "")
138  anima::writeImage <InputImageType> (b0OutArg.getValue(),mainFilter->GetEstimatedB0Image());
139 
140  if (varOutArg.getValue() != "")
141  anima::writeImage <InputImageType> (varOutArg.getValue(),mainFilter->GetEstimatedVarianceImage());
142 
143  return EXIT_SUCCESS;
144 }
int main(int argc, char **argv)
void SetGradientFileName(std::string fName)
void eventCallback(itk::Object *caller, const itk::EventObject &event, void *clientData)