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 clitkAffineTransformGenericFilter_txx
19 #define clitkAffineTransformGenericFilter_txx
24 #include <itkCenteredEuler3DTransform.h>
25 #include <itkRecursiveGaussianImageFilter.h>
26 #include "clitkElastix.h"
27 #include "clitkResampleImageWithOptionsFilter.h"
32 //-----------------------------------------------------------
34 //-----------------------------------------------------------
35 template<class args_info_type>
36 AffineTransformGenericFilter<args_info_type>::AffineTransformGenericFilter()
41 //-------------------------------------------------------------------
44 //-----------------------------------------------------------
46 //-----------------------------------------------------------
47 template<class args_info_type>
48 void AffineTransformGenericFilter<args_info_type>::Update()
50 // Read the Dimension and PixelType
51 int Dimension, Components;
52 std::string PixelType;
53 ReadImageDimensionAndPixelType(m_InputFileName, Dimension, PixelType, Components);
56 if(Dimension==2) UpdateWithDim<2>(PixelType, Components);
58 if(Dimension==3) UpdateWithDim<3>(PixelType, Components);
59 else if (Dimension==4)UpdateWithDim<4>(PixelType, Components);
61 std::cout<<"Error, Only for 2, 3 or 4 Dimensions!!!"<<std::endl ;
65 //-------------------------------------------------------------------
68 //-------------------------------------------------------------------
69 // Update with the number of dimensions
70 //-------------------------------------------------------------------
71 template<class args_info_type>
72 template<unsigned int Dimension>
74 AffineTransformGenericFilter<args_info_type>::UpdateWithDim(std::string PixelType, int Components)
76 if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D and "<<Components<<" component(s) of "<< PixelType<<"..."<<std::endl;
79 if(PixelType == "short") {
80 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed short..." << std::endl;
81 UpdateWithDimAndPixelType<Dimension, signed short>();
83 else if(PixelType == "unsigned_short"){
84 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_short..." << std::endl;
85 UpdateWithDimAndPixelType<Dimension, unsigned short>();
88 else if (PixelType == "unsigned_char") {
89 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl;
90 UpdateWithDimAndPixelType<Dimension, unsigned char>();
93 // else if (PixelType == "char"){
94 // if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed_char..." << std::endl;
95 // UpdateWithDimAndPixelType<Dimension, signed char>();
97 else if(PixelType == "double"){
98 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and double..." << std::endl;
99 UpdateWithDimAndPixelType<Dimension, double>();
102 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
103 UpdateWithDimAndPixelType<Dimension, float>();
107 else if (Components==3) {
108 if (m_Verbose) std::cout << "Launching transform in "<< Dimension <<"D and 3D float (DVF)" << std::endl;
109 UpdateWithDimAndVectorType<Dimension, itk::Vector<float, Dimension> >();
112 else std::cerr<<"Number of components is "<<Components<<", not supported!"<<std::endl;
115 //-------------------------------------------------------------------
118 //-------------------------------------------------------------------
119 // Compute updated bounding box
120 //-------------------------------------------------------------------
121 template<class args_info_type>
123 AffineTransformGenericFilter<args_info_type>::ComputeSize(vnl_vector<double> inputSize, vnl_matrix<double> transformationMatrix, bool returnMin)
125 //Compute input corners
126 int Dimension = inputSize.size();
127 vnl_matrix<double> vnlOutputSize(std::pow(2, Dimension), Dimension);
128 vnlOutputSize.fill(0);
129 if (Dimension == 2) {
130 for(unsigned int i=0; i< Dimension; i++)
131 vnlOutputSize[3][i] = inputSize[i];
132 vnlOutputSize[1][0] = inputSize[0];
133 vnlOutputSize[2][1] = inputSize[1];
134 } else if (Dimension == 3) {
135 for(unsigned int i=0; i< Dimension; i++)
136 vnlOutputSize[7][i] = inputSize[i];
137 vnlOutputSize[1][0] = inputSize[0];
138 vnlOutputSize[2][1] = inputSize[1];
139 vnlOutputSize[3][2] = inputSize[2];
140 vnlOutputSize[4][0] = inputSize[0];
141 vnlOutputSize[4][1] = inputSize[1];
142 vnlOutputSize[5][1] = inputSize[1];
143 vnlOutputSize[5][2] = inputSize[2];
144 vnlOutputSize[6][0] = inputSize[0];
145 vnlOutputSize[6][2] = inputSize[2];
146 } else { //Dimension ==4
147 for(unsigned int i=0; i< Dimension; i++)
148 vnlOutputSize[15][i] = inputSize[i];
149 vnlOutputSize[1][0] = inputSize[0];
150 vnlOutputSize[2][1] = inputSize[1];
151 vnlOutputSize[3][2] = inputSize[2];
152 vnlOutputSize[4][3] = inputSize[3];
153 vnlOutputSize[5][0] = inputSize[0];
154 vnlOutputSize[5][1] = inputSize[1];
155 vnlOutputSize[6][0] = inputSize[0];
156 vnlOutputSize[6][2] = inputSize[2];
157 vnlOutputSize[7][0] = inputSize[0];
158 vnlOutputSize[7][3] = inputSize[3];
159 vnlOutputSize[8][1] = inputSize[1];
160 vnlOutputSize[8][2] = inputSize[2];
161 vnlOutputSize[9][1] = inputSize[1];
162 vnlOutputSize[9][3] = inputSize[3];
163 vnlOutputSize[10][2] = inputSize[2];
164 vnlOutputSize[10][3] = inputSize[3];
165 vnlOutputSize[11][0] = inputSize[0];
166 vnlOutputSize[11][1] = inputSize[1];
167 vnlOutputSize[11][2] = inputSize[2];
168 vnlOutputSize[12][0] = inputSize[0];
169 vnlOutputSize[12][1] = inputSize[1];
170 vnlOutputSize[12][3] = inputSize[3];
171 vnlOutputSize[13][0] = inputSize[0];
172 vnlOutputSize[13][2] = inputSize[2];
173 vnlOutputSize[13][3] = inputSize[3];
174 vnlOutputSize[14][1] = inputSize[1];
175 vnlOutputSize[14][2] = inputSize[2];
176 vnlOutputSize[14][3] = inputSize[3];
179 //Compute the transformation of all corner
180 for (unsigned int i=0; i< std::pow(2, Dimension); ++i)
181 vnlOutputSize.set_row(i, transformationMatrix*vnlOutputSize.get_row(i));
183 //Compute the bounding box taking the max and the min
184 vnl_vector<double> minBB(vnlOutputSize.get_row(0)), maxBB(vnlOutputSize.get_row(0));
185 for (unsigned int i=0; i< std::pow(2, Dimension); ++i) {
186 for (unsigned int j=0; j< Dimension; ++j) {
187 if (vnlOutputSize[i][j] < minBB[j])
188 minBB[j] = vnlOutputSize[i][j];
189 if (vnlOutputSize[i][j] > maxBB[j])
190 maxBB[j] = vnlOutputSize[i][j];
198 vnl_vector<double> size;
199 size = maxBB - minBB;
204 //-------------------------------------------------------------------
207 //-------------------------------------------------------------------
208 // Update with the number of dimensions and the pixeltype
209 //-------------------------------------------------------------------
210 template<class args_info_type>
211 template <unsigned int Dimension, class PixelType>
213 AffineTransformGenericFilter<args_info_type>::UpdateWithDimAndPixelType()
217 typedef itk::Image<PixelType, Dimension> InputImageType;
218 typedef itk::Image<PixelType, Dimension> OutputImageType;
221 typedef itk::ImageFileReader<InputImageType> InputReaderType;
222 typename InputReaderType::Pointer reader = InputReaderType::New();
223 reader->SetFileName( m_InputFileName);
225 typename InputImageType::Pointer input= reader->GetOutput();
227 //Adaptative size, spacing origin (use previous clitkResampleImage)
228 if (m_ArgsInfo.adaptive_given) {
230 typedef clitk::ResampleImageWithOptionsFilter<InputImageType, OutputImageType> ResampleImageFilterType;
231 typename ResampleImageFilterType::Pointer filter = ResampleImageFilterType::New();
232 filter->SetInput(input);
235 filter->SetVerboseOptions(m_ArgsInfo.verbose_flag);
237 // Set size / spacing
238 static const unsigned int dim = OutputImageType::ImageDimension;
239 typename OutputImageType::SpacingType spacing;
240 typename OutputImageType::SizeType size;
241 typename OutputImageType::PointType origin;
242 typename OutputImageType::DirectionType direction;
244 if (m_ArgsInfo.like_given) {
245 itk::ImageIOBase::Pointer header = clitk::readImageHeader(m_ArgsInfo.like_arg);
247 for(unsigned int i=0; i<dim; i++){
248 spacing[i] = header->GetSpacing(i);
249 size[i] = header->GetDimensions(i);
250 origin[i] = header->GetOrigin(i);
252 for(unsigned int i=0; i<dim; i++) {
253 for(unsigned int j=0;j<dim;j++) {
254 direction(i,j) = header->GetDirection(i)[j];
257 filter->SetOutputSpacing(spacing);
258 filter->SetOutputSize(size);
259 filter->SetOutputOrigin(origin);
260 filter->SetOutputDirection(direction);
263 std::cerr << "*** Warning : I could not read '" << m_ArgsInfo.like_arg << "' ***" << std::endl;
268 if (m_ArgsInfo.spacing_given == 1) {
269 filter->SetOutputIsoSpacing(m_ArgsInfo.spacing_arg[0]);
271 else if ((m_ArgsInfo.spacing_given != 0) && (m_ArgsInfo.size_given != 0)) {
272 std::cerr << "Error: use spacing or size, not both." << std::endl;
275 else if (m_ArgsInfo.spacing_given) {
276 if ((m_ArgsInfo.spacing_given != 0) && (m_ArgsInfo.spacing_given != dim)) {
277 std::cerr << "Error: spacing should have one or " << dim << " values." << std::endl;
280 for(unsigned int i=0; i<dim; i++)
281 spacing[i] = m_ArgsInfo.spacing_arg[i];
282 filter->SetOutputSpacing(spacing);
284 else if (m_ArgsInfo.size_given) {
285 if ((m_ArgsInfo.size_given != 0) && (m_ArgsInfo.size_given != dim)) {
286 std::cerr << "Error: size should have " << dim << " values." << std::endl;
289 for(unsigned int i=0; i<dim; i++)
290 size[i] = m_ArgsInfo.size_arg[i];
291 filter->SetOutputSize(size);
293 for(unsigned int i=0; i<dim; i++){
294 origin[i] = input->GetOrigin()[i];
296 for(unsigned int i=0; i<dim; i++) {
297 for(unsigned int j=0;j<dim;j++) {
298 direction(i,j) = input->GetDirection()[i][j];
301 filter->SetOutputOrigin(origin);
302 filter->SetOutputDirection(direction);
305 // Set temporal dimension
306 //filter->SetLastDimensionIsTime(m_ArgsInfo.time_flag);
309 filter->SetGaussianFilteringEnabled(m_ArgsInfo.autogauss_flag);
310 if (m_ArgsInfo.gauss_given != 0) {
311 typename ResampleImageFilterType::GaussianSigmaType g;
312 for(unsigned int i=0; i<dim; i++) {
313 g[i] = m_ArgsInfo.gauss_arg[i];
315 filter->SetGaussianSigma(g);
319 int interp = m_ArgsInfo.interp_arg;
321 filter->SetInterpolationType(ResampleImageFilterType::NearestNeighbor);
324 filter->SetInterpolationType(ResampleImageFilterType::Linear);
327 filter->SetInterpolationType(ResampleImageFilterType::BSpline);
330 filter->SetInterpolationType(ResampleImageFilterType::B_LUT);
332 std::cerr << "Error. I do not know interpolation '" << m_ArgsInfo.interp_arg
333 << "'. Choose among: nn, linear, bspline, blut, windowed sinc" << std::endl;
340 // Set default pixel value
341 filter->SetDefaultPixelValue(m_ArgsInfo.pad_arg);
344 //if (m_ArgsInfo.thread_given) {
345 // filter->SetNumberOfThreads(m_ArgsInfo.thread_arg);
350 typename OutputImageType::Pointer output = filter->GetOutput();
351 //this->template SetNextOutput<OutputImageType>(outputImage);
354 typedef itk::ImageFileWriter<OutputImageType> WriterType;
355 typename WriterType::Pointer writer = WriterType::New();
356 writer->SetFileName(m_ArgsInfo.output_arg);
357 writer->SetInput(output);
363 //Gaussian pre-filtering
364 typename itk::Vector<double, Dimension> gaussianSigma;
365 gaussianSigma.Fill(0);
366 bool gaussianFilteringEnabled(false);
367 bool autoGaussEnabled(false);
368 if (m_ArgsInfo.autogauss_given) { // Gaussian filter auto
369 autoGaussEnabled = m_ArgsInfo.autogauss_flag;
371 if (m_ArgsInfo.gauss_given) { // Gaussian filter set by user
372 gaussianFilteringEnabled = true;
373 if (m_ArgsInfo.gauss_given == 1)
375 for (unsigned int i=0; i<Dimension; i++)
377 gaussianSigma[i] = m_ArgsInfo.gauss_arg[0];
380 else if (m_ArgsInfo.gauss_given == Dimension)
382 for (unsigned int i=0; i<Dimension; i++)
384 gaussianSigma[i] = m_ArgsInfo.gauss_arg[i];
389 std::cerr << "Gaussian sigma dimension is incorrect" << std::endl;
395 typedef itk::ResampleImageFilter< InputImageType,OutputImageType > ResampleFilterType;
396 typename ResampleFilterType::Pointer resampler = ResampleFilterType::New();
399 typename itk::Matrix<double, Dimension+1, Dimension+1> matrix;
400 if (m_ArgsInfo.rotate_given || m_ArgsInfo.translate_given)
402 if (m_ArgsInfo.matrix_given)
404 std::cerr << "You must use either rotate/translate or matrix options" << std::endl;
407 itk::Array<double> transformParameters(2 * Dimension);
408 transformParameters.Fill(0.0);
409 if (m_ArgsInfo.rotate_given)
412 transformParameters[0] = m_ArgsInfo.rotate_arg[0];
414 for (unsigned int i = 0; i < 3; i++)
415 transformParameters[i] = m_ArgsInfo.rotate_arg[i];
417 if (m_ArgsInfo.translate_given)
422 for (unsigned int i = 0; i < Dimension && i < 3; i++)
423 transformParameters[pos++] = m_ArgsInfo.translate_arg[i];
427 matrix.SetIdentity();
428 itk::Matrix<double, 4, 4> tmp = GetForwardAffineMatrix3D(transformParameters);
429 for (unsigned int i = 0; i < 3; ++i)
430 for (unsigned int j = 0; j < 3; ++j)
431 matrix[i][j] = tmp[i][j];
432 for (unsigned int i = 0; i < 3; ++i)
433 matrix[i][4] = tmp[i][3];
436 matrix = GetForwardAffineMatrix<Dimension>(transformParameters);
440 if (m_ArgsInfo.matrix_given)
442 matrix= clitk::ReadMatrix<Dimension>(m_ArgsInfo.matrix_arg);
443 if (m_Verbose) std::cout << "Reading the matrix..." << std::endl;
446 if (m_ArgsInfo.elastix_given) {
447 std::string filename(m_ArgsInfo.elastix_arg);
448 matrix = createMatrixFromElastixFile<Dimension>(filename, m_Verbose);
451 matrix.SetIdentity();
455 std::cout << "Using the following matrix:" << std::endl
456 << matrix << std::endl;
457 typename itk::Matrix<double, Dimension, Dimension> rotationMatrix = clitk::GetRotationalPartMatrix(matrix);
458 typename itk::Vector<double, Dimension> translationPart = clitk::GetTranslationPartMatrix(matrix);
461 typedef itk::AffineTransform<double, Dimension> AffineTransformType;
462 typename AffineTransformType::Pointer affineTransform=AffineTransformType::New();
463 affineTransform->SetMatrix(rotationMatrix);
464 affineTransform->SetTranslation(translationPart);
467 typedef clitk::GenericInterpolator<args_info_type, InputImageType, double> GenericInterpolatorType;
468 typename GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New();
469 genericInterpolator->SetArgsInfo(m_ArgsInfo);
472 if (m_ArgsInfo.like_given) {
473 typename InputReaderType::Pointer likeReader=InputReaderType::New();
474 likeReader->SetFileName(m_ArgsInfo.like_arg);
475 likeReader->Update();
476 resampler->SetOutputParametersFromImage(likeReader->GetOutput());
477 resampler->SetOutputDirection(likeReader->GetOutput()->GetDirection());
478 if (autoGaussEnabled) { // Automated sigma when downsample
479 for(unsigned int i=0; i<Dimension; i++) {
480 if (likeReader->GetOutput()->GetSpacing()[i] > input->GetSpacing()[i]) { // downsample
481 gaussianSigma[i] = 0.5*likeReader->GetOutput()->GetSpacing()[i];// / inputSpacing[i]);
483 else gaussianSigma[i] = 0; // will be ignore after
486 } else if(m_ArgsInfo.transform_grid_flag) {
487 typename itk::Matrix<double, Dimension+1, Dimension+1> invMatrix( matrix.GetInverse() );
488 typename itk::Matrix<double, Dimension, Dimension> invRotMatrix( clitk::GetRotationalPartMatrix(invMatrix) );
489 typename itk::Vector<double,Dimension> invTrans = clitk::GetTranslationPartMatrix(invMatrix);
492 if (m_ArgsInfo.spacing_given)
493 std::cout << "Warning --spacing ignored (because --transform_grid_flag)" << std::endl;
494 if (m_ArgsInfo.origin_given)
495 std::cout << "Warning --origin ignored (because --transform_grid_flag)" << std::endl;
497 // Origin is influenced by translation but not by input direction
498 typename InputImageType::PointType outputOrigin;
499 outputOrigin = invRotMatrix *
503 // Size is influenced by affine transform matrix and input direction
504 // Size is converted to double, transformed and converted back to size type.
505 // Determine the bounding box tranforming all corners
506 vnl_vector<double> vnlOutputSize(Dimension), vnlOutputmmSize(Dimension), vnlOutputOffset(Dimension);
507 typename InputImageType::SpacingType outputSpacing;
508 for(unsigned int i=0; i< Dimension; i++) {
509 vnlOutputSize[i] = input->GetLargestPossibleRegion().GetSize()[i];
510 vnlOutputmmSize[i] = input->GetLargestPossibleRegion().GetSize()[i]*input->GetSpacing()[i];
511 vnlOutputOffset[i] = input->GetLargestPossibleRegion().GetSize()[i]*input->GetSpacing()[i];
513 vnlOutputSize = ComputeSize(vnlOutputSize, invRotMatrix.GetVnlMatrix() * input->GetDirection().GetVnlMatrix(), 0);
514 vnlOutputmmSize = ComputeSize(vnlOutputmmSize, invRotMatrix.GetVnlMatrix() * input->GetDirection().GetVnlMatrix(), 0);
515 vnlOutputOffset = ComputeSize(vnlOutputOffset, invRotMatrix.GetVnlMatrix() * input->GetDirection().GetVnlMatrix(), 1);
516 for(unsigned int i=0; i< Dimension; i++) {
517 outputSpacing[i] = vnlOutputmmSize[i]/lrint(vnlOutputSize[i]);
518 outputOrigin[i] += vnlOutputOffset[i];
520 if (autoGaussEnabled) { // Automated sigma when downsample
521 for(unsigned int i=0; i<Dimension; i++) {
522 if (outputSpacing[i] > input->GetSpacing()[i]) { // downsample
523 gaussianSigma[i] = 0.5*outputSpacing[i];// / inputSpacing[i]);
525 else gaussianSigma[i] = 0; // will be ignore after
529 typename OutputImageType::SizeType outputSize;
530 for(unsigned int i=0; i< Dimension; i++) {
531 // If the size is negative, we have a flip and we must modify
532 // the origin and the spacing accordingly.
533 if(vnlOutputSize[i]<0.) {
534 vnlOutputSize[i] *= -1.;
535 outputOrigin[i] = outputOrigin[i] + outputSpacing[i] * (vnlOutputSize[i]-1);
536 outputSpacing[i] *= -1.;
538 outputSize[i] = lrint(vnlOutputSize[i]);
540 resampler->SetSize( outputSize );
541 resampler->SetOutputSpacing( outputSpacing );
542 resampler->SetOutputOrigin( outputOrigin );
545 typename OutputImageType::SizeType outputSize;
546 if (m_ArgsInfo.size_given) {
547 for(unsigned int i=0; i< Dimension; i++)
548 outputSize[i]=m_ArgsInfo.size_arg[i];
549 } else outputSize=input->GetLargestPossibleRegion().GetSize();
552 typename OutputImageType::SpacingType outputSpacing;
553 if (m_ArgsInfo.spacing_given) {
554 for(unsigned int i=0; i< Dimension; i++)
555 outputSpacing[i]=m_ArgsInfo.spacing_arg[i];
556 } else outputSpacing=input->GetSpacing();
557 if (autoGaussEnabled) { // Automated sigma when downsample
558 for(unsigned int i=0; i<Dimension; i++) {
559 if (outputSpacing[i] > input->GetSpacing()[i]) { // downsample
560 gaussianSigma[i] = 0.5*outputSpacing[i];// / inputSpacing[i]);
562 else gaussianSigma[i] = 0; // will be ignore after
567 typename OutputImageType::PointType outputOrigin;
568 if (m_ArgsInfo.origin_given) {
569 for(unsigned int i=0; i< Dimension; i++)
570 outputOrigin[i]=m_ArgsInfo.origin_arg[i];
571 } else outputOrigin=input->GetOrigin();
574 typename OutputImageType::DirectionType outputDirection;
575 if (m_ArgsInfo.direction_given) {
576 for(unsigned int j=0; j< Dimension; j++)
577 for(unsigned int i=0; i< Dimension; i++)
578 outputDirection[j][i]=m_ArgsInfo.direction_arg[i+Dimension*j];
579 } else outputDirection=input->GetDirection();
582 resampler->SetSize( outputSize );
583 resampler->SetOutputSpacing( outputSpacing );
584 resampler->SetOutputOrigin( outputOrigin );
585 resampler->SetOutputDirection( outputDirection );
589 if (m_ArgsInfo.spacinglike_given) {
590 typename InputReaderType::Pointer likeReader=InputReaderType::New();
591 likeReader->SetFileName(m_ArgsInfo.spacinglike_arg);
592 likeReader->Update();
594 // set the support like the image
595 if (m_ArgsInfo.like_given) {
596 typename OutputImageType::SizeType outputSize;
597 outputSize[0] = ceil(resampler->GetSize()[0]*resampler->GetOutputSpacing()[0]
598 /likeReader->GetOutput()->GetSpacing()[0]);
599 outputSize[1] = ceil(resampler->GetSize()[1]*resampler->GetOutputSpacing()[1]
600 /likeReader->GetOutput()->GetSpacing()[1]);
601 outputSize[2] = ceil(resampler->GetSize()[2]*resampler->GetOutputSpacing()[2]
602 /likeReader->GetOutput()->GetSpacing()[2]);
603 if (m_ArgsInfo.verbose_flag) {
604 std::cout << "Compute the number of pixels such as the support is like " << m_ArgsInfo.like_arg << std::endl;
606 resampler->SetSize( outputSize );
609 resampler->SetOutputSpacing( likeReader->GetOutput()->GetSpacing() );
612 if (m_ArgsInfo.verbose_flag) {
613 std::cout << "Setting the output size to " << resampler->GetSize() << "..." << std::endl;
614 std::cout << "Setting the output spacing to " << resampler->GetOutputSpacing() << "..." << std::endl;
615 std::cout << "Setting the output origin to " << resampler->GetOutputOrigin() << "..." << std::endl;
616 std::cout << "Setting the output direction to " << resampler->GetOutputDirection() << "..." << std::endl;
619 typedef itk::RecursiveGaussianImageFilter<InputImageType, InputImageType> GaussianFilterType;
620 std::vector<typename GaussianFilterType::Pointer> gaussianFilters;
621 if (gaussianFilteringEnabled || autoGaussEnabled) {
622 for(unsigned int i=0; i<Dimension; i++) {
623 if (gaussianSigma[i] != 0) {
624 gaussianFilters.push_back(GaussianFilterType::New());
625 gaussianFilters[i]->SetDirection(i);
626 gaussianFilters[i]->SetOrder(GaussianFilterType::ZeroOrder);
627 gaussianFilters[i]->SetNormalizeAcrossScale(false);
628 gaussianFilters[i]->SetSigma(gaussianSigma[i]); // in millimeter !
629 if (gaussianFilters.size() == 1) { // first
630 gaussianFilters[0]->SetInput(input);
632 gaussianFilters[i]->SetInput(gaussianFilters[i-1]->GetOutput());
636 if (gaussianFilters.size() > 0) {
637 resampler->SetInput(gaussianFilters[gaussianFilters.size()-1]->GetOutput());
638 } else resampler->SetInput(input);
639 } else resampler->SetInput(input);
641 resampler->SetTransform( affineTransform );
642 resampler->SetInterpolator( genericInterpolator->GetInterpolatorPointer());
643 resampler->SetDefaultPixelValue( static_cast<PixelType>(m_ArgsInfo.pad_arg) );
647 } catch(itk::ExceptionObject) {
648 std::cerr<<"Error resampling the image"<<std::endl;
651 typename OutputImageType::Pointer output = resampler->GetOutput();
654 typedef itk::ImageFileWriter<OutputImageType> WriterType;
655 typename WriterType::Pointer writer = WriterType::New();
656 writer->SetFileName(m_ArgsInfo.output_arg);
657 writer->SetInput(output);
661 //-------------------------------------------------------------------
664 //-------------------------------------------------------------------
665 // Update with the number of dimensions and the pixeltype (components)
666 //-------------------------------------------------------------------
667 template<class args_info_type>
668 template<unsigned int Dimension, class PixelType>
669 void AffineTransformGenericFilter<args_info_type>::UpdateWithDimAndVectorType()
672 typedef itk::Image<PixelType, Dimension> InputImageType;
673 typedef itk::Image<PixelType, Dimension> OutputImageType;
676 typedef itk::ImageFileReader<InputImageType> InputReaderType;
677 typename InputReaderType::Pointer reader = InputReaderType::New();
678 reader->SetFileName( m_InputFileName);
680 typename InputImageType::Pointer input= reader->GetOutput();
682 //Gaussian pre-filtering
683 typename itk::Vector<double, Dimension> gaussianSigma;
684 gaussianSigma.Fill(0);
685 bool gaussianFilteringEnabled(false);
686 bool autoGaussEnabled(false);
687 if (m_ArgsInfo.autogauss_given) { // Gaussian filter auto
688 autoGaussEnabled = m_ArgsInfo.autogauss_flag;
690 if (m_ArgsInfo.gauss_given) { // Gaussian filter set by user
691 gaussianFilteringEnabled = true;
692 if (m_ArgsInfo.gauss_given == 1)
694 for (unsigned int i=0; i<Dimension; i++)
696 gaussianSigma[i] = m_ArgsInfo.gauss_arg[0];
699 else if (m_ArgsInfo.gauss_given == Dimension)
701 for (unsigned int i=0; i<Dimension; i++)
703 gaussianSigma[i] = m_ArgsInfo.gauss_arg[i];
708 std::cerr << "Gaussian sigma dimension is incorrect" << std::endl;
714 typedef itk::VectorResampleImageFilter< InputImageType,OutputImageType, double > ResampleFilterType;
715 typename ResampleFilterType::Pointer resampler = ResampleFilterType::New();
718 typename itk::Matrix<double, Dimension+1, Dimension+1> matrix;
719 if (m_ArgsInfo.rotate_given || m_ArgsInfo.translate_given)
721 if (m_ArgsInfo.matrix_given)
723 std::cerr << "You must use either rotate/translate or matrix options" << std::endl;
726 itk::Array<double> transformParameters(2 * Dimension);
727 transformParameters.Fill(0.0);
728 if (m_ArgsInfo.rotate_given)
731 transformParameters[0] = m_ArgsInfo.rotate_arg[0];
733 for (unsigned int i = 0; i < 3; i++)
734 transformParameters[i] = m_ArgsInfo.rotate_arg[i];
736 if (m_ArgsInfo.translate_given)
741 for (unsigned int i = 0; i < Dimension && i < 3; i++)
742 transformParameters[pos++] = m_ArgsInfo.translate_arg[i];
746 matrix.SetIdentity();
747 itk::Matrix<double, 4, 4> tmp = GetForwardAffineMatrix3D(transformParameters);
748 for (unsigned int i = 0; i < 3; ++i)
749 for (unsigned int j = 0; j < 3; ++j)
750 matrix[i][j] = tmp[i][j];
751 for (unsigned int i = 0; i < 3; ++i)
752 matrix[i][4] = tmp[i][3];
755 matrix = GetForwardAffineMatrix<Dimension>(transformParameters);
759 if (m_ArgsInfo.matrix_given)
761 matrix= clitk::ReadMatrix<Dimension>(m_ArgsInfo.matrix_arg);
762 if (m_Verbose) std::cout << "Reading the matrix..." << std::endl;
765 matrix.SetIdentity();
768 std::cout << "Using the following matrix:" << std::endl
769 << matrix << std::endl;
770 typename itk::Matrix<double, Dimension, Dimension> rotationMatrix = clitk::GetRotationalPartMatrix(matrix);
771 typename itk::Vector<double, Dimension> translationPart = clitk::GetTranslationPartMatrix(matrix);
774 typedef itk::AffineTransform<double, Dimension> AffineTransformType;
775 typename AffineTransformType::Pointer affineTransform=AffineTransformType::New();
776 affineTransform->SetMatrix(rotationMatrix);
777 affineTransform->SetTranslation(translationPart);
780 typedef clitk::GenericVectorInterpolator<args_info_type, InputImageType, double> GenericInterpolatorType;
781 typename GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New();
782 genericInterpolator->SetArgsInfo(m_ArgsInfo);
785 if (m_ArgsInfo.like_given) {
786 typename InputReaderType::Pointer likeReader=InputReaderType::New();
787 likeReader->SetFileName(m_ArgsInfo.like_arg);
788 likeReader->Update();
789 resampler->SetSize( likeReader->GetOutput()->GetLargestPossibleRegion().GetSize() );
790 resampler->SetOutputSpacing( likeReader->GetOutput()->GetSpacing() );
791 resampler->SetOutputOrigin( likeReader->GetOutput()->GetOrigin() );
792 resampler->SetOutputDirection( likeReader->GetOutput()->GetDirection() );
793 if (autoGaussEnabled) { // Automated sigma when downsample
794 for(unsigned int i=0; i<Dimension; i++) {
795 if (likeReader->GetOutput()->GetSpacing()[i] > input->GetSpacing()[i]) { // downsample
796 gaussianSigma[i] = 0.5*likeReader->GetOutput()->GetSpacing()[i];// / inputSpacing[i]);
798 else gaussianSigma[i] = 0; // will be ignore after
803 typename OutputImageType::SizeType outputSize;
804 if (m_ArgsInfo.size_given) {
805 for(unsigned int i=0; i< Dimension; i++)
806 outputSize[i]=m_ArgsInfo.size_arg[i];
807 } else outputSize=input->GetLargestPossibleRegion().GetSize();
808 std::cout<<"Setting the size to "<<outputSize<<"..."<<std::endl;
811 typename OutputImageType::SpacingType outputSpacing;
812 if (m_ArgsInfo.spacing_given) {
813 for(unsigned int i=0; i< Dimension; i++)
814 outputSpacing[i]=m_ArgsInfo.spacing_arg[i];
815 } else outputSpacing=input->GetSpacing();
816 if (autoGaussEnabled) { // Automated sigma when downsample
817 for(unsigned int i=0; i<Dimension; i++) {
818 if (outputSpacing[i] > input->GetSpacing()[i]) { // downsample
819 gaussianSigma[i] = 0.5*outputSpacing[i];// / inputSpacing[i]);
821 else gaussianSigma[i] = 0; // will be ignore after
824 std::cout<<"Setting the spacing to "<<outputSpacing<<"..."<<std::endl;
827 typename OutputImageType::PointType outputOrigin;
828 if (m_ArgsInfo.origin_given) {
829 for(unsigned int i=0; i< Dimension; i++)
830 outputOrigin[i]=m_ArgsInfo.origin_arg[i];
831 } else outputOrigin=input->GetOrigin();
832 std::cout<<"Setting the origin to "<<outputOrigin<<"..."<<std::endl;
835 typename OutputImageType::DirectionType outputDirection;
836 if (m_ArgsInfo.direction_given) {
837 for(unsigned int j=0; j< Dimension; j++)
838 for(unsigned int i=0; i< Dimension; i++)
839 outputDirection[j][i]=m_ArgsInfo.direction_arg[i+Dimension*j];
840 } else outputDirection=input->GetDirection();
841 std::cout<<"Setting the direction to "<<outputDirection<<"..."<<std::endl;
844 resampler->SetSize( outputSize );
845 resampler->SetOutputSpacing( outputSpacing );
846 resampler->SetOutputOrigin( outputOrigin );
847 resampler->SetOutputDirection( outputDirection );
851 typedef itk::RecursiveGaussianImageFilter<InputImageType, InputImageType> GaussianFilterType;
852 std::vector<typename GaussianFilterType::Pointer> gaussianFilters;
853 if (gaussianFilteringEnabled || autoGaussEnabled) {
854 for(unsigned int i=0; i<Dimension; i++) {
855 if (gaussianSigma[i] != 0) {
856 gaussianFilters.push_back(GaussianFilterType::New());
857 gaussianFilters[i]->SetDirection(i);
858 gaussianFilters[i]->SetOrder(GaussianFilterType::ZeroOrder);
859 gaussianFilters[i]->SetNormalizeAcrossScale(false);
860 gaussianFilters[i]->SetSigma(gaussianSigma[i]); // in millimeter !
861 if (gaussianFilters.size() == 1) { // first
862 gaussianFilters[0]->SetInput(input);
864 gaussianFilters[i]->SetInput(gaussianFilters[i-1]->GetOutput());
868 if (gaussianFilters.size() > 0) {
869 resampler->SetInput(gaussianFilters[gaussianFilters.size()-1]->GetOutput());
870 } else resampler->SetInput(input);
871 } else resampler->SetInput(input);
873 resampler->SetInput( input );
874 resampler->SetTransform( affineTransform );
875 resampler->SetInterpolator( genericInterpolator->GetInterpolatorPointer());
876 resampler->SetDefaultPixelValue( static_cast<PixelType>(m_ArgsInfo.pad_arg) );
880 } catch(itk::ExceptionObject) {
881 std::cerr<<"Error resampling the image"<<std::endl;
884 typename OutputImageType::Pointer output = resampler->GetOutput();
887 typedef itk::ImageFileWriter<OutputImageType> WriterType;
888 typename WriterType::Pointer writer = WriterType::New();
889 writer->SetFileName(m_ArgsInfo.output_arg);
890 writer->SetInput(output);
894 //-------------------------------------------------------------------
898 #endif //#define clitkAffineTransformGenericFilter_txx