X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=filters%2FclitkVFResampleGenericFilter.cxx;h=f335415a2d88ccc8913dcb0d059e157cfb966785;hb=7f78d7b6b6bca61ca0aae964a6747cd769c0ba94;hp=9eb995b836254593acf1457eaa0ea38b1f5d06ce;hpb=931a42358442f4ee4f314613c991c838d4b4e3b7;p=clitk.git

diff --git a/filters/clitkVFResampleGenericFilter.cxx b/filters/clitkVFResampleGenericFilter.cxx
index 9eb995b..f335415 100644
--- a/filters/clitkVFResampleGenericFilter.cxx
+++ b/filters/clitkVFResampleGenericFilter.cxx
@@ -14,7 +14,11 @@
 #include "clitkVFResampleGenericFilter.h"
 
 //--------------------------------------------------------------------
-clitk::VFResampleGenericFilter::VFResampleGenericFilter() {
+clitk::VFResampleGenericFilter::VFResampleGenericFilter():
+  clitk::ImageToImageGenericFilter<Self>("VFResample") {
+  InitializeImageType<2>();
+  InitializeImageType<3>();
+  InitializeImageType<4>();
   mApplyGaussianFilterBefore = false;
   mDefaultPixelValue = 0.0;
   mInterpolatorName = "NN";
@@ -22,26 +26,158 @@ clitk::VFResampleGenericFilter::VFResampleGenericFilter() {
 }
 //--------------------------------------------------------------------
 
+
+//--------------------------------------------------------------------
+template<unsigned int Dim>
+void clitk::VFResampleGenericFilter::InitializeImageType() {
+  ADD_IMAGE_TYPE(Dim, float);
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template<class ImageType>
+void clitk::VFResampleGenericFilter::UpdateWithInputImageType() {
+
+  if (mNbOfComponents == 1) { 
+    std::cerr << "Error, only one components ? Use clitkImageResample instead." << std::endl;
+    exit(0);
+  }
+  typedef typename ImageType::PixelType PixelType;
+  if (mNbOfComponents == 2) Update_WithDimAndPixelTypeAndComponent<ImageType::ImageDimension,PixelType,2>();
+  if (mNbOfComponents == 3) Update_WithDimAndPixelTypeAndComponent<ImageType::ImageDimension,PixelType,3>();  
+  if (mNbOfComponents == 4) Update_WithDimAndPixelTypeAndComponent<ImageType::ImageDimension,PixelType,4>();  
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template<unsigned int Dim, class PixelType, unsigned int DimCompo>
+void clitk::VFResampleGenericFilter::Update_WithDimAndPixelTypeAndComponent() {
+  // Reading input
+  typedef itk::Vector<PixelType, DimCompo> DisplacementType;
+  typedef itk::Image< DisplacementType, Dim > ImageType;
+
+  typename ImageType::Pointer input = clitk::readImage<ImageType>(mInputFilenames, mIOVerbose);
+
+  // Main filter
+  typename ImageType::Pointer outputImage = ComputeImage<ImageType>(input);
+
+  // Write results
+  SetNextOutput<ImageType>(outputImage);
+}
+//--------------------------------------------------------------------
+
 //--------------------------------------------------------------------
-void clitk::VFResampleGenericFilter::Update() {
-  // Load image header
-  itk::ImageIOBase::Pointer header = clitk::readImageHeader(mInputFilenames[0]);
+template<class ImageType>
+typename ImageType::Pointer 
+clitk::VFResampleGenericFilter::ComputeImage(typename ImageType::Pointer inputImage) {
+
+  // Check options
+  static unsigned int dim = ImageType::ImageDimension;
+  if (mOutputSize.size() != dim) {
+    std::cerr << "Please set size with " << dim << " dimensions." << std::endl;
+    return NULL;
+  }
+  if (mOutputSpacing.size() != dim) {
+    std::cerr << "Please set spacing with " << dim << " dimensions." << std::endl;
+    return NULL;
+  }
+  mOutputOrigin.resize(dim);
+
+  if (mApplyGaussianFilterBefore && mSigma.size() != dim) {
+    std::cerr << "Please set sigma with " << dim << " dimensions." << std::endl;
+    return NULL;
+  }
+
+  // Some typedefs
+  typedef typename ImageType::SizeType SizeType;
+  typedef typename ImageType::SpacingType SpacingType;
+  typedef typename ImageType::PointType PointType;
+
+  // Create Image Filter
+  typedef itk::VectorResampleImageFilter<ImageType,ImageType> FilterType;
+  typename FilterType::Pointer filter = FilterType::New();
+    
+  // Instance of the transform object to be passed to the resample
+  // filter. By default, identity transform is applied
+  typedef itk::AffineTransform<double, ImageType::ImageDimension> TransformType;
+  typename TransformType::Pointer transform =  TransformType::New();
+  filter->SetTransform(transform);
+
+  // Set filter's parameters
+  SizeType outputSize;
+  SpacingType outputSpacing;
+  PointType outputOrigin;
+  for(unsigned int i=0; i<ImageType::ImageDimension; i++) {
+    outputSize[i] = mOutputSize[i];
+    outputSpacing[i] = mOutputSpacing[i];
+    outputOrigin[i] = inputImage->GetOrigin()[i];
+  }
+
+  filter->SetSize(outputSize);
+  filter->SetOutputSpacing(outputSpacing);
+  filter->SetOutputOrigin(outputOrigin);
+  filter->SetDefaultPixelValue(static_cast<typename ImageType::PixelType>(mDefaultPixelValue));
+
+  // Select interpolator
+  if (mInterpolatorName == "nn") {
+    typedef itk::VectorNearestNeighborInterpolateImageFunction<ImageType, double> InterpolatorType;     
+    typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
+    filter->SetInterpolator(interpolator);
+  }
+  else { 
+    if (mInterpolatorName == "linear") {
+      typedef itk::VectorLinearInterpolateImageFunction<ImageType, double> 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);
+    }
+  }
+
+  // Build initial Gaussian bluring (if needed)
+  typedef itk::RecursiveGaussianImageFilter<ImageType, ImageType> GaussianFilterType;
+  std::vector<typename GaussianFilterType::Pointer> gaussianFilters;
+  if (mApplyGaussianFilterBefore) {
+    for(unsigned int i=0; i<ImageType::ImageDimension; i++) {
+      // Create filter
+      gaussianFilters.push_back(GaussianFilterType::New());
+      // Set options
+      gaussianFilters[i]->SetDirection(i);
+      gaussianFilters[i]->SetOrder(GaussianFilterType::ZeroOrder);
+      gaussianFilters[i]->SetNormalizeAcrossScale(false);
+      gaussianFilters[i]->SetSigma(mSigma[i]); // in millimeter !
+      // Set input
+      if (i==0) gaussianFilters[i]->SetInput(inputImage);
+      else gaussianFilters[i]->SetInput(gaussianFilters[i-1]->GetOutput());
+    }
+    filter->SetInput(gaussianFilters[ImageType::ImageDimension-1]->GetOutput());
+  }
+  else {
+    filter->SetInput(inputImage);
+  }
+
+  // Go !
+  try { 
+    filter->Update();
+  }
+  catch( itk::ExceptionObject & err ) {
+    std::cerr << "Error while filtering " << mInputFilenames[0].c_str() 
+	      << " " << err << std::endl;
+    exit(0);
+  }
+
+  // Return result
+  return filter->GetOutput();
   
-  // Determine dim, pixel type, number of components
-  mDim = header->GetNumberOfDimensions();
-  mPixelTypeName = header->GetComponentTypeAsString(header->GetComponentType());  
-  mNbOfComponents = header->GetNumberOfComponents();
-
-  // Switch by dimension
-  if (mDim == 2) { Update_WithDim<2>(); return; }
-  if (mDim == 3) { Update_WithDim<3>(); return; }
-  if (mDim == 4) { Update_WithDim<4>(); return; }
-
-  std::cerr << "Error, dimension of input image is " << mDim << ", but I only work with 2,3,4." << std::endl;
-  exit(0);
 }
 //--------------------------------------------------------------------
 
+
 //--------------------------------------------------------------------
 void clitk::VFResampleGenericFilter::SetOutputSize(const std::vector<int> & size) {
   mOutputSize.resize(size.size());