/*========================================================================= Program: vv http://www.creatis.insa-lyon.fr/rio/vv Authors belong to: - University of LYON http://www.universite-lyon.fr/ - Léon Bérard cancer center http://www.centreleonberard.fr - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the copyright notices for more information. It is distributed under dual licence - BSD See included LICENSE.txt file - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html ===========================================================================**/ #ifndef clitkAffineTransformGenericFilter_txx #define clitkAffineTransformGenericFilter_txx #include #include #include #include namespace clitk { //----------------------------------------------------------- // Constructor //----------------------------------------------------------- template AffineTransformGenericFilter::AffineTransformGenericFilter() { m_Verbose=false; m_InputFileName=""; } //------------------------------------------------------------------- //----------------------------------------------------------- // Update //----------------------------------------------------------- template void AffineTransformGenericFilter::Update() { // Read the Dimension and PixelType int Dimension, Components; std::string PixelType; ReadImageDimensionAndPixelType(m_InputFileName, Dimension, PixelType, Components); // Call UpdateWithDim if(Dimension==2) UpdateWithDim<2>(PixelType, Components); else if(Dimension==3) UpdateWithDim<3>(PixelType, Components); else if (Dimension==4)UpdateWithDim<4>(PixelType, Components); else { std::cout<<"Error, Only for 2, 3 or 4 Dimensions!!!"< template void AffineTransformGenericFilter::UpdateWithDim(std::string PixelType, int Components) { if (m_Verbose) std::cout << "Image was detected to be "<(); } else if(PixelType == "unsigned_short"){ if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_short..." << std::endl; UpdateWithDimAndPixelType(); } else if (PixelType == "unsigned_char") { if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl; UpdateWithDimAndPixelType(); } // else if (PixelType == "char"){ // if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed_char..." << std::endl; // UpdateWithDimAndPixelType(); // } else { if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and float..." << std::endl; UpdateWithDimAndPixelType(); } } else if (Components==3) { if (m_Verbose) std::cout << "Launching transform in "<< Dimension <<"D and 3D float (DVF)" << std::endl; UpdateWithDimAndVectorType >(); } else std::cerr<<"Number of components is "< template void AffineTransformGenericFilter::UpdateWithDimAndPixelType() { // ImageTypes typedef itk::Image InputImageType; typedef itk::Image OutputImageType; // Read the input typedef itk::ImageFileReader InputReaderType; typename InputReaderType::Pointer reader = InputReaderType::New(); reader->SetFileName( m_InputFileName); reader->Update(); typename InputImageType::Pointer input= reader->GetOutput(); //Filter typedef itk::ResampleImageFilter< InputImageType,OutputImageType > ResampleFilterType; typename ResampleFilterType::Pointer resampler = ResampleFilterType::New(); // Matrix typename itk::Matrix matrix; if (m_ArgsInfo.rotate_given || m_ArgsInfo.translate_given) { if (m_ArgsInfo.matrix_given) { std::cerr << "You must use either rotate/translate or matrix options" << std::cout; return; } itk::Array transformParameters(2 * Dimension); transformParameters.Fill(0.0); if (m_ArgsInfo.rotate_given) { if (Dimension == 2) transformParameters[0] = m_ArgsInfo.rotate_arg[0]; else for (unsigned int i = 0; i < 3; i++) transformParameters[i] = m_ArgsInfo.rotate_arg[i]; } if (m_ArgsInfo.translate_given) { int pos = 3; if (Dimension == 2) pos = 1; for (unsigned int i = 0; i < Dimension && i < 3; i++) transformParameters[pos++] = m_ArgsInfo.translate_arg[i]; } if (Dimension == 4) { matrix.SetIdentity(); itk::Matrix tmp = GetForwardAffineMatrix3D(transformParameters); for (unsigned int i = 0; i < 3; ++i) for (unsigned int j = 0; j < 3; ++j) matrix[i][j] = tmp[i][j]; for (unsigned int i = 0; i < 3; ++i) matrix[i][4] = tmp[i][3]; } else matrix = GetForwardAffineMatrix(transformParameters); } else { if (m_ArgsInfo.matrix_given) { matrix= clitk::ReadMatrix(m_ArgsInfo.matrix_arg); if (m_Verbose) std::cout << "Reading the matrix..." << std::endl; } else { if (m_ArgsInfo.elastix_given) { std::vector s; for(uint i=0; i(s); } else matrix.SetIdentity(); } } if (m_Verbose) std::cout << "Using the following matrix:" << std::endl << matrix << std::endl; typename itk::Matrix rotationMatrix = clitk::GetRotationalPartMatrix(matrix); typename itk::Vector translationPart = clitk::GetTranslationPartMatrix(matrix); // Transform typedef itk::AffineTransform AffineTransformType; typename AffineTransformType::Pointer affineTransform=AffineTransformType::New(); affineTransform->SetMatrix(rotationMatrix); affineTransform->SetTranslation(translationPart); // Interp typedef clitk::GenericInterpolator GenericInterpolatorType; typename GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New(); genericInterpolator->SetArgsInfo(m_ArgsInfo); // Properties if (m_ArgsInfo.like_given) { typename InputReaderType::Pointer likeReader=InputReaderType::New(); likeReader->SetFileName(m_ArgsInfo.like_arg); likeReader->Update(); resampler->SetOutputParametersFromImage(likeReader->GetOutput()); } else if(m_ArgsInfo.transform_grid_flag) { typename itk::Matrix invMatrix( matrix.GetInverse() ); typename itk::Matrix invRotMatrix( clitk::GetRotationalPartMatrix(invMatrix) ); typename itk::Vector invTrans = clitk::GetTranslationPartMatrix(invMatrix); // Display warning if (m_ArgsInfo.spacing_given) std::cout << "Warning --spacing ignored (because --transform_grid_flag)" << std::endl; if (m_ArgsInfo.origin_given) std::cout << "Warning --origin ignored (because --transform_grid_flag)" << std::endl; // Spacing is influenced by affine transform matrix and input direction typename InputImageType::SpacingType outputSpacing; outputSpacing = invRotMatrix * input->GetDirection() * input->GetSpacing(); // Origin is influenced by translation but not by input direction typename InputImageType::PointType outputOrigin; outputOrigin = invRotMatrix * input->GetOrigin() + invTrans; // Size is influenced by affine transform matrix and input direction // Size is converted to double, transformed and converted back to size type. vnl_vector vnlOutputSize(Dimension); for(unsigned int i=0; i< Dimension; i++) { vnlOutputSize[i] = input->GetLargestPossibleRegion().GetSize()[i]; } vnlOutputSize = invRotMatrix * input->GetDirection().GetVnlMatrix() * vnlOutputSize; typename OutputImageType::SizeType outputSize; for(unsigned int i=0; i< Dimension; i++) { // If the size is negative, we have a flip and we must modify // the origin and the spacing accordingly. if(vnlOutputSize[i]<0.) { vnlOutputSize[i] *= -1.; outputOrigin[i] = outputOrigin[i] + outputSpacing[i] * (vnlOutputSize[i]-1); outputSpacing[i] *= -1.; } outputSize[i] = lrint(vnlOutputSize[i]); } resampler->SetSize( outputSize ); resampler->SetOutputSpacing( outputSpacing ); resampler->SetOutputOrigin( outputOrigin ); } else { //Size typename OutputImageType::SizeType outputSize; if (m_ArgsInfo.size_given) { for(unsigned int i=0; i< Dimension; i++) outputSize[i]=m_ArgsInfo.size_arg[i]; } else outputSize=input->GetLargestPossibleRegion().GetSize(); //Spacing typename OutputImageType::SpacingType outputSpacing; if (m_ArgsInfo.spacing_given) { for(unsigned int i=0; i< Dimension; i++) outputSpacing[i]=m_ArgsInfo.spacing_arg[i]; } else outputSpacing=input->GetSpacing(); //Origin typename OutputImageType::PointType outputOrigin; if (m_ArgsInfo.origin_given) { for(unsigned int i=0; i< Dimension; i++) outputOrigin[i]=m_ArgsInfo.origin_arg[i]; } else outputOrigin=input->GetOrigin(); // Set resampler->SetSize( outputSize ); resampler->SetOutputSpacing( outputSpacing ); resampler->SetOutputOrigin( outputOrigin ); } if (m_ArgsInfo.spacinglike_given) { typename InputReaderType::Pointer likeReader=InputReaderType::New(); likeReader->SetFileName(m_ArgsInfo.spacinglike_arg); likeReader->Update(); // set the support like the image if (m_ArgsInfo.like_given) { typename OutputImageType::SizeType outputSize; outputSize[0] = ceil(resampler->GetSize()[0]*resampler->GetOutputSpacing()[0] /likeReader->GetOutput()->GetSpacing()[0]); outputSize[1] = ceil(resampler->GetSize()[1]*resampler->GetOutputSpacing()[1] /likeReader->GetOutput()->GetSpacing()[1]); outputSize[2] = ceil(resampler->GetSize()[2]*resampler->GetOutputSpacing()[2] /likeReader->GetOutput()->GetSpacing()[2]); if (m_ArgsInfo.verbose_flag) { std::cout << "Compute the number of pixels such as the support is like " << m_ArgsInfo.like_arg << std::endl; } resampler->SetSize( outputSize ); } resampler->SetOutputSpacing( likeReader->GetOutput()->GetSpacing() ); } if (m_ArgsInfo.verbose_flag) { std::cout << "Setting the output size to " << resampler->GetSize() << "..." << std::endl; std::cout << "Setting the output spacing to " << resampler->GetOutputSpacing() << "..." << std::endl; std::cout << "Setting the output origin to " << resampler->GetOutputOrigin() << "..." << std::endl; } resampler->SetInput( input ); resampler->SetTransform( affineTransform ); resampler->SetInterpolator( genericInterpolator->GetInterpolatorPointer()); resampler->SetDefaultPixelValue( static_cast(m_ArgsInfo.pad_arg) ); try { resampler->Update(); } catch(itk::ExceptionObject) { std::cerr<<"Error resampling the image"<GetOutput(); // Output typedef itk::ImageFileWriter WriterType; typename WriterType::Pointer writer = WriterType::New(); writer->SetFileName(m_ArgsInfo.output_arg); writer->SetInput(output); writer->Update(); } //------------------------------------------------------------------- //------------------------------------------------------------------- // Update with the number of dimensions and the pixeltype (components) //------------------------------------------------------------------- template template void AffineTransformGenericFilter::UpdateWithDimAndVectorType() { // ImageTypes typedef itk::Image InputImageType; typedef itk::Image OutputImageType; // Read the input typedef itk::ImageFileReader InputReaderType; typename InputReaderType::Pointer reader = InputReaderType::New(); reader->SetFileName( m_InputFileName); reader->Update(); typename InputImageType::Pointer input= reader->GetOutput(); //Filter typedef itk::VectorResampleImageFilter< InputImageType,OutputImageType, double > ResampleFilterType; typename ResampleFilterType::Pointer resampler = ResampleFilterType::New(); // Matrix typename itk::Matrix matrix; if (m_ArgsInfo.rotate_given || m_ArgsInfo.translate_given) { if (m_ArgsInfo.matrix_given) { std::cerr << "You must use either rotate/translate or matrix options" << std::cout; return; } itk::Array transformParameters(2 * Dimension); transformParameters.Fill(0.0); if (m_ArgsInfo.rotate_given) { if (Dimension == 2) transformParameters[0] = m_ArgsInfo.rotate_arg[0]; else for (unsigned int i = 0; i < 3; i++) transformParameters[i] = m_ArgsInfo.rotate_arg[i]; } if (m_ArgsInfo.translate_given) { int pos = 3; if (Dimension == 2) pos = 1; for (unsigned int i = 0; i < Dimension && i < 3; i++) transformParameters[pos++] = m_ArgsInfo.translate_arg[i]; } if (Dimension == 4) { matrix.SetIdentity(); itk::Matrix tmp = GetForwardAffineMatrix3D(transformParameters); for (unsigned int i = 0; i < 3; ++i) for (unsigned int j = 0; j < 3; ++j) matrix[i][j] = tmp[i][j]; for (unsigned int i = 0; i < 3; ++i) matrix[i][4] = tmp[i][3]; } else matrix = GetForwardAffineMatrix(transformParameters); } else { if (m_ArgsInfo.matrix_given) { matrix= clitk::ReadMatrix(m_ArgsInfo.matrix_arg); if (m_Verbose) std::cout << "Reading the matrix..." << std::endl; } else matrix.SetIdentity(); } if (m_Verbose) std::cout << "Using the following matrix:" << std::endl << matrix << std::endl; typename itk::Matrix rotationMatrix = clitk::GetRotationalPartMatrix(matrix); typename itk::Vector translationPart = clitk::GetTranslationPartMatrix(matrix); // Transform typedef itk::AffineTransform AffineTransformType; typename AffineTransformType::Pointer affineTransform=AffineTransformType::New(); affineTransform->SetMatrix(rotationMatrix); affineTransform->SetTranslation(translationPart); // Interp typedef clitk::GenericVectorInterpolator GenericInterpolatorType; typename GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New(); genericInterpolator->SetArgsInfo(m_ArgsInfo); // Properties if (m_ArgsInfo.like_given) { typename InputReaderType::Pointer likeReader=InputReaderType::New(); likeReader->SetFileName(m_ArgsInfo.like_arg); likeReader->Update(); resampler->SetSize( likeReader->GetOutput()->GetLargestPossibleRegion().GetSize() ); resampler->SetOutputSpacing( likeReader->GetOutput()->GetSpacing() ); resampler->SetOutputOrigin( likeReader->GetOutput()->GetOrigin() ); } else { //Size typename OutputImageType::SizeType outputSize; if (m_ArgsInfo.size_given) { for(unsigned int i=0; i< Dimension; i++) outputSize[i]=m_ArgsInfo.size_arg[i]; } else outputSize=input->GetLargestPossibleRegion().GetSize(); std::cout<<"Setting the size to "<GetSpacing(); std::cout<<"Setting the spacing to "<GetOrigin(); std::cout<<"Setting the origin to "<SetSize( outputSize ); resampler->SetOutputSpacing( outputSpacing ); resampler->SetOutputOrigin( outputOrigin ); } resampler->SetInput( input ); resampler->SetTransform( affineTransform ); resampler->SetInterpolator( genericInterpolator->GetInterpolatorPointer()); resampler->SetDefaultPixelValue( static_cast(m_ArgsInfo.pad_arg) ); try { resampler->Update(); } catch(itk::ExceptionObject) { std::cerr<<"Error resampling the image"<GetOutput(); // Output typedef itk::ImageFileWriter WriterType; typename WriterType::Pointer writer = WriterType::New(); writer->SetFileName(m_ArgsInfo.output_arg); writer->SetInput(output); writer->Update(); } //------------------------------------------------------------------- //------------------------------------------------------------------- template template typename itk::Matrix AffineTransformGenericFilter::createMatrixFromElastixFile(std::vector & filename) { if (Dimension != 3) { FATAL("Only 3D yet" << std::endl); } typename itk::Matrix matrix; itk::CenteredEuler3DTransform::Pointer mat = itk::CenteredEuler3DTransform::New(); itk::CenteredEuler3DTransform::Pointer previous; for(uint j=0; j cor; GetValuesFromValue(s, cor); // Get Transformparameters GetElastixValueFromTag(is, "TransformParameters ", s); // space is needed if (!b) { FATAL("Error must read 'TransformParameters' in " << filename[j] << std::endl); } std::vector results; GetValuesFromValue(s, results); // construct a stream from the string itk::CenteredEuler3DTransform::ParametersType p; p.SetSize(9); for(uint i=0; i<3; i++) p[i] = atof(results[i].c_str()); // Rotation for(uint i=0; i<3; i++) p[i+3] = atof(cor[i].c_str()); // Centre of rotation for(uint i=0; i<3; i++) p[i+6] = atof(results[i+3].c_str()); // Translation mat->SetParameters(p); if (m_Verbose) { std::cout << "Rotation (deg) : " << rad2deg(p[0]) << " " << rad2deg(p[1]) << " " << rad2deg(p[2]) << std::endl; std::cout << "Center of rot (phy) : " << p[3] << " " << p[4] << " " << p[5] << std::endl; std::cout << "Translation (phy) : " << p[6] << " " << p[7] << " " << p[8] << std::endl; } // Compose with previous if needed if (j!=0) { mat->Compose(previous); if (m_Verbose) { std::cout << "Composed rotation (deg) : " << rad2deg(mat->GetAngleX()) << " " << rad2deg(mat->GetAngleY()) << " " << rad2deg(mat->GetAngleZ()) << std::endl; std::cout << "Composed center of rot (phy) : " << mat->GetCenter() << std::endl; std::cout << "Compsoed translation (phy) : " << mat->GetTranslation() << std::endl; } } previous = mat->Clone(); } mat = previous; for(uint i=0; i<3; i++) for(uint j=0; j<3; j++) matrix[i][j] = mat->GetMatrix()[i][j]; // Offset is -Rc + t + c matrix[0][3] = mat->GetOffset()[0]; matrix[1][3] = mat->GetOffset()[1]; matrix[2][3] = mat->GetOffset()[2]; matrix[3][3] = 1; return matrix; } //------------------------------------------------------------------- template bool AffineTransformGenericFilter::GetElastixValueFromTag(std::ifstream & is, std::string tag, std::string & value) { std::string line; is.seekg (0, is.beg); while(std::getline(is, line)) { unsigned pos = line.find(tag); if (pos void AffineTransformGenericFilter::GetValuesFromValue(const std::string & s, std::vector & values) { std::stringstream strstr(s); std::istream_iterator it(strstr); std::istream_iterator end; std::vector results(it, end); values.clear(); values.resize(results.size()); for(uint i=0; i