tools/clitkInvertVFGenericFilter.txx

   1 /*=========================================================================
   2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
   3
   4   Authors belong to:
   5   - University of LYON              http://www.universite-lyon.fr/
   6   - Léon Bérard cancer center       http://www.centreleonberard.fr
   7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
   8
   9   This software is distributed WITHOUT ANY WARRANTY; without even
  10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  11   PURPOSE.  See the copyright notices for more information.
  12
  13   It is distributed under dual licence
  14
  15   - BSD        See included LICENSE.txt file
  16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  17 ===========================================================================**/
  18 #ifndef clitkInvertVFGenericFilter_txx
  19 #define clitkInvertVFGenericFilter_txx
  20
  21 #include "itkVectorResampleImageFilter.h"
  22 #include "clitkBSplineDeformableTransform.h"
  23 #if ITK_VERSION_MAJOR >= 4
  24 #include "itkTransformToDisplacementFieldSource.h"
  25 #else
  26 #include "itkTransformToDeformationFieldSource.h"
  27 #endif
  28
  29 /* =================================================
  30  * @file   clitkInvertVFGenericFilter.txx
  31  * @author
  32  * @date
  33  *
  34  * @brief
  35  *
  36  ===================================================*/
  37
  38
  39 namespace clitk
  40 {
  41
  42 //-----------------------------------------------------------
  43 // Constructor
  44 //-----------------------------------------------------------
  45 template<class args_info_type>
  46 InvertVFGenericFilter<args_info_type>::InvertVFGenericFilter()
  47 {
  48   m_Verbose=false;
  49   m_InputFileName="";
  50 }
  51
  52
  53 //-----------------------------------------------------------
  54 // Update
  55 //-----------------------------------------------------------
  56 template<class args_info_type>
  57 void InvertVFGenericFilter<args_info_type>::Update()
  58 {
  59   // Read the Dimension and PixelType
  60   int Dimension;
  61   std::string PixelType;
  62   ReadImageDimensionAndPixelType(m_InputFileName, Dimension, PixelType);
  63
  64
  65   // Call UpdateWithDim
  66   if(Dimension==2) UpdateWithDim<2>(PixelType);
  67   else if(Dimension==3) UpdateWithDim<3>(PixelType);
  68   // else if (Dimension==4)UpdateWithDim<4>(PixelType);
  69   else {
  70     std::cout<<"Error, Only for 2 or 3  Dimensions!!!"<<std::endl ;
  71     return;
  72   }
  73 }
  74
  75 //-------------------------------------------------------------------
  76 // Update with the number of dimensions
  77 //-------------------------------------------------------------------
  78 template<class args_info_type>
  79 template<unsigned int Dimension>
  80 void
  81 InvertVFGenericFilter<args_info_type>::UpdateWithDim(std::string PixelType)
  82 {
  83   if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D and "<< PixelType<<"..."<<std::endl;
  84
  85   //    if(PixelType == "short"){
  86   //       if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed short..." << std::endl;
  87   //       UpdateWithDimAndPixelType<Dimension, signed short>();
  88   //     }
  89   //    else if(PixelType == "unsigned_short"){
  90   //       if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and unsigned_short..." << std::endl;
  91   //       UpdateWithDimAndPixelType<Dimension, unsigned short>();
  92   //     }
  93
  94   //     else if (PixelType == "unsigned_char"){
  95   //       if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl;
  96   //       UpdateWithDimAndPixelType<Dimension, unsigned char>();
  97   //     }
  98
  99   //     else if (PixelType == "char"){
 100   //       if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and signed_char..." << std::endl;
 101   //       UpdateWithDimAndPixelType<Dimension, signed char>();
 102   //     }
 103   //  else {
 104   if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
 105   UpdateWithDimAndPixelType<Dimension, itk::Vector<float, Dimension> >();
 106   // }
 107 }
 108
 109 //-------------------------------------------------------------------
 110 // Convert Coefficient image to DVF
 111 //-------------------------------------------------------------------
 112 template<class args_info_type>
 113 template<class DisplacementFieldType>
 114 typename DisplacementFieldType::Pointer
 115 InvertVFGenericFilter<args_info_type>::CoeffsToDVF(std::string fileName, std::string likeFileName)
 116 {
 117   typedef clitk::BSplineDeformableTransform<double, DisplacementFieldType::ImageDimension, DisplacementFieldType::ImageDimension> TransformType;
 118   typedef typename TransformType::CoefficientImageType CoefficientImageType;
 119
 120   typedef itk::ImageFileReader<CoefficientImageType> CoeffReaderType;
 121   typename CoeffReaderType::Pointer reader = CoeffReaderType::New();
 122   reader->SetFileName(fileName);
 123   reader->Update();
 124
 125   typename TransformType::Pointer transform = TransformType::New();
 126   transform->SetCoefficientImage(reader->GetOutput());
 127
 128 #if ITK_VERSION_MAJOR >= 4
 129       typedef itk::TransformToDisplacementFieldSource<DisplacementFieldType, double> ConvertorType;
 130 #else
 131       typedef itk::TransformToDeformationFieldSource<DisplacementFieldType, double> ConvertorType;
 132 #endif
 133
 134   typedef itk::ImageIOBase ImageIOType;
 135   typename ImageIOType::Pointer imageIO = itk::ImageIOFactory::CreateImageIO(likeFileName.c_str(), itk::ImageIOFactory::ReadMode);
 136   imageIO->SetFileName(likeFileName);
 137   imageIO->ReadImageInformation();
 138
 139   typename ConvertorType::Pointer convertor= ConvertorType::New();
 140   typename ConvertorType::SizeType output_size;
 141   typename ConvertorType::SpacingType output_spacing;
 142   typename ConvertorType::OriginType output_origin;
 143   typename ConvertorType::DirectionType output_direction;
 144   for (unsigned int i = 0; i < DisplacementFieldType::ImageDimension; i++) {
 145     output_size[i] = imageIO->GetDimensions(i);
 146     output_spacing[i] = imageIO->GetSpacing(i);
 147     output_origin[i] = imageIO->GetOrigin(i);
 148     for (unsigned int j = 0; j < DisplacementFieldType::ImageDimension; j++)
 149       output_direction[i][j] = imageIO->GetDirection(i)[j];
 150   }
 151
 152   if (m_Verbose) {
 153     std::cout << "Interpolating coefficients with grid:" << std::endl;
 154     std::cout << output_size << output_spacing << std::endl;
 155   }
 156
 157   convertor->SetNumberOfThreads(1);
 158   convertor->SetTransform(transform);
 159   convertor->SetOutputOrigin(output_origin);
 160   convertor->SetOutputSpacing(output_spacing);
 161   convertor->SetOutputSize(output_size);
 162   convertor->SetOutputDirection(output_direction);
 163   convertor->Update();
 164
 165   return convertor->GetOutput();
 166 }
 167
 168
 169 //-------------------------------------------------------------------
 170 // Update with the number of dimensions and the pixeltype
 171 //-------------------------------------------------------------------
 172 template<class args_info_type>
 173 template <unsigned int Dimension, class  PixelType>
 174 void
 175 InvertVFGenericFilter<args_info_type>::UpdateWithDimAndPixelType()
 176 {
 177
 178   // ImageTypes
 179   typedef itk::Image<PixelType, Dimension> InputImageType;
 180   typedef itk::Image<PixelType, Dimension> OutputImageType;
 181
 182   // Read the input
 183   typedef itk::ImageFileReader<InputImageType> InputReaderType;
 184   typename InputReaderType::Pointer reader = InputReaderType::New();
 185   reader->SetFileName( m_InputFileName);
 186   reader->Update();
 187   typename InputImageType::Pointer input = reader->GetOutput();
 188
 189   // Filter
 190   typename OutputImageType::Pointer output;
 191   switch (m_ArgsInfo.type_arg) {
 192
 193     // clitk filter
 194   case 0: {
 195     // Create the InvertVFFilter
 196     typedef clitk::InvertVFFilter<InputImageType,OutputImageType> FilterType;
 197     typename FilterType::Pointer filter =FilterType::New();
 198     if (m_ArgsInfo.like_given) {
 199       typename FilterType::SpacingType spacing;
 200       typename FilterType::SizeType size;
 201       itk::ImageIOBase::Pointer header = readImageHeader(m_ArgsInfo.like_arg);
 202       for(unsigned int i=0; i<InputImageType::ImageDimension; i++) {
 203         size[i] = header->GetDimensions(i);
 204         spacing[i] = header->GetSpacing(i);
 205       }
 206
 207       typedef itk::VectorResampleImageFilter<InputImageType, OutputImageType> ResampleFilterType;
 208       typename ResampleFilterType::Pointer resampler = ResampleFilterType::New();
 209       resampler->SetInput(input);
 210       resampler->SetOutputOrigin(input->GetOrigin());
 211       resampler->SetOutputDirection(input->GetDirection());
 212       resampler->SetOutputSpacing(spacing);
 213       resampler->SetSize(size);
 214       resampler->Update();
 215       spacing = resampler->GetOutput()->GetSpacing();
 216       size = resampler->GetOutput()->GetLargestPossibleRegion().GetSize();
 217       filter->SetInput(resampler->GetOutput());
 218     }
 219     else
 220       filter->SetInput(input);
 221
 222     filter->SetVerbose(m_Verbose);
 223     if (m_ArgsInfo.threads_given) filter->SetNumberOfThreads(m_ArgsInfo.threads_arg);
 224     if (m_ArgsInfo.pad_given) {
 225       PixelType pad;
 226       if (m_ArgsInfo.pad_given !=  (pad.GetNumberOfComponents()) )
 227         pad.Fill(m_ArgsInfo.pad_arg[0]);
 228       else
 229         for(unsigned int i=0; i<Dimension; i++)
 230           pad[i]=m_ArgsInfo.pad_arg[i];
 231     }
 232     filter->SetThreadSafe(m_ArgsInfo.threadSafe_flag);
 233     filter->Update();
 234     output=filter->GetOutput();
 235
 236     break;
 237   }
 238
 239   case 1: {
 240     // Create the InvertVFFilter
 241     typedef clitk::InvertVFFilter<InputImageType,OutputImageType> FilterType;
 242     typename FilterType::Pointer filter =FilterType::New();
 243     if (m_ArgsInfo.like_given) {
 244       filter->SetInput(CoeffsToDVF<OutputImageType>(m_InputFileName, m_ArgsInfo.like_arg));
 245     }
 246
 247     filter->SetVerbose(m_Verbose);
 248     if (m_ArgsInfo.threads_given) filter->SetNumberOfThreads(m_ArgsInfo.threads_arg);
 249     if (m_ArgsInfo.pad_given) {
 250       PixelType pad;
 251       if (m_ArgsInfo.pad_given !=  (pad.GetNumberOfComponents()) )
 252         pad.Fill(m_ArgsInfo.pad_arg[0]);
 253       else
 254         for(unsigned int i=0; i<Dimension; i++)
 255           pad[i]=m_ArgsInfo.pad_arg[i];
 256     }
 257     filter->SetThreadSafe(m_ArgsInfo.threadSafe_flag);
 258     filter->Update();
 259     output=filter->GetOutput();
 260
 261     break;
 262   }
 263
 264   case 2: {
 265     // Create the InverseDeformationFieldFilter
 266 #if ITK_VERSION_MAJOR >= 4
 267     typedef itk::InverseDisplacementFieldImageFilter<InputImageType,OutputImageType> FilterType;
 268 #else
 269     typedef itk::InverseDeformationFieldImageFilter<InputImageType,OutputImageType> FilterType;
 270 #endif
 271     typename FilterType::Pointer filter =FilterType::New();
 272     filter->SetInput(input);
 273     filter->SetOutputOrigin(input->GetOrigin());
 274     filter->SetOutputSpacing(input->GetSpacing());
 275     filter->SetSize(input->GetLargestPossibleRegion().GetSize());
 276     filter->SetSubsamplingFactor(m_ArgsInfo.sampling_arg);
 277     filter->Update();
 278     output=filter->GetOutput();
 279
 280     break;
 281   }
 282
 283
 284   }
 285
 286   // Output
 287   typedef itk::ImageFileWriter<OutputImageType> WriterType;
 288   typename WriterType::Pointer writer = WriterType::New();
 289   writer->SetFileName(m_ArgsInfo.output_arg);
 290   writer->SetInput(output);
 291   writer->Update();
 292
 293 }
 294
 295
 296 }//end clitk
 297
 298 #endif //#define clitkInvertVFGenericFilter_txx