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 #if ITK_VERSION_MAJOR >= 4
  11 #include "itkTransformToDisplacementFieldSource.h"
  12 #else
  13 #include "itkTransformToDeformationFieldSource.h"
  14 #endif
  15 #include "itkBSplineDeformableTransform.h"
  16
  17 namespace clitk
  18 {
  19   template <class TDVFType>
  20   ConvertBLUTCoeffsToVFFilter<TDVFType>::ConvertBLUTCoeffsToVFFilter()
  21   {
  22     m_Verbose = false;
  23     m_TransformType = 0;
  24     m_OutputOrigin.Fill(0);
  25     m_OutputSpacing.Fill(1);
  26     m_OutputSize.Fill(1);
  27     m_BLUTSplineOrders.Fill(3);
  28
  29     m_Filter = ConvertorType::New();
  30     m_BLUTTransform = BLUTTransformType::New();
  31     m_ITKTransform = ITKTransformType::New();
  32   }
  33
  34   template <class TDVFType>
  35   void ConvertBLUTCoeffsToVFFilter<TDVFType>::GenerateOutputInformation()
  36   {
  37     if (m_Verbose)
  38       std::cout << "ConvertBLUTCoeffsToVFFilter<TDVFType>::GenerateOutputInformation()" << std::endl;
  39
  40     Superclass::GenerateOutputInformation();
  41     OutputImagePointer output = this->GetOutput();
  42     if (!m_LikeFileName.empty())
  43     {
  44         typename LikeImageType::Pointer imageIO = itk::ImageIOFactory::CreateImageIO(m_LikeFileName.c_str(), itk::ImageIOFactory::ReadMode);
  45         imageIO->SetFileName(m_LikeFileName.c_str());
  46         imageIO->ReadImageInformation();
  47
  48         typename ConvertorType::SizeType output_size;
  49         typename ConvertorType::SpacingType output_spacing;
  50         typename ConvertorType::OriginType output_origin;
  51         typename ConvertorType::DirectionType output_direction;
  52         for (unsigned int i = 0; i < OutputImageType::ImageDimension; i++) {
  53           output_size[i] = imageIO->GetDimensions(i);
  54           output_spacing[i] = imageIO->GetSpacing(i);
  55           output_origin[i] = imageIO->GetOrigin(i);
  56           //for (unsigned int j = 0; j < Dimension; j++)
  57           //  output_direction[i][j] = imageIO->GetDirection(i)[j];
  58         }
  59
  60         if (m_Verbose)
  61           std::cout << output_origin << output_size << output_spacing << std::endl;
  62
  63         output->SetOrigin(output_origin);
  64         output->SetSpacing(output_spacing);
  65         //output->SetDirection(output_direction);
  66         OutputImageRegionType output_region;
  67         output_region.SetSize(output_size);
  68         output->SetRegions(output_region);
  69     }
  70     else
  71     {
  72         if (m_Verbose)
  73           std::cout << m_OutputOrigin << m_OutputSize << m_OutputSpacing << std::endl;
  74
  75         output->SetOrigin(m_OutputOrigin);
  76         output->SetSpacing(m_OutputSpacing);
  77         OutputImageRegionType output_region;
  78         output_region.SetSize(m_OutputSize);
  79         output->SetRegions(output_region);
  80     }
  81
  82   }
  83
  84   template <class TDVFType>
  85   void ConvertBLUTCoeffsToVFFilter<TDVFType>::GenerateData()
  86   {
  87     // Read the input
  88     typedef itk::ImageFileReader<BLUTCoefficientImageType> InputReaderType;
  89     typename InputReaderType::Pointer reader = InputReaderType::New();
  90     reader->SetFileName(m_InputFileName.c_str());
  91     reader->Update();
  92     typename BLUTCoefficientImageType::Pointer input = reader->GetOutput();
  93     OutputImagePointer output = this->GetOutput();
  94
  95     if (m_TransformType != 0 ) { // using BLUT
  96       // Spline orders:  Default is cubic splines
  97         if (m_Verbose) {
  98           std::cout << "Using clitk::BLUT." << std::endl;
  99           std::cout << "Setting spline orders and sampling factors." << std::endl;
 100         }
 101
 102       m_BLUTTransform->SetSplineOrders(m_BLUTSplineOrders);
 103
 104       typename BLUTCoefficientImageType::SizeType samplingFactors;
 105       for (unsigned int i=0; i< OutputImageType::ImageDimension; i++)
 106       {
 107           samplingFactors[i]= (int) ( input->GetSpacing()[i]/ output->GetSpacing()[i]);
 108           if (m_Verbose) std::cout<<"Setting sampling factor "<<i<<" to "<<samplingFactors[i]<<"..."<<std::endl;
 109       }
 110       m_BLUTTransform->SetLUTSamplingFactors(samplingFactors);
 111       m_BLUTTransform->SetCoefficientImage(input);
 112
 113       // Mask
 114       typedef itk::ImageMaskSpatialObject<BLUTCoefficientImageType::ImageDimension >   MaskType;
 115       typename MaskType::Pointer  spatialObjectMask=NULL;
 116       if (!m_MaskFileName.empty())
 117       {
 118           typedef itk::Image< unsigned char, BLUTCoefficientImageType::ImageDimension >   ImageMaskType;
 119           typedef itk::ImageFileReader< ImageMaskType >    MaskReaderType;
 120           typename MaskReaderType::Pointer  maskReader = MaskReaderType::New();
 121           maskReader->SetFileName(m_MaskFileName.c_str());
 122
 123           try
 124           {
 125               maskReader->Update();
 126           }
 127           catch ( itk::ExceptionObject & err )
 128           {
 129               std::cerr << "ExceptionObject caught while reading mask !" << std::endl;
 130               std::cerr << err << std::endl;
 131               return;
 132           }
 133           if (m_Verbose)std::cout <<"Mask was read..." <<std::endl;
 134
 135           // Set the image to the spatialObject
 136           spatialObjectMask = MaskType::New();
 137           spatialObjectMask->SetImage( maskReader->GetOutput() );
 138           m_BLUTTransform->SetMask(spatialObjectMask);
 139       }
 140
 141       m_GenericTransform = m_BLUTTransform;
 142     }
 143     else { // using ITK transform
 144       if (m_Verbose) {
 145         std::cout << "Using itk::BSpline" << std::endl;
 146         std::cout << "Extracting components from input coefficient image and creating one coefficient image per-component" << std::endl;
 147       }
 148
 149       typedef clitk::VectorImageToImageFilter<BLUTCoefficientImageType, typename ITKTransformType::ImageType> FilterType;
 150       typename FilterType::Pointer component_filter[BLUTCoefficientImageType::ImageDimension];
 151
 152 #if ITK_VERSION_MAJOR >= 4
 153       typename ITKTransformType::CoefficientImageArray coefficient_images;
 154 #else
 155       typename ITKTransformType::ImagePointer coefficient_images[BLUTCoefficientImageType::ImageDimension];
 156 #endif
 157
 158       for (unsigned int i=0; i < BLUTCoefficientImageType::ImageDimension; i++) {
 159           component_filter[i] = FilterType::New();
 160           component_filter[i]->SetInput(input);
 161           component_filter[i]->SetComponentIndex(i);
 162           component_filter[i]->Update();
 163           coefficient_images[i] = component_filter[i]->GetOutput();
 164       }
 165 #if ITK_VERSION_MAJOR >= 4
 166       m_ITKTransform->SetCoefficientImages(coefficient_images);
 167 #else
 168       m_ITKTransform->SetCoefficientImage(coefficient_images);
 169 #endif
 170
 171       m_GenericTransform = m_ITKTransform;
 172     }
 173
 174     m_Filter->SetOutputOrigin(output->GetOrigin());
 175     m_Filter->SetOutputSpacing(output->GetSpacing());
 176     m_Filter->SetOutputSize(output->GetLargestPossibleRegion().GetSize());
 177     m_Filter->SetTransform(m_GenericTransform);
 178
 179     m_Filter->Update();
 180
 181     SetNthOutput(0, m_Filter->GetOutput());
 182   }
 183
 184
 185 }