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>();
98 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
99 UpdateWithDimAndPixelType<Dimension, float>();
103 else if (Components==3) {
104 if (m_Verbose) std::cout << "Launching transform in "<< Dimension <<"D and 3D float (DVF)" << std::endl;
105 UpdateWithDimAndVectorType<Dimension, itk::Vector<float, Dimension> >();
108 else std::cerr<<"Number of components is "<<Components<<", not supported!"<<std::endl;
111 //-------------------------------------------------------------------
114 //-------------------------------------------------------------------
115 // Update with the number of dimensions and the pixeltype
116 //-------------------------------------------------------------------
117 template<class args_info_type>
118 template <unsigned int Dimension, class PixelType>
120 AffineTransformGenericFilter<args_info_type>::UpdateWithDimAndPixelType()
124 typedef itk::Image<PixelType, Dimension> InputImageType;
125 typedef itk::Image<PixelType, Dimension> OutputImageType;
128 typedef itk::ImageFileReader<InputImageType> InputReaderType;
129 typename InputReaderType::Pointer reader = InputReaderType::New();
130 reader->SetFileName( m_InputFileName);
132 typename InputImageType::Pointer input= reader->GetOutput();
134 //Adaptative size, spacing origin (use previous clitkResampleImage)
135 if (m_ArgsInfo.adaptive_given) {
137 typedef clitk::ResampleImageWithOptionsFilter<InputImageType, OutputImageType> ResampleImageFilterType;
138 typename ResampleImageFilterType::Pointer filter = ResampleImageFilterType::New();
139 filter->SetInput(input);
142 filter->SetVerboseOptions(m_ArgsInfo.verbose_flag);
144 // Set size / spacing
145 static const unsigned int dim = OutputImageType::ImageDimension;
146 typename OutputImageType::SpacingType spacing;
147 typename OutputImageType::SizeType size;
148 typename OutputImageType::PointType origin;
149 typename OutputImageType::DirectionType direction;
151 if (m_ArgsInfo.like_given) {
152 itk::ImageIOBase::Pointer header = clitk::readImageHeader(m_ArgsInfo.like_arg);
154 for(unsigned int i=0; i<dim; i++){
155 spacing[i] = header->GetSpacing(i);
156 size[i] = header->GetDimensions(i);
157 origin[i] = header->GetOrigin(i);
159 for(unsigned int i=0; i<dim; i++) {
160 for(unsigned int j=0;j<dim;j++) {
161 direction(i,j) = header->GetDirection(i)[j];
164 filter->SetOutputSpacing(spacing);
165 filter->SetOutputSize(size);
166 filter->SetOutputOrigin(origin);
167 filter->SetOutputDirection(direction);
170 std::cerr << "*** Warning : I could not read '" << m_ArgsInfo.like_arg << "' ***" << std::endl;
175 if (m_ArgsInfo.spacing_given == 1) {
176 filter->SetOutputIsoSpacing(m_ArgsInfo.spacing_arg[0]);
178 else if ((m_ArgsInfo.spacing_given != 0) && (m_ArgsInfo.size_given != 0)) {
179 std::cerr << "Error: use spacing or size, not both." << std::endl;
182 else if (m_ArgsInfo.spacing_given) {
183 if ((m_ArgsInfo.spacing_given != 0) && (m_ArgsInfo.spacing_given != dim)) {
184 std::cerr << "Error: spacing should have one or " << dim << " values." << std::endl;
187 for(unsigned int i=0; i<dim; i++)
188 spacing[i] = m_ArgsInfo.spacing_arg[i];
189 filter->SetOutputSpacing(spacing);
191 else if (m_ArgsInfo.size_given) {
192 if ((m_ArgsInfo.size_given != 0) && (m_ArgsInfo.size_given != dim)) {
193 std::cerr << "Error: size should have " << dim << " values." << std::endl;
196 for(unsigned int i=0; i<dim; i++)
197 size[i] = m_ArgsInfo.size_arg[i];
198 filter->SetOutputSize(size);
200 for(unsigned int i=0; i<dim; i++){
201 origin[i] = input->GetOrigin()[i];
203 for(unsigned int i=0; i<dim; i++) {
204 for(unsigned int j=0;j<dim;j++) {
205 direction(i,j) = input->GetDirection()[i][j];
208 filter->SetOutputOrigin(origin);
209 filter->SetOutputDirection(direction);
212 // Set temporal dimension
213 //filter->SetLastDimensionIsTime(m_ArgsInfo.time_flag);
216 filter->SetGaussianFilteringEnabled(m_ArgsInfo.autogauss_flag);
217 if (m_ArgsInfo.gauss_given != 0) {
218 typename ResampleImageFilterType::GaussianSigmaType g;
219 for(unsigned int i=0; i<dim; i++) {
220 g[i] = m_ArgsInfo.gauss_arg[i];
222 filter->SetGaussianSigma(g);
226 int interp = m_ArgsInfo.interp_arg;
228 filter->SetInterpolationType(ResampleImageFilterType::NearestNeighbor);
231 filter->SetInterpolationType(ResampleImageFilterType::Linear);
234 filter->SetInterpolationType(ResampleImageFilterType::BSpline);
237 filter->SetInterpolationType(ResampleImageFilterType::B_LUT);
239 std::cerr << "Error. I do not know interpolation '" << m_ArgsInfo.interp_arg
240 << "'. Choose among: nn, linear, bspline, blut, windowed sinc" << std::endl;
247 // Set default pixel value
248 filter->SetDefaultPixelValue(m_ArgsInfo.pad_arg);
251 //if (m_ArgsInfo.thread_given) {
252 // filter->SetNumberOfThreads(m_ArgsInfo.thread_arg);
257 typename OutputImageType::Pointer output = filter->GetOutput();
258 //this->template SetNextOutput<OutputImageType>(outputImage);
261 typedef itk::ImageFileWriter<OutputImageType> WriterType;
262 typename WriterType::Pointer writer = WriterType::New();
263 writer->SetFileName(m_ArgsInfo.output_arg);
264 writer->SetInput(output);
270 //Gaussian pre-filtering
271 typename itk::Vector<double, Dimension> gaussianSigma;
272 gaussianSigma.Fill(0);
273 bool gaussianFilteringEnabled(false);
274 bool autoGaussEnabled(false);
275 if (m_ArgsInfo.autogauss_given) { // Gaussian filter auto
276 autoGaussEnabled = m_ArgsInfo.autogauss_flag;
278 if (m_ArgsInfo.gauss_given) { // Gaussian filter set by user
279 gaussianFilteringEnabled = true;
280 if (m_ArgsInfo.gauss_given == 1)
282 for (unsigned int i=0; i<Dimension; i++)
284 gaussianSigma[i] = m_ArgsInfo.gauss_arg[0];
287 else if (m_ArgsInfo.gauss_given == Dimension)
289 for (unsigned int i=0; i<Dimension; i++)
291 gaussianSigma[i] = m_ArgsInfo.gauss_arg[i];
296 std::cerr << "Gaussian sigma dimension is incorrect" << std::endl;
302 typedef itk::ResampleImageFilter< InputImageType,OutputImageType > ResampleFilterType;
303 typename ResampleFilterType::Pointer resampler = ResampleFilterType::New();
306 typename itk::Matrix<double, Dimension+1, Dimension+1> matrix;
307 if (m_ArgsInfo.rotate_given || m_ArgsInfo.translate_given)
309 if (m_ArgsInfo.matrix_given)
311 std::cerr << "You must use either rotate/translate or matrix options" << std::endl;
314 itk::Array<double> transformParameters(2 * Dimension);
315 transformParameters.Fill(0.0);
316 if (m_ArgsInfo.rotate_given)
319 transformParameters[0] = m_ArgsInfo.rotate_arg[0];
321 for (unsigned int i = 0; i < 3; i++)
322 transformParameters[i] = m_ArgsInfo.rotate_arg[i];
324 if (m_ArgsInfo.translate_given)
329 for (unsigned int i = 0; i < Dimension && i < 3; i++)
330 transformParameters[pos++] = m_ArgsInfo.translate_arg[i];
334 matrix.SetIdentity();
335 itk::Matrix<double, 4, 4> tmp = GetForwardAffineMatrix3D(transformParameters);
336 for (unsigned int i = 0; i < 3; ++i)
337 for (unsigned int j = 0; j < 3; ++j)
338 matrix[i][j] = tmp[i][j];
339 for (unsigned int i = 0; i < 3; ++i)
340 matrix[i][4] = tmp[i][3];
343 matrix = GetForwardAffineMatrix<Dimension>(transformParameters);
347 if (m_ArgsInfo.matrix_given)
349 matrix= clitk::ReadMatrix<Dimension>(m_ArgsInfo.matrix_arg);
350 if (m_Verbose) std::cout << "Reading the matrix..." << std::endl;
353 if (m_ArgsInfo.elastix_given) {
354 std::string filename(m_ArgsInfo.elastix_arg);
355 matrix = createMatrixFromElastixFile<Dimension>(filename, m_Verbose);
358 matrix.SetIdentity();
362 std::cout << "Using the following matrix:" << std::endl
363 << matrix << std::endl;
364 typename itk::Matrix<double, Dimension, Dimension> rotationMatrix = clitk::GetRotationalPartMatrix(matrix);
365 typename itk::Vector<double, Dimension> translationPart = clitk::GetTranslationPartMatrix(matrix);
368 typedef itk::AffineTransform<double, Dimension> AffineTransformType;
369 typename AffineTransformType::Pointer affineTransform=AffineTransformType::New();
370 affineTransform->SetMatrix(rotationMatrix);
371 affineTransform->SetTranslation(translationPart);
374 typedef clitk::GenericInterpolator<args_info_type, InputImageType, double> GenericInterpolatorType;
375 typename GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New();
376 genericInterpolator->SetArgsInfo(m_ArgsInfo);
379 if (m_ArgsInfo.like_given) {
380 typename InputReaderType::Pointer likeReader=InputReaderType::New();
381 likeReader->SetFileName(m_ArgsInfo.like_arg);
382 likeReader->Update();
383 resampler->SetOutputParametersFromImage(likeReader->GetOutput());
384 resampler->SetOutputDirection(likeReader->GetOutput()->GetDirection());
385 if (autoGaussEnabled) { // Automated sigma when downsample
386 for(unsigned int i=0; i<Dimension; i++) {
387 if (likeReader->GetOutput()->GetSpacing()[i] > input->GetSpacing()[i]) { // downsample
388 gaussianSigma[i] = 0.5*likeReader->GetOutput()->GetSpacing()[i];// / inputSpacing[i]);
390 else gaussianSigma[i] = 0; // will be ignore after
393 } else if(m_ArgsInfo.transform_grid_flag) {
394 typename itk::Matrix<double, Dimension+1, Dimension+1> invMatrix( matrix.GetInverse() );
395 typename itk::Matrix<double, Dimension, Dimension> invRotMatrix( clitk::GetRotationalPartMatrix(invMatrix) );
396 typename itk::Vector<double,Dimension> invTrans = clitk::GetTranslationPartMatrix(invMatrix);
399 if (m_ArgsInfo.spacing_given)
400 std::cout << "Warning --spacing ignored (because --transform_grid_flag)" << std::endl;
401 if (m_ArgsInfo.origin_given)
402 std::cout << "Warning --origin ignored (because --transform_grid_flag)" << std::endl;
404 // Spacing is influenced by affine transform matrix and input direction
405 typename InputImageType::SpacingType outputSpacing;
406 outputSpacing = invRotMatrix *
407 input->GetDirection() *
409 if (autoGaussEnabled) { // Automated sigma when downsample
410 for(unsigned int i=0; i<Dimension; i++) {
411 if (outputSpacing[i] > input->GetSpacing()[i]) { // downsample
412 gaussianSigma[i] = 0.5*outputSpacing[i];// / inputSpacing[i]);
414 else gaussianSigma[i] = 0; // will be ignore after
418 // Origin is influenced by translation but not by input direction
419 typename InputImageType::PointType outputOrigin;
420 outputOrigin = invRotMatrix *
424 // Size is influenced by affine transform matrix and input direction
425 // Size is converted to double, transformed and converted back to size type.
426 vnl_vector<double> vnlOutputSize(Dimension);
427 for(unsigned int i=0; i< Dimension; i++) {
428 vnlOutputSize[i] = input->GetLargestPossibleRegion().GetSize()[i];
430 vnlOutputSize = invRotMatrix *
431 input->GetDirection().GetVnlMatrix() *
433 typename OutputImageType::SizeType outputSize;
434 for(unsigned int i=0; i< Dimension; i++) {
435 // If the size is negative, we have a flip and we must modify
436 // the origin and the spacing accordingly.
437 if(vnlOutputSize[i]<0.) {
438 vnlOutputSize[i] *= -1.;
439 outputOrigin[i] = outputOrigin[i] + outputSpacing[i] * (vnlOutputSize[i]-1);
440 outputSpacing[i] *= -1.;
442 outputSize[i] = lrint(vnlOutputSize[i]);
444 resampler->SetSize( outputSize );
445 resampler->SetOutputSpacing( outputSpacing );
446 resampler->SetOutputOrigin( outputOrigin );
449 typename OutputImageType::SizeType outputSize;
450 if (m_ArgsInfo.size_given) {
451 for(unsigned int i=0; i< Dimension; i++)
452 outputSize[i]=m_ArgsInfo.size_arg[i];
453 } else outputSize=input->GetLargestPossibleRegion().GetSize();
456 typename OutputImageType::SpacingType outputSpacing;
457 if (m_ArgsInfo.spacing_given) {
458 for(unsigned int i=0; i< Dimension; i++)
459 outputSpacing[i]=m_ArgsInfo.spacing_arg[i];
460 } else outputSpacing=input->GetSpacing();
461 if (autoGaussEnabled) { // Automated sigma when downsample
462 for(unsigned int i=0; i<Dimension; i++) {
463 if (outputSpacing[i] > input->GetSpacing()[i]) { // downsample
464 gaussianSigma[i] = 0.5*outputSpacing[i];// / inputSpacing[i]);
466 else gaussianSigma[i] = 0; // will be ignore after
471 typename OutputImageType::PointType outputOrigin;
472 if (m_ArgsInfo.origin_given) {
473 for(unsigned int i=0; i< Dimension; i++)
474 outputOrigin[i]=m_ArgsInfo.origin_arg[i];
475 } else outputOrigin=input->GetOrigin();
478 typename OutputImageType::DirectionType outputDirection;
479 if (m_ArgsInfo.direction_given) {
480 for(unsigned int j=0; j< Dimension; j++)
481 for(unsigned int i=0; i< Dimension; i++)
482 outputDirection[j][i]=m_ArgsInfo.direction_arg[i+Dimension*j];
483 } else outputDirection=input->GetDirection();
486 resampler->SetSize( outputSize );
487 resampler->SetOutputSpacing( outputSpacing );
488 resampler->SetOutputOrigin( outputOrigin );
489 resampler->SetOutputDirection( outputDirection );
493 if (m_ArgsInfo.spacinglike_given) {
494 typename InputReaderType::Pointer likeReader=InputReaderType::New();
495 likeReader->SetFileName(m_ArgsInfo.spacinglike_arg);
496 likeReader->Update();
498 // set the support like the image
499 if (m_ArgsInfo.like_given) {
500 typename OutputImageType::SizeType outputSize;
501 outputSize[0] = ceil(resampler->GetSize()[0]*resampler->GetOutputSpacing()[0]
502 /likeReader->GetOutput()->GetSpacing()[0]);
503 outputSize[1] = ceil(resampler->GetSize()[1]*resampler->GetOutputSpacing()[1]
504 /likeReader->GetOutput()->GetSpacing()[1]);
505 outputSize[2] = ceil(resampler->GetSize()[2]*resampler->GetOutputSpacing()[2]
506 /likeReader->GetOutput()->GetSpacing()[2]);
507 if (m_ArgsInfo.verbose_flag) {
508 std::cout << "Compute the number of pixels such as the support is like " << m_ArgsInfo.like_arg << std::endl;
510 resampler->SetSize( outputSize );
513 resampler->SetOutputSpacing( likeReader->GetOutput()->GetSpacing() );
516 if (m_ArgsInfo.verbose_flag) {
517 std::cout << "Setting the output size to " << resampler->GetSize() << "..." << std::endl;
518 std::cout << "Setting the output spacing to " << resampler->GetOutputSpacing() << "..." << std::endl;
519 std::cout << "Setting the output origin to " << resampler->GetOutputOrigin() << "..." << std::endl;
520 std::cout << "Setting the output direction to " << resampler->GetOutputDirection() << "..." << std::endl;
523 typedef itk::RecursiveGaussianImageFilter<InputImageType, InputImageType> GaussianFilterType;
524 std::vector<typename GaussianFilterType::Pointer> gaussianFilters;
525 if (gaussianFilteringEnabled || autoGaussEnabled) {
526 for(unsigned int i=0; i<Dimension; i++) {
527 if (gaussianSigma[i] != 0) {
528 gaussianFilters.push_back(GaussianFilterType::New());
529 gaussianFilters[i]->SetDirection(i);
530 gaussianFilters[i]->SetOrder(GaussianFilterType::ZeroOrder);
531 gaussianFilters[i]->SetNormalizeAcrossScale(false);
532 gaussianFilters[i]->SetSigma(gaussianSigma[i]); // in millimeter !
533 if (gaussianFilters.size() == 1) { // first
534 gaussianFilters[0]->SetInput(input);
536 gaussianFilters[i]->SetInput(gaussianFilters[i-1]->GetOutput());
540 if (gaussianFilters.size() > 0) {
541 resampler->SetInput(gaussianFilters[gaussianFilters.size()-1]->GetOutput());
542 } else resampler->SetInput(input);
543 } else resampler->SetInput(input);
545 resampler->SetTransform( affineTransform );
546 resampler->SetInterpolator( genericInterpolator->GetInterpolatorPointer());
547 resampler->SetDefaultPixelValue( static_cast<PixelType>(m_ArgsInfo.pad_arg) );
551 } catch(itk::ExceptionObject) {
552 std::cerr<<"Error resampling the image"<<std::endl;
555 typename OutputImageType::Pointer output = resampler->GetOutput();
558 typedef itk::ImageFileWriter<OutputImageType> WriterType;
559 typename WriterType::Pointer writer = WriterType::New();
560 writer->SetFileName(m_ArgsInfo.output_arg);
561 writer->SetInput(output);
565 //-------------------------------------------------------------------
568 //-------------------------------------------------------------------
569 // Update with the number of dimensions and the pixeltype (components)
570 //-------------------------------------------------------------------
571 template<class args_info_type>
572 template<unsigned int Dimension, class PixelType>
573 void AffineTransformGenericFilter<args_info_type>::UpdateWithDimAndVectorType()
576 typedef itk::Image<PixelType, Dimension> InputImageType;
577 typedef itk::Image<PixelType, Dimension> OutputImageType;
580 typedef itk::ImageFileReader<InputImageType> InputReaderType;
581 typename InputReaderType::Pointer reader = InputReaderType::New();
582 reader->SetFileName( m_InputFileName);
584 typename InputImageType::Pointer input= reader->GetOutput();
586 //Gaussian pre-filtering
587 typename itk::Vector<double, Dimension> gaussianSigma;
588 gaussianSigma.Fill(0);
589 bool gaussianFilteringEnabled(false);
590 bool autoGaussEnabled(false);
591 if (m_ArgsInfo.autogauss_given) { // Gaussian filter auto
592 autoGaussEnabled = m_ArgsInfo.autogauss_flag;
594 if (m_ArgsInfo.gauss_given) { // Gaussian filter set by user
595 gaussianFilteringEnabled = true;
596 if (m_ArgsInfo.gauss_given == 1)
598 for (unsigned int i=0; i<Dimension; i++)
600 gaussianSigma[i] = m_ArgsInfo.gauss_arg[0];
603 else if (m_ArgsInfo.gauss_given == Dimension)
605 for (unsigned int i=0; i<Dimension; i++)
607 gaussianSigma[i] = m_ArgsInfo.gauss_arg[i];
612 std::cerr << "Gaussian sigma dimension is incorrect" << std::endl;
618 typedef itk::VectorResampleImageFilter< InputImageType,OutputImageType, double > ResampleFilterType;
619 typename ResampleFilterType::Pointer resampler = ResampleFilterType::New();
622 typename itk::Matrix<double, Dimension+1, Dimension+1> matrix;
623 if (m_ArgsInfo.rotate_given || m_ArgsInfo.translate_given)
625 if (m_ArgsInfo.matrix_given)
627 std::cerr << "You must use either rotate/translate or matrix options" << std::endl;
630 itk::Array<double> transformParameters(2 * Dimension);
631 transformParameters.Fill(0.0);
632 if (m_ArgsInfo.rotate_given)
635 transformParameters[0] = m_ArgsInfo.rotate_arg[0];
637 for (unsigned int i = 0; i < 3; i++)
638 transformParameters[i] = m_ArgsInfo.rotate_arg[i];
640 if (m_ArgsInfo.translate_given)
645 for (unsigned int i = 0; i < Dimension && i < 3; i++)
646 transformParameters[pos++] = m_ArgsInfo.translate_arg[i];
650 matrix.SetIdentity();
651 itk::Matrix<double, 4, 4> tmp = GetForwardAffineMatrix3D(transformParameters);
652 for (unsigned int i = 0; i < 3; ++i)
653 for (unsigned int j = 0; j < 3; ++j)
654 matrix[i][j] = tmp[i][j];
655 for (unsigned int i = 0; i < 3; ++i)
656 matrix[i][4] = tmp[i][3];
659 matrix = GetForwardAffineMatrix<Dimension>(transformParameters);
663 if (m_ArgsInfo.matrix_given)
665 matrix= clitk::ReadMatrix<Dimension>(m_ArgsInfo.matrix_arg);
666 if (m_Verbose) std::cout << "Reading the matrix..." << std::endl;
669 matrix.SetIdentity();
672 std::cout << "Using the following matrix:" << std::endl
673 << matrix << std::endl;
674 typename itk::Matrix<double, Dimension, Dimension> rotationMatrix = clitk::GetRotationalPartMatrix(matrix);
675 typename itk::Vector<double, Dimension> translationPart = clitk::GetTranslationPartMatrix(matrix);
678 typedef itk::AffineTransform<double, Dimension> AffineTransformType;
679 typename AffineTransformType::Pointer affineTransform=AffineTransformType::New();
680 affineTransform->SetMatrix(rotationMatrix);
681 affineTransform->SetTranslation(translationPart);
684 typedef clitk::GenericVectorInterpolator<args_info_type, InputImageType, double> GenericInterpolatorType;
685 typename GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New();
686 genericInterpolator->SetArgsInfo(m_ArgsInfo);
689 if (m_ArgsInfo.like_given) {
690 typename InputReaderType::Pointer likeReader=InputReaderType::New();
691 likeReader->SetFileName(m_ArgsInfo.like_arg);
692 likeReader->Update();
693 resampler->SetSize( likeReader->GetOutput()->GetLargestPossibleRegion().GetSize() );
694 resampler->SetOutputSpacing( likeReader->GetOutput()->GetSpacing() );
695 resampler->SetOutputOrigin( likeReader->GetOutput()->GetOrigin() );
696 resampler->SetOutputDirection( likeReader->GetOutput()->GetDirection() );
697 if (autoGaussEnabled) { // Automated sigma when downsample
698 for(unsigned int i=0; i<Dimension; i++) {
699 if (likeReader->GetOutput()->GetSpacing()[i] > input->GetSpacing()[i]) { // downsample
700 gaussianSigma[i] = 0.5*likeReader->GetOutput()->GetSpacing()[i];// / inputSpacing[i]);
702 else gaussianSigma[i] = 0; // will be ignore after
707 typename OutputImageType::SizeType outputSize;
708 if (m_ArgsInfo.size_given) {
709 for(unsigned int i=0; i< Dimension; i++)
710 outputSize[i]=m_ArgsInfo.size_arg[i];
711 } else outputSize=input->GetLargestPossibleRegion().GetSize();
712 std::cout<<"Setting the size to "<<outputSize<<"..."<<std::endl;
715 typename OutputImageType::SpacingType outputSpacing;
716 if (m_ArgsInfo.spacing_given) {
717 for(unsigned int i=0; i< Dimension; i++)
718 outputSpacing[i]=m_ArgsInfo.spacing_arg[i];
719 } else outputSpacing=input->GetSpacing();
720 if (autoGaussEnabled) { // Automated sigma when downsample
721 for(unsigned int i=0; i<Dimension; i++) {
722 if (outputSpacing[i] > input->GetSpacing()[i]) { // downsample
723 gaussianSigma[i] = 0.5*outputSpacing[i];// / inputSpacing[i]);
725 else gaussianSigma[i] = 0; // will be ignore after
728 std::cout<<"Setting the spacing to "<<outputSpacing<<"..."<<std::endl;
731 typename OutputImageType::PointType outputOrigin;
732 if (m_ArgsInfo.origin_given) {
733 for(unsigned int i=0; i< Dimension; i++)
734 outputOrigin[i]=m_ArgsInfo.origin_arg[i];
735 } else outputOrigin=input->GetOrigin();
736 std::cout<<"Setting the origin to "<<outputOrigin<<"..."<<std::endl;
739 typename OutputImageType::DirectionType outputDirection;
740 if (m_ArgsInfo.direction_given) {
741 for(unsigned int j=0; j< Dimension; j++)
742 for(unsigned int i=0; i< Dimension; i++)
743 outputDirection[j][i]=m_ArgsInfo.direction_arg[i+Dimension*j];
744 } else outputDirection=input->GetDirection();
745 std::cout<<"Setting the direction to "<<outputDirection<<"..."<<std::endl;
748 resampler->SetSize( outputSize );
749 resampler->SetOutputSpacing( outputSpacing );
750 resampler->SetOutputOrigin( outputOrigin );
751 resampler->SetOutputDirection( outputDirection );
755 typedef itk::RecursiveGaussianImageFilter<InputImageType, InputImageType> GaussianFilterType;
756 std::vector<typename GaussianFilterType::Pointer> gaussianFilters;
757 if (gaussianFilteringEnabled || autoGaussEnabled) {
758 for(unsigned int i=0; i<Dimension; i++) {
759 if (gaussianSigma[i] != 0) {
760 gaussianFilters.push_back(GaussianFilterType::New());
761 gaussianFilters[i]->SetDirection(i);
762 gaussianFilters[i]->SetOrder(GaussianFilterType::ZeroOrder);
763 gaussianFilters[i]->SetNormalizeAcrossScale(false);
764 gaussianFilters[i]->SetSigma(gaussianSigma[i]); // in millimeter !
765 if (gaussianFilters.size() == 1) { // first
766 gaussianFilters[0]->SetInput(input);
768 gaussianFilters[i]->SetInput(gaussianFilters[i-1]->GetOutput());
772 if (gaussianFilters.size() > 0) {
773 resampler->SetInput(gaussianFilters[gaussianFilters.size()-1]->GetOutput());
774 } else resampler->SetInput(input);
775 } else resampler->SetInput(input);
777 resampler->SetInput( input );
778 resampler->SetTransform( affineTransform );
779 resampler->SetInterpolator( genericInterpolator->GetInterpolatorPointer());
780 resampler->SetDefaultPixelValue( static_cast<PixelType>(m_ArgsInfo.pad_arg) );
784 } catch(itk::ExceptionObject) {
785 std::cerr<<"Error resampling the image"<<std::endl;
788 typename OutputImageType::Pointer output = resampler->GetOutput();
791 typedef itk::ImageFileWriter<OutputImageType> WriterType;
792 typename WriterType::Pointer writer = WriterType::New();
793 writer->SetFileName(m_ArgsInfo.output_arg);
794 writer->SetInput(output);
798 //-------------------------------------------------------------------
802 #endif //#define clitkAffineTransformGenericFilter_txx