added tool to interpolate VFs

[clitk.git] / tools / clitkVFInterpolateGenericFilter.cxx

diff --git a/tools/clitkVFInterpolateGenericFilter.cxx b/tools/clitkVFInterpolateGenericFilter.cxx
new file mode 100644 (file)
index 0000000..777bb13
--- /dev/null
+++ b/tools/clitkVFInterpolateGenericFilter.cxx
@@ -0,0 +1,181 @@
+/*=========================================================================
+  Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
+
+  Authors belong to:
+  - University of LYON              http://www.universite-lyon.fr/
+  - Léon Bérard cancer center       http://www.centreleonberard.fr
+  - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
+
+  This software is distributed WITHOUT ANY WARRANTY; without even
+  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+  PURPOSE.  See the copyright notices for more information.
+
+  It is distributed under dual licence
+
+  - BSD        See included LICENSE.txt file
+  - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
+===========================================================================**/
+#ifndef CLITKVFRESAMPLEGENERICFILTER_CXX
+#define CLITKVFRESAMPLEGENERICFILTER_CXX
+
+#include "clitkVFInterpolateGenericFilter.h"
+#include "itkInterpolateImageFilter.h"
+#include "itkLinearInterpolateImageFunction.h"
+#include "itkNearestNeighborInterpolateImageFunction.h"
+#include "itkNthElementImageAdaptor.h"
+
+//--------------------------------------------------------------------
+clitk::VFInterpolateGenericFilter::VFInterpolateGenericFilter():
+  clitk::ImageToImageGenericFilter<Self>("VFInterpolate")
+{
+  //InitializeImageType<2>();
+  InitializeImageType<3>();
+  //  InitializeImageType<4>();
+  mInterpolatorName = "nn";
+  mBSplineOrder=3;
+  mDistance=0;
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template<unsigned int Dim>
+void clitk::VFInterpolateGenericFilter::InitializeImageType()
+{
+  //typedef itk::Vector<float,Dim> v3f;
+  //ADD_IMAGE_TYPE(Dim, v3f);
+//   ADD_VEC_IMAGE_TYPE(Dim,2,double)
+  ADD_VEC_IMAGE_TYPE(Dim,3,float)
+  ADD_VEC_IMAGE_TYPE(Dim,3,double)
+  
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template<class ImageType>
+void clitk::VFInterpolateGenericFilter::UpdateWithInputImageType()
+{
+
+  if (m_NbOfComponents == 1) {
+    std::cerr << "Error, only one components ? Use clitkImageInterpolate instead." << std::endl;
+    exit(0);
+  }
+  typedef typename ImageType::PixelType PixelType;
+//   if (m_NbOfComponents == 2) Update_WithDimAndPixelTypeAndComponent<ImageType::ImageDimension,PixelType,2>();
+  if (m_NbOfComponents == 3) Update_WithDimAndPixelTypeAndComponent<ImageType::ImageDimension,PixelType,3>();
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template<unsigned int Dim, class PixelType, unsigned int DimCompo>
+void clitk::VFInterpolateGenericFilter::Update_WithDimAndPixelTypeAndComponent()
+{
+  // Reading input
+  //  typedef itk::Vector<PixelType, DimCompo> DisplacementType;
+  typedef PixelType DisplacementType;
+  typedef itk::Image< DisplacementType, Dim > ImageType;
+
+  typename ImageType::Pointer input1 = clitk::readImage<ImageType>(m_InputFilenames[0], m_IOVerbose);
+  typename ImageType::Pointer input2 = clitk::readImage<ImageType>(mInputFilename2, m_IOVerbose);
+  
+  // Main filter
+  typename ImageType::Pointer outputImage = ComputeImage<ImageType>(input1, input2);
+
+  // Write results
+  SetNextOutput<ImageType>(outputImage);
+}
+//--------------------------------------------------------------------
+
+//--------------------------------------------------------------------
+template<class ImageType>
+typename ImageType::Pointer
+clitk::VFInterpolateGenericFilter::ComputeImage(typename ImageType::Pointer inputImage1, typename ImageType::Pointer inputImage2)
+{
+
+  // Some typedefs
+  typedef itk::Image<typename ImageType::PixelType::ValueType, ImageType::ImageDimension> ScalarImageType;
+  typedef itk::Image<typename ImageType::PixelType::ValueType, ImageType::ImageDimension + 1> InterpolationImageType;
+  typedef itk::NthElementImageAdaptor<ImageType, typename ImageType::PixelType::ValueType> ImageAdaptorType;
+  typename ImageAdaptorType::Pointer adaptor1 = ImageAdaptorType::New();
+  typename ImageAdaptorType::Pointer adaptor2 = ImageAdaptorType::New();
+  adaptor1->SetImage(inputImage1);
+  adaptor2->SetImage(inputImage2);
+  
+  // Create Image Filter
+  typedef itk::InterpolateImageFilter<ImageAdaptorType, ScalarImageType> FilterType;
+  typename FilterType::Pointer filter = FilterType::New();
+  filter->SetInput1(adaptor1);
+  filter->SetInput2(adaptor2);
+  filter->SetDistance(mDistance);
+
+  // Select interpolator
+  if (mInterpolatorName == "nn") {
+    typedef itk::NearestNeighborInterpolateImageFunction<InterpolationImageType> InterpolatorType;
+    typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
+    filter->SetInterpolator(interpolator);
+  } else {
+    if (mInterpolatorName == "linear") {
+      typedef itk::LinearInterpolateImageFunction<InterpolationImageType> InterpolatorType;
+      typename InterpolatorType::Pointer interpolator =  InterpolatorType::New();
+      filter->SetInterpolator(interpolator);
+    } else {
+      std::cerr << "Sorry, I do not know the interpolator (for vector field) '" << mInterpolatorName
+                << "'. Known interpolators are :  nn, linear" << std::endl;
+      exit(0);
+    }
+  }
+
+  typename ImageType::Pointer output = ImageType::New();
+  typename ImageAdaptorType::Pointer adaptorOutput = ImageAdaptorType::New();
+  output->CopyInformation(inputImage1);
+  output->SetRegions(inputImage1->GetLargestPossibleRegion());
+  output->Allocate();
+
+  typedef itk::ImageRegionIterator<ScalarImageType> IteratorType1;
+  typedef itk::ImageRegionIterator<ImageAdaptorType> IteratorType2;
+  
+  for (unsigned int i = 0; i < ImageType::PixelType::Dimension; i++) {
+    adaptor1->SelectNthElement(i);
+    adaptor2->SelectNthElement(i);
+
+    // Go !
+    try {
+      filter->Update();
+    } catch( itk::ExceptionObject & err ) {
+      std::cerr << "Error while filtering " << m_InputFilenames[0].c_str()
+                << " " << err << std::endl;
+      exit(0);
+    }
+
+    adaptorOutput->SelectNthElement(i);
+    adaptorOutput->SetImage(output);
+    
+    IteratorType1 it1(filter->GetOutput(), filter->GetOutput()->GetLargestPossibleRegion());
+    IteratorType2 it2(adaptorOutput, adaptorOutput->GetLargestPossibleRegion());
+    
+    it1.GoToBegin();
+    it2.GoToBegin();
+    while ( ! it1.IsAtEnd() ) {
+      it2.Set(it1.Get());
+      ++it1;
+      ++it2;
+    }
+  }
+  
+  // Return result
+  return output;
+
+}
+//--------------------------------------------------------------------
+
+//--------------------------------------------------------------------
+void clitk::VFInterpolateGenericFilter::SetInterpolationName(const std::string & inter)
+{
+  mInterpolatorName = inter;
+}
+//--------------------------------------------------------------------
+
+#endif
+