tools/clitkVFInterpolateGenericFilter.cxx

   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 CLITKVFRESAMPLEGENERICFILTER_CXX
  19 #define CLITKVFRESAMPLEGENERICFILTER_CXX
  20
  21 #include "clitkVFInterpolateGenericFilter.h"
  22 #include "itkInterpolateImageFilter.h"
  23 #include "itkLinearInterpolateImageFunction.h"
  24 #include "itkNearestNeighborInterpolateImageFunction.h"
  25 #include "itkNthElementImageAdaptor.h"
  26
  27 //--------------------------------------------------------------------
  28 clitk::VFInterpolateGenericFilter::VFInterpolateGenericFilter():
  29   clitk::ImageToImageGenericFilter<Self>("VFInterpolate")
  30 {
  31   //InitializeImageType<2>();
  32   InitializeImageType<3>();
  33   //  InitializeImageType<4>();
  34   mInterpolatorName = "nn";
  35   mBSplineOrder=3;
  36   mDistance=0;
  37 }
  38 //--------------------------------------------------------------------
  39
  40
  41 //--------------------------------------------------------------------
  42 template<unsigned int Dim>
  43 void clitk::VFInterpolateGenericFilter::InitializeImageType()
  44 {
  45   //typedef itk::Vector<float,Dim> v3f;
  46   //ADD_IMAGE_TYPE(Dim, v3f);
  47 //   ADD_VEC_IMAGE_TYPE(Dim,2,double)
  48   ADD_VEC_IMAGE_TYPE(Dim,3,float)
  49   ADD_VEC_IMAGE_TYPE(Dim,3,double)
  50
  51 }
  52 //--------------------------------------------------------------------
  53
  54
  55 //--------------------------------------------------------------------
  56 template<class ImageType>
  57 void clitk::VFInterpolateGenericFilter::UpdateWithInputImageType()
  58 {
  59
  60   if (m_NbOfComponents == 1) {
  61     std::cerr << "Error, only one components ? Use clitkImageInterpolate instead." << std::endl;
  62     exit(0);
  63   }
  64   typedef typename ImageType::PixelType PixelType;
  65 //   if (m_NbOfComponents == 2) Update_WithDimAndPixelTypeAndComponent<ImageType::ImageDimension,PixelType,2>();
  66   if (m_NbOfComponents == 3) Update_WithDimAndPixelTypeAndComponent<ImageType::ImageDimension,PixelType,3>();
  67 }
  68 //--------------------------------------------------------------------
  69
  70
  71 //--------------------------------------------------------------------
  72 template<unsigned int Dim, class PixelType, unsigned int DimCompo>
  73 void clitk::VFInterpolateGenericFilter::Update_WithDimAndPixelTypeAndComponent()
  74 {
  75   // Reading input
  76   //  typedef itk::Vector<PixelType, DimCompo> DisplacementType;
  77   typedef PixelType DisplacementType;
  78   typedef itk::Image< DisplacementType, Dim > ImageType;
  79
  80   typename ImageType::Pointer input1 = clitk::readImage<ImageType>(m_InputFilenames[0], m_IOVerbose);
  81   typename ImageType::Pointer input2 = clitk::readImage<ImageType>(mInputFilename2, m_IOVerbose);
  82
  83   // Main filter
  84   typename ImageType::Pointer outputImage = ComputeImage<ImageType>(input1, input2);
  85
  86   // Write results
  87   SetNextOutput<ImageType>(outputImage);
  88 }
  89 //--------------------------------------------------------------------
  90
  91 //--------------------------------------------------------------------
  92 template<class ImageType>
  93 typename ImageType::Pointer
  94 clitk::VFInterpolateGenericFilter::ComputeImage(typename ImageType::Pointer inputImage1, typename ImageType::Pointer inputImage2)
  95 {
  96
  97   // Some typedefs
  98   typedef itk::Image<typename ImageType::PixelType::ValueType, ImageType::ImageDimension> ScalarImageType;
  99   typedef itk::Image<typename ImageType::PixelType::ValueType, ImageType::ImageDimension + 1> InterpolationImageType;
 100   typedef itk::NthElementImageAdaptor<ImageType, typename ImageType::PixelType::ValueType> ImageAdaptorType;
 101   typename ImageAdaptorType::Pointer adaptor1 = ImageAdaptorType::New();
 102   typename ImageAdaptorType::Pointer adaptor2 = ImageAdaptorType::New();
 103   adaptor1->SetImage(inputImage1);
 104   adaptor2->SetImage(inputImage2);
 105
 106   // Create Image Filter
 107   typedef itk::InterpolateImageFilter<ImageAdaptorType, ScalarImageType> FilterType;
 108   typename FilterType::Pointer filter = FilterType::New();
 109   filter->SetInput1(adaptor1);
 110   filter->SetInput2(adaptor2);
 111   filter->SetDistance(mDistance);
 112
 113   // Select interpolator
 114   if (mInterpolatorName == "nn") {
 115     typedef itk::NearestNeighborInterpolateImageFunction<InterpolationImageType> InterpolatorType;
 116     typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
 117     filter->SetInterpolator(interpolator);
 118   } else {
 119     if (mInterpolatorName == "linear") {
 120       typedef itk::LinearInterpolateImageFunction<InterpolationImageType> InterpolatorType;
 121       typename InterpolatorType::Pointer interpolator =  InterpolatorType::New();
 122       filter->SetInterpolator(interpolator);
 123     } else {
 124       std::cerr << "Sorry, I do not know the interpolator (for vector field) '" << mInterpolatorName
 125                 << "'. Known interpolators are :  nn, linear" << std::endl;
 126       exit(0);
 127     }
 128   }
 129
 130   typename ImageType::Pointer output = ImageType::New();
 131   typename ImageAdaptorType::Pointer adaptorOutput = ImageAdaptorType::New();
 132   output->CopyInformation(inputImage1);
 133   output->SetRegions(inputImage1->GetLargestPossibleRegion());
 134   output->Allocate();
 135
 136   typedef itk::ImageRegionIterator<ScalarImageType> IteratorType1;
 137   typedef itk::ImageRegionIterator<ImageAdaptorType> IteratorType2;
 138
 139   for (unsigned int i = 0; i < ImageType::PixelType::Dimension; i++) {
 140     adaptor1->SelectNthElement(i);
 141     adaptor2->SelectNthElement(i);
 142
 143     // Go !
 144     try {
 145       filter->Update();
 146     } catch( itk::ExceptionObject & err ) {
 147       std::cerr << "Error while filtering " << m_InputFilenames[0].c_str()
 148                 << " " << err << std::endl;
 149       exit(0);
 150     }
 151
 152     adaptorOutput->SelectNthElement(i);
 153     adaptorOutput->SetImage(output);
 154
 155     IteratorType1 it1(filter->GetOutput(), filter->GetOutput()->GetLargestPossibleRegion());
 156     IteratorType2 it2(adaptorOutput, adaptorOutput->GetLargestPossibleRegion());
 157
 158     it1.GoToBegin();
 159     it2.GoToBegin();
 160     while ( ! it1.IsAtEnd() ) {
 161       it2.Set(it1.Get());
 162       ++it1;
 163       ++it2;
 164     }
 165   }
 166
 167   // Return result
 168   return output;
 169
 170 }
 171 //--------------------------------------------------------------------
 172
 173 //--------------------------------------------------------------------
 174 void clitk::VFInterpolateGenericFilter::SetInterpolationName(const std::string & inter)
 175 {
 176   mInterpolatorName = inter;
 177 }
 178 //--------------------------------------------------------------------
 179
 180 #endif
 181