2 #include "clitkConvertBLUTCoeffsToVFFilter.h"
3 #include "clitkBSplineDeformableTransform.h"
4 //#include "clitkTransformToDeformationFieldSource.h"
5 //#include "clitkShapedBLUTSpatioTemporalDeformableTransform.h"
6 #include "itkImageMaskSpatialObject.h"
8 //#include "clitkConvertBSplineDeformableTransformToVFGenericFilter.h"
9 #include "clitkVectorImageToImageFilter.h"
10 #include "itkBSplineDeformableTransform.h"
14 template <class TDVFType>
15 ConvertBLUTCoeffsToVFFilter<TDVFType>::ConvertBLUTCoeffsToVFFilter()
19 m_OutputOrigin.Fill(0);
20 m_OutputSpacing.Fill(1);
22 m_BLUTSplineOrders.Fill(3);
24 m_Filter = ConvertorType::New();
25 m_BLUTTransform = BLUTTransformType::New();
26 m_ITKTransform = ITKTransformType::New();
29 template <class TDVFType>
30 void ConvertBLUTCoeffsToVFFilter<TDVFType>::GenerateOutputInformation()
33 std::cout << "ConvertBLUTCoeffsToVFFilter<TDVFType>::GenerateOutputInformation()" << std::endl;
35 Superclass::GenerateOutputInformation();
36 OutputImagePointer output = this->GetOutput();
37 if (!m_LikeFileName.empty())
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();
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];
56 std::cout << output_origin << output_size << output_spacing << std::endl;
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);
68 std::cout << m_OutputOrigin << m_OutputSize << m_OutputSpacing << std::endl;
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);
79 template <class TDVFType>
80 void ConvertBLUTCoeffsToVFFilter<TDVFType>::GenerateData()
83 typedef itk::ImageFileReader<BLUTCoefficientImageType> InputReaderType;
84 typename InputReaderType::Pointer reader = InputReaderType::New();
85 reader->SetFileName(m_InputFileName.c_str());
87 typename BLUTCoefficientImageType::Pointer input = reader->GetOutput();
88 OutputImagePointer output = this->GetOutput();
90 if (m_TransformType != 0 ) { // using BLUT
91 // Spline orders: Default is cubic splines
93 std::cout << "Using clitk::BLUT." << std::endl;
94 std::cout << "Setting spline orders and sampling factors." << std::endl;
97 m_BLUTTransform->SetSplineOrders(m_BLUTSplineOrders);
99 typename BLUTCoefficientImageType::SizeType samplingFactors;
100 for (unsigned int i=0; i< OutputImageType::ImageDimension; i++)
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;
105 m_BLUTTransform->SetLUTSamplingFactors(samplingFactors);
106 m_BLUTTransform->SetCoefficientImage(input);
109 typedef itk::ImageMaskSpatialObject<BLUTCoefficientImageType::ImageDimension > MaskType;
110 typename MaskType::Pointer spatialObjectMask=ITK_NULLPTR;
111 if (!m_MaskFileName.empty())
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());
120 maskReader->Update();
122 catch ( itk::ExceptionObject & err )
124 std::cerr << "ExceptionObject caught while reading mask !" << std::endl;
125 std::cerr << err << std::endl;
128 if (m_Verbose)std::cout <<"Mask was read..." <<std::endl;
130 // Set the image to the spatialObject
131 spatialObjectMask = MaskType::New();
132 spatialObjectMask->SetImage( maskReader->GetOutput() );
133 m_BLUTTransform->SetMask(spatialObjectMask);
136 m_GenericTransform = m_BLUTTransform;
138 else { // using ITK transform
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;
144 typedef clitk::VectorImageToImageFilter<BLUTCoefficientImageType, typename ITKTransformType::ImageType> FilterType;
145 typename FilterType::Pointer component_filter[BLUTCoefficientImageType::ImageDimension];
146 typename ITKTransformType::CoefficientImageArray coefficient_images;
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();
157 // ATTENTION: Apparently, there's a bug in the SetCoefficientImages function of ITK 4.x
158 // I needed to use SetParametersByValue instead.
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];
165 typename ITKTransformTypeZero::Pointer zero;
166 typename ITKTransformTypeOne::Pointer one;
167 typename ITKTransformTypeTwo::Pointer two;
168 typename ITKTransformTypeFour::Pointer four;
170 switch(m_BLUTSplineOrders[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;
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;
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;
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;
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;
213 m_GenericTransform->SetParametersByValue(params);
216 #if ITK_VERSION_MAJOR > 4 || (ITK_VERSION_MAJOR == 4 && ITK_VERSION_MINOR >= 6)
217 m_Filter->SetSize(output->GetLargestPossibleRegion().GetSize());
219 m_Filter->SetOutputSize(output->GetLargestPossibleRegion().GetSize());
221 m_Filter->SetOutputOrigin(output->GetOrigin());
222 m_Filter->SetOutputSpacing(output->GetSpacing());
223 m_Filter->SetTransform(m_GenericTransform);
227 this->SetNthOutput(0, m_Filter->GetOutput());