1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
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
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.
13 It is distributed under dual licence
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
21 #include "itkImageFileReader.h"
22 #include "itkImageIOBase.h"
24 #include "clitkBSplineDeformableTransform.h"
25 #include "clitkResampleBSplineDeformableTransformImageFilter.h"
26 #if ITK_VERSION_MAJOR >= 4
27 #include "itkTransformToDisplacementFieldSource.h"
29 #include "itkTransformToDeformationFieldSource.h"
31 #include "itkBSplineDeformableTransform.h"
35 //-------------------------------------------------------------------
36 // Initialize transform from coefficient images
37 //-------------------------------------------------------------------
38 template <class TransformType>
40 SetInitialTransformParameters(typename TransformType::Pointer transform, const typename TransformType::CoefficientImageType::Pointer coefficientImage, typename TransformType::CoefficientImageType::SpacingType outputSpacing)
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]);
52 transform->SetLUTSamplingFactors(samplingFactors);
54 typedef typename TransformType::ParametersType ParametersType;
55 const unsigned int numberOfParameters = transform->GetNumberOfParameters();
56 ParametersType params(numberOfParameters);
59 typedef itk::ImageRegionConstIterator<typename TransformType::CoefficientImageType> Iterator;
60 Iterator it (coefficientImage, coefficientImage->GetLargestPossibleRegion() );
63 while (! it.IsAtEnd()) {
64 for (unsigned int j = 0; j < dim; j++)
65 params[i+j]=it.Get()[j];
71 transform->SetParameters(params);
72 transform->SetBulkTransform(NULL);
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)
82 const unsigned int dim = DisplacementFieldType::ImageDimension;
83 typedef clitk::BSplineDeformableTransform<double, DisplacementFieldType::ImageDimension, DisplacementFieldType::ImageDimension> TransformType;
84 typedef typename TransformType::CoefficientImageType CoefficientImageType;
86 typedef itk::ImageFileReader<CoefficientImageType> CoeffReaderType;
87 typename CoeffReaderType::Pointer reader = CoeffReaderType::New();
88 reader->SetFileName(fileName);
91 #if ITK_VERSION_MAJOR >= 4
92 typedef itk::TransformToDisplacementFieldSource<DisplacementFieldType, double> ConvertorType;
94 typedef itk::TransformToDeformationFieldSource<DisplacementFieldType, double> ConvertorType;
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();
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];
115 typename CoefficientImageType::Pointer coeffs = reader->GetOutput();
116 typename TransformType::Pointer transform = TransformType::New();
117 SetInitialTransformParameters<TransformType>(transform, coeffs, output_spacing);
120 std::cout << "Interpolating coefficients with grid:" << std::endl;
121 std::cout << output_size << output_spacing << std::endl;
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);
132 return convertor->GetOutput();