ANIMA  4.0
animaODFApplyTransformSerie.cxx
Go to the documentation of this file.
1 #include <tclap/CmdLine.h>
2 
3 #include <itkImageFileReader.h>
4 #include <itkImageFileWriter.h>
5 
7 #include <itkNearestNeighborInterpolateImageFunction.h>
8 
10 
11 int main(int ac, const char** av)
12 {
13  std::string descriptionMessage;
14  descriptionMessage += "Resampler tool to apply a series of transformations to an image of ODF models in Max Descoteaux's basis. ";
15  descriptionMessage += "Input transform is an XML file describing all transforms to apply. ";
16  descriptionMessage += "Such a file should look like this:\n";
17  descriptionMessage += "<TransformationList>\n";
18  descriptionMessage += "<Transformation>\n";
19  descriptionMessage += "<Type>linear</Type> (it can be svf or dense too)\n";
20  descriptionMessage += "<Path>FileName</Path>\n";
21  descriptionMessage += "<Inversion>0</Inversion>\n";
22  descriptionMessage += "</Transformation>\n";
23  descriptionMessage += "...\n";
24  descriptionMessage += "</TransformationList>\n\n";
25  descriptionMessage += "INRIA / IRISA - VisAGeS/Empenn Team";
26 
27  TCLAP::CmdLine cmd(descriptionMessage, ' ',ANIMA_VERSION);
28 
29  TCLAP::ValueArg<std::string> inArg("i","input","Input image",true,"","input image",cmd);
30  TCLAP::ValueArg<std::string> trArg("t","trsf","Transformations XML list",true,"","transformations list",cmd);
31  TCLAP::ValueArg<std::string> outArg("o","output","Output resampled image",true,"","output image",cmd);
32  TCLAP::ValueArg<std::string> geomArg("g","geometry","Geometry image",true,"","geometry image",cmd);
33 
34  TCLAP::ValueArg<unsigned int> expOrderArg("e","exp-order","Order of field exponentiation approximation (in between 0 and 1, default: 0)",false,0,"exponentiation order",cmd);
35  TCLAP::SwitchArg invertArg("I","invert","Invert the transformation series",cmd,false);
36  TCLAP::SwitchArg nearestArg("N","nearest","Use nearest neighbor interpolation",cmd,false);
37 
38  TCLAP::ValueArg<unsigned int> nbpArg("p","numberofthreads","Number of threads to run on (default: all cores)",false,itk::MultiThreaderBase::GetGlobalDefaultNumberOfThreads(),"number of threads",cmd);
39 
40  try
41  {
42  cmd.parse(ac,av);
43  }
44  catch (TCLAP::ArgException& e)
45  {
46  std::cerr << "Error: " << e.error() << "for argument " << e.argId() << std::endl;
47  return EXIT_FAILURE;
48  }
49 
50  const unsigned int Dimension = 3;
51  typedef double PixelType;
52 
53  typedef itk::VectorImage< PixelType, Dimension > ImageType;
54  typedef anima::TransformSeriesReader <double, Dimension> TransformSeriesReaderType;
55  typedef TransformSeriesReaderType::OutputTransformType TransformType;
56 
57  typedef itk::ImageFileReader <ImageType> ReaderType;
58  typedef itk::ImageFileWriter <ImageType> WriterType;
59 
60  ReaderType::Pointer reader = ReaderType::New();
61  reader->SetFileName(inArg.getValue());
62  reader->Update();
63 
64  itk::ImageIOBase::Pointer imageIO = itk::ImageIOFactory::CreateImageIO(geomArg.getValue().c_str(),
65  itk::ImageIOFactory::ReadMode);
66 
67  if (!imageIO)
68  {
69  std::cerr << "Itk could not find suitable IO factory for the input" << std::endl;
70  return EXIT_FAILURE;
71  }
72 
73  // Now that we found the appropriate ImageIO class, ask it to read the meta data from the image file.
74  imageIO->SetFileName(geomArg.getValue());
75  imageIO->ReadImageInformation();
76 
77  TransformSeriesReaderType *trReader = new TransformSeriesReaderType;
78  trReader->SetInput(trArg.getValue());
79  trReader->SetInvertTransform(invertArg.isSet());
80  trReader->SetExponentiationOrder(expOrderArg.getValue());
81  trReader->SetNumberOfWorkUnits(nbpArg.getValue());
82 
83  try
84  {
85  trReader->Update();
86  }
87  catch (itk::ExceptionObject &e)
88  {
89  std::cerr << e << std::endl;
90  return EXIT_FAILURE;
91  }
92 
93  TransformType::Pointer trsf = trReader->GetOutputTransform();
94 
95  itk::InterpolateImageFunction <ImageType>::Pointer interpolator;
96 
97  if (nearestArg.isSet())
98  interpolator = itk::NearestNeighborInterpolateImageFunction<ImageType>::New();
99  else
101 
102  typedef anima::ODFResampleImageFilter <ImageType, double> ResampleFilterType;
103 
104  ResampleFilterType::Pointer resample = ResampleFilterType::New();
105 
106  resample->SetTransform(trsf);
107  resample->SetFiniteStrainReorientation(true);
108  resample->SetInterpolator(interpolator.GetPointer());
109  resample->SetNumberOfWorkUnits(nbpArg.getValue());
110 
111  ImageType::DirectionType directionMatrix;
112  ImageType::PointType origin;
113  ImageType::SpacingType spacing;
114  ImageType::RegionType largestRegion;
115 
116  for (unsigned int i = 0;i < Dimension;++i)
117  {
118  origin[i] = imageIO->GetOrigin(i);
119  spacing[i] = imageIO->GetSpacing(i);
120  largestRegion.SetIndex(i,0);
121  largestRegion.SetSize(i,imageIO->GetDimensions(i));
122 
123  for (unsigned int j = 0;j < Dimension;++j)
124  directionMatrix(i,j) = imageIO->GetDirection(j)[i];
125  }
126 
127  resample->SetOutputLargestPossibleRegion(largestRegion);
128  resample->SetOutputOrigin(origin);
129  resample->SetOutputSpacing(spacing);
130  resample->SetOutputDirection(directionMatrix);
131 
132  resample->SetInput(reader->GetOutput());
133 
134  std::cout << "Applying transform... " << std::flush;
135  resample->Update();
136  std::cout << "Done..." << std::endl;
137 
138  WriterType::Pointer writer = WriterType::New();
139 
140  writer->SetUseCompression(true);
141 
142  writer->SetInput(resample->GetOutput());
143 
144  writer->SetFileName(outArg.getValue());
145 
146  writer->Update();
147 
148  return EXIT_SUCCESS;
149 }
int main(int ac, const char **av)