registration/clitkConvertBLUTCoeffsToVFFilter.txx

   1
   2 #include "clitkConvertBLUTCoeffsToVFFilter.h"
   3 #include "clitkBSplineDeformableTransform.h"
   4 //#include "clitkTransformToDeformationFieldSource.h"
   5 //#include "clitkShapedBLUTSpatioTemporalDeformableTransform.h"
   6 #include "itkImageMaskSpatialObject.h"
   7
   8 //#include "clitkConvertBSplineDeformableTransformToVFGenericFilter.h"
   9 #include "clitkVectorImageToImageFilter.h"
  10 #include "itkBSplineDeformableTransform.h"
  11
  12 namespace clitk
  13 {
  14   template <class TDVFType>
  15   ConvertBLUTCoeffsToVFFilter<TDVFType>::ConvertBLUTCoeffsToVFFilter()
  16   {
  17     m_Verbose = false;
  18     m_TransformType = 0;
  19     m_OutputOrigin.Fill(0);
  20     m_OutputSpacing.Fill(1);
  21     m_OutputSize.Fill(1);
  22     m_BLUTSplineOrders.Fill(3);
  23
  24     m_Filter = ConvertorType::New();
  25     m_BLUTTransform = BLUTTransformType::New();
  26     m_ITKTransform = ITKTransformType::New();
  27   }
  28
  29   template <class TDVFType>
  30   void ConvertBLUTCoeffsToVFFilter<TDVFType>::GenerateOutputInformation()
  31   {
  32     if (m_Verbose)
  33       std::cout << "ConvertBLUTCoeffsToVFFilter<TDVFType>::GenerateOutputInformation()" << std::endl;
  34
  35     Superclass::GenerateOutputInformation();
  36     OutputImagePointer output = this->GetOutput();
  37     if (!m_LikeFileName.empty())
  38     {
  39         typename LikeImageType::Pointer imageIO = itk::ImageIOFactory::CreateImageIO(m_LikeFileName.c_str(), itk::ImageIOFactory::ReadMode);
  40         imageIO->SetFileName(m_LikeFileName.c_str());
  41         imageIO->ReadImageInformation();
  42
  43         typename ConvertorType::SizeType output_size;
  44         typename ConvertorType::SpacingType output_spacing;
  45         typename ConvertorType::OriginType output_origin;
  46         typename ConvertorType::DirectionType output_direction;
  47         for (unsigned int i = 0; i < OutputImageType::ImageDimension; i++) {
  48           output_size[i] = imageIO->GetDimensions(i);
  49           output_spacing[i] = imageIO->GetSpacing(i);
  50           output_origin[i] = imageIO->GetOrigin(i);
  51           //for (unsigned int j = 0; j < Dimension; j++)
  52           //  output_direction[i][j] = imageIO->GetDirection(i)[j];
  53         }
  54
  55         if (m_Verbose)
  56           std::cout << output_origin << output_size << output_spacing << std::endl;
  57
  58         output->SetOrigin(output_origin);
  59         output->SetSpacing(output_spacing);
  60         //output->SetDirection(output_direction);
  61         OutputImageRegionType output_region;
  62         output_region.SetSize(output_size);
  63         output->SetRegions(output_region);
  64     }
  65     else
  66     {
  67         if (m_Verbose)
  68           std::cout << m_OutputOrigin << m_OutputSize << m_OutputSpacing << std::endl;
  69
  70         output->SetOrigin(m_OutputOrigin);
  71         output->SetSpacing(m_OutputSpacing);
  72         OutputImageRegionType output_region;
  73         output_region.SetSize(m_OutputSize);
  74         output->SetRegions(output_region);
  75     }
  76
  77   }
  78
  79   template <class TDVFType>
  80   void ConvertBLUTCoeffsToVFFilter<TDVFType>::GenerateData()
  81   {
  82     // Read the input
  83     typedef itk::ImageFileReader<BLUTCoefficientImageType> InputReaderType;
  84     typename InputReaderType::Pointer reader = InputReaderType::New();
  85     reader->SetFileName(m_InputFileName.c_str());
  86     reader->Update();
  87     typename BLUTCoefficientImageType::Pointer input = reader->GetOutput();
  88     OutputImagePointer output = this->GetOutput();
  89
  90     if (m_TransformType != 0 ) { // using BLUT
  91       // Spline orders:  Default is cubic splines
  92         if (m_Verbose) {
  93           std::cout << "Using clitk::BLUT." << std::endl;
  94           std::cout << "Setting spline orders and sampling factors." << std::endl;
  95         }
  96
  97       m_BLUTTransform->SetSplineOrders(m_BLUTSplineOrders);
  98
  99       typename BLUTCoefficientImageType::SizeType samplingFactors;
 100       for (unsigned int i=0; i< OutputImageType::ImageDimension; i++)
 101       {
 102           samplingFactors[i]= (int) ( input->GetSpacing()[i]/ output->GetSpacing()[i]);
 103           if (m_Verbose) std::cout<<"Setting sampling factor "<<i<<" to "<<samplingFactors[i]<<"..."<<std::endl;
 104       }
 105       m_BLUTTransform->SetLUTSamplingFactors(samplingFactors);
 106       m_BLUTTransform->SetCoefficientImage(input);
 107
 108       // Mask
 109       typedef itk::ImageMaskSpatialObject<BLUTCoefficientImageType::ImageDimension >   MaskType;
 110       typename MaskType::Pointer  spatialObjectMask=ITK_NULLPTR;
 111       if (!m_MaskFileName.empty())
 112       {
 113           typedef itk::Image< unsigned char, BLUTCoefficientImageType::ImageDimension >   ImageMaskType;
 114           typedef itk::ImageFileReader< ImageMaskType >    MaskReaderType;
 115           typename MaskReaderType::Pointer  maskReader = MaskReaderType::New();
 116           maskReader->SetFileName(m_MaskFileName.c_str());
 117
 118           try
 119           {
 120               maskReader->Update();
 121           }
 122           catch ( itk::ExceptionObject & err )
 123           {
 124               std::cerr << "ExceptionObject caught while reading mask !" << std::endl;
 125               std::cerr << err << std::endl;
 126               return;
 127           }
 128           if (m_Verbose)std::cout <<"Mask was read..." <<std::endl;
 129
 130           // Set the image to the spatialObject
 131           spatialObjectMask = MaskType::New();
 132           spatialObjectMask->SetImage( maskReader->GetOutput() );
 133           m_BLUTTransform->SetMask(spatialObjectMask);
 134       }
 135
 136       m_GenericTransform = m_BLUTTransform;
 137     }
 138     else { // using ITK transform
 139       if (m_Verbose) {
 140         std::cout << "Using itk::BSpline" << std::endl;
 141         std::cout << "Extracting components from input coefficient image and creating one coefficient image per-component" << std::endl;
 142       }
 143
 144       typedef clitk::VectorImageToImageFilter<BLUTCoefficientImageType, typename ITKTransformType::ImageType> FilterType;
 145       typename FilterType::Pointer component_filter[BLUTCoefficientImageType::ImageDimension];
 146       typename ITKTransformType::CoefficientImageArray coefficient_images;
 147
 148       for (unsigned int i=0; i < BLUTCoefficientImageType::ImageDimension; i++) {
 149           component_filter[i] = FilterType::New();
 150           component_filter[i]->SetInput(input);
 151           component_filter[i]->SetComponentIndex(i);
 152           component_filter[i]->Update();
 153           coefficient_images[i] = component_filter[i]->GetOutput();
 154       }
 155
 156       // RP: 16/01/2013
 157       // ATTENTION: Apparently, there's a bug in the SetCoefficientImages function of ITK 4.x
 158       // I needed to use SetParametersByValue instead.
 159       //
 160       typename ITKTransformType::ParametersType params(input->GetPixelContainer()->Size() * BLUTCoefficientImageType::ImageDimension);
 161       for (unsigned int i=0; i < BLUTCoefficientImageType::ImageDimension; i++) {
 162         for (unsigned int j=0; j < coefficient_images[i]->GetPixelContainer()->Size(); j++)
 163           params[input->GetPixelContainer()->Size() * i + j] = coefficient_images[i]->GetPixelContainer()->GetBufferPointer()[j];
 164       }
 165       typename ITKTransformTypeZero::Pointer zero;
 166       typename ITKTransformTypeOne::Pointer one;
 167       typename ITKTransformTypeTwo::Pointer two;
 168       typename ITKTransformTypeFour::Pointer four;
 169
 170       switch(m_BLUTSplineOrders[0])
 171       {
 172         case 0:
 173           zero = ITKTransformTypeZero::New();
 174           zero->SetGridOrigin(input->GetOrigin());
 175           zero->SetGridDirection(input->GetDirection());
 176           zero->SetGridRegion(input->GetLargestPossibleRegion());
 177           zero->SetGridSpacing(input->GetSpacing());
 178           m_GenericTransform = zero;
 179           break;
 180         case 1:
 181           one = ITKTransformTypeOne::New();
 182           one->SetGridOrigin(input->GetOrigin());
 183           one->SetGridDirection(input->GetDirection());
 184           one->SetGridRegion(input->GetLargestPossibleRegion());
 185           one->SetGridSpacing(input->GetSpacing());
 186           m_GenericTransform = one;
 187           break;
 188         case 2:
 189           two = ITKTransformTypeTwo::New();
 190           two->SetGridOrigin(input->GetOrigin());
 191           two->SetGridDirection(input->GetDirection());
 192           two->SetGridRegion(input->GetLargestPossibleRegion());
 193           two->SetGridSpacing(input->GetSpacing());
 194           m_GenericTransform = two;
 195           break;
 196         case 3:
 197           m_ITKTransform = ITKTransformType::New();
 198           m_ITKTransform->SetGridOrigin(input->GetOrigin());
 199           m_ITKTransform->SetGridDirection(input->GetDirection());
 200           m_ITKTransform->SetGridRegion(input->GetLargestPossibleRegion());
 201           m_ITKTransform->SetGridSpacing(input->GetSpacing());
 202           m_GenericTransform = m_ITKTransform;
 203           break;
 204         case 4:
 205           four = ITKTransformTypeFour::New();
 206           four->SetGridOrigin(input->GetOrigin());
 207           four->SetGridDirection(input->GetDirection());
 208           four->SetGridRegion(input->GetLargestPossibleRegion());
 209           four->SetGridSpacing(input->GetSpacing());
 210           m_GenericTransform = four;
 211           break;
 212       }
 213       m_GenericTransform->SetParametersByValue(params);
 214     }
 215
 216 #if ITK_VERSION_MAJOR > 4 || (ITK_VERSION_MAJOR == 4 && ITK_VERSION_MINOR >= 6)
 217     m_Filter->SetSize(output->GetLargestPossibleRegion().GetSize());
 218 #else
 219     m_Filter->SetOutputSize(output->GetLargestPossibleRegion().GetSize());
 220 #endif
 221     m_Filter->SetOutputOrigin(output->GetOrigin());
 222     m_Filter->SetOutputSpacing(output->GetSpacing());
 223     m_Filter->SetTransform(m_GenericTransform);
 224
 225     m_Filter->Update();
 226
 227     this->SetNthOutput(0, m_Filter->GetOutput());
 228   }
 229
 230
 231 }