From 80552a943483c27a0aa5dc295965d621645ace9d Mon Sep 17 00:00:00 2001 From: Gang Song Date: Thu, 30 Jun 2011 23:06:14 +0000 Subject: [PATCH] ENH: Average affine transform: AverageAffineTransform, support both 2D and 3D, and it's weighted. --- CMakeLists.txt | 1 + Examples/AverageAffineTransform.cxx | 420 ++++++++++++++++++ Examples/CMakeLists.txt | 34 +- .../itkANTSAffine3DTransform.txx | 243 ++++++---- Utilities/itkAverageAffineTransformFunction.h | 203 +++++++++ .../itkAverageAffineTransformFunction.txx | 333 ++++++++++++++ 6 files changed, 1140 insertions(+), 94 deletions(-) create mode 100644 Examples/AverageAffineTransform.cxx create mode 100644 Utilities/itkAverageAffineTransformFunction.h create mode 100644 Utilities/itkAverageAffineTransformFunction.txx diff --git a/CMakeLists.txt b/CMakeLists.txt index ae1b30e6b..739c99212 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,4 +1,5 @@ project(ANTS) cmake_minimum_required(VERSION 2.8.2) set(CMAKE_BUILD_TYPE "Release") + add_subdirectory(Examples) diff --git a/Examples/AverageAffineTransform.cxx b/Examples/AverageAffineTransform.cxx new file mode 100644 index 000000000..4fb81dcaf --- /dev/null +++ b/Examples/AverageAffineTransform.cxx @@ -0,0 +1,420 @@ +// compute the average of a list of affine transform + +#include +#include +#include "itkImageFileReader.h" +#include "itkVector.h" +// #include "itkVectorImageFileReader.h" +// #include "itkVectorImageFileWriter.h" +#include "itkImageFileWriter.h" +#include "itkMatrixOffsetTransformBase.h" +#include "itkTransformFactory.h" +// #include "itkWarpImageMultiTransformFilter.h" +// #include "itkDeformationFieldFromMultiTransformFilter.h" + +#include "itkAverageAffineTransformFunction.h" + +#include "itkTransformFileReader.h" +#include "itkTransformFileWriter.h" +#include +#include + +typedef enum + { + INVALID_FILE = 1, AFFINE_FILE, DEFORMATION_FILE + } TRAN_FILE_TYPE; +typedef struct + { + char * filename; + TRAN_FILE_TYPE file_type; + bool do_affine_inv; + + double weight; // for average + } TRAN_OPT; + +typedef std::vector TRAN_OPT_QUEUE; + +TRAN_FILE_TYPE CheckFileType(char *str) +{ + std::string filename = str; + std::string::size_type pos = filename.rfind("."); + std::string filepre = std::string(filename, 0, pos); + + if( pos != std::string::npos ) + { + std::string extension = std::string(filename, pos, + filename.length() - 1); + if( extension == std::string(".gz") ) + { + pos = filepre.rfind("."); + extension = std::string(filepre, pos, filepre.length() - 1); + } + if( extension == ".txt" ) + { + return AFFINE_FILE; + } + else + { + return DEFORMATION_FILE; + } + } + else + { + return INVALID_FILE; + } + return AFFINE_FILE; +} + +// adapted from http://stackoverflow.com/questions/194465/how-to-parse-a-string-to-an-int-in-c +bool get_a_double_number(const char *str, double & v) +{ + char *end; + + errno = 0; + v = strtod(str, &end); + if( (errno == ERANGE && v == HUGE_VAL) ) + { + return false; // OVERFLOW + } + if( (errno == ERANGE && v == 0.0) ) + { + return false; // UNDERFLOW + } + if( *str == '\0' || *end != '\0' ) + { + return false; // INCONVERTIBLE; + } + return true; +} + +bool ParseInput(int argc, char * *argv, char *& output_transform_filename, + char *& reference_transform_filename, TRAN_OPT_QUEUE & opt_queue) +{ + opt_queue.clear(); + opt_queue.reserve(argc); + + output_transform_filename = argv[0]; + + reference_transform_filename = NULL; + + int ind = 1; + while( ind < argc ) + { + if( strcmp(argv[ind], "-R") == 0 ) + { + ind++; + if( ind >= argc ) + { + return false; + } + reference_transform_filename = argv[ind]; + } + else if( strcmp(argv[ind], "-i") == 0 ) + { + ind++; + if( ind >= argc ) + { + return false; + } + TRAN_OPT opt; + opt.filename = argv[ind]; + if( CheckFileType(opt.filename) != AFFINE_FILE ) + { + std::cout << "file: " << opt.filename + << " is not an affine .txt file. Invalid to use '-i' " + << std::endl; + return false; + } + opt.file_type = AFFINE_FILE; + opt.do_affine_inv = true; + + opt.weight = 1.0; // default value + if( ind < argc - 1 ) // test if still has extra parameters + { + double weight; + if( get_a_double_number(argv[ind + 1], weight) ) + { + ind++; + opt.weight = weight; + } + } + opt_queue.push_back(opt); + } + else + { + TRAN_OPT opt; + opt.filename = argv[ind]; + if( CheckFileType(opt.filename) != AFFINE_FILE ) + { + std::cout << "file: " << opt.filename + << " is not an affine .txt file." << std::endl; + return false; + } + opt.file_type = CheckFileType(opt.filename); + opt.do_affine_inv = false; + + opt.weight = 1.0; // default value + if( ind < argc - 1 ) // test if still has extra parameters + { + double weight; + if( get_a_double_number(argv[ind + 1], weight) ) + { + ind++; + opt.weight = weight; + } + } + + opt_queue.push_back(opt); + } + ind++; + } + +// if (reference_image_filename == NULL) { +// std::cout << "the reference image file (-R) must be given!!!" +// << std::endl; +// return false; +// } + + return true; +} + +void DisplayOptQueue(const TRAN_OPT_QUEUE & opt_queue) +{ + const int kQueueSize = opt_queue.size(); + + for( int i = 0; i < kQueueSize; i++ ) + { + std::cout << "[" << i << "/" << kQueueSize << "]: "; + + switch( opt_queue[i].file_type ) + { + case AFFINE_FILE: + { + std::cout << "AFFINE"; + if( opt_queue[i].do_affine_inv ) + { + std::cout << "-INV"; + } + std::cout << " (" << opt_queue[i].weight << ") "; + } + break; +// case DEFORMATION_FILE: +// std::cout << "FIELD"; +// break; + default: + { + std::cout << "Invalid Format!!!"; + } + break; + } + std::cout << ": " << opt_queue[i].filename << std::endl; + } +} + +template +void AverageAffineTransform(char *output_affine_txt, char *reference_affine_txt, + TRAN_OPT_QUEUE & opt_queue) +{ +// typedef itk::Image ImageType; +// typedef itk::Vector VectorType; +// typedef itk::Image DeformationFieldType; + + typedef itk::MatrixOffsetTransformBase AffineTransformType; +// typedef itk::WarpImageMultiTransformFilter WarperType; + + typedef itk::AverageAffineTransformFunction WarperType; + + itk::TransformFactory::RegisterTransform(); + + // typedef itk::ImageFileReader ImageFileReaderType; + // typename ImageFileReaderType::Pointer reader_img = ImageFileReaderType::New(); + // typename ImageType::Pointer img_ref = ImageType::New(); + + // typename ImageFileReaderType::Pointer reader_img_ref = ImageFileReaderType::New(); + + WarperType average_func; + // warper->SetInput(img_mov); + // warper->SetEdgePaddingValue( 0); +// VectorType pad; +// pad.Fill(0); + // warper->SetEdgePaddingValue(pad); + + typedef itk::TransformFileReader TranReaderType; + +// typedef itk::ImageFileReader FieldReaderType; + + int cnt_affine = 0; + const int kOptQueueSize = opt_queue.size(); + for( int i = 0; i < kOptQueueSize; i++ ) + { + const TRAN_OPT & opt = opt_queue[i]; + + switch( opt.file_type ) + { + case AFFINE_FILE: + { + typename TranReaderType::Pointer tran_reader = + TranReaderType::New(); + tran_reader->SetFileName(opt.filename); + tran_reader->Update(); + typename AffineTransformType::Pointer aff = + dynamic_cast( (tran_reader->GetTransformList() )->front().GetPointer() ); + + if( opt_queue[i].do_affine_inv ) + { + aff->GetInverse(aff); + } + // std::cout << aff << std::endl; + + double weight = opt.weight; + average_func.PushBackAffineTransform(aff, weight); + cnt_affine++; + break; + } + case DEFORMATION_FILE: + { + std::cout << "Average affine only files: ignore " << opt.filename + << std::endl; + break; + } + default: + std::cout << "Unknown file type!" << std::endl; + } + } + + typedef typename WarperType::PointType PointType; + PointType aff_center; + + typename AffineTransformType::Pointer aff_ref_tmp; + if( reference_affine_txt ) + { + typename TranReaderType::Pointer tran_reader = TranReaderType::New(); + tran_reader->SetFileName(reference_affine_txt); + tran_reader->Update(); + aff_ref_tmp = + dynamic_cast( (tran_reader->GetTransformList() )->front().GetPointer() ); + } + else + { + if( cnt_affine > 0 ) + { + std::cout << "the reference affine file for center is selected as the first affine!" << std::endl; + aff_ref_tmp = average_func.GetTransformList().begin()->aff; + } + else + { + std::cout << "No affine input is given. nothing to do ......" << std::endl; + return; + } + } + + aff_center = aff_ref_tmp->GetCenter(); + std::cout << "new center is : " << aff_center << std::endl; + + // warper->PrintTransformList(); + + // typename AffineTransformType::Pointer aff_output = warper->ComposeAffineOnlySequence(aff_center); + typename AffineTransformType::Pointer aff_output = AffineTransformType::New(); + + average_func.AverageMultipleAffineTransform(aff_center, aff_output); + + typedef itk::TransformFileWriter TranWriterType; + typename TranWriterType::Pointer tran_writer = TranWriterType::New(); + tran_writer->SetFileName(output_affine_txt); + tran_writer->SetInput(aff_output); + tran_writer->Update(); + + std::cout << "wrote file to : " << output_affine_txt << std::endl; +} + +int main(int argc, char * *argv) +{ + if( argc <= 3 ) + { + std::cout + << "AverageAffineTransform ImageDimension output_affine_transform [-R reference_affine_transform] " + << "{[-i] affine_transform_txt [weight(=1)] ]}" + << std::endl + << std::endl + << " Usage: Compute weighted average of input affine transforms. " + << std::endl + << "For 2D and 3D transform, the affine transform is first decomposed into " + "scale x shearing x rotation. Then these parameters are averaged, using the weights if they provided. " + "For 3D transform, the rotation component is the quaternion. After averaging, the quaternion will also " + "be normalized to have unit norm. For 2D transform, the rotation component is the rotation angle. " + "The weight for each transform is a non-negative number. The sum of all weights will be normalized to 1 " + "before averaging. The default value for each weight is 1.0. " + << std::endl + << std::endl + << "All affine transforms is a \"centerd\" transform, following ITK convention. A reference_affine_transform" + " defines the center for the output transform. The first provided transform is the default reference " + "transform" + << std::endl + << "Output affine transform is a MatrixOffsetBaseTransform." + << std::endl + << " -i option takes the inverse of the affine mapping." + << std::endl + << " For example: " + << std::endl + << " 2 output_affine.txt -R A.txt A1.txt 1.0 -i A2.txt 2.0 A3.txt A4.txt 6.0 A5.txt" + << std::endl + << "This computes: (1*A1 + 2*(A2)^-1 + A3 + A4*6 + A5 ) / (1+2+1+6+5)" + << std::endl; + return EXIT_SUCCESS; + } + + TRAN_OPT_QUEUE opt_queue; + + char *output_transform_filename = NULL; + char *reference_transform_filename = NULL; + + bool is_parsing_ok = false; + int kImageDim = atoi(argv[1]); + + is_parsing_ok = ParseInput(argc - 2, argv + 2, output_transform_filename, + reference_transform_filename, opt_queue); + + if( is_parsing_ok ) + { + std::cout << "output_transform_filename: " << output_transform_filename + << std::endl; + std::cout << "reference_transform_filename: "; + + if( reference_transform_filename ) + { + std::cout << reference_transform_filename << std::endl; + } + else + { + std::cout << "NULL" << std::endl; + } + + DisplayOptQueue(opt_queue); + + switch( kImageDim ) + { + case 2: + { + AverageAffineTransform<2>(output_transform_filename, + reference_transform_filename, opt_queue); + break; + } + case 3: + { + AverageAffineTransform<3>(output_transform_filename, + reference_transform_filename, opt_queue); + break; + } + } + } + + else + { + std::cout << "Input error!" << std::endl; + return EXIT_FAILURE; + } + + return EXIT_SUCCESS; +} diff --git a/Examples/CMakeLists.txt b/Examples/CMakeLists.txt index e9d1bdf93..63e3d8ef7 100644 --- a/Examples/CMakeLists.txt +++ b/Examples/CMakeLists.txt @@ -1,7 +1,7 @@ project(ANTS) cmake_minimum_required(VERSION 2.8.2) -set(CMAKE_BUILD_TYPE "Release") +# set(CMAKE_BUILD_TYPE "Release") #Change PROJECT_NAME to the name of your project @@ -10,7 +10,6 @@ option(USE_VTK "Use VTK Libraries" OFF) set (CMAKE_INCLUDE_DIRECTORIES_BEFORE ON) - # Set up ITK if(USE_ITK) find_package(ITK) @@ -42,11 +41,36 @@ if(USE_VTK) endif(VTK_FOUND) endif(USE_VTK) +option(USE_FFTWD "Use double precision fftw if found" OFF) +option(USE_FFTWF "Use single precision fftw if found" OFF) +option(USE_SYSTEM_FFTW "Use an installed version of fftw" OFF) +if (USE_FFTWD OR USE_FFTWF) + if(USE_SYSTEM_FFTW) + find_package( FFTW ) + link_directories(${FFTW_LIBDIR}) + else(USE_SYSTEM_FFTW) + link_directories(${ITK_DIR}/fftw/lib) + include_directories(${ITK_DIR}/fftw/include) + endif(USE_SYSTEM_FFTW) +endif(USE_FFTWD OR USE_FFTWF) + + + #The following lines are required to use Dart include(CTest) enable_testing() + + + + + + + + + + #INCLUDE (${CMAKE_ROOT}/Modules/FindITK.cmake) #IF (USE_ITK_FILE) # include(${USE_ITK_FILE}) @@ -262,6 +286,10 @@ target_link_libraries(KellySlater ${ITK_LIBRARIES}) add_executable(KellyKapowski KellyKapowski.cxx ${UI_SOURCES}) target_link_libraries(KellyKapowski ${ITK_LIBRARIES}) +add_executable(AverageAffineTransform AverageAffineTransform.cxx) + target_link_libraries(AverageAffineTransform ${ITK_LIBRARIES}) + + IF (USE_VTK) INCLUDE (${CMAKE_ROOT}/Modules/FindVTK.cmake) IF (USE_VTK_FILE) @@ -295,6 +323,8 @@ target_link_libraries(WarpVTKPolyDataMultiTransform ${ITK_LIBRARIES} vtksys vtkI add_executable(ANTSConformalMapping ANTSConformalMapping.cxx ${UI_SOURCES}) target_link_libraries(ANTSConformalMapping ${ITK_LIBRARIES} FLATFEM ${ITK_LIBRARIES} ITKMetaIO vtkRendering vtkIO) + + ENDIF (USE_VTK_FILE) endif(USE_VTK) diff --git a/ImageRegistration/itkANTSAffine3DTransform.txx b/ImageRegistration/itkANTSAffine3DTransform.txx index 3486d00da..60daeb9c2 100644 --- a/ImageRegistration/itkANTSAffine3DTransform.txx +++ b/ImageRegistration/itkANTSAffine3DTransform.txx @@ -8,11 +8,10 @@ namespace itk { // Constructor with default arguments template -ANTSAffine3DTransform -::ANTSAffine3DTransform() : +ANTSAffine3DTransform::ANTSAffine3DTransform() : Superclass(SpaceDimension, ParametersDimension) { - m_Rotation = VnlQuaternionType(0, 0, 0, 1); // axis * vcl_sin(t/2), vcl_cos(t/2) + m_Rotation = VnlQuaternionType(0, 0, 0, 1); // axis * vcl_sin(t/2), vcl_cos(t/2) m_S1 = NumericTraits::One; m_S2 = NumericTraits::One; m_S3 = NumericTraits::One; @@ -23,11 +22,11 @@ ANTSAffine3DTransform // Constructor with default arguments template -ANTSAffine3DTransform::ANTSAffine3DTransform( unsigned int outputSpaceDimension, - unsigned int parametersDimension ) : +ANTSAffine3DTransform::ANTSAffine3DTransform( + unsigned int outputSpaceDimension, unsigned int parametersDimension) : Superclass(outputSpaceDimension, parametersDimension) { - m_Rotation = VnlQuaternionType(0, 0, 0, 1); // axis * vcl_sin(t/2), vcl_cos(t/2) + m_Rotation = VnlQuaternionType(0, 0, 0, 1); // axis * vcl_sin(t/2), vcl_cos(t/2) m_S1 = NumericTraits::One; m_S2 = NumericTraits::One; m_S3 = NumericTraits::One; @@ -48,22 +47,24 @@ ANTSAffine3DTransform::ANTSAffine3DTransform( unsigned int outputSp // Print self template -void -ANTSAffine3DTransform::PrintSelf(std::ostream & os, Indent indent ) const +void ANTSAffine3DTransform::PrintSelf(std::ostream & os, + Indent indent) const { Superclass::PrintSelf(os, indent); - os << indent << "Rotation: " << m_Rotation << std::endl; - os << indent << "S1, S2, S3: " << m_S1 << ", " << m_S2 << ", " << m_S3 << std::endl; - os << indent << "K1, K2, K3: " << m_K1 << ", " << m_K2 << ", " << m_K3 << std::endl; + os << indent << "Rotation: " << m_Rotation << std::endl; + os << indent << "S1, S2, S3: " << m_S1 << ", " << m_S2 << ", " << m_S3 + << std::endl; + os << indent << "K1, K2, K3: " << m_K1 << ", " << m_K2 << ", " << m_K3 + << std::endl; } // Set rotation template -void -ANTSAffine3DTransform::SetRotation(const VnlQuaternionType & rotation ) +void ANTSAffine3DTransform::SetRotation( + const VnlQuaternionType & rotation) { - m_Rotation = rotation; + m_Rotation = rotation; this->ComputeMatrix(); this->ComputeOffset(); @@ -73,8 +74,7 @@ ANTSAffine3DTransform::SetRotation(const VnlQuaternionType & rotati } template -void -ANTSAffine3DTransform::SetS1(const TScalarType S1) +void ANTSAffine3DTransform::SetS1(const TScalarType S1) { m_S1 = S1; @@ -86,8 +86,7 @@ ANTSAffine3DTransform::SetS1(const TScalarType S1) } template -void -ANTSAffine3DTransform::SetS2(const TScalarType S2) +void ANTSAffine3DTransform::SetS2(const TScalarType S2) { m_S2 = S2; @@ -99,8 +98,7 @@ ANTSAffine3DTransform::SetS2(const TScalarType S2) } template -void -ANTSAffine3DTransform::SetS3(const TScalarType S3) +void ANTSAffine3DTransform::SetS3(const TScalarType S3) { m_S3 = S3; @@ -112,8 +110,7 @@ ANTSAffine3DTransform::SetS3(const TScalarType S3) } template -void -ANTSAffine3DTransform::SetK1(const TScalarType K1) +void ANTSAffine3DTransform::SetK1(const TScalarType K1) { m_K1 = K1; @@ -125,8 +122,7 @@ ANTSAffine3DTransform::SetK1(const TScalarType K1) } template -void -ANTSAffine3DTransform::SetK2(const TScalarType K2) +void ANTSAffine3DTransform::SetK2(const TScalarType K2) { m_K2 = K2; @@ -138,8 +134,7 @@ ANTSAffine3DTransform::SetK2(const TScalarType K2) } template -void -ANTSAffine3DTransform::SetK3(const TScalarType K3) +void ANTSAffine3DTransform::SetK3(const TScalarType K3) { m_K3 = K3; @@ -152,8 +147,7 @@ ANTSAffine3DTransform::SetK3(const TScalarType K3) // Set the parameters in order to fit an Identity transform template -void -ANTSAffine3DTransform::SetIdentity( void ) +void ANTSAffine3DTransform::SetIdentity(void) { m_Rotation = VnlQuaternionType(0, 0, 0, 1); m_S1 = NumericTraits::One; @@ -167,9 +161,8 @@ ANTSAffine3DTransform::SetIdentity( void ) // Set Parameters template -void -ANTSAffine3DTransform -::SetParameters( const ParametersType & parameters ) +void ANTSAffine3DTransform::SetParameters( + const ParametersType & parameters) { OutputVectorType translation; @@ -196,7 +189,7 @@ ANTSAffine3DTransform translation[i] = parameters[par]; ++par; } - this->SetVarTranslation( translation ); + this->SetVarTranslation(translation); this->ComputeOffset(); @@ -207,12 +200,9 @@ ANTSAffine3DTransform // Set Parameters template -const -typename ANTSAffine3DTransform::ParametersType -& ANTSAffine3DTransform -::GetParameters() const - { - VnlQuaternionType quaternion = this->GetRotation(); +const typename ANTSAffine3DTransform::ParametersType +& ANTSAffine3DTransform::GetParameters() const { + VnlQuaternionType quaternion = this->GetRotation(); OutputVectorType translation = this->GetTranslation(); // Transfer the quaternion part @@ -242,9 +232,8 @@ typename ANTSAffine3DTransform::ParametersType // Get parameters template const typename ANTSAffine3DTransform::JacobianType -& ANTSAffine3DTransform:: -GetJacobian( const InputPointType &p ) const - { +& ANTSAffine3DTransform::GetJacobian( + const InputPointType &p) const { this->m_Jacobian.Fill(0.0); TScalarType c1 = this->GetCenter()[0]; @@ -270,24 +259,30 @@ GetJacobian( const InputPointType &p ) const TScalarType z3 = s3 * w3; // compute Jacobian with respect to quaternion parameters - this->m_Jacobian[0][0] = 2.0 * ( m_Rotation.x() * z1 + m_Rotation.y() * z2 - + m_Rotation.z() * z3 ); - this->m_Jacobian[0][1] = 2.0 * (-m_Rotation.y() * z1 + m_Rotation.x() * z2 - + m_Rotation.r() * z3 ); - this->m_Jacobian[0][2] = 2.0 * (-m_Rotation.z() * z1 - m_Rotation.r() * z2 - + m_Rotation.x() * z3 ); - this->m_Jacobian[0][3] = -2.0 * (-m_Rotation.r() * z1 + m_Rotation.z() * z2 - - m_Rotation.y() * z3 ); + this->m_Jacobian[0][0] = 2.0 + * (m_Rotation.x() * z1 + m_Rotation.y() * z2 + m_Rotation.z() * z3); + this->m_Jacobian[0][1] = + 2.0 + * (-m_Rotation.y() * z1 + m_Rotation.x() * z2 + + m_Rotation.r() * z3); + this->m_Jacobian[0][2] = + 2.0 + * (-m_Rotation.z() * z1 - m_Rotation.r() * z2 + + m_Rotation.x() * z3); + this->m_Jacobian[0][3] = + -2.0 + * (-m_Rotation.r() * z1 + m_Rotation.z() * z2 + - m_Rotation.y() * z3); this->m_Jacobian[1][0] = -this->m_Jacobian[0][1]; - this->m_Jacobian[1][1] = this->m_Jacobian[0][0]; - this->m_Jacobian[1][2] = this->m_Jacobian[0][3]; + this->m_Jacobian[1][1] = this->m_Jacobian[0][0]; + this->m_Jacobian[1][2] = this->m_Jacobian[0][3]; this->m_Jacobian[1][3] = -this->m_Jacobian[0][2]; this->m_Jacobian[2][0] = -this->m_Jacobian[0][2]; this->m_Jacobian[2][1] = -this->m_Jacobian[0][3]; - this->m_Jacobian[2][2] = this->m_Jacobian[0][0]; - this->m_Jacobian[2][3] = this->m_Jacobian[0][1]; + this->m_Jacobian[2][2] = this->m_Jacobian[0][0]; + this->m_Jacobian[2][3] = this->m_Jacobian[0][1]; // get rotation matrix first // this is done to compensate for the transposed representation @@ -403,7 +398,8 @@ GetJacobian( const InputPointType &p ) const // } template -typename ANTSAffine3DTransform::MatrixType ANTSAffine3DTransform::ComputeMyRotationMatrix() +typename ANTSAffine3DTransform::MatrixType ANTSAffine3DTransform< + TScalarType>::ComputeMyRotationMatrix() { VnlQuaternionType conjugateRotation = m_Rotation.conjugate(); // this is done to compensate for the transposed representation @@ -415,8 +411,7 @@ typename ANTSAffine3DTransform::MatrixType ANTSAffine3DTransform -void -ANTSAffine3DTransform::ComputeMatrix() +void ANTSAffine3DTransform::ComputeMatrix() { VnlQuaternionType conjugateRotation = m_Rotation.conjugate(); // this is done to compensate for the transposed representation @@ -448,8 +443,7 @@ ANTSAffine3DTransform::ComputeMatrix() } template -void -ANTSAffine3DTransform::ComputeMatrixParameters() +void ANTSAffine3DTransform::ComputeMatrixParameters() { // VnlQuaternionType quat(this->GetMatrix().GetVnlMatrix()); // m_Rotation = quat; @@ -460,67 +454,132 @@ ANTSAffine3DTransform::ComputeMatrixParameters() typedef vnl_matrix TMatrix; - TMatrix A, Q, R; + TMatrix A, R, Q; A = this->GetMatrix().GetVnlMatrix(); vnl_qr myqr(A); - R = myqr.Q(); // Q() is the rotation - Q = myqr.R(); // R() is the upper triangluar + Q = myqr.Q(); // Q() is the rotation + R = myqr.R(); // R() is the upper triangluar - // this is not necessary, for the mirror case - // the scale factor could not negative - // normalize Q + // songgang: anyone added this??? + // this is not necessary, for the mirror case + // the scale factor could not negative + // normalize R + // need to run this, otherwise, identity matrix have negative scale!!! + +// // force diagnoal of rotation matrix to be positive // TMatrix dq(3,3,0); // for(unsigned i=0;i<3;i++){ -// dq(i,i) = (Q(i,i)>=0)? 1 : -1; +// dq(i,i) = (R(i,i)>=0)? 1 : -1; // } -// R = R * dq; -// Q = dq * Q; +// Q = Q * dq; +// R = dq * R; + + // force trace of rotation matrix be maximum possible by multiplying + // a diagonal (+/-1 +/-1 +/-1) whose determinant is always positive +1. + + double trA = Q(0, 0) + Q(1, 1) + Q(2, 2); // 1, 1, 1 + double trB = Q(0, 0) - Q(1, 1) - Q(2, 2); // 1, -1, -1 + double trC = -Q(0, 0) + Q(1, 1) - Q(2, 2); // -1, 1, -1 + double trD = -Q(0, 0) - Q(1, 1) + Q(2, 2); // -1, -1, 1 + double maxTr = trA; // find the maximum of all terms; + if( trB > maxTr ) + { + maxTr = trB; // dividing by the maximum makes + } + if( trC > maxTr ) + { + maxTr = trC; // the computations more stable + } + if( trD > maxTr ) + { + maxTr = trD; // and avoid division by zero + } + if( maxTr == trB ) + { + TMatrix dq(3, 3, 0); + dq(0, 0) = 1; + dq(1, 1) = -1; + dq(2, 2) = -1; + Q = Q * dq; + R = dq * R; + } + + if( maxTr == trC ) + { + TMatrix dq(3, 3, 0); + dq(0, 0) = -1; + dq(1, 1) = 1; + dq(2, 2) = -1; + Q = Q * dq; + R = dq * R; + } + + if( maxTr == trD ) + { + TMatrix dq(3, 3, 0); + dq(0, 0) = -1; + dq(1, 1) = -1; + dq(2, 2) = 1; + Q = Q * dq; + R = dq * R; + } - double tr = 1 + R(0, 0) + R(1, 1) + R(2, 2); + double tr = 1 + Q(0, 0) + Q(1, 1) + Q(2, 2); double s, r, u, v, w; if( tr > 0 ) { s = 0.5 / sqrt(tr); r = 0.25 / s; - u = (R(2, 1) - R(1, 2) ) * s; - v = (R(0, 2) - R(2, 0) ) * s; - w = (R(1, 0) - R(0, 1) ) * s; + u = (Q(2, 1) - Q(1, 2) ) * s; + v = (Q(0, 2) - Q(2, 0) ) * s; + w = (Q(1, 0) - Q(0, 1) ) * s; } - else if( R(0, 0) > R(1, 1) && R(0, 0) > R(2, 2) ) + else if( Q(0, 0) > Q(1, 1) && Q(0, 0) > Q(2, 2) ) { - s = 2 * sqrt(1 + R(0, 0) - R(1, 1) - R(2, 2) ); + s = 2 * sqrt(1 + Q(0, 0) - Q(1, 1) - Q(2, 2) ); u = 0.25 * s; - v = (R(0, 1) + R(1, 0) ) / s; - w = (R(0, 2) + R(2, 0) ) / s; - r = (R(1, 2) - R(2, 1) ) / s; + v = (Q(0, 1) + Q(1, 0) ) / s; + w = (Q(0, 2) + Q(2, 0) ) / s; + r = (Q(1, 2) - Q(2, 1) ) / s; } - else if( R(0, 0) > R(1, 1) ) + else if( Q(0, 0) > Q(1, 1) ) { - s = 2 * sqrt(1 + R(1, 1) - R(0, 0) - R(2, 2) ); - u = (R(0, 1) + R(1, 0) ) / s; + s = 2 * sqrt(1 + Q(1, 1) - Q(0, 0) - Q(2, 2) ); + u = (Q(0, 1) + Q(1, 0) ) / s; v = 0.25 * s; - w = (R(1, 2) + R(2, 1) ) / s; - r = (R(0, 2) - R(2, 0) ) / s; + w = (Q(1, 2) + Q(2, 1) ) / s; + r = (Q(0, 2) - Q(2, 0) ) / s; } else { - s = 2 * sqrt(1 + R(2, 2) - R(0, 0) - R(1, 1) ); - u = (R(0, 2) + R(2, 0) ) / s; - v = (R(1, 2) + R(2, 1) ) / s; + s = 2 * sqrt(1 + Q(2, 2) - Q(0, 0) - Q(1, 1) ); + u = (Q(0, 2) + Q(2, 0) ) / s; + v = (Q(1, 2) + Q(2, 1) ) / s; w = 0.25 * s; - r = (R(0, 1) - R(1, 0) ) / s; + r = (Q(0, 1) - Q(1, 0) ) / s; } + std::cout << "A=" << A << std::endl; + std::cout << "rotation R" << Q << std::endl; + std::cout << "upper R" << R << std::endl; + std::cout << "s=" << s << " u=" << u << " v=" << v << " w" << w << " r=" + << r << std::endl; + m_Rotation = VnlQuaternionType(u, v, w, r); - m_S1 = Q(0, 0); - m_S2 = Q(1, 1); - m_S3 = Q(2, 2); - m_K1 = Q(0, 1) / Q(0, 0); - m_K2 = Q(0, 2) / Q(0, 0); - m_K3 = Q(1, 2) / Q(1, 1); + + std::cout << "m_Rotation from vnl" << VnlQuaternionType(u, v, w, r) + << std::endl; + + m_S1 = R(0, 0); + m_S2 = R(1, 1); + m_S3 = R(2, 2); + + m_K1 = R(0, 1) / R(0, 0); + m_K2 = R(0, 2) / R(0, 0); + m_K3 = R(1, 2) / R(1, 1); // std::cout << "before: this->GetMatrix(): " << this->GetMatrix(); @@ -529,7 +588,7 @@ ANTSAffine3DTransform::ComputeMatrixParameters() // std::cout << "after: this->GetMatrix(): " << this->GetMatrix(); // std::cout << "A=" << A << std::endl; -// std::cout << "Q=" << Q << std::endl; +// std::cout << "R=" << R << std::endl; // std::cout << "R=" << R << std::endl; // std::cout << "dq=" << dq << std::endl; } diff --git a/Utilities/itkAverageAffineTransformFunction.h b/Utilities/itkAverageAffineTransformFunction.h new file mode 100644 index 000000000..f10d5a678 --- /dev/null +++ b/Utilities/itkAverageAffineTransformFunction.h @@ -0,0 +1,203 @@ +/*========================================================================= + + Program: Advanced Normalization Tools + Module: $RCSfile: itkWarpImageMultiTransformFilter.h,v $ + Language: C++ + Date: $Date: 2009/01/08 21:36:48 $ + Version: $Revision: 1.18 $ + + Copyright (c) ConsortiumOfANTS. All rights reserved. + See accompanying COPYING.txt or + http://sourceforge.net/projects/advants/files/ANTS/ANTSCopyright.txt for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + + =========================================================================*/ +#ifndef __itkAverageAffineTransformFunction_h +#define __itkAverageAffineTransformFunction_h + +#include "itkMacro.h" +#include "itkConceptChecking.h" + +#include "itkANTSCenteredAffine2DTransform.h" +#include "itkANTSAffine3DTransform.h" + +#include + +namespace itk +{ +/* + */ + +namespace AverageAffineTransformFunctionHelperNameSpace +{ +template +struct Dispatcher + { + }; + +template +struct HelperCommonType + { + typedef TAffine InternalAffineTransformType; + + typedef typename InternalAffineTransformType::Pointer InternalAffineTransformPointerType; + typedef struct + { + InternalAffineTransformPointerType aff; + double weight; + } SingleInternalTransformItemType; + + typedef std::list InternalTransformListType; + typedef typename InternalAffineTransformType::ParametersType ParametersType; + + static void ComputeAveragePartialParameters(InternalTransformListType & transform_list, + ParametersType & average_parameters, unsigned int iStart, + unsigned int iEnd); + }; + +template +class HelperType; + +// { +//// purposely not include any types +// }; + +// explicit specialization for 2D affine transform +template <> +struct HelperType > + { +// public: + typedef ANTSCenteredAffine2DTransform InternalAffineTransformType; + + typedef HelperCommonType::InternalAffineTransformPointerType + InternalAffineTransformPointerType; + typedef HelperCommonType::SingleInternalTransformItemType + SingleInternalTransformItemType; + typedef HelperCommonType::InternalTransformListType InternalTransformListType; + typedef HelperCommonType::ParametersType ParametersType; + + static void ComputeAverageScaleParameters(InternalTransformListType & transform_list, + ParametersType & average_parameters); + + static void ComputeAverageShearingParameters(InternalTransformListType & transform_list, + ParametersType & average_parameters); + + static void ComputeAverageRotationParameters(InternalTransformListType & transform_list, + ParametersType & average_parameters); + + static void ComputeAverageTranslationParameters(InternalTransformListType & transform_list, + ParametersType & average_parameters); + }; + +// explicit specialization for 3D affine transform +template <> +struct HelperType > + { + typedef ANTSAffine3DTransform InternalAffineTransformType; + + typedef HelperCommonType::InternalAffineTransformPointerType + InternalAffineTransformPointerType; + typedef HelperCommonType::SingleInternalTransformItemType + SingleInternalTransformItemType; + typedef HelperCommonType::InternalTransformListType InternalTransformListType; + typedef HelperCommonType::ParametersType ParametersType; + + static void ComputeAverageScaleParameters(InternalTransformListType & transform_list, + ParametersType & average_parameters); + + static void ComputeAverageShearingParameters(InternalTransformListType & transform_list, + ParametersType & average_parameters); + + static void ComputeAverageRotationParameters(InternalTransformListType & transform_list, + ParametersType & average_parameters); + + static void ComputeAverageTranslationParameters(InternalTransformListType & transform_list, + ParametersType & average_parameters); + }; +} + +template +class ITK_EXPORT AverageAffineTransformFunction +{ +public: + + typedef TTransform GenericAffineTransformType; + + itkStaticConstMacro(InputSpaceDimension, unsigned int, GenericAffineTransformType::InputSpaceDimension); + itkStaticConstMacro(OutputSpaceDimension, unsigned int, GenericAffineTransformType::OutputSpaceDimension); + itkStaticConstMacro(SpaceDimension, unsigned int, InputSpaceDimension); + + typedef double InternalScalarType; + typedef Point PointType; + + AverageAffineTransformFunction(); + ~AverageAffineTransformFunction() + { + } + + ; + + // void PrintSelf(std::ostream& os, Indent indent) const; + + typedef typename GenericAffineTransformType::Pointer GenericAffineTransformPointerType; + + typedef struct + { + GenericAffineTransformPointerType aff; + double weight; + } SingleTransformItemType; + + typedef std::list TransformListType; + + TransformListType & GetTransformList() + { + return m_TransformList; + } + + void PrintTransformList(); + + void PushBackAffineTransform(const GenericAffineTransformType* t, double weight); + + void AverageMultipleAffineTransform(const PointType & center_output, + GenericAffineTransformPointerType & affine_output); + +#ifdef ITK_USE_CONCEPT_CHECKING + /** Begin concept checking */ + itkConceptMacro(SameDimensionCheck1, + (Concept::SameDimension ) ); + /** End concept checking */ +#endif +protected: + TransformListType m_TransformList; + +// type declaration to include support both 2D and 3D affine transform +protected: + + typedef typename::itk::AverageAffineTransformFunctionHelperNameSpace::HelperType< + ::itk::AverageAffineTransformFunctionHelperNameSpace::Dispatcher< + SpaceDimension> > HelperType; + + typedef typename HelperType::InternalAffineTransformType InternalAffineTransformType; + typedef typename HelperType::InternalAffineTransformPointerType InternalAffineTransformPointerType; + typedef typename HelperType::SingleInternalTransformItemType SingleInternalTransformItemType; + typedef typename HelperType::InternalTransformListType InternalTransformListType; + + InternalTransformListType m_InternalTransformList; + + void ConvertGenericAffineToInternalAffineByFixingCenter(GenericAffineTransformPointerType & aff, + InternalAffineTransformPointerType & iaff, + const PointType & center); + + void ConvertInternalAffineToGenericAffine(InternalAffineTransformPointerType & iaff, + GenericAffineTransformPointerType & aff); +}; +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkAverageAffineTransformFunction.txx" +#endif + +#endif /*__itkAverageAffineTransformFunction_h*/ diff --git a/Utilities/itkAverageAffineTransformFunction.txx b/Utilities/itkAverageAffineTransformFunction.txx new file mode 100644 index 000000000..6fe7eb96e --- /dev/null +++ b/Utilities/itkAverageAffineTransformFunction.txx @@ -0,0 +1,333 @@ +/*========================================================================= + + Program: Advanced Normalization Tools + Module: $RCSfile: itkWarpImageMultiTransformFilter.txx,v $ + Language: C++ + Date: $Date: 2009/01/08 21:36:48 $ + Version: $Revision: 1.18 $ + + Copyright (c) ConsortiumOfANTS. All rights reserved. + See accompanying COPYING.txt or + http://sourceforge.net/projects/advants/files/ANTS/ANTSCopyright.txt for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + + =========================================================================*/ +#ifndef __itkAverageAffineTransformFunction_txx +#define __itkAverageAffineTransformFunction_txx +#include "itkAverageAffineTransformFunction.h" + +#include "itkNumericTraits.h" +#include + +namespace itk +{ +/** + * Default constructor. + */ +template +AverageAffineTransformFunction::AverageAffineTransformFunction() +{ +} + +template +void AverageAffineTransformFunction::PrintTransformList() +{ + std::cout << "transform list: " << std::endl; + + typename TransformListType::iterator it = (m_TransformList.begin() ); + for( int ii = 0; it != m_TransformList.end(); it++, ii++ ) + { + std::cout << '[' << ii << ":" << it->weight << "]:" << it->aff + << std::endl; + } +} + +///** +// * Standard PrintSelf method. +// */ +// template +// void +// WarpImageMultiTransformFilter +// ::PrintSelf(std::ostream& os, Indent indent) const +// { +// +// Superclass::PrintSelf(os, indent); +// } + +template +void AverageAffineTransformFunction::PushBackAffineTransform( + const GenericAffineTransformType* t, double weight) +{ + if( t ) + { + SingleTransformItemType item; + item.aff = const_cast(t); + item.weight = weight; + m_TransformList.push_back(SingleTransformItemType(item) ); + } +} + +template +void AverageAffineTransformFunction::AverageMultipleAffineTransform( + const PointType & reference_center, + GenericAffineTransformPointerType & affine_output) +{ +// std::cout << "test " ; +// TransformTypePointer affine_output = TransformType::New(); + + affine_output->SetIdentity(); + affine_output->SetCenter(reference_center); + + unsigned int number_of_affine = m_TransformList.size(); + + number_of_affine--; + +// std::cout << affine_output; + + typename TransformListType::iterator it = m_TransformList.begin(); + +// typename InternalAffineTransformType::InputPointType center_ref = m_ReferenceCenter; + typename InternalAffineTransformType::Pointer average_iaff = + InternalAffineTransformType::New(); + + typename InternalAffineTransformType::ParametersType average_parameters = + average_iaff->GetParameters(); + for( ; it != m_TransformList.end(); it++ ) + { + SingleInternalTransformItemType internal_item; + internal_item.aff = InternalAffineTransformType::New(); + ConvertGenericAffineToInternalAffineByFixingCenter(it->aff, + internal_item.aff, reference_center); + internal_item.weight = it->weight; + m_InternalTransformList.push_back(internal_item); + + std::cout << "internal_transform: " << internal_item.aff << std::endl; + } + + HelperType::ComputeAverageScaleParameters(m_InternalTransformList, + average_parameters); + HelperType::ComputeAverageShearingParameters(m_InternalTransformList, + average_parameters); + HelperType::ComputeAverageRotationParameters(m_InternalTransformList, + average_parameters); + HelperType::ComputeAverageTranslationParameters(m_InternalTransformList, + average_parameters); + + average_iaff->SetParameters(average_parameters); + average_iaff->SetCenter(reference_center); + + std::cout << "average_iaff" << average_iaff << std::endl; + + ConvertInternalAffineToGenericAffine(average_iaff, affine_output); + + std::cout << "affine_output" << affine_output << std::endl; + return; +} + +template +void AverageAffineTransformFunction::ConvertGenericAffineToInternalAffineByFixingCenter( + GenericAffineTransformPointerType & aff, + InternalAffineTransformPointerType & iaff, const PointType & center) +{ + iaff->SetCenter(center); + iaff->SetMatrix(aff->GetMatrix() ); + iaff->SetTranslation(aff->GetTranslation() ); + + return; +} + +template +void AverageAffineTransformFunction::ConvertInternalAffineToGenericAffine( + InternalAffineTransformPointerType & iaff, + GenericAffineTransformPointerType & aff) +{ + aff->SetCenter(iaff->GetCenter() ); + aff->SetTranslation(iaff->GetTranslation() ); + aff->SetMatrix(iaff->GetMatrix() ); + + return; +} + +namespace AverageAffineTransformFunctionHelperNameSpace +{ +template +void HelperCommonType::ComputeAveragePartialParameters( + InternalTransformListType & transform_list, + ParametersType & average_parameters, unsigned int istart, + unsigned int iend) +{ + double w = 0.0; + + // initialize partial parameters to zero + for( unsigned int k = istart; k <= iend; k++ ) + { + average_parameters[k] = 0.0; + } + + typename InternalTransformListType::iterator it = transform_list.begin(); + unsigned int cnt = 0; + for( ; it != transform_list.end(); it++ ) + { + ParametersType current_parameters = it->aff->GetParameters(); + w += it->weight; + + std::cout << "[" << cnt++ << "]:" << it->weight << "\t"; + for( unsigned int k = istart; k <= iend; k++ ) + { + average_parameters[k] += it->weight * current_parameters[k]; + + std::cout << current_parameters[k] << " "; + } + + std::cout << std::endl; + } + + if( w <= 0.0 ) + { + std::cout << "Total weight smaller than 0!!!" << std::endl; + exit(-1); + } + + // normalize by weight + std::cout << "sum:w=" << w << "\t"; + for( unsigned int k = istart; k <= iend; k++ ) + { + std::cout << average_parameters[k] << " "; + } + std::cout << std::endl; + + // normalize by weight + std::cout << "average" << "\t"; + for( unsigned int k = istart; k <= iend; k++ ) + { + average_parameters[k] /= w; + std::cout << average_parameters[k] << " "; + } + + std::cout << std::endl; + return; +} + +void HelperType >::ComputeAverageScaleParameters( + InternalTransformListType & transform_list, + ParametersType & average_parameters) +{ + unsigned int istart = 1; + unsigned int iend = 2; + + std::cout << "average 2D scale parameter " << std::endl; + + HelperCommonType::ComputeAveragePartialParameters( + transform_list, average_parameters, istart, iend); +} + +void HelperType >::ComputeAverageShearingParameters( + InternalTransformListType & transform_list, + ParametersType & average_parameters) +{ + unsigned int istart = 3; + unsigned int iend = 3; + + std::cout << "average 2D shearing parameter " << std::endl; + + HelperCommonType::ComputeAveragePartialParameters( + transform_list, average_parameters, istart, iend); +} + +void HelperType >::ComputeAverageRotationParameters( + InternalTransformListType & transform_list, + ParametersType & average_parameters) +{ + unsigned int istart = 0; + unsigned int iend = 0; + + std::cout << "average 2D rotation parameter " << std::endl; + + HelperCommonType::ComputeAveragePartialParameters( + transform_list, average_parameters, istart, iend); +} + +void HelperType >::ComputeAverageTranslationParameters( + InternalTransformListType & transform_list, + ParametersType & average_parameters) +{ + unsigned int istart = 6; + unsigned int iend = 7; + + std::cout << "average 2D translation parameter " << std::endl; + + HelperCommonType::ComputeAveragePartialParameters( + transform_list, average_parameters, istart, iend); +} + +void HelperType >::ComputeAverageScaleParameters( + InternalTransformListType & transform_list, + ParametersType & average_parameters) +{ + unsigned int istart = 4; + unsigned int iend = 6; + + std::cout << "average 3D scale parameter " << std::endl; + + HelperCommonType::ComputeAveragePartialParameters( + transform_list, average_parameters, istart, iend); +} + +void HelperType >::ComputeAverageShearingParameters( + InternalTransformListType & transform_list, + ParametersType & average_parameters) +{ + unsigned int istart = 7; + unsigned int iend = 9; + + std::cout << "average 3D shearing parameter " << std::endl; + + HelperCommonType::ComputeAveragePartialParameters( + transform_list, average_parameters, istart, iend); +} + +void HelperType >::ComputeAverageRotationParameters( + InternalTransformListType & transform_list, + ParametersType & average_parameters) +{ + unsigned int istart = 0; + unsigned int iend = 3; + + std::cout << "average 3D rotation parameter " << std::endl; + + HelperCommonType::ComputeAveragePartialParameters( + transform_list, average_parameters, istart, iend); + + // extra normalization for quaternion + + double quat_mag = 0.0; + for( unsigned int j = istart; j <= iend; j++ ) + { + quat_mag += average_parameters[j] * average_parameters[j]; + } + quat_mag = sqrt(quat_mag); + for( unsigned int j = 0; j < 4; j++ ) + { + average_parameters[j] /= quat_mag; + } +} + +void HelperType >::ComputeAverageTranslationParameters( + InternalTransformListType & transform_list, + ParametersType & average_parameters) +{ + unsigned int istart = 10; + unsigned int iend = 12; + + std::cout << "average 3D translation parameter " << std::endl; + + HelperCommonType::ComputeAveragePartialParameters( + transform_list, average_parameters, istart, iend); +} +} // end namespace AverageAffineTransformFunctionHelperNameSpace +} // end namespace itk + +#endif // __itkAverageAffineTransformFunction_txx