common/clitkCoeffsToDVF.h

   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 clitkCoeffsToDVF_h
  19 #define clitkCoeffsToDVF_h
  20
  21 #include "itkImageFileReader.h"
  22 #include "itkImageIOBase.h"
  23
  24 #include "clitkBSplineDeformableTransform.h"
  25 #include "clitkResampleBSplineDeformableTransformImageFilter.h"
  26 #if ITK_VERSION_MAJOR >= 4
  27 #include "itkTransformToDisplacementFieldSource.h"
  28 #else
  29 #include "itkTransformToDeformationFieldSource.h"
  30 #endif
  31 #include "itkBSplineDeformableTransform.h"
  32
  33 namespace clitk
  34 {
  35   //-------------------------------------------------------------------
  36   // Initialize transform from coefficient images
  37   //-------------------------------------------------------------------
  38   template <class TransformType>
  39   void
  40   SetInitialTransformParameters(typename TransformType::Pointer transform, const typename TransformType::CoefficientImageType::Pointer coefficientImage, typename TransformType::CoefficientImageType::SpacingType outputSpacing)
  41   {
  42     unsigned int dim = TransformType::CoefficientImageType::ImageDimension;
  43     transform->SetSplineOrder(3);
  44     transform->SetGridRegion( coefficientImage->GetLargestPossibleRegion() );
  45     transform->SetGridOrigin( coefficientImage->GetOrigin() );
  46     transform->SetGridSpacing( coefficientImage->GetSpacing() );
  47     transform->SetGridDirection( coefficientImage->GetDirection() );
  48     typename TransformType::RegionType::SizeType samplingFactors;
  49     for (unsigned int i=0; i< dim; i++) {
  50         samplingFactors[i]= (int) ( coefficientImage->GetSpacing()[i]/ outputSpacing[i]);
  51     }
  52     transform->SetLUTSamplingFactors(samplingFactors);
  53
  54     typedef typename TransformType::ParametersType     ParametersType;
  55     const unsigned int numberOfParameters = transform->GetNumberOfParameters();
  56     ParametersType params(numberOfParameters);
  57     params.Fill( 0.0 );
  58
  59     typedef itk::ImageRegionConstIterator<typename TransformType::CoefficientImageType> Iterator;
  60     Iterator it (coefficientImage, coefficientImage->GetLargestPossibleRegion() );
  61     it.GoToBegin();
  62     unsigned int k = 0;
  63     while (! it.IsAtEnd()) {
  64         for (unsigned int j = 0; j < dim; j++)
  65             params[k+j]=it.Get()[j];
  66
  67         ++it;
  68         k += dim;
  69     }
  70
  71     transform->SetParameters(params);
  72     transform->SetBulkTransform(NULL);
  73   }
  74
  75   //-------------------------------------------------------------------
  76   // Convert Coefficient image to DVF
  77   //-------------------------------------------------------------------
  78   template<class DisplacementFieldType>
  79   typename DisplacementFieldType::Pointer
  80   BLUTCoeffsToDVF(std::string fileName, std::string likeFileName, bool verbose = false)
  81   {
  82     const unsigned int dim = DisplacementFieldType::ImageDimension;
  83     typedef clitk::BSplineDeformableTransform<double, DisplacementFieldType::ImageDimension, DisplacementFieldType::ImageDimension> TransformType;
  84     typedef typename TransformType::CoefficientImageType CoefficientImageType;
  85
  86     typedef itk::ImageFileReader<CoefficientImageType> CoeffReaderType;
  87     typename CoeffReaderType::Pointer reader = CoeffReaderType::New();
  88     reader->SetFileName(fileName);
  89     reader->Update();
  90
  91   #if ITK_VERSION_MAJOR >= 4
  92         typedef itk::TransformToDisplacementFieldSource<DisplacementFieldType, double> ConvertorType;
  93   #else
  94         typedef itk::TransformToDeformationFieldSource<DisplacementFieldType, double> ConvertorType;
  95   #endif
  96
  97     typedef itk::ImageIOBase ImageIOType;
  98     typename ImageIOType::Pointer imageIO = itk::ImageIOFactory::CreateImageIO(likeFileName.c_str(), itk::ImageIOFactory::ReadMode);
  99     imageIO->SetFileName(likeFileName);
 100     imageIO->ReadImageInformation();
 101
 102     typename ConvertorType::Pointer convertor= ConvertorType::New();
 103     typename ConvertorType::SizeType output_size;
 104     typename ConvertorType::SpacingType output_spacing;
 105     typename ConvertorType::OriginType output_origin;
 106     typename ConvertorType::DirectionType output_direction;
 107     for (unsigned int i = 0; i < dim; i++) {
 108       output_size[i] = imageIO->GetDimensions(i);
 109       output_spacing[i] = imageIO->GetSpacing(i);
 110       output_origin[i] = imageIO->GetOrigin(i);
 111       for (unsigned int j = 0; j < DisplacementFieldType::ImageDimension; j++)
 112         output_direction[i][j] = imageIO->GetDirection(i)[j];
 113     }
 114
 115     typename CoefficientImageType::Pointer coeffs = reader->GetOutput();
 116     typename TransformType::Pointer transform = TransformType::New();
 117     SetInitialTransformParameters<TransformType>(transform, coeffs, output_spacing);
 118
 119     if (verbose) {
 120       std::cout << "Interpolating coefficients with grid:" << std::endl;
 121       std::cout << output_size << output_spacing << std::endl;
 122     }
 123
 124     convertor->SetNumberOfThreads(1);
 125     convertor->SetTransform(transform);
 126     convertor->SetOutputOrigin(output_origin);
 127     convertor->SetOutputSpacing(output_spacing);
 128     convertor->SetOutputSize(output_size);
 129     convertor->SetOutputDirection(output_direction);
 130     convertor->Update();
 131
 132     return convertor->GetOutput();
 133   }
 134 }
 135
 136 #endif