diff --git a/CMake/elxModules.cmake b/CMake/elxModules.cmake index 894204aebbf5c5098a57bfd09d967c535fee639b..5cc11cd58dd02129e7b9591c53a2f7221917bc4e 100644 --- a/CMake/elxModules.cmake +++ b/CMake/elxModules.cmake @@ -35,7 +35,7 @@ macro( _elxmodule_enable MODULE_NAME ) endmacro() macro( _elxmodule_disable MODULE_NAME ) - # TODO + # TODO: elxmodule_disable endmacro() macro( _elxmodules_initialize ) diff --git a/CMake/elxITKRequiredModules.cmake b/CMake/elxRequiredITKModules.cmake similarity index 70% rename from CMake/elxITKRequiredModules.cmake rename to CMake/elxRequiredITKModules.cmake index c9cce4a03c670dab607b646114c1f44a241e44d8..61e56e78903db275dbfc4be42cad8266853a7a62 100644 --- a/CMake/elxITKRequiredModules.cmake +++ b/CMake/elxRequiredITKModules.cmake @@ -1,8 +1,8 @@ -list( APPEND ElastixRequiredITKModules +list( APPEND RequiredITKModules ITKReview ) -foreach( ITKModule ${ElastixRequiredITKModules} ) +foreach( ITKModule ${RequiredITKModules} ) if( NOT ${ITKModule}_LOADED ) message( FATAL_ERROR "Elastix requires that ITK is build with ${ITKModule}. Please rebuild ITK with Module_${ITKModule} set to ON." ) endif() diff --git a/CMake/elxWinConfig.cmake b/CMake/elxWinConfig.cmake deleted file mode 100644 index e0b1b38c129f9632c7c919475551c72a54c8bdf7..0000000000000000000000000000000000000000 --- a/CMake/elxWinConfig.cmake +++ /dev/null @@ -1,18 +0,0 @@ -# Visual Studio complains if paths are too long -string( LENGTH "${CMAKE_CURRENT_SOURCE_DIR}" n ) -if( n GREATER 50 ) -message( - FATAL_ERROR - "Source code directory path length is too long for MSVC (${n} > 50)." - "Please move the source code directory to a directory with a shorter path." - ) -endif() - -string( LENGTH "${CMAKE_CURRENT_BINARY_DIR}" n ) -if( n GREATER 50 ) -message( - FATAL_ERROR - "Build directory path length is too long for MSVC (${n} > 50)." - "Please move the build directory to a directory with a shorter path." - ) -endif() \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt index 9215420345da586c296658a2a617844fb37b2bd4..4f3daa923ad33969697cacfe13c932ea64c1fb53 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,4 +1,4 @@ -cmake_minimum_required( VERSION 2.8 ) +cmake_minimum_required( VERSION 3.0.2 ) # Explicitly add INCREMENTAL linking option to command lines. # http://www.cmake.org/pipermail/cmake/2010-February/035174.html @@ -17,9 +17,36 @@ set( CMAKE_RUNTIME_OUTPUT_DIRECTORY # Include SuperElastix CMake scripts list( APPEND CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/CMake" ) +# ----------------------------------------------------------------- +# Compiler-dependent settings +# GCC +if( ${CMAKE_CXX_COMPILER_ID} STREQUAL GNU ) + add_definitions( + -DVCL_CAN_STATIC_CONST_INIT_FLOAT=0 + -DVCL_CAN_STATIC_CONST_INIT_INT=0 + ) +endif() + +# MSVC if( ${CMAKE_CXX_COMPILER_ID} STREQUAL MSVC ) - include( elxWinConfig ) + string( LENGTH "${CMAKE_CURRENT_SOURCE_DIR}" n ) + if( n GREATER 50 ) + message( + FATAL_ERROR + "Source code directory path length is too long for MSVC (${n} > 50)." + "Please move the source code directory to a directory with a shorter path." + ) + endif() + + string( LENGTH "${CMAKE_CURRENT_BINARY_DIR}" n ) + if( n GREATER 50 ) + message( + FATAL_ERROR + "Build directory path length is too long for MSVC (${n} > 50)." + "Please move the build directory to a directory with a shorter path." + ) + endif() endif() # --------------------------------------------------------------------- @@ -30,11 +57,11 @@ endif() find_package( ITK REQUIRED ) include( ${ITK_USE_FILE} ) -include( "${CMAKE_CURRENT_SOURCE_DIR}/CMake/elxITKRequiredModules.cmake" ) +include("${CMAKE_CURRENT_SOURCE_DIR}/CMake/elxRequiredITKModules.cmake") # --------------------------------------------------------------------- # Boost Graph Library -find_package( Boost ) +find_package( Boost REQUIRED ) include_directories( ${Boost_INCLUDE_DIRS} ) # --------------------------------------------------------------------- @@ -49,24 +76,20 @@ endif() # --------------------------------------------------------------------- # Build SuperElastix + # For now we just enable all modules include( elxModules ) elxmodule_enable( ModuleCore ) elxmodule_enable( ModuleElastix ) +# TODO: Build tests depending on enabled modules + # --------------------------------------------------------------------- # Testing -# Testing requires CMake version 2.8.11 to download test data -if( CMAKE_VERSION VERSION_LESS 2.8.11 ) - set( SUPERELASTIX_BUILD_TESTING_DEFAULT OFF ) - message( STATUS "ELASTIX_BUILD_TESTING is set to OFF because CMake version is less than 2.8.11" ) -else() - set( SUPERELASTIX_BUILD_TESTING_DEFAULT ON ) -endif() - -option( SUPERELASTIX_BUILD_TESTING "Enable building tests." ${SUPERELASTIX_BUILD_TESTING_DEFAULT} ) +option( SUPERELASTIX_BUILD_TESTING "Enable building tests." ON ) if( ${SUPERELASTIX_BUILD_TESTING} ) + option( SUPERELASTIX_BUILD_LONG_TESTS OFF ) enable_testing() add_subdirectory( Testing ) endif() diff --git a/Modules/Components/Elastix/ModuleElastix.cmake b/Modules/Components/Elastix/ModuleElastix.cmake index 62e52cbe820ec5017440d40d59e32e9de35f8529..4d54550cce7922242e43706778cd745a47d29829 100644 --- a/Modules/Components/Elastix/ModuleElastix.cmake +++ b/Modules/Components/Elastix/ModuleElastix.cmake @@ -1,4 +1,4 @@ -# Module that exposes the old elastix as a component +# Module that exposes old elastix as an ITK filter set( MODULE ModuleElastix ) if( NOT EXISTS ${ELASTIX_USE_FILE} ) @@ -8,6 +8,14 @@ endif() include( ${ELASTIX_USE_FILE} ) +# If OpenMP is supported by this machine, elastix will be compiled with +# OpenMP flags, and we need to add them here as well +find_package( OpenMP ) +if (OPENMP_FOUND) + set( CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}" ) + set( CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}" ) +endif() + # Export include files set( ${MODULE}_INCLUDE_DIRS ${${MODULE}_SOURCE_DIR}/include diff --git a/Modules/Components/Elastix/include/elxElastixFilter.h b/Modules/Components/Elastix/include/elxElastixFilter.h index 3ae5102b1b0e56e06e9c8bd75c45b76cec9c0cf2..c333e0d1d167de899a20bdcac3241f5a14a240b4 100644 --- a/Modules/Components/Elastix/include/elxElastixFilter.h +++ b/Modules/Components/Elastix/include/elxElastixFilter.h @@ -1,8 +1,7 @@ -#ifndef ElastixComponent_h -#define ElastixComponent_h +#ifndef ElastixFilter_h +#define ElastixFilter_h // ITK -#include "itkProcessObject.h" #include "itkImageSource.h" // Elastix @@ -15,9 +14,8 @@ namespace selx { template< typename TFixedImage, - typename TMovingImage, - typename TOutputImage > -class ElastixFilter : public itk::ImageSource< TOutputImage > + typename TMovingImage > +class ElastixFilter : public itk::ImageSource< TFixedImage > { public: @@ -37,8 +35,9 @@ public: typedef ParameterObject::ParameterMapListType ParameterMapListType; typedef ParameterObject::ParameterMapType ParameterMapType; - typedef ParameterObject::ParameterValuesType ParameterValuesType; + typedef ParameterObject::ParameterVectorType ParameterVectorType; typedef typename ParameterObject::Pointer ParameterObjectPointer; + typedef typename ParameterObject::ConstPointer ParameterObjectConstPointer; typedef typename TFixedImage::Pointer FixedImagePointer; typedef typename TMovingImage::Pointer MovingImagePointer; @@ -47,18 +46,31 @@ public: void SetFixedImage( DataObjectContainerPointer fixedImage ); void SetMovingImage( MovingImagePointer movingImage ); void SetMovingImage( DataObjectContainerPointer movingImage ); - void SetParameterObject( ParameterObjectPointer parameterObject ); void SetFixedMask( FixedImagePointer fixedMask ); void SetMovingMask( MovingImagePointer movingMask ); + void SetParameterObject( ParameterObjectPointer parameterObject ); ParameterObjectPointer GetTransformParameters( void ); + itkSetMacro( FixedMeshFileName, std::string ); + itkGetConstMacro( FixedMeshFileName, std::string ); + void DeleteFixedMeshFileName( void ) { this->SetFixedMeshFileName( std::string() ); }; + + itkSetMacro( MovingMeshFileName, std::string ); + itkGetConstMacro( MovingMeshFileName, std::string ); + void DeleteMovingMeshFileName( void ) { this->SetMovingMeshFileName( std::string() ); }; + + itkSetMacro( OutputDirectory, std::string ); + itkGetConstMacro( OutputDirectory, std::string ); + void DeleteOutputDirectory() { this->m_OutputDirectory = std::string(); }; + itkSetMacro( LogToConsole, bool ); - itkGetMacro( LogToConsole, bool ); + itkGetConstMacro( LogToConsole, bool ); itkBooleanMacro( LogToConsole ); - itkSetMacro( LogToFile, std::string ); - itkGetMacro( LogToFile, std::string ); + itkSetMacro( LogToFile, bool ); + itkGetConstMacro( LogToFile, bool ); + itkBooleanMacro( LogToFile ); protected: @@ -68,13 +80,18 @@ private: ElastixFilter(); + // TODO: When set to true, ReleaseDataFlag should also touch these input containers DataObjectContainerPointer m_FixedImageContainer; DataObjectContainerPointer m_MovingImageContainer; DataObjectContainerPointer m_FixedMaskContainer; DataObjectContainerPointer m_MovingMaskContainer; + std::string m_FixedMeshFileName; + std::string m_MovingMeshFileName; + + std::string m_OutputDirectory; bool m_LogToConsole; - std::string m_LogToFile; + bool m_LogToFile; }; diff --git a/Modules/Components/Elastix/include/elxElastixFilter.hxx b/Modules/Components/Elastix/include/elxElastixFilter.hxx index 1cacc6e9bac66f1b06d9ac37de4c46aca5142a51..a213d50c6cf86453f347c8d7aa2892a8bc74989e 100644 --- a/Modules/Components/Elastix/include/elxElastixFilter.hxx +++ b/Modules/Components/Elastix/include/elxElastixFilter.hxx @@ -1,10 +1,10 @@ -#ifndef ElastixComponent_hxx -#define ElastixComponent_hxx +#ifndef ElastixFilter_hxx +#define ElastixFilter_hxx namespace selx { -template< typename TFixedImage, typename TMovingImage, typename TOutputImage > -ElastixFilter< TFixedImage, TMovingImage, TOutputImage > +template< typename TFixedImage, typename TMovingImage > +ElastixFilter< TFixedImage, TMovingImage > ::ElastixFilter( void ) { this->AddRequiredInputName( "FixedImage" ); @@ -16,60 +16,99 @@ ElastixFilter< TFixedImage, TMovingImage, TOutputImage > this->m_FixedImageContainer = DataObjectContainerType::New(); this->m_MovingImageContainer = DataObjectContainerType::New(); + + this->m_FixedMeshFileName = std::string(); + this->m_MovingMeshFileName = std::string(); + + this->m_OutputDirectory = std::string(); + this->LogToConsoleOff(); + this->LogToFileOff(); } -template< typename TFixedImage, typename TMovingImage, typename TOutputImage > +template< typename TFixedImage, typename TMovingImage > void -ElastixFilter< TFixedImage, TMovingImage, TOutputImage > +ElastixFilter< TFixedImage, TMovingImage > ::GenerateData( void ) { + // Initialize variables here so they don't go out of scope between iterations of the main loop ElastixMainObjectPointer transform = 0; DataObjectContainerPointer fixedImageContainer = this->m_FixedImageContainer; DataObjectContainerPointer movingImageContainer = this->m_MovingImageContainer; DataObjectContainerPointer fixedMaskContainer = 0; DataObjectContainerPointer movingMaskContainer = 0; DataObjectContainerPointer resultImageContainer = 0; - ParameterMapListType TransformParameterMapList; + FlatDirectionCosinesType fixedImageOriginalDirection; // Fixed mask (optional) if( this->HasInput( "FixedMask" ) ) { - fixedMaskContainer = DataObjectContainerType::New(); + fixedMaskContainer = DataObjectContainerType::New(); fixedMaskContainer->CreateElementAt( 0 ) = this->GetInput( "FixedMask" ); } // Moving mask (optional) if( this->HasInput( "MovingMask" ) ) { - movingMaskContainer = DataObjectContainerType::New(); + movingMaskContainer = DataObjectContainerType::New(); movingMaskContainer->CreateElementAt( 0 ) = this->GetInput( "MovingMask" ); } - // Get ParameterMap - ParameterObjectPointer parameterObject = static_cast< ParameterObject* >( this->GetInput( "ParameterObject" ) ); - ParameterMapListType parameterMapList = parameterObject->GetParameterMapList(); - - // Emulate command line parameters. There must be an "-out", this is checked later in code + // There must be an "-out", this is checked later in the code ArgumentMapType argumentMap; - argumentMap.insert( ArgumentMapEntryType( std::string( "-out" ), std::string( "output_path_not_set" ) ) ); + if( this->GetOutputDirectory().empty() ) { + // There must be an "-out", this is checked later in the code + argumentMap.insert( ArgumentMapEntryType( "-out", "output_path_not_set" ) ); + } + else + { + argumentMap.insert( ArgumentMapEntryType( "-out", this->GetOutputDirectory() ) ); + } + + // Fixed mesh (optional) + if( !this->m_FixedMeshFileName.empty() ) + { + argumentMap.insert( ArgumentMapEntryType( "-fp", std::string( this->m_FixedMeshFileName ) ) ); + } - // Direction Cosines - FlatDirectionCosinesType fixedImageOriginalDirection; + // Moving mesh (optional) + if( !this->m_MovingMeshFileName.empty() ) + { + argumentMap.insert( ArgumentMapEntryType( "-mp", std::string( this->m_MovingMeshFileName ) ) ); + } // Setup xout - if( elx::xoutSetup( this->GetLogToFile().c_str(), this->GetLogToFile() != std::string(), this->GetLogToConsole() ) ) + std::string logFileName; + if( this->GetLogToFile() ) + { + if( this->GetOutputDirectory().empty() ) + { + itkExceptionMacro( "LogToFileOn() requires an output directory to be specified. Use SetOutputDirectory().") + } + logFileName = this->GetOutputDirectory() + "transformix.log"; + } + + if( elx::xoutSetup( logFileName.c_str(), this->GetLogToFile(), this->GetLogToConsole() ) ) { - itkExceptionMacro( "ERROR while setting up xout." ); + // TODO: The following exception thrown if two different filters are initialized by the same process + itkExceptionMacro( "ERROR while setting up xout" ); } + + // Get ParameterMap + ParameterObjectConstPointer parameterObject = static_cast< const ParameterObject* >( this->GetInput( "ParameterObject" ) ); + ParameterMapListType parameterMapList = parameterObject->GetParameterMapList(); + // Run the (possibly multiple) registration(s) - unsigned int isError; for( unsigned int i = 0; i < parameterMapList.size(); ++i ) { // Create another instance of ElastixMain ElastixMainPointer elastix = ElastixMainType::New(); + // Set the current elastix-level + elastix->SetElastixLevel( i ); + elastix->SetTotalNumberOfElastixLevels( parameterMapList.size() ); + // Set stuff we get from a previous registration elastix->SetInitialTransform( transform ); elastix->SetFixedImageContainer( fixedImageContainer ); @@ -79,20 +118,26 @@ ElastixFilter< TFixedImage, TMovingImage, TOutputImage > elastix->SetResultImageContainer( resultImageContainer ); elastix->SetOriginalFixedImageDirectionFlat( fixedImageOriginalDirection ); - // Set the current elastix-level - elastix->SetElastixLevel( i ); - elastix->SetTotalNumberOfElastixLevels( 1 ); - - // ITK pipe-line mechanism need a result image - parameterMapList[ i ][ "WriteResultImage"] = ParameterValuesType( 1, std::string( "true" ) ); - // Start registration - isError = elastix->Run( argumentMap, parameterMapList[ i ] ); + unsigned int isError = 0; + try + { + unsigned int isError = elastix->Run( argumentMap, parameterMapList[ i ] ); + } + catch( itk::ExceptionObject &e ) + { + itkExceptionMacro( << "Errors occurred during registration: " << e.what() ); + } + + if( isError == -2 ) + { + itkExceptionMacro( << "Errors occurred during registration: Output directory does not exist." ); + } - // Assert success if( isError != 0 ) { - itkExceptionMacro( "Errors occured" ); + std::cout << std::endl << isError << std::endl; + itkExceptionMacro( << "Uncought errors occurred during registration." ); } // Get the transform, the fixedImage and the movingImage @@ -125,30 +170,30 @@ ElastixFilter< TFixedImage, TMovingImage, TOutputImage > // Save parameter map ParameterObject::Pointer TransformParameters = ParameterObject::New(); TransformParameters->SetParameterMapList( TransformParameterMapList ); - this->SetOutput( "TransformParametersObject", static_cast< itk::DataObject* >( TransformParameters ) ); + this->SetOutput( "TransformParameterObject", static_cast< itk::DataObject* >( TransformParameters ) ); // Clean up - transform = 0; - fixedImageContainer = 0; - movingImageContainer = 0; - fixedMaskContainer = 0; - movingMaskContainer = 0; - resultImageContainer = 0; + transform = ITK_NULLPTR; + fixedImageContainer = ITK_NULLPTR; + movingImageContainer = ITK_NULLPTR; + fixedMaskContainer = ITK_NULLPTR; + movingMaskContainer = ITK_NULLPTR; + resultImageContainer = ITK_NULLPTR; // Close the modules ElastixMainType::UnloadComponents(); } -template< typename TFixedImage, typename TMovingImage, typename TOutputImage > +template< typename TFixedImage, typename TMovingImage > void -ElastixFilter< TFixedImage, TMovingImage, TOutputImage > +ElastixFilter< TFixedImage, TMovingImage > ::SetFixedImage( FixedImagePointer fixedImage ) { + // Free references to fixed images that has already been set if( this->m_FixedImageContainer->Size() > 0 ) { - // Free images that has already been given - this->m_FixedImageContainer = 0; + this->m_FixedImageContainer = ITK_NULLPTR; } // Input for elastix @@ -158,9 +203,9 @@ ElastixFilter< TFixedImage, TMovingImage, TOutputImage > this->SetInput( DataObjectIdentifierType( "FixedImage" ), this->m_FixedImageContainer->ElementAt( 0 ) ); } -template< typename TFixedImage, typename TMovingImage, typename TOutputImage > +template< typename TFixedImage, typename TMovingImage > void -ElastixFilter< TFixedImage, TMovingImage, TOutputImage > +ElastixFilter< TFixedImage, TMovingImage > ::SetFixedImage( DataObjectContainerPointer fixedImages ) { // Input for elastix @@ -170,68 +215,69 @@ ElastixFilter< TFixedImage, TMovingImage, TOutputImage > this->SetInput( DataObjectIdentifierType( "FixedImage" ), this->m_FixedImageContainer->ElementAt( 0 ) ); } -template< typename TFixedImage, typename TMovingImage, typename TOutputImage > +template< typename TFixedImage, typename TMovingImage > void -ElastixFilter< TFixedImage, TMovingImage, TOutputImage > +ElastixFilter< TFixedImage, TMovingImage > ::SetMovingImage( MovingImagePointer movingImage ) { + // Free references to moving images that has already been set if( this->m_MovingImageContainer->Size() > 0 ) { - // Free images that has already been given - this->m_MovingImageContainer = 0; + + this->m_MovingImageContainer = ITK_NULLPTR; } // Input for elastix this->m_MovingImageContainer->CreateElementAt( 0 ) = static_cast< itk::DataObject* >( movingImage ) ; - // Primary input for ITK pipeline + // Input for ITK pipeline this->SetInput( DataObjectIdentifierType( "MovingImage" ), this->m_MovingImageContainer->ElementAt( 0 ) ); } -template< typename TFixedImage, typename TMovingImage, typename TOutputImage > +template< typename TFixedImage, typename TMovingImage > void -ElastixFilter< TFixedImage, TMovingImage, TOutputImage > +ElastixFilter< TFixedImage, TMovingImage > ::SetMovingImage( DataObjectContainerPointer movingImages ) { // Input for elastix this->m_MovingImageContainer = movingImages; - // Primary input for ITK pipeline + // Input for ITK pipeline this->SetInput( DataObjectIdentifierType( "MovingImage" ), this->m_MovingImageContainer->ElementAt( 0 ) ); } -template< typename TFixedImage, typename TMovingImage, typename TOutputImage > +template< typename TFixedImage, typename TMovingImage > void -ElastixFilter< TFixedImage, TMovingImage, TOutputImage > +ElastixFilter< TFixedImage, TMovingImage > ::SetParameterObject( ParameterObjectPointer parameterObject ) { this->SetInput( DataObjectIdentifierType( "ParameterObject" ), static_cast< itk::DataObject* >( parameterObject ) ); } -template< typename TFixedImage, typename TMovingImage, typename TOutputImage > +template< typename TFixedImage, typename TMovingImage > void -ElastixFilter< TFixedImage, TMovingImage, TOutputImage > +ElastixFilter< TFixedImage, TMovingImage > ::SetFixedMask( FixedImagePointer fixedMask ) { this->SetInput( DataObjectIdentifierType( "FixedMask" ), static_cast< itk::DataObject* >( fixedMask ) ); } -template< typename TFixedImage, typename TMovingImage, typename TOutputImage > +template< typename TFixedImage, typename TMovingImage > void -ElastixFilter< TFixedImage, TMovingImage, TOutputImage > +ElastixFilter< TFixedImage, TMovingImage > ::SetMovingMask( MovingImagePointer movingMask ) { this->SetInput( DataObjectIdentifierType( "MovingMask" ), static_cast< itk::DataObject* >( movingMask ) ); } -template< typename TFixedImage, typename TMovingImage, typename TOutputImage > -typename selx::ElastixFilter< TFixedImage, TMovingImage, TOutputImage >::ParameterObjectPointer -ElastixFilter< TFixedImage, TMovingImage, TOutputImage > +template< typename TFixedImage, typename TMovingImage > +typename selx::ElastixFilter< TFixedImage, TMovingImage >::ParameterObjectPointer +ElastixFilter< TFixedImage, TMovingImage > ::GetTransformParameters( void ) { - return static_cast< ParameterObject* >( itk::ProcessObject::GetOutput( DataObjectIdentifierType( "TransformParametersObject" ) ) ); + return static_cast< ParameterObject* >( itk::ProcessObject::GetOutput( DataObjectIdentifierType( "TransformParameterObject" ) ) ); } } // namespace selx -#endif // ElastixComponent_hxx \ No newline at end of file +#endif // ElastixFilter_hxx \ No newline at end of file diff --git a/Modules/Components/Elastix/include/elxTransformixFilter.h b/Modules/Components/Elastix/include/elxTransformixFilter.h new file mode 100644 index 0000000000000000000000000000000000000000..a981680a92257b5d20103fc2771ad2d7239aaa69 --- /dev/null +++ b/Modules/Components/Elastix/include/elxTransformixFilter.h @@ -0,0 +1,99 @@ +#ifndef TransformixFilter_h +#define TransformixFilter_h + +// ITK +#include "itkImageSource.h" + +// Transformix +#include "elxTransformixMain.h" + +// SuperElastix +#include "elxMacro.h" +#include "elxParameterObject.h" + +namespace selx { + +template< typename TInputImage > +class TransformixFilter : public itk::ImageSource< TInputImage > +{ +public: + + elxNewMacro( TransformixFilter, itk::ImageSource ); + + typedef elastix::TransformixMain TransformixMainType; + typedef TransformixMainType::Pointer TransformixMainPointer; + typedef TransformixMainType::ArgumentMapType ArgumentMapType; + typedef ArgumentMapType::value_type ArgumentMapEntryType; + + typedef itk::ProcessObject::DataObjectIdentifierType DataObjectIdentifierType; + typedef TransformixMainType::DataObjectContainerType DataObjectContainerType; + typedef TransformixMainType::DataObjectContainerPointer DataObjectContainerPointer; + + typedef ParameterObject::ParameterMapListType ParameterMapListType; + typedef ParameterObject::ParameterMapType ParameterMapType; + typedef ParameterObject::ParameterVectorType ParameterVectorType; + typedef typename ParameterObject::Pointer ParameterObjectPointer; + typedef typename ParameterObject::ConstPointer ParameterObjectConstPointer; + + typedef typename TInputImage::Pointer InputImagePointer; + + void SetInputImage( InputImagePointer inputImage ); + + void SetTransformParameters( ParameterObjectPointer parameterObject ); + ParameterObjectPointer GetTransformParameters( void ); + + itkSetMacro( ComputeSpatialJacobian, bool ); + itkGetConstMacro( ComputeSpatialJacobian, bool ); + itkBooleanMacro( ComputeSpatialJacobian ); + + itkSetMacro( ComputeDeterminantOfSpatialJacobian, bool ); + itkGetConstMacro( ComputeDeterminantOfSpatialJacobian, bool ); + itkBooleanMacro( ComputeDeterminantOfSpatialJacobian ); + + itkSetMacro( ComputeDeformationField, bool ); + itkGetConstMacro( ComputeDeformationField, bool ); + itkBooleanMacro( ComputeDeformationField ); + + itkSetMacro( PointSetFileName, std::string ); + itkGetConstMacro( PointSetFileName, std::string ); + void DeletePointSetFileName() { this->m_PointSetFileName = std::string(); }; + + itkSetMacro( OutputDirectory, std::string ); + itkGetConstMacro( OutputDirectory, std::string ); + void DeleteOutputDirectory() { this->m_OutputDirectory = std::string(); }; + + itkSetMacro( LogToConsole, bool ); + itkGetConstMacro( LogToConsole, bool ); + itkBooleanMacro( LogToConsole ); + + itkSetMacro( LogToFile, bool ); + itkGetConstMacro( LogToFile, bool ); + itkBooleanMacro( LogToFile ); + +protected: + + void GenerateData( void ) ITK_OVERRIDE; + +private: + + TransformixFilter(); + + bool m_ComputeSpatialJacobian; + bool m_ComputeDeterminantOfSpatialJacobian; + bool m_ComputeDeformationField; + std::string m_PointSetFileName; + + std::string m_OutputDirectory; + bool m_LogToConsole; + bool m_LogToFile; + + +}; + +} // namespace selx + +#ifndef ITK_MANUAL_INSTANTIATION +#include "elxTransformixFilter.hxx" +#endif + +#endif // TransformixFilter_h \ No newline at end of file diff --git a/Modules/Components/Elastix/include/elxTransformixFilter.hxx b/Modules/Components/Elastix/include/elxTransformixFilter.hxx new file mode 100644 index 0000000000000000000000000000000000000000..29136c2c0cf9af0696c0c8149dd9ece9c3d70616 --- /dev/null +++ b/Modules/Components/Elastix/include/elxTransformixFilter.hxx @@ -0,0 +1,177 @@ +#ifndef TransformixFilter_hxx +#define TransformixFilter_hxx + +namespace selx { + +template< typename TInputImage > +TransformixFilter< TInputImage > +::TransformixFilter( void ) +{ + this->AddRequiredInputName( "TransformParameterObject"); + + this->SetPrimaryInputName( "InputImage" ); + this->SetPrimaryOutputName( "ResultImage" ); + + this->ComputeSpatialJacobianOff(); + this->ComputeDeterminantOfSpatialJacobianOff(); + this->m_PointSetFileName = std::string(); + + this->LogToConsoleOff(); + this->LogToFileOff(); +} + +template< typename TInputImage > +void +TransformixFilter< TInputImage > +::GenerateData( void ) +{ + // Assert that at least one output has been requested + if( !this->HasInput( "InputImage" ) && + !this->GetComputeSpatialJacobian() && + !this->GetComputeDeterminantOfSpatialJacobian() && + !this->GetComputeDeformationField() && + this->GetPointSetFileName().empty() ) + { + itkExceptionMacro( << "Expected at least one of SetInputImage(), ComputeSpatialJacobianOn(), " + << "ComputeDeterminantOfSpatialJacobianOn(), SetPointSetFileName() or " ); + } + + // Check if an output directory is needed + // TODO: Change behaviour upstream to have transformix save all output to its resultImageContainer + if( ( this->GetComputeSpatialJacobian() || + this->GetComputeDeterminantOfSpatialJacobian() || + this->GetComputeDeformationField() || + !this->GetPointSetFileName().empty() || + this->GetLogToFile() ) && + this->GetOutputDirectory().empty() ) + { + itkExceptionMacro( << "The requested outputs require an output directory to be specified." + << "Use SetOutputDirectory()." ) + } + + // Transformix uses "-def" for path to point sets AND as flag for writing deformatin field + // TODO: Change behaviour upstream: Split into seperate arguments + if( this->GetComputeDeformationField() && !this->GetPointSetFileName().empty() ) + { + itkExceptionMacro( << "Only one of ComputeDeformationFieldOn() or SetPointSetFileName() can be " + << "active at any one time for backwards compatibility." ) + } + + // Setup arguments that transformix uses to figure out what needs to be done + ArgumentMapType argumentMap; + if( this->GetOutputDirectory().empty() ) { + // There must be an "-out", this is checked later in the code + argumentMap.insert( ArgumentMapEntryType( "-out", "output_path_not_set" ) ); + } + else + { + argumentMap.insert( ArgumentMapEntryType( "-out", this->GetOutputDirectory() ) ); + } + + if( this->GetComputeSpatialJacobian() ) + { + argumentMap.insert( ArgumentMapEntryType( "-jacmat", "all" ) ); + } + + if( this->GetComputeDeterminantOfSpatialJacobian() ) + { + argumentMap.insert( ArgumentMapEntryType( "-jac", "all" ) ); + } + + if( this->GetComputeDeformationField() ) + { + argumentMap.insert( ArgumentMapEntryType( "-def" , "all" ) ); + } + + if( !this->GetPointSetFileName().empty() ) + { + argumentMap.insert( ArgumentMapEntryType( "-def", this->GetPointSetFileName() ) ); + } + + // Setup xout + std::string logFileName; + if( this->GetLogToFile() ) + { + if( this->GetOutputDirectory().empty() ) + { + itkExceptionMacro( "LogToFileOn() requires an output directory to be specified. Use SetOutputDirectory().") + } + logFileName = this->GetOutputDirectory() + "transformix.log"; + } + + if( elx::xoutSetup( logFileName.c_str(), this->GetLogToFile(), this->GetLogToConsole() ) ) + { + // TODO: The following exception is thrown if two different filters are initialized by the same process + itkExceptionMacro( "ERROR while setting up xout" ); + } + + TransformixMainPointer transformix = TransformixMainType::New(); + + // Setup image containers + DataObjectContainerPointer inputImageContainer = 0; + if( this->HasInput( "InputImage" ) ) { + std::cout << "Input image is set" << std::endl; + DataObjectContainerPointer inputImageContainer = DataObjectContainerType::New(); + inputImageContainer->CreateElementAt(0) = this->GetInput("InputImage"); + } + transformix->SetMovingImageContainer( inputImageContainer ); + + DataObjectContainerPointer resultImageContainer = 0; + transformix->SetResultImageContainer( resultImageContainer ); + + // Get ParameterMap + ParameterObjectConstPointer transformParameterObject = static_cast< const ParameterObject* >( this->GetInput( "TransformParameterObject" ) ); + ParameterMapListType transformParameterMapList = transformParameterObject->GetParameterMapList(); + + // Run transformix + unsigned int isError = 0; + try + { + isError = transformix->Run( argumentMap, transformParameterMapList ); + } + catch( itk::ExceptionObject &e ) + { + itkExceptionMacro( << "Errors occured during registration: " << e.what() ); + } + + if( isError != 0 ) + { + std::cout << std::endl << isError << std::endl; + itkExceptionMacro( << "Uncought errors occured during registration." ); + } + + // Save result image + resultImageContainer = transformix->GetResultImageContainer(); + if( resultImageContainer.IsNotNull() && resultImageContainer->Size() > 0 ) + { + this->SetOutput( "ResultImage", resultImageContainer->ElementAt( 0 ) ); + } +} + +template< typename TInputImage > +void +TransformixFilter< TInputImage > +::SetInputImage( InputImagePointer inputImage ) +{ + this->SetInput( DataObjectIdentifierType( "InputImage" ), static_cast< itk::DataObject* >( inputImage ) ); +} + +template< typename TInputImage > +void +TransformixFilter< TInputImage > +::SetTransformParameters( ParameterObjectPointer parameterObject ) +{ + this->SetInput( DataObjectIdentifierType( "TransformParameterObject" ), static_cast< itk::DataObject* >( parameterObject ) ); +} + +template< typename TInputImage > +typename selx::TransformixFilter< TInputImage >::ParameterObjectPointer +TransformixFilter< TInputImage > +::GetTransformParameters( void ) +{ + return static_cast< ParameterObject* >( itk::ProcessObject::GetInput( DataObjectIdentifierType( "TransformParameterObject" ) ) ); +} + +} // namespace selx + +#endif // TransformixFilter_hxx \ No newline at end of file diff --git a/Modules/Core/Common/include/elxMacro.h b/Modules/Core/Common/include/elxMacro.h index 15f63f9ee33b2af82f4136af788092f07e9e3095..b2d2de1c8351e4d5bdb508b0437412f483497e19 100644 --- a/Modules/Core/Common/include/elxMacro.h +++ b/Modules/Core/Common/include/elxMacro.h @@ -1,5 +1,5 @@ -#ifndef __elxMacro_h -#define __elxMacro_h +#ifndef elxMacro_h +#define elxMacro_h /** * Register class with the object factory and set RTTI (Run-Time Type @@ -12,4 +12,4 @@ itkNewMacro( Self ); \ itkTypeMacro( Self, superClassName ); \ -#endif // __elxMacro_h \ No newline at end of file +#endif // elxMacro_h \ No newline at end of file diff --git a/Modules/Core/ParameterObject/include/elxParameterObject.h b/Modules/Core/ParameterObject/include/elxParameterObject.h index 3bfbe1f49cbd15f614cd49c7daa0adb7e698701e..9cb838f17127694bfdf312f8f9412578e57c40a8 100644 --- a/Modules/Core/ParameterObject/include/elxParameterObject.h +++ b/Modules/Core/ParameterObject/include/elxParameterObject.h @@ -14,9 +14,11 @@ public: elxNewMacro( ParameterObject, itk::DataObject ); - typedef std::vector< std::string > ParameterValuesType; - typedef std::map< std::string, ParameterValuesType > ParameterMapType; - typedef std::vector< ParameterMapType > ParameterMapListType; + typedef std::string ParameterKeyType; + typedef std::string ParameterValueType; + typedef std::vector< ParameterKeyType > ParameterVectorType; + typedef std::map< ParameterKeyType, ParameterVectorType > ParameterMapType; + typedef std::vector< ParameterMapType > ParameterMapListType; void SetParameterMap( ParameterMapType parameterMap ) { @@ -30,20 +32,16 @@ public: this->m_ParameterMapList = parameterMapList; }; - ParameterMapListType GetParameterMapList( void ) + ParameterMapListType& GetParameterMapList( void ) { + this->Modified(); return this->m_ParameterMapList; }; - // TODO: - // itkSetMacro( ParameterMap, ParameterMapType ) - // itkGetMacro( ParameterMap, ParameterMapType ) - - // friend ITKCommon_EXPORT std::ostream& operator<<( std::ostream& os, const ParameterObject& parameterObject ) - // { - // os << parameterObject.m_ParameterMapList; - // return os; - // } + const ParameterMapListType& GetParameterMapList( void ) const + { + return this->m_ParameterMapList; + }; private: diff --git a/SuperBuild/CMakeLists.txt b/SuperBuild/CMakeLists.txt index 0b216a8174184a8c62d742f48b41f3cb0824b267..3b25ce094423203ccf017218f4033fec70ee2c87 100644 --- a/SuperBuild/CMakeLists.txt +++ b/SuperBuild/CMakeLists.txt @@ -1,4 +1,4 @@ -cmake_minimum_required( VERSION 2.8 ) +cmake_minimum_required( VERSION 3.0.2 ) # --------------------------------------------------------------------- project( SuperElastixSuperBuild ) @@ -56,6 +56,28 @@ else() include( ExternalITK ) endif() +# --------------------------------------------------------------------- +# Build Old Elastix + +set( ELASTIX_VERSION_MAJOR "4" ) +set( ELASTIX_VERSION_MINOR "8" ) +set( ELASTIX_VERSION_STRING "${ELASTIX_VERSION_MAJOR}.${ELASTIX_VERSION_MINOR}" ) + +mark_as_advanced( SUPERELASTIX_BUILD_ELASTIX ) +option( SUPERELASTIX_BUILD_ELASTIX "Build SuperElastix support Elastix ${ELASTIX_VERSION_STRING}." ON ) +if( SUPERELASTIX_BUILD_ELASTIX ) + if( USE_SYSTEM_ELASTIX ) + if( NOT EXISTS ${ELASTIX_USE_FILE} ) + # Expose CMake variable to user + set( ELASTIX_USE_FILE ) + # Stop the build + message( FATAL_ERROR "Please point ELASTIX_USE_FILE to your systems UseElastix.cmake file." ) + endif() + else() + include( ExternalElastix ) + endif() +endif() + # --------------------------------------------------------------------- # Boost Graph Library @@ -81,24 +103,6 @@ if ( SUPERELASTIX_BUILD_TESTING ) endif() endif() -# --------------------------------------------------------------------- -# Build Elastix - -mark_as_advanced( SUPERELASTIX_BUILD_ELASTIX ) -option( SUPERELASTIX_BUILD_ELASTIX ON ) -if( SUPERELASTIX_BUILD_ELASTIX ) - if( USE_SYSTEM_ELASTIX ) - if( NOT EXISTS ${ELASTIX_USE_FILE} ) - # Expose CMake variable to user - set( ELASTIX_USE_FILE ) - # Stop the build - message( FATAL_ERROR "Please point ELASTIX_USE_FILE to your systems UseElastix.cmake file." ) - endif() - else() - include( ExternalElastix ) - endif() -endif() - # --------------------------------------------------------------------- # Build SuperElastix diff --git a/SuperBuild/ExternalElastix.cmake b/SuperBuild/ExternalElastix.cmake index c1fe6b59df0d56977f82eead1b340b63d644d1ff..8a922598183686109bb27a7e9b283e18fb369d56 100644 --- a/SuperBuild/ExternalElastix.cmake +++ b/SuperBuild/ExternalElastix.cmake @@ -1,6 +1,6 @@ set( proj Elastix ) -set( ELASTIX_REPOSITORY git://github.com/mstaring/elastix.git ) +set( ELASTIX_REPOSITORY https://github.com/mstaring/elastix.git ) set( ELASTIX_TAG c35842cd0152d8fd2894e7c6706d2dd8396f0fed ) ExternalProject_Add( ${proj} diff --git a/Testing/Unit/CMakeLists.txt b/Testing/Unit/CMakeLists.txt index 6270de443aa16e688de9ea1d79e562e6ce98df6f..0e9d8951907de462a75b703d71b6cd1756cd6fef 100644 --- a/Testing/Unit/CMakeLists.txt +++ b/Testing/Unit/CMakeLists.txt @@ -6,10 +6,16 @@ set( ElastixUnitTestFilenames elxBlueprintTest.cxx elxElastixFilterTest.cxx + elxTransformixFilterTest.cxx elxComponentFactoryTest.cxx elxComponentInterfaceTest.cxx ) +# --------------------------------------------------------------------- +# Options +if( SUPERELASTIX_BUILD_LONG_TESTS ) + add_definitions( -DSUPERELASTIX_RUN_LONG_TESTS ) +endif() # --------------------------------------------------------------------- # Set data directories @@ -44,6 +50,7 @@ list( APPEND TEST_LIBRARIES # --------------------------------------------------------------------- # Build tests +message( STATUS "ITK_LIBRARIES: ${ITK_LIBRARIES}" ) foreach( ElastixUnitTestFilename ${ElastixUnitTestFilenames} ) # Build tests executables string( REPLACE ".cxx" "" ElastixUnitTest ${ElastixUnitTestFilename} ) diff --git a/Testing/Unit/Data/Input/4D.nii.gz.md5 b/Testing/Unit/Data/Input/4D.nii.gz.md5 new file mode 100644 index 0000000000000000000000000000000000000000..6a3350182109aafc0de6a8976e5ccf192bd3b48b --- /dev/null +++ b/Testing/Unit/Data/Input/4D.nii.gz.md5 @@ -0,0 +1 @@ +58e10b6898a5d5271a890eecf25093c2 diff --git a/Testing/Unit/Data/Input/BrainProtonDensitySlice.png.md5 b/Testing/Unit/Data/Input/BrainProtonDensitySlice.png.md5 index 317f8a125578f6b03363ddc5c3ff4ce81209960e..274d1daaee9552f0b254c0e66d8e740c54d37aa7 100644 --- a/Testing/Unit/Data/Input/BrainProtonDensitySlice.png.md5 +++ b/Testing/Unit/Data/Input/BrainProtonDensitySlice.png.md5 @@ -1 +1 @@ -073df8eb397d1746d2343c78dd4bd964 \ No newline at end of file +073df8eb397d1746d2343c78dd4bd964 diff --git a/Testing/Unit/Data/Input/BrainProtonDensitySliceBorder20.png.md5 b/Testing/Unit/Data/Input/BrainProtonDensitySliceBorder20.png.md5 new file mode 100644 index 0000000000000000000000000000000000000000..6302135e8a2fc9ca1e3abff8507d6765ea90ad40 --- /dev/null +++ b/Testing/Unit/Data/Input/BrainProtonDensitySliceBorder20.png.md5 @@ -0,0 +1 @@ +9967b9c9bae3af9c3aac9676663ee05d diff --git a/Testing/Unit/Data/Input/BrainProtonDensitySliceBorder20Mask.png.md5 b/Testing/Unit/Data/Input/BrainProtonDensitySliceBorder20Mask.png.md5 new file mode 100644 index 0000000000000000000000000000000000000000..3e5cf3a48eb5af617335075c95d615cbdab909b4 --- /dev/null +++ b/Testing/Unit/Data/Input/BrainProtonDensitySliceBorder20Mask.png.md5 @@ -0,0 +1 @@ +646be5659d98e7cef40b63b3dff86726 diff --git a/Testing/Unit/Data/Input/BrainProtonDensitySliceR10X13Y17.png.md5 b/Testing/Unit/Data/Input/BrainProtonDensitySliceR10X13Y17.png.md5 index 93f3cd126feb7e33c9878ecf57b35692788d44e1..8058585db4a233381a44807b77344690a92a41f6 100644 --- a/Testing/Unit/Data/Input/BrainProtonDensitySliceR10X13Y17.png.md5 +++ b/Testing/Unit/Data/Input/BrainProtonDensitySliceR10X13Y17.png.md5 @@ -1 +1 @@ -610392a128986d934dfc0a1b0dc27e91 \ No newline at end of file +610392a128986d934dfc0a1b0dc27e91 diff --git a/Testing/Unit/Data/Input/BrainProtonDensitySliceR10X13Y17S12.png.md5 b/Testing/Unit/Data/Input/BrainProtonDensitySliceR10X13Y17S12.png.md5 new file mode 100644 index 0000000000000000000000000000000000000000..42e5aa4981dbb5804312f4fa60344d470cd2f977 --- /dev/null +++ b/Testing/Unit/Data/Input/BrainProtonDensitySliceR10X13Y17S12.png.md5 @@ -0,0 +1 @@ +eefa617b3220c1cea47bdc4672541b2d diff --git a/Testing/Unit/Data/Input/BrainProtonDensitySliceShifted13x17y.png.md5 b/Testing/Unit/Data/Input/BrainProtonDensitySliceShifted13x17y.png.md5 new file mode 100644 index 0000000000000000000000000000000000000000..98dda290136268cbbe59b6847fc7262db576f320 --- /dev/null +++ b/Testing/Unit/Data/Input/BrainProtonDensitySliceShifted13x17y.png.md5 @@ -0,0 +1 @@ +d3568c0d5c44afedb78305f620c0a79b diff --git a/Testing/Unit/Data/Input/OAS1_0001_MR1_mpr-1_anon.nrrd.md5 b/Testing/Unit/Data/Input/OAS1_0001_MR1_mpr-1_anon.nrrd.md5 new file mode 100644 index 0000000000000000000000000000000000000000..52d8753b024de51c98f35a06c0b28e5d36ec314f --- /dev/null +++ b/Testing/Unit/Data/Input/OAS1_0001_MR1_mpr-1_anon.nrrd.md5 @@ -0,0 +1 @@ +034e616e71ade1535cba4a77180ad79d \ No newline at end of file diff --git a/Testing/Unit/Data/Input/OAS1_0002_MR1_mpr-1_anon.nrrd.md5 b/Testing/Unit/Data/Input/OAS1_0002_MR1_mpr-1_anon.nrrd.md5 new file mode 100644 index 0000000000000000000000000000000000000000..151510adabe0f63adc8671662a149bf68cdf1eb1 --- /dev/null +++ b/Testing/Unit/Data/Input/OAS1_0002_MR1_mpr-1_anon.nrrd.md5 @@ -0,0 +1 @@ +49129770e7152ea5d203a519ac17be43 \ No newline at end of file diff --git a/Testing/Unit/Data/Input/cthead1-Float.mha.md5 b/Testing/Unit/Data/Input/cthead1-Float.mha.md5 deleted file mode 100644 index 815300480e124b9a3eca98655a895642354483c2..0000000000000000000000000000000000000000 --- a/Testing/Unit/Data/Input/cthead1-Float.mha.md5 +++ /dev/null @@ -1 +0,0 @@ -25de5707b18c0c684fd5fa30351bf787 \ No newline at end of file diff --git a/Testing/Unit/elxElastixFilterTest.cxx b/Testing/Unit/elxElastixFilterTest.cxx index 6dac45f0457d4ea648045eefedbf9c59e265a677..29fdc2e3cc79cc1f9957c44016a70faf382d7695 100644 --- a/Testing/Unit/elxElastixFilterTest.cxx +++ b/Testing/Unit/elxElastixFilterTest.cxx @@ -7,6 +7,15 @@ #include "elxDataManager.h" #include "gtest/gtest.h" +// TODO: +// - In the following examples, why do we need to call update on both reader and elastixFilter +// before using the writer? Should the call to update on the writer not be propagated upstream? +// (See http://itk.org/Wiki/ITK/Examples/Segmentation/OtsuThresholdImageFilter which defeats the whole +// purpose of having a pipeline, no?) +// - Compare result images against baseline +// - SetUp() runs before every test. Does GoogleTest have a function that is called once before tests are run? + + using namespace selx; class ElastixFilterTest : public ::testing::Test @@ -15,141 +24,368 @@ protected: typedef DataManager DataManagerType; - typedef itk::Image< float, 2u > ImageType; - ImageType::Pointer fixedImage; - ImageType::Pointer movingImage; - ImageType::Pointer resultImage; - - typedef itk::ImageFileReader< ImageType > ImageReaderType; - typedef itk::ImageFileWriter< ImageType > ImageWriterType; - - typedef ParameterObject::ParameterValuesType ParameterValuesType; - typedef ParameterObject::ParameterMapType ParameterMapType; - ParameterObject::Pointer parameterObject; - ParameterObject::Pointer TransformParameterObject; - - typedef ElastixFilter< ImageType, ImageType, ImageType > ElastixFilterType; - typedef ElastixFilterType::DataObjectContainerType DataObjectContainerType; - typedef ElastixFilterType::DataObjectContainerPointer DataObjectContainerPointer; + // Parameter typedefs + typedef ParameterObject::ParameterVectorType ParameterVectorType; + typedef ParameterObject::ParameterMapType ParameterMapType; + ParameterMapType parameterMap; + ParameterObject::Pointer eulerTransformParameterObject; + ParameterObject::Pointer affineTransformParameterObject; + ParameterObject::Pointer correspondingPointsParameterObject; + ParameterObject::Pointer nonRigidParameterObject; + ParameterObject::Pointer nonRigidParameterObject3D; + ParameterObject::Pointer groupwiseParameterObject; + ParameterObject::ConstPointer transformParameterObject; virtual void SetUp() { - DataManagerType::Pointer dataManager = DataManagerType::New(); - - ImageReaderType::Pointer reader = ImageReaderType::New(); - - reader->SetFileName( dataManager->GetInputFile( "cthead1-Float.mha" ) ); - fixedImage = reader->GetOutput(); - - // TODO: Only call update on writer - reader->Update(); - - reader->SetFileName( dataManager->GetInputFile( "cthead1-Float.mha" ) ); - movingImage = reader->GetOutput(); - - // TODO: Only call update on writer - reader->Update(); - - // ParameterMap + // Nonrigid ParameterMap ParameterMapType parameterMap = ParameterMapType(); // Images - parameterMap[ "FixedInternalImagePixelType" ] = ParameterValuesType( 1, "float" ); - parameterMap[ "FixedImageDimension" ] = ParameterValuesType( 1, "2" ); - parameterMap[ "MovingInternalImagePixelType" ] = ParameterValuesType( 1, "float" ); - parameterMap[ "MovingImageDimension" ] = ParameterValuesType( 1, "2" ); - parameterMap[ "ResultImagePixelType" ] = ParameterValuesType( 1, "float" ); + parameterMap[ "FixedInternalImagePixelType" ] = ParameterVectorType( 1, "float" ); + parameterMap[ "FixedImageDimension" ] = ParameterVectorType( 1, "2" ); + parameterMap[ "MovingInternalImagePixelType" ] = ParameterVectorType( 1, "float" ); + parameterMap[ "MovingImageDimension" ] = ParameterVectorType( 1, "2" ); + parameterMap[ "ResultImagePixelType" ] = ParameterVectorType( 1, "float" ); // Common components - parameterMap[ "FixedImagePyramid" ] = ParameterValuesType( 1, "FixedSmoothingImagePyramid" ); - parameterMap[ "MovingImagePyramid" ] = ParameterValuesType( 1, "MovingSmoothingImagePyramid" ); - parameterMap[ "Interpolator"] = ParameterValuesType( 1, "LinearInterpolator"); - parameterMap[ "Optimizer" ] = ParameterValuesType( 1, "AdaptiveStochasticGradientDescent" ); - parameterMap[ "Resampler"] = ParameterValuesType( 1, "DefaultResampler" ); - parameterMap[ "ResampleInterpolator"] = ParameterValuesType( 1, "FinalLinearInterpolator" ); - parameterMap[ "FinalBSplineInterpolationOrder" ] = ParameterValuesType( 1, "2" ); - parameterMap[ "NumberOfResolutions" ] = ParameterValuesType( 1, "2" ); + parameterMap[ "FixedImagePyramid" ] = ParameterVectorType( 1, "FixedSmoothingImagePyramid" ); + parameterMap[ "MovingImagePyramid" ] = ParameterVectorType( 1, "MovingSmoothingImagePyramid" ); + parameterMap[ "Interpolator"] = ParameterVectorType( 1, "LinearInterpolator"); + parameterMap[ "Optimizer" ] = ParameterVectorType( 1, "AdaptiveStochasticGradientDescent" ); + parameterMap[ "Resampler"] = ParameterVectorType( 1, "DefaultResampler" ); + parameterMap[ "ResampleInterpolator"] = ParameterVectorType( 1, "FinalLinearInterpolator" ); + parameterMap[ "FinalBSplineInterpolationOrder" ] = ParameterVectorType( 1, "2" ); + parameterMap[ "NumberOfResolutions" ] = ParameterVectorType( 1, "2" ); // Image Sampler - parameterMap[ "ImageSampler" ] = ParameterValuesType( 1, "RandomCoordinate" ); - parameterMap[ "NumberOfSpatialSamples"] = ParameterValuesType( 1, "2048" ); - parameterMap[ "CheckNumberOfSamples" ] = ParameterValuesType( 1, "true" ); - parameterMap[ "MaximumNumberOfSamplingAttempts" ] = ParameterValuesType( 1, "8" ); - parameterMap[ "NewSamplesEveryIteration" ] = ParameterValuesType( 1, "true"); + parameterMap[ "ImageSampler" ] = ParameterVectorType( 1, "RandomCoordinate" ); + parameterMap[ "NumberOfSpatialSamples"] = ParameterVectorType( 1, "2048" ); + parameterMap[ "CheckNumberOfSamples" ] = ParameterVectorType( 1, "true" ); + parameterMap[ "MaximumNumberOfSamplingAttempts" ] = ParameterVectorType( 1, "8" ); + parameterMap[ "NewSamplesEveryIteration" ] = ParameterVectorType( 1, "true"); // Optimizer - parameterMap[ "NumberOfSamplesForExactGradient" ] = ParameterValuesType( 1, "4096" ); - parameterMap[ "DefaultPixelValue" ] = ParameterValuesType( 1, "0" ); - parameterMap[ "AutomaticParameterEstimation" ] = ParameterValuesType( 1, "true" ); + parameterMap[ "NumberOfSamplesForExactGradient" ] = ParameterVectorType( 1, "4096" ); + parameterMap[ "DefaultPixelValue" ] = ParameterVectorType( 1, "0" ); + parameterMap[ "AutomaticParameterEstimation" ] = ParameterVectorType( 1, "true" ); // Output - parameterMap[ "WriteResultImage" ] = ParameterValuesType( 1, "false" ); - parameterMap[ "ResultImageFormat" ] = ParameterValuesType( 1, "nii" ); + parameterMap[ "WriteResultImage" ] = ParameterVectorType( 1, "true" ); + parameterMap[ "ResultImageFormat" ] = ParameterVectorType( 1, "nii" ); // Registration - parameterMap[ "Registration" ] = ParameterValuesType( 1, "MultiResolutionRegistration" ); - parameterMap[ "Transform" ] = ParameterValuesType( 1, "EulerTransform" ); - parameterMap[ "Metric" ] = ParameterValuesType( 1, "AdvancedMattesMutualInformation" ); - parameterMap[ "MaximumNumberOfIterations" ] = ParameterValuesType( 1, "128" ); - - parameterObject = ParameterObject::New(); - parameterObject->SetParameterMap( parameterMap ); + parameterMap[ "Registration" ] = ParameterVectorType( 1, "MultiResolutionRegistration" ); + parameterMap[ "Transform" ] = ParameterVectorType( 1, "EulerTransform" ); + parameterMap[ "Metric" ] = ParameterVectorType( 1, "AdvancedMattesMutualInformation" ); + parameterMap[ "MaximumNumberOfIterations" ] = ParameterVectorType( 1, "128" ); + + eulerTransformParameterObject = ParameterObject::New(); + eulerTransformParameterObject->SetParameterMap( parameterMap ); + + ParameterMapType affineTransformParameterMap = parameterMap; + affineTransformParameterMap[ "Transform" ] = ParameterVectorType( 1, "AffineTransform" ); + affineTransformParameterObject = ParameterObject::New(); + affineTransformParameterObject->SetParameterMap( affineTransformParameterMap ); + + ParameterMapType correspondingPointsParameterMap = parameterMap; + correspondingPointsParameterMap[ "Registration" ] = ParameterVectorType( 1, "MultiMetricMultiResolutionRegistration" ); + correspondingPointsParameterMap[ "Transform" ] = ParameterVectorType( 1, "TranslationTransform" ); + correspondingPointsParameterMap[ "Metric" ].push_back( "CorrespondingPointsEuclideanDistanceMetric" ); + correspondingPointsParameterMap[ "Metric0Weight"] = ParameterVectorType( 1, "0.0" ); + correspondingPointsParameterObject = ParameterObject::New(); + correspondingPointsParameterObject->SetParameterMap( correspondingPointsParameterMap ); + + ParameterMapType nonRigidParameterMap = parameterMap; + nonRigidParameterMap[ "Transform" ] = ParameterVectorType( 1, "BSplineTransform" ); + nonRigidParameterObject = ParameterObject::New(); + nonRigidParameterObject->SetParameterMap( nonRigidParameterMap ); + + ParameterMapType nonRigidParameterMap3D = parameterMap; + nonRigidParameterMap3D[ "FixedInternalImagePixelType" ] = ParameterVectorType( 1, "float" ); + nonRigidParameterMap3D[ "FixedImageDimension" ] = ParameterVectorType( 1, "3" ); + nonRigidParameterMap3D[ "MovingInternalImagePixelType" ] = ParameterVectorType( 1, "float" ); + nonRigidParameterMap3D[ "MovingImageDimension" ] = ParameterVectorType( 1, "3" ); + nonRigidParameterMap3D[ "ResultImagePixelType" ] = ParameterVectorType( 1, "float" ); + nonRigidParameterMap3D[ "Transform" ] = ParameterVectorType( 1, "BSplineTransform" ); + nonRigidParameterMap3D[ "MaximumNumberOfIterations" ] = ParameterVectorType( 1, "8" ); + nonRigidParameterObject3D = ParameterObject::New(); + nonRigidParameterObject3D->SetParameterMap( nonRigidParameterMap3D ); + + ParameterMapType groupwiseParameterMap = parameterMap; + groupwiseParameterMap[ "FixedInternalImagePixelType" ] = ParameterVectorType( 1, "float" ); + groupwiseParameterMap[ "FixedImageDimension" ] = ParameterVectorType( 1, "4" ); + groupwiseParameterMap[ "MovingInternalImagePixelType" ] = ParameterVectorType( 1, "float" ); + groupwiseParameterMap[ "MovingImageDimension" ] = ParameterVectorType( 1, "4" ); + groupwiseParameterMap[ "ResultImagePixelType" ] = ParameterVectorType( 1, "float" ); + groupwiseParameterMap[ "Transform" ] = ParameterVectorType( 1, "BSplineStackTransform" ); + groupwiseParameterMap[ "Metric" ] = ParameterVectorType( 1, "VarianceOverLastDimensionMetric" ); + groupwiseParameterMap[ "MaximumNumberOfIterations" ] = ParameterVectorType( 1, "512" ); + groupwiseParameterMap[ "Interpolator"] = ParameterVectorType( 1, "ReducedDimensionBSplineInterpolator" ); + groupwiseParameterMap[ "ResampleInterpolator" ] = ParameterVectorType( 1, "FinalReducedDimensionBSplineInterpolator" ); + groupwiseParameterMap[ "MaximumNumberOfIterations" ] = ParameterVectorType( 1, "8" ); + groupwiseParameterObject = ParameterObject::New(); + groupwiseParameterObject->SetParameterMap( nonRigidParameterMap ); } }; TEST_F( ElastixFilterTest, Instantiation ) { + typedef itk::Image< float, 2 > ImageType; + typedef ElastixFilter< ImageType, ImageType > ElastixFilterType; EXPECT_NO_THROW( ElastixFilterType::Pointer elastixFilter = ElastixFilterType::New() ); } -TEST_F( ElastixFilterTest, PairwiseRegistration ) +TEST_F( ElastixFilterTest, Euler2D ) { - ElastixFilterType::Pointer elastixFilter = ElastixFilterType::New(); + typedef itk::Image< float, 2 > ImageType; + typedef itk::ImageFileReader< ImageType > ImageFileReaderType; + typedef itk::ImageFileWriter< ImageType > ImageFileWriterType; + typedef ElastixFilter< ImageType, ImageType > ElastixFilterType; + + DataManagerType::Pointer dataManager = DataManagerType::New(); - elastixFilter->LogToConsoleOn(); - elastixFilter->SetFixedImage( fixedImage ); - elastixFilter->SetMovingImage( movingImage ); - elastixFilter->SetParameterObject( parameterObject ); + ImageFileReaderType::Pointer fixedImageReader = ImageFileReaderType::New(); + fixedImageReader->SetFileName( dataManager->GetInputFile( "BrainProtonDensitySliceBorder20.png" ) ); + fixedImageReader->Update(); - // TODO: This update should not be needed - elastixFilter->Update(); + ImageFileReaderType::Pointer movingImageReader = ImageFileReaderType::New(); + movingImageReader->SetFileName( dataManager->GetInputFile( "BrainProtonDensitySliceR10X13Y17.png" ) ); + movingImageReader->Update(); - EXPECT_NO_THROW( resultImage = elastixFilter->GetOutput() ); - EXPECT_NO_THROW( TransformParameterObject = elastixFilter->GetTransformParameters() ); + ElastixFilterType::Pointer elastixFilter; + EXPECT_NO_THROW( elastixFilter = ElastixFilterType::New() ); + EXPECT_NO_THROW( elastixFilter->LogToConsoleOn() ); + EXPECT_NO_THROW( elastixFilter->SetFixedImage( fixedImageReader->GetOutput() ) ); + EXPECT_NO_THROW( elastixFilter->SetMovingImage( movingImageReader->GetOutput() ) ); + EXPECT_NO_THROW( elastixFilter->SetParameterObject( eulerTransformParameterObject ) ); - ImageWriterType::Pointer writer = ImageWriterType::New(); - writer->SetFileName( "ElastixResultImage.nii" ); + // TODO: This update should not be needed (see description above) + EXPECT_NO_THROW( elastixFilter->Update() ); + + ImageFileWriterType::Pointer writer = ImageFileWriterType::New(); + writer->SetFileName( dataManager->GetOutputFile( "Euler2DResultImage.nii" ) ); writer->SetInput( elastixFilter->GetOutput() ); EXPECT_NO_THROW( writer->Update() ); + EXPECT_NO_THROW( transformParameterObject = elastixFilter->GetTransformParameters() ); } -TEST_F( ElastixFilterTest, MultiPairwiseRegistration ) +TEST_F( ElastixFilterTest, AffineWithMultipleFixedAndMovingImages2D ) + { + typedef itk::Image< float, 2 > ImageType; + typedef itk::ImageFileReader< ImageType > ImageFileReaderType; + typedef itk::ImageFileWriter< ImageType > ImageFileWriterType; + typedef ElastixFilter< ImageType, ImageType > ElastixFilterType; + typedef ElastixFilterType::DataObjectContainerType DataObjectContainerType; + typedef ElastixFilterType::DataObjectContainerPointer DataObjectContainerPointer; + + DataManagerType::Pointer dataManager = DataManagerType::New(); + + ImageFileReaderType::Pointer fixedImageReader0 = ImageFileReaderType::New(); + fixedImageReader0->SetFileName( dataManager->GetInputFile( "BrainProtonDensitySliceBorder20.png" ) ); + fixedImageReader0->Update(); + + ImageFileReaderType::Pointer fixedImageReader1 = ImageFileReaderType::New(); + fixedImageReader1->SetFileName( dataManager->GetInputFile( "BrainProtonDensitySliceBorder20.png" ) ); + fixedImageReader1->Update(); + + ImageFileReaderType::Pointer movingImageReader0 = ImageFileReaderType::New(); + movingImageReader0->SetFileName( dataManager->GetInputFile( "BrainProtonDensitySliceR10X13Y17S12.png" ) ); + movingImageReader0->Update(); + + ImageFileReaderType::Pointer movingImageReader1 = ImageFileReaderType::New(); + movingImageReader1->SetFileName( dataManager->GetInputFile( "BrainProtonDensitySliceR10X13Y17S12.png" ) ); + movingImageReader1->Update(); + + DataObjectContainerPointer fixedImages = DataObjectContainerType::New(); + fixedImages->CreateElementAt( 0 ) = static_cast< itk::DataObject* >( fixedImageReader0->GetOutput() ); + fixedImages->CreateElementAt( 1 ) = static_cast< itk::DataObject* >( fixedImageReader1->GetOutput() ); + + DataObjectContainerPointer movingImages = DataObjectContainerType::New(); + movingImages->CreateElementAt( 0 ) = static_cast< itk::DataObject* >( movingImageReader0->GetOutput() ); + movingImages->CreateElementAt( 1 ) = static_cast< itk::DataObject* >( movingImageReader1->GetOutput() ); + + ElastixFilterType::Pointer elastixFilter = ElastixFilterType::New(); + EXPECT_NO_THROW( elastixFilter->LogToConsoleOn() ); + EXPECT_NO_THROW( elastixFilter->SetFixedImage( fixedImages ) ); + EXPECT_NO_THROW( elastixFilter->SetMovingImage( movingImages ) ); + EXPECT_NO_THROW( elastixFilter->SetParameterObject( affineTransformParameterObject ) ); + + // TODO: This update should not be needed (see description above) + EXPECT_NO_THROW( elastixFilter->Update() ); + + ImageFileWriterType::Pointer writer = ImageFileWriterType::New(); + EXPECT_NO_THROW( writer->SetFileName( dataManager->GetOutputFile( "AffineWithMultipleFixedAndMovingImages2DResultImage.nii" ) ) ); + EXPECT_NO_THROW( writer->SetInput( elastixFilter->GetOutput() ) ); + EXPECT_NO_THROW( writer->Update() ); + + EXPECT_NO_THROW( transformParameterObject = elastixFilter->GetTransformParameters() ); + } + +TEST_F( ElastixFilterTest, TranslationWithPointSets2D ) { + typedef itk::Image< float, 2 > ImageType; + typedef itk::ImageFileReader< ImageType > ImageFileReaderType; + typedef itk::ImageFileWriter< ImageType > ImageFileWriterType; + typedef ElastixFilter< ImageType, ImageType > ElastixFilterType; + + DataManagerType::Pointer dataManager = DataManagerType::New(); + + // We generate the point sets manually + std::ofstream fixedMeshFile; + fixedMeshFile.open( dataManager->GetInputFile( "FixedMesh.pts" )); + fixedMeshFile << "point\n"; + fixedMeshFile << "1\n"; + fixedMeshFile << "128.0 128.0\n"; + fixedMeshFile.close(); + + std::ofstream movingMeshFile; + movingMeshFile.open( dataManager->GetInputFile( "MovingMesh.pts" )); + movingMeshFile << "point\n"; + movingMeshFile << "1\n"; + movingMeshFile << "115.0 111.0\n"; + movingMeshFile.close(); + + ImageFileReaderType::Pointer fixedImageReader = ImageFileReaderType::New(); + fixedImageReader->SetFileName( dataManager->GetInputFile( "BrainProtonDensitySliceBorder20.png" ) ); + fixedImageReader->Update(); + + ImageFileReaderType::Pointer movingImageReader = ImageFileReaderType::New(); + movingImageReader->SetFileName( dataManager->GetInputFile( "BrainProtonDensitySliceShifted13x17y.png" ) ); + movingImageReader->Update(); + ElastixFilterType::Pointer elastixFilter = ElastixFilterType::New(); + EXPECT_NO_THROW( elastixFilter->LogToConsoleOn() ); + EXPECT_NO_THROW( elastixFilter->SetFixedImage( fixedImageReader->GetOutput() ) ); + EXPECT_NO_THROW( elastixFilter->SetFixedMeshFileName( dataManager->GetInputFile( "FixedMesh.pts" ) ) ); + EXPECT_NO_THROW( elastixFilter->SetMovingImage( movingImageReader->GetOutput() ) ); + EXPECT_NO_THROW( elastixFilter->SetMovingMeshFileName( dataManager->GetInputFile( "MovingMesh.pts" ) ) ); + EXPECT_NO_THROW( elastixFilter->SetParameterObject( correspondingPointsParameterObject ) ); + + // TODO: This update should not be needed (see description above) + EXPECT_NO_THROW( elastixFilter->Update() ); + + ImageFileWriterType::Pointer writer = ImageFileWriterType::New(); + writer->SetFileName( "TranslationWithPointSets2DResultImage.nii" ); + writer->SetInput( elastixFilter->GetOutput() ); + EXPECT_NO_THROW( writer->Update() ); - DataObjectContainerPointer fixedImages = DataObjectContainerType::New(); - fixedImages->CreateElementAt( 0 ) = static_cast< itk::DataObject* >( fixedImage ); - fixedImages->CreateElementAt( 1 ) = static_cast< itk::DataObject* >( fixedImage ); + EXPECT_NO_THROW( transformParameterObject = elastixFilter->GetTransformParameters() ); + } - DataObjectContainerPointer movingImages = DataObjectContainerType::New(); - movingImages->CreateElementAt( 0 ) = static_cast< itk::DataObject* >( movingImage ); - movingImages->CreateElementAt( 1 ) = static_cast< itk::DataObject* >( movingImage ); +TEST_F( ElastixFilterTest, BSplineWithFixedMask2D ) +{ + typedef itk::Image< float, 2 > ImageType; + typedef itk::ImageFileReader< ImageType > ImageFileReaderType; + typedef itk::ImageFileWriter< ImageType > ImageFileWriterType; + typedef ElastixFilter< ImageType, ImageType > ElastixFilterType; + + DataManagerType::Pointer dataManager = DataManagerType::New(); + + ImageFileReaderType::Pointer fixedImageReader = ImageFileReaderType::New(); + fixedImageReader->SetFileName( dataManager->GetInputFile( "BrainProtonDensitySliceBorder20.png" ) ); + fixedImageReader->Update(); + + ImageFileReaderType::Pointer fixedMaskReader = ImageFileReaderType::New(); + fixedMaskReader->SetFileName( dataManager->GetInputFile( "BrainProtonDensitySliceBorder20Mask.png" ) ); + fixedMaskReader->Update(); + + ImageFileReaderType::Pointer movingImageReader = ImageFileReaderType::New(); + movingImageReader->SetFileName( dataManager->GetInputFile( "BrainProtonDensitySliceR10X13Y17.png" ) ); + movingImageReader->Update(); + + ElastixFilterType::Pointer elastixFilter; + EXPECT_NO_THROW( elastixFilter = ElastixFilterType::New() ); + EXPECT_NO_THROW( elastixFilter->LogToConsoleOn() ); + EXPECT_NO_THROW( elastixFilter->SetFixedImage( fixedImageReader->GetOutput() ) ); + EXPECT_NO_THROW( elastixFilter->SetFixedMask( fixedMaskReader->GetOutput() ) ); + EXPECT_NO_THROW( elastixFilter->SetMovingImage( movingImageReader->GetOutput() ) ); + EXPECT_NO_THROW( elastixFilter->SetParameterObject( nonRigidParameterObject ) ); + + // TODO: This update should not be needed (see description above) + EXPECT_NO_THROW( elastixFilter->Update() ); + + ImageFileWriterType::Pointer writer = ImageFileWriterType::New(); + writer->SetFileName( dataManager->GetOutputFile( "BSplineWithFixedMask2DResultImage.nii" ) ); + writer->SetInput( elastixFilter->GetOutput() ); + EXPECT_NO_THROW( writer->Update() ); + + EXPECT_NO_THROW( transformParameterObject = elastixFilter->GetTransformParameters() ); +} + +#ifdef SUPERELASTIX_RUN_LONG_TESTS + +TEST_F( ElastixFilterTest, BSpline3D ) +{ + typedef itk::Image< float, 3 > ImageType; + typedef itk::ImageFileReader< ImageType > ImageFileReaderType; + typedef itk::ImageFileWriter< ImageType > ImageFileWriterType; + typedef ElastixFilter< ImageType, ImageType > ElastixFilterType; + + DataManagerType::Pointer dataManager = DataManagerType::New(); + + ImageFileReaderType::Pointer fixedImageReader = ImageFileReaderType::New(); + fixedImageReader->SetFileName( dataManager->GetInputFile( "OAS1_0001_MR1_mpr-1_anon.nrrd" ) ); + fixedImageReader->Update(); - elastixFilter->LogToConsoleOn(); - elastixFilter->SetFixedImage( fixedImages ); - elastixFilter->SetMovingImage( movingImages ); - elastixFilter->SetParameterObject( parameterObject ); + ImageFileReaderType::Pointer movingImageReader = ImageFileReaderType::New(); + movingImageReader->SetFileName( dataManager->GetInputFile( "OAS1_0002_MR1_mpr-1_anon.nrrd" ) ); + movingImageReader->Update(); - // TODO: This update should not be needed - elastixFilter->Update(); + ElastixFilterType::Pointer elastixFilter; + EXPECT_NO_THROW( elastixFilter = ElastixFilterType::New() ); + EXPECT_NO_THROW( elastixFilter->LogToConsoleOn() ); + EXPECT_NO_THROW( elastixFilter->SetFixedImage( fixedImageReader->GetOutput() ) ); + EXPECT_NO_THROW( elastixFilter->SetMovingImage( movingImageReader->GetOutput() ) ); + EXPECT_NO_THROW( elastixFilter->SetParameterObject( nonRigidParameterObject3D ) ); - EXPECT_NO_THROW( resultImage = elastixFilter->GetOutput() ); - EXPECT_NO_THROW( TransformParameterObject = elastixFilter->GetTransformParameters() ); + // TODO: This update should not be needed (see description above) + EXPECT_NO_THROW( elastixFilter->Update() ); - ImageWriterType::Pointer writer = ImageWriterType::New(); - writer->SetFileName( "ElastixResultImage.nii" ); + ImageFileWriterType::Pointer writer = ImageFileWriterType::New(); + writer->SetFileName( dataManager->GetOutputFile( "BSpline3DResultImage.nii" ) ); writer->SetInput( elastixFilter->GetOutput() ); EXPECT_NO_THROW( writer->Update() ); + + EXPECT_NO_THROW( transformParameterObject = elastixFilter->GetTransformParameters() ); } + +// TODO: Find out why this segfaults +//TEST_F( ElastixFilterTest, BSpline4D ) +//{ +// typedef itk::Image< float, 4 > ImageType; +// typedef itk::ImageFileReader< ImageType > ImageFileReaderType; +// typedef itk::ImageFileWriter< ImageType > ImageFileWriterType; +// typedef ElastixFilter< ImageType, ImageType > ElastixFilterType; +// +// DataManagerType::Pointer dataManager = DataManagerType::New(); +// +// ImageFileReaderType::Pointer fixedImageReader = ImageFileReaderType::New(); +// fixedImageReader->SetFileName( dataManager->GetInputFile( "4D.nii.gz" ) ); +// fixedImageReader->Update(); +// +// ImageFileReaderType::Pointer movingImageReader = ImageFileReaderType::New(); +// movingImageReader->SetFileName( dataManager->GetInputFile( "4D.nii.gz" ) ); +// movingImageReader->Update(); +// +// ElastixFilterType::Pointer elastixFilter; +// EXPECT_NO_THROW( elastixFilter = ElastixFilterType::New() ); +// EXPECT_NO_THROW( elastixFilter->LogToConsoleOn() ); +// EXPECT_NO_THROW( elastixFilter->SetFixedImage( fixedImageReader->GetOutput() ) ); +// EXPECT_NO_THROW( elastixFilter->SetMovingImage( movingImageReader->GetOutput() ) ); +// EXPECT_NO_THROW( elastixFilter->SetParameterObject( groupwiseParameterObject ) ); +// +// // TODO: This update should not be needed (see description above) +// EXPECT_NO_THROW( elastixFilter->Update() ); +// +// ImageFileWriterType::Pointer writer = ImageFileWriterType::New(); +// writer->SetFileName( dataManager->GetOutputFile( "BSpline4DResultImage.nii" ) ); +// writer->SetInput( elastixFilter->GetOutput() ); +// EXPECT_NO_THROW( writer->Update() ); +// +// EXPECT_NO_THROW( transformParameterObject = elastixFilter->GetTransformParameters() ); +//} + +#endif // SUPERELASTIX_RUN_LONG_TESTS diff --git a/Testing/Unit/elxTransformixFilterTest.cxx b/Testing/Unit/elxTransformixFilterTest.cxx new file mode 100644 index 0000000000000000000000000000000000000000..cfe47b1462947cca8ca15244220f158420c3b038 --- /dev/null +++ b/Testing/Unit/elxTransformixFilterTest.cxx @@ -0,0 +1,134 @@ +#include "elxElastixFilter.h" +#include "elxTransformixFilter.h" +#include "elxParameterObject.h" + +#include "itkImageFileReader.h" +#include "itkImageFileWriter.h" + +#include "elxDataManager.h" +#include "gtest/gtest.h" + +using namespace selx; + +class TransformixFilterTest : public ::testing::Test +{ +protected: + + typedef DataManager DataManagerType; + + typedef ParameterObject::ParameterVectorType ParameterVectorType; + typedef ParameterObject::ParameterMapType ParameterMapType; + + typedef itk::Image< float, 2 > ImageType; + typedef itk::ImageFileReader< ImageType > ImageFileReaderType; + typedef itk::ImageFileWriter< ImageType > ImageFileWriterType; + typedef ElastixFilter< ImageType, ImageType > ElastixFilterType; + typedef TransformixFilter< ImageType > TransformixFilterType; + + ElastixFilterType::Pointer elastixFilter; + ParameterObject::Pointer eulerTransformParameterObject; + + virtual void SetUp() + { + // Nonrigid ParameterMap + ParameterMapType parameterMap = ParameterMapType(); + + // Images + parameterMap[ "FixedInternalImagePixelType" ] = ParameterVectorType( 1, "float" ); + parameterMap[ "FixedImageDimension" ] = ParameterVectorType( 1, "2" ); + parameterMap[ "MovingInternalImagePixelType" ] = ParameterVectorType( 1, "float" ); + parameterMap[ "MovingImageDimension" ] = ParameterVectorType( 1, "2" ); + parameterMap[ "ResultImagePixelType" ] = ParameterVectorType( 1, "float" ); + + // Common components + parameterMap[ "FixedImagePyramid" ] = ParameterVectorType( 1, "FixedSmoothingImagePyramid" ); + parameterMap[ "MovingImagePyramid" ] = ParameterVectorType( 1, "MovingSmoothingImagePyramid" ); + parameterMap[ "Interpolator"] = ParameterVectorType( 1, "LinearInterpolator"); + parameterMap[ "Optimizer" ] = ParameterVectorType( 1, "AdaptiveStochasticGradientDescent" ); + parameterMap[ "Resampler"] = ParameterVectorType( 1, "DefaultResampler" ); + parameterMap[ "ResampleInterpolator"] = ParameterVectorType( 1, "FinalLinearInterpolator" ); + parameterMap[ "FinalBSplineInterpolationOrder" ] = ParameterVectorType( 1, "2" ); + parameterMap[ "NumberOfResolutions" ] = ParameterVectorType( 1, "2" ); + + // Image Sampler + parameterMap[ "ImageSampler" ] = ParameterVectorType( 1, "RandomCoordinate" ); + parameterMap[ "NumberOfSpatialSamples"] = ParameterVectorType( 1, "2048" ); + parameterMap[ "CheckNumberOfSamples" ] = ParameterVectorType( 1, "true" ); + parameterMap[ "MaximumNumberOfSamplingAttempts" ] = ParameterVectorType( 1, "8" ); + parameterMap[ "NewSamplesEveryIteration" ] = ParameterVectorType( 1, "true"); + + // Optimizer + parameterMap[ "NumberOfSamplesForExactGradient" ] = ParameterVectorType( 1, "4096" ); + parameterMap[ "DefaultPixelValue" ] = ParameterVectorType( 1, "0" ); + parameterMap[ "AutomaticParameterEstimation" ] = ParameterVectorType( 1, "true" ); + + // Output + parameterMap[ "WriteResultImage" ] = ParameterVectorType( 1, "true" ); + parameterMap[ "ResultImageFormat" ] = ParameterVectorType( 1, "nii" ); + + // Registration + parameterMap[ "Registration" ] = ParameterVectorType( 1, "MultiResolutionRegistration" ); + parameterMap[ "Transform" ] = ParameterVectorType( 1, "EulerTransform" ); + parameterMap[ "Metric" ] = ParameterVectorType( 1, "AdvancedMattesMutualInformation" ); + parameterMap[ "MaximumNumberOfIterations" ] = ParameterVectorType( 1, "128" ); + + eulerTransformParameterObject = ParameterObject::New(); + eulerTransformParameterObject->SetParameterMap( parameterMap ); + + // We generate the point sets manually + DataManagerType::Pointer dataManager = DataManagerType::New(); + std::ofstream inputMeshFile; + inputMeshFile.open( dataManager->GetInputFile( "InputMesh.pts" )); + inputMeshFile << "point\n"; + inputMeshFile << "1\n"; + inputMeshFile << "115.0 111.0\n"; + inputMeshFile.close(); + } +}; + +TEST_F( TransformixFilterTest, Instantiation ) +{ + typedef itk::Image< float, 2 > ImageType; + typedef ElastixFilter< ImageType, ImageType > ElastixFilterType; + EXPECT_NO_THROW( ElastixFilterType::Pointer elastixFilter = ElastixFilterType::New() ); +} + +TEST_F( TransformixFilterTest, Euler2D ) +{ + DataManagerType::Pointer dataManager = DataManagerType::New(); + + ImageFileReaderType::Pointer fixedImageReader = ImageFileReaderType::New(); + fixedImageReader->SetFileName( dataManager->GetInputFile( "BrainProtonDensitySliceBorder20.png" ) ); + fixedImageReader->Update(); + + ImageFileReaderType::Pointer movingImageReader = ImageFileReaderType::New(); + movingImageReader->SetFileName( dataManager->GetInputFile( "BrainProtonDensitySliceR10X13Y17.png" ) ); + movingImageReader->Update(); + + ElastixFilterType::Pointer elastixFilter; + EXPECT_NO_THROW( elastixFilter = ElastixFilterType::New() ); + EXPECT_NO_THROW( elastixFilter->SetFixedImage( fixedImageReader->GetOutput() ) ); + EXPECT_NO_THROW( elastixFilter->SetMovingImage( movingImageReader->GetOutput() ) ); + EXPECT_NO_THROW( elastixFilter->SetParameterObject( eulerTransformParameterObject ) ); + + // TODO: This update should not be needed (see description in elxElastixFilterTest.cxx) + elastixFilter->Update(); + + TransformixFilterType::Pointer transformixFilter; + EXPECT_NO_THROW( transformixFilter = TransformixFilterType::New() ); + EXPECT_NO_THROW( transformixFilter->SetInputImage( movingImageReader->GetOutput() ) ); + EXPECT_NO_THROW( transformixFilter->SetTransformParameters( elastixFilter->GetTransformParameters() ) ); + EXPECT_NO_THROW( transformixFilter->ComputeSpatialJacobianOn() ); + EXPECT_NO_THROW( transformixFilter->ComputeDeterminantOfSpatialJacobianOn() ); + EXPECT_NO_THROW( transformixFilter->ComputeDeformationFieldOn() ); + // EXPECT_NO_THROW( transformixFilter->SetPointSetFileName( dataManager->GetInputFile( "InputMesh.pts" ) ) ); + EXPECT_NO_THROW( transformixFilter->SetOutputDirectory( dataManager->GetOutputDirectory() ) ); + + transformixFilter->Update(); + + ImageFileWriterType::Pointer writer = ImageFileWriterType::New(); + writer->SetFileName( dataManager->GetOutputFile( "Euler2DTransformixResultImage.nii" ) ); + writer->SetInput( transformixFilter->GetOutput() ); + writer->Update(); + EXPECT_NO_THROW( writer->Update() ); +}