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=NULL;
 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
 147 #if ITK_VERSION_MAJOR >= 4
 148       typename ITKTransformType::CoefficientImageArray coefficient_images;
 149 #else
 150       typename ITKTransformType::ImagePointer coefficient_images[BLUTCoefficientImageType::ImageDimension];
 151 #endif
 152
 153       for (unsigned int i=0; i < BLUTCoefficientImageType::ImageDimension; i++) {
 154           component_filter[i] = FilterType::New();
 155           component_filter[i]->SetInput(input);
 156           component_filter[i]->SetComponentIndex(i);
 157           component_filter[i]->Update();
 158           coefficient_images[i] = component_filter[i]->GetOutput();
 159       }
 160
 161 #if ITK_VERSION_MAJOR >= 4
 162       // RP: 16/01/2013
 163       // ATTENTION: Apparently, there's a bug in the SetCoefficientImages function of ITK 4.x
 164       // I needed to use SetParametersByValue instead.
 165       //
 166       typename ITKTransformType::ParametersType params(input->GetPixelContainer()->Size() * BLUTCoefficientImageType::ImageDimension);
 167       for (unsigned int i=0; i < BLUTCoefficientImageType::ImageDimension; i++) {
 168         for (unsigned int j=0; j < coefficient_images[i]->GetPixelContainer()->Size(); j++)
 169           params[input->GetPixelContainer()->Size() * i + j] = coefficient_images[i]->GetPixelContainer()->GetBufferPointer()[j];
 170       }
 171
 172       m_ITKTransform->SetGridOrigin(input->GetOrigin());
 173       m_ITKTransform->SetGridDirection(input->GetDirection());
 174       m_ITKTransform->SetGridRegion(input->GetLargestPossibleRegion());
 175       m_ITKTransform->SetGridSpacing(input->GetSpacing());
 176       m_ITKTransform->SetParametersByValue(params);
 177 #else
 178       m_ITKTransform->SetCoefficientImage(coefficient_images);
 179 #endif
 180
 181       m_GenericTransform = m_ITKTransform;
 182     }
 183
 184 #if ITK_VERSION_MAJOR > 4 || (ITK_VERSION_MAJOR == 4 && ITK_VERSION_MINOR >= 6)
 185     m_Filter->SetReferenceImage(output);
 186 #else
 187     m_Filter->SetOutputOrigin(output->GetOrigin());
 188     m_Filter->SetOutputSpacing(output->GetSpacing());
 189     m_Filter->SetOutputSize(output->GetLargestPossibleRegion().GetSize());
 190 #endif
 191     m_Filter->SetTransform(m_GenericTransform);
 192
 193     m_Filter->Update();
 194
 195     this->SetNthOutput(0, m_Filter->GetOutput());
 196   }
 197
 198
 199 }