ANIMA  4.0
animaMCMApplyTransformSerie.cxx
Go to the documentation of this file.
1 #include <tclap/CmdLine.h>
2 
3 #include <animaMCMFileReader.h>
4 #include <animaMCMFileWriter.h>
6 #include <itkImageFileReader.h>
7 
9 #include <itkNearestNeighborInterpolateImageFunction.h>
11 
13 
14 int main(int ac, const char** av)
15 {
16  std::string descriptionMessage;
17  descriptionMessage += "Resampler tool to apply a series of transformations to a multi-tensor image. ";
18  descriptionMessage += "Input transform is an XML file describing all transforms to apply. ";
19  descriptionMessage += "Such a file should look like this:\n";
20  descriptionMessage += "<TransformationList>\n";
21  descriptionMessage += "<Transformation>\n";
22  descriptionMessage += "<Type>linear</Type> (it can be svf or dense too)\n";
23  descriptionMessage += "<Path>FileName</Path>\n";
24  descriptionMessage += "<Inversion>0</Inversion>\n";
25  descriptionMessage += "</Transformation>\n";
26  descriptionMessage += "...\n";
27  descriptionMessage += "</TransformationList>\n\n";
28  descriptionMessage += "INRIA / IRISA - VisAGeS/Empenn Team";
29 
30  TCLAP::CmdLine cmd(descriptionMessage, ' ',ANIMA_VERSION);
31 
32  TCLAP::ValueArg<std::string> inArg("i","input","Input image",true,"","input image",cmd);
33  TCLAP::ValueArg<std::string> trArg("t","trsf","Transformations XML list",true,"","transformations list",cmd);
34  TCLAP::ValueArg<std::string> outArg("o","output","Output resampled image",true,"","output image",cmd);
35  TCLAP::ValueArg<std::string> geomArg("g","geometry","Geometry image",true,"","geometry image",cmd);
36 
37  TCLAP::ValueArg<int> numFasciclesOut("n","num-fasc","Number of MCM output fascicles (default: same as input)",false,-1,"number of MCM output fascicles",cmd);
38  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);
39 
40  TCLAP::SwitchArg ppdArg("P","ppd","Use PPD re-orientation scheme (default: no)",cmd,false);
41  TCLAP::SwitchArg invertArg("I","invert","Invert the transformation series",cmd,false);
42  TCLAP::SwitchArg nearestArg("N","nearest","Use nearest neighbor interpolation",cmd,false);
43 
44  TCLAP::ValueArg<unsigned int> nbpArg("p","numberofthreads","Number of threads to run on (default: all cores)",false,itk::MultiThreaderBase::GetGlobalDefaultNumberOfThreads(),"number of threads",cmd);
45 
46  try
47  {
48  cmd.parse(ac,av);
49  }
50  catch (TCLAP::ArgException& e)
51  {
52  std::cerr << "Error: " << e.error() << "for argument " << e.argId() << std::endl;
53  return EXIT_FAILURE;
54  }
55 
56  const unsigned int Dimension = 3;
57  typedef double PixelType;
58 
60  typedef anima::TransformSeriesReader <double, Dimension> TransformSeriesReaderType;
61  typedef TransformSeriesReaderType::OutputTransformType TransformType;
62 
63  typedef itk::ImageFileReader <ImageType> ReaderType;
64  typedef anima::MCMFileReader <double,3> MCMReaderType;
65  typedef anima::MCMFileWriter <double, 3> MCMWriterType;
66 
67  MCMReaderType reader;
68  reader.SetFileName(inArg.getValue());
69  reader.Update();
70 
71  itk::ImageIOBase::Pointer imageIO = itk::ImageIOFactory::CreateImageIO(geomArg.getValue().c_str(),
72  itk::ImageIOFactory::ReadMode);
73 
74  if( !imageIO )
75  {
76  std::cerr << "Itk could not find suitable IO factory for the input" << std::endl;
77  return EXIT_FAILURE;
78  }
79 
80  // Now that we found the appropriate ImageIO class, ask it to read the meta data from the image file.
81  imageIO->SetFileName(geomArg.getValue());
82  imageIO->ReadImageInformation();
83 
84  TransformSeriesReaderType *trReader = new TransformSeriesReaderType;
85  trReader->SetInput(trArg.getValue());
86  trReader->SetExponentiationOrder(expOrderArg.getValue());
87  trReader->SetNumberOfWorkUnits(nbpArg.getValue());
88  trReader->SetInvertTransform(invertArg.isSet());
89 
90  try
91  {
92  trReader->Update();
93  }
94  catch (itk::ExceptionObject &e)
95  {
96  std::cerr << e << std::endl;
97  return EXIT_FAILURE;
98  }
99 
100  TransformType::Pointer trsf = trReader->GetOutputTransform();
101 
102  itk::InterpolateImageFunction <ImageType>::Pointer interpolator;
103  unsigned int numIsotropicCompartments = reader.GetModelVectorImage()->GetDescriptionModel()->GetNumberOfIsotropicCompartments();
104  unsigned int numFasciclesInput = reader.GetModelVectorImage()->GetDescriptionModel()->GetNumberOfCompartments() - numIsotropicCompartments;
105  unsigned int numFasciclesOutput = numFasciclesInput;
106 
107  if (numFasciclesOut.getValue() >= 0)
108  numFasciclesOutput = std::min((unsigned int)numFasciclesOut.getValue(),numFasciclesInput);
109 
111 
112  mcmCreator.SetModelWithFreeWaterComponent(false);
113  mcmCreator.SetModelWithRestrictedWaterComponent(false);
114  mcmCreator.SetModelWithStaniszComponent(false);
115  mcmCreator.SetModelWithStationaryWaterComponent(false);
116 
117  for (unsigned int i = 0;i < numIsotropicCompartments;++i)
118  {
119  switch (reader.GetModelVectorImage()->GetDescriptionModel()->GetCompartment(i)->GetCompartmentType())
120  {
121  case anima::FreeWater:
122  mcmCreator.SetModelWithFreeWaterComponent(true);
123  break;
124 
126  mcmCreator.SetModelWithRestrictedWaterComponent(true);
127  break;
128 
129  case anima::Stanisz:
130  mcmCreator.SetModelWithStaniszComponent(true);
131  break;
132 
134  default:
135  mcmCreator.SetModelWithStationaryWaterComponent(true);
136  break;
137  }
138  }
139 
140  if (numFasciclesInput > 0)
141  mcmCreator.SetCompartmentType(reader.GetModelVectorImage()->GetDescriptionModel()->GetCompartment(numIsotropicCompartments)->GetCompartmentType());
142  mcmCreator.SetNumberOfCompartments(numFasciclesOutput);
143 
144  anima::MultiCompartmentModel::Pointer outputReferenceModel = mcmCreator.GetNewMultiCompartmentModel();
145 
146  if (nearestArg.isSet())
147  interpolator = itk::NearestNeighborInterpolateImageFunction<ImageType>::New();
148  else
149  {
151  tmpInterpolator->SetReferenceOutputModel(outputReferenceModel);
152  interpolator = tmpInterpolator;
153  }
154 
155  typedef anima::MCMResampleImageFilter <ImageType, double> ResampleFilterType;
156  ResampleFilterType::Pointer resample = ResampleFilterType::New();
157 
158  resample->SetTransform(trsf);
159  resample->SetFiniteStrainReorientation(!ppdArg.isSet());
160  if (nearestArg.isSet())
161  resample->SetReferenceOutputModel(reader.GetModelVectorImage()->GetDescriptionModel());
162  else
163  resample->SetReferenceOutputModel(outputReferenceModel);
164 
165  resample->SetInterpolator(interpolator.GetPointer());
166  resample->SetNumberOfWorkUnits(nbpArg.getValue());
167 
168  ImageType::DirectionType directionMatrix;
169  ImageType::PointType origin;
170  ImageType::SpacingType spacing;
171  ImageType::RegionType largestRegion;
172 
173  for (unsigned int i = 0;i < Dimension;++i)
174  {
175  origin[i] = imageIO->GetOrigin(i);
176  spacing[i] = imageIO->GetSpacing(i);
177  largestRegion.SetIndex(i,0);
178  largestRegion.SetSize(i,imageIO->GetDimensions(i));
179 
180  for (unsigned int j = 0;j < Dimension;++j)
181  directionMatrix(i,j) = imageIO->GetDirection(j)[i];
182  }
183 
184  resample->SetOutputLargestPossibleRegion(largestRegion);
185  resample->SetOutputOrigin(origin);
186  resample->SetOutputSpacing(spacing);
187  resample->SetOutputDirection(directionMatrix);
188 
189  resample->SetInput(reader.GetModelVectorImage());
190 
191  std::cout << "Applying transform... " << std::flush;
192  resample->Update();
193  std::cout << "Done..." << std::endl;
194 
195  MCMWriterType mcmWriter;
196  mcmWriter.SetInputImage(resample->GetOutput());
197 
198  mcmWriter.SetFileName(outArg.getValue());
199  mcmWriter.Update();
200 
201  return EXIT_SUCCESS;
202 }
int main(int ac, const char **av)
void SetFileName(std::string fileName)
Really this class is some simplified factory that creates the MCM that it knows.