Commit f67b3a3b authored by Floris Berendsen's avatar Floris Berendsen
Browse files

Merge branch 'ELASTIX-147-Create-SYNImageRegistrationMethod-Component' into develop

parents 8d07d8a3 f7aa64c8
......@@ -92,6 +92,7 @@ selxmodule_enable( ModuleExamples )
selxmodule_enable( ModuleSinksAndSources )
selxmodule_enable( ModuleItkSmoothingRecursiveGaussianImageFilter )
selxmodule_enable( ModuleItkImageRegistrationMethodv4 )
selxmodule_enable( ModuleItkSyNImageRegistrationMethod )
selxmodule_enable( ModuleElastix )
selxmodule_enable( ModuleController )
#selxmodule_enable( ModuleCommandLine )
......
......@@ -39,14 +39,11 @@ namespace po = boost::program_options;
#include <string>
#include <stdexcept>
using namespace std;
using namespace selx;
template< class T >
ostream &
operator<<( ostream & os, const vector< T > & v )
std::ostream &
operator<<( std::ostream & os, const std::vector< T > & v )
{
copy( v.begin(), v.end(), ostream_iterator< T >( os, " " ) );
std::copy( v.begin(), v.end(), std::ostream_iterator< T >( os, " " ) );
return os;
}
......@@ -56,26 +53,26 @@ main( int ac, char * av[] )
{
try
{
typedef vector< string > split_vector_type;
typedef std::vector< std::string > VectorOfStringsType;
SuperElastixFilter< DefaultComponents >::Pointer superElastixFilter = SuperElastixFilter< DefaultComponents >::New();
selx::SuperElastixFilter< selx::DefaultComponents >::Pointer superElastixFilter = selx::SuperElastixFilter< selx::DefaultComponents >::New();
fs::path configurationPath;
vector< string > inputPairs;
vector< string > outputPairs;
fs::path configurationPath;
VectorOfStringsType inputPairs;
VectorOfStringsType outputPairs;
// Store the reader so that they will not be destroyed before the pipeline is executed.
vector< AnyFileReader::Pointer > fileReaders;
std::vector< selx::AnyFileReader::Pointer > fileReaders;
// Store the writers for the update call
//vector<ImageWriter2DType::Pointer> fileWriters;
vector< AnyFileWriter::Pointer > fileWriters;
std::vector< selx::AnyFileWriter::Pointer > fileWriters;
po::options_description desc( "Allowed options" );
desc.add_options()
( "help", "produce help message" )
( "conf", po::value< fs::path >( &configurationPath )->required(), "Configuration file" )
( "in", po::value< vector< string >>( &inputPairs )->multitoken(), "Input data: images, labels, meshes, etc. Usage <name>=<path>" )
( "out", po::value< vector< string >>( &outputPairs )->multitoken(), "Output data: images, labels, meshes, etc. Usage <name>=<path>" )
( "in", po::value< VectorOfStringsType >( &inputPairs )->multitoken(), "Input data: images, labels, meshes, etc. Usage <name>=<path>" )
( "out", po::value< VectorOfStringsType >( &outputPairs )->multitoken(), "Output data: images, labels, meshes, etc. Usage <name>=<path>" )
( "graphout", po::value< fs::path >(), "Output Graphviz dot file" )
;
......@@ -85,24 +82,24 @@ main( int ac, char * av[] )
if( vm.count( "help" ) )
{
cout << desc << "\n";
std::cout << desc << "\n";
return 0;
}
Blueprint::Pointer blueprint;
selx::Blueprint::Pointer blueprint;
if( configurationPath.extension() == ".xml" )
{
// TODO: open file here and pass a stream to the ConfigurationReader
blueprint = ConfigurationReader::FromXML( configurationPath.string() );
blueprint = selx::ConfigurationReader::FromXML( configurationPath.string() );
}
else if( configurationPath.extension() == ".json" )
{
// TODO: open file here and pass a stream to the ConfigurationReader
blueprint = ConfigurationReader::FromJson( configurationPath.string() );
blueprint = selx::ConfigurationReader::FromJson( configurationPath.string() );
}
else
{
throw invalid_argument( "Configuration file requires extension .xml or .json" );
throw std::invalid_argument( "Configuration file requires extension .xml or .json" );
}
if( vm.count( "graphout" ) )
......@@ -114,22 +111,22 @@ main( int ac, char * av[] )
if( vm.count( "in" ) )
{
cout << "Number of input data: " << inputPairs.size() << "\n";
std::cout << "Number of input data: " << inputPairs.size() << "\n";
int index = 0;
for( const auto & inputPair : inputPairs )
{
split_vector_type nameAndPath;
VectorOfStringsType nameAndPath;
boost::split( nameAndPath, inputPair, boost::is_any_of( "=" ) ); // NameAndPath == { "name","path" }
const string & name = nameAndPath[ 0 ];
const string & path = nameAndPath[ 1 ];
const std::string & name = nameAndPath[ 0 ];
const std::string & path = nameAndPath[ 1 ];
cout << " " << index << " " << name << " : " << path << "\n";
std::cout << " " << index << " " << name << " : " << path << "\n";
++index;
// since we do not know which reader type we should instantiate for input "name",
// we ask SuperElastix for a reader that matches the type of the source component "name"
AnyFileReader::Pointer reader = superElastixFilter->GetInputFileReader( name );
selx::AnyFileReader::Pointer reader = superElastixFilter->GetInputFileReader( name );
reader->SetFileName( path );
superElastixFilter->SetInput( name, reader->GetOutput() );
fileReaders.push_back( reader );
......@@ -138,21 +135,21 @@ main( int ac, char * av[] )
if( vm.count( "out" ) )
{
cout << "Number of output data: " << outputPairs.size() << "\n";
std::cout << "Number of output data: " << outputPairs.size() << "\n";
int index = 0;
for( const auto & outputPair : outputPairs )
{
split_vector_type nameAndPath;
VectorOfStringsType nameAndPath;
boost::split( nameAndPath, outputPair, boost::is_any_of( "=" ) ); // NameAndPath == { "name","path" }
const string & name = nameAndPath[ 0 ];
const string & path = nameAndPath[ 1 ];
const std::string & name = nameAndPath[ 0 ];
const std::string & path = nameAndPath[ 1 ];
cout << " " << index << " " << name << " : " << path << "\n";
std::cout << " " << index << " " << name << " : " << path << "\n";
++index;
// since we do not know which writer type we should instantiate for output "name",
// we ask SuperElastix for a writer that matches the type of the sink component "name"
AnyFileWriter::Pointer writer = superElastixFilter->GetOutputFileWriter( name );
selx::AnyFileWriter::Pointer writer = superElastixFilter->GetOutputFileWriter( name );
//ImageWriter2DType::Pointer writer = ImageWriter2DType::New();
writer->SetFileName( path );
......@@ -169,14 +166,14 @@ main( int ac, char * av[] )
writer->Update();
}
}
catch( exception & e )
catch( std::exception & e )
{
cerr << "error: " << e.what() << "\n";
std::cerr << "error: " << e.what() << "\n";
return 1;
}
catch( ... )
{
cerr << "Exception of unknown type!\n";
std::cerr << "Exception of unknown type!\n";
}
return 0;
......
......@@ -40,7 +40,7 @@ set( ${MODULE}_TESTS
# Module source files
set( ${MODULE}_SOURCE_FILES
${${MODULE}_SOURCE_DIR}/src/selxItkImageRegistrationMethodv4.cxx
${${MODULE}_SOURCE_DIR}/src/selxItkImageRegistrationMethodv4Component.cxx
)
# Compile library
......
......@@ -389,14 +389,14 @@ TEST_F( WBIRDemoTest, elastix_BS_NCC )
blueprint->AddConnection( "ResultImageSink", "Controller", { {} } ); // { { "NameOfInterface", { "AfterRegistrationInterface" } } } ;
blueprint->WriteBlueprint( "elastix_BS_NCC.dot" );
// Data manager provides the paths to the input and output data for unit tests
DataManagerType::Pointer dataManager = DataManagerType::New();
blueprint->WriteBlueprint(dataManager->GetOutputFile("elastix_BS_NCC.dot"));
// Instantiate SuperElastix
EXPECT_NO_THROW( superElastixFilter = SuperElastixFilterType::New() );
// Data manager provides the paths to the input and output data for unit tests
DataManagerType::Pointer dataManager = DataManagerType::New();
// Set up the readers and writers
ImageReader2DType::Pointer fixedImageReader = ImageReader2DType::New();
fixedImageReader->SetFileName( dataManager->GetInputFile( "coneA2d64.mhd" ) );
......@@ -483,14 +483,14 @@ TEST_F( WBIRDemoTest, elastix_BS_MSD )
blueprint->AddConnection( "ResultImageSink", "Controller", { {} } ); // { { "NameOfInterface", { "AfterRegistrationInterface" } } } ;
blueprint->WriteBlueprint( "elastix_BS_MSD.dot" );
// Data manager provides the paths to the input and output data for unit tests
DataManagerType::Pointer dataManager = DataManagerType::New();
blueprint->WriteBlueprint(dataManager->GetOutputFile("elastix_BS_MSD.dot"));
// Instantiate SuperElastix
EXPECT_NO_THROW( superElastixFilter = SuperElastixFilterType::New() );
// Data manager provides the paths to the input and output data for unit tests
DataManagerType::Pointer dataManager = DataManagerType::New();
// Set up the readers and writers
ImageReader2DType::Pointer fixedImageReader = ImageReader2DType::New();
fixedImageReader->SetFileName( dataManager->GetInputFile( "coneA2d64.mhd" ) );
......
#=========================================================================
#
# Copyright Leiden University Medical Center, Erasmus University Medical
# Center and contributors
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0.txt
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
#=========================================================================
set( MODULE ModuleItkSyNImageRegistrationMethod )
# Export include files
set( ${MODULE}_INCLUDE_DIRS
${${MODULE}_SOURCE_DIR}/include
)
# Collect header files for Visual Studio Project
file(GLOB ${MODULE}_HEADER_FILES "${${MODULE}_SOURCE_DIR}/include/*.*")
# Export libraries
set( ${MODULE}_LIBRARIES
${MODULE}
)
# Export tests
set( ${MODULE}_TESTS
${${MODULE}_SOURCE_DIR}/test/selxSyNRegistrationItkv4Test.cxx
)
# Module source files
set( ${MODULE}_SOURCE_FILES
${${MODULE}_SOURCE_DIR}/src/selxItkSyNImageRegistrationMethodComponent.cxx
)
# Compile library
add_library( ${MODULE} STATIC ${${MODULE}_SOURCE_FILES} ${${MODULE}_HEADER_FILES})
target_link_libraries( ${MODULE} ${SUPERELASTIX_LIBRARIES} )
/*=========================================================================
*
* Copyright Leiden University Medical Center, Erasmus University Medical
* Center and contributors
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
#ifndef selxItkSyNImageRegistrationMethodComponent_h
#define selxItkSyNImageRegistrationMethodComponent_h
#include "selxSuperElastixComponent.h"
#include "selxInterfaces.h"
#include "itkSyNImageRegistrationMethod.h"
#include "itkGradientDescentOptimizerv4.h"
#include "itkImageSource.h"
#include <itkTransformToDisplacementFieldFilter.h>
#include <string.h>
#include "selxMacro.h"
#include "itkComposeDisplacementFieldsImageFilter.h"
#include "itkGaussianExponentialDiffeomorphicTransform.h"
#include "itkGaussianExponentialDiffeomorphicTransformParametersAdaptor.h"
namespace selx
{
template< int Dimensionality, class TPixel >
class ItkSyNImageRegistrationMethodComponent :
public SuperElastixComponent<
Accepting< itkImageFixedInterface< Dimensionality, TPixel >,
itkImageMovingInterface< Dimensionality, TPixel >,
itkMetricv4Interface< Dimensionality, TPixel >
>,
Providing< itkTransformInterface< double, Dimensionality >,
RunRegistrationInterface
>
>
{
public:
selxNewMacro( ItkSyNImageRegistrationMethodComponent, ComponentBase );
//itkStaticConstMacro(Dimensionality, unsigned int, Dimensionality);
ItkSyNImageRegistrationMethodComponent();
virtual ~ItkSyNImageRegistrationMethodComponent();
typedef TPixel PixelType;
// Get the type definitions from the interfaces
typedef typename itkImageFixedInterface< Dimensionality, TPixel >::ItkImageType FixedImageType;
typedef typename itkImageMovingInterface< Dimensionality, TPixel >::ItkImageType MovingImageType;
typedef typename itkTransformInterface< double, Dimensionality >::TransformType TransformType;
typedef typename itkTransformInterface< double, Dimensionality >::TransformPointer TransformPointer;
typedef typename itkTransformInterface< double, Dimensionality >::InternalComputationValueType TransformInternalComputationValueType; //should be from class template
typedef itk::SyNImageRegistrationMethod< FixedImageType, MovingImageType > TheItkFilterType;
typedef typename TheItkFilterType::ImageMetricType ImageMetricType;
typedef itk::RegistrationParameterScalesFromPhysicalShift< ImageMetricType > ScalesEstimatorType;
//Accepting Interfaces:
virtual int Set( itkImageFixedInterface< Dimensionality, TPixel > * ) override;
virtual int Set( itkImageMovingInterface< Dimensionality, TPixel > * ) override;
virtual int Set( itkMetricv4Interface< Dimensionality, TPixel > * ) override;
//Providing Interfaces:
virtual TransformPointer GetItkTransform() override;
virtual void RunRegistration() override;
//BaseClass methods
virtual bool MeetsCriterion( const ComponentBase::CriterionType & criterion ) override;
//static const char * GetName() { return "ItkSyNImageRegistrationMethod"; } ;
static const char * GetDescription() { return "ItkSyNImageRegistrationMethod Component"; }
private:
typename TheItkFilterType::Pointer m_theItkFilter;
protected:
/* The following struct returns the string name of computation type */
/* default implementation */
static inline const std::string GetTypeNameString()
{
itkGenericExceptionMacro( << "Unknown ScalarType" << typeid( TPixel ).name() );
// TODO: provide the user instructions how to enable the compilation of the component with the required template types (if desired)
// We might define an exception object that can communicate various error messages: for simple user, for developer user, etc
}
static inline const std::string GetPixelTypeNameString()
{
itkGenericExceptionMacro( << "Unknown PixelType" << typeid( TPixel ).name() );
// TODO: provide the user instructions how to enable the compilation of the component with the required template types (if desired)
// We might define an exception object that can communicate various error messages: for simple user, for developer user, etc
}
};
template< >
inline const std::string
ItkSyNImageRegistrationMethodComponent< 2, float >
::GetPixelTypeNameString()
{
return std::string( "float" );
}
template< >
inline const std::string
ItkSyNImageRegistrationMethodComponent< 2, double >
::GetPixelTypeNameString()
{
return std::string( "double" );
}
template< >
inline const std::string
ItkSyNImageRegistrationMethodComponent< 3, float >
::GetPixelTypeNameString()
{
return std::string( "float" );
}
template< >
inline const std::string
ItkSyNImageRegistrationMethodComponent< 3, double >
::GetPixelTypeNameString()
{
return std::string( "double" );
}
template< >
inline const std::string
ItkSyNImageRegistrationMethodComponent< 2, float >
::GetTypeNameString()
{
return std::string( "2_float" );
}
template< >
inline const std::string
ItkSyNImageRegistrationMethodComponent< 2, double >
::GetTypeNameString()
{
return std::string( "2_double" );
}
template< >
inline const std::string
ItkSyNImageRegistrationMethodComponent< 3, float >
::GetTypeNameString()
{
return std::string( "3_float" );
}
template< >
inline const std::string
ItkSyNImageRegistrationMethodComponent< 3, double >
::GetTypeNameString()
{
return std::string( "3_double" );
}
} //end namespace selx
#ifndef ITK_MANUAL_INSTANTIATION
#include "selxItkSyNImageRegistrationMethodComponent.hxx"
#endif
#endif // #define selxItkSyNImageRegistrationMethodComponent_h
/*=========================================================================
*
* Copyright Leiden University Medical Center, Erasmus University Medical
* Center and contributors
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
#include "selxItkSyNImageRegistrationMethodComponent.h"
#include "selxItkImageRegistrationMethodv4Component.h"
#include "itkDisplacementFieldTransformParametersAdaptor.h"
//TODO: get rid of these
#include "itkGradientDescentOptimizerv4.h"
namespace selx
{
template< int Dimensionality, class TPixel >
ItkSyNImageRegistrationMethodComponent< Dimensionality, TPixel >::ItkSyNImageRegistrationMethodComponent()
{
m_theItkFilter = TheItkFilterType::New();
m_theItkFilter->InPlaceOn();
//TODO: instantiating the filter in the constructor might be heavy for the use in component selector factory, since all components of the database are created during the selection process.
// we could choose to keep the component light weighted (for checking criteria such as names and connections) until the settings are passed to the filter, but this requires an additional initialization step.
}
template< int Dimensionality, class TPixel >
ItkSyNImageRegistrationMethodComponent< Dimensionality, TPixel >::~ItkSyNImageRegistrationMethodComponent()
{
}
template< int Dimensionality, class TPixel >
int
ItkSyNImageRegistrationMethodComponent< Dimensionality, TPixel >
::Set( itkImageFixedInterface< Dimensionality, TPixel > * component )
{
auto fixedImage = component->GetItkImageFixed();
// connect the itk pipeline
this->m_theItkFilter->SetFixedImage( fixedImage );
return 0;
}
template< int Dimensionality, class TPixel >
int
ItkSyNImageRegistrationMethodComponent< Dimensionality, TPixel >
::Set( itkImageMovingInterface< Dimensionality, TPixel > * component )
{
auto movingImage = component->GetItkImageMoving();
// connect the itk pipeline
this->m_theItkFilter->SetMovingImage( movingImage );
return 0;
}
template< int Dimensionality, class TPixel >
int
ItkSyNImageRegistrationMethodComponent< Dimensionality, TPixel >::Set( itkMetricv4Interface< Dimensionality, TPixel > * component )
{
this->m_theItkFilter->SetMetric( component->GetItkMetricv4() );
return 0;
}
template< int Dimensionality, class TPixel >
void
ItkSyNImageRegistrationMethodComponent< Dimensionality, TPixel >::RunRegistration( void )
{
typename FixedImageType::ConstPointer fixedImage = this->m_theItkFilter->GetFixedImage();
typename MovingImageType::ConstPointer movingImage = this->m_theItkFilter->GetMovingImage();
// Scale estimator is not used in current implementation yet
typename ScalesEstimatorType::Pointer scalesEstimator = ScalesEstimatorType::New();
ImageMetricType * theMetric = dynamic_cast< ImageMetricType * >( this->m_theItkFilter->GetModifiableMetric() );
auto optimizer = dynamic_cast< itk::GradientDescentOptimizerv4 * >( this->m_theItkFilter->GetModifiableOptimizer() );
//auto optimizer = dynamic_cast<itk::ObjectToObjectOptimizerBaseTemplate< InternalComputationValueType > *>(this->m_theItkFilter->GetModifiableOptimizer());
typedef itk::Vector< TransformInternalComputationValueType, Dimensionality > VectorType;
VectorType zeroVector( 0.0 );
typedef itk::Image< VectorType, Dimensionality > DisplacementFieldType;
typename DisplacementFieldType::Pointer displacementField = DisplacementFieldType::New();
displacementField->CopyInformation( fixedImage );
displacementField->SetRegions( fixedImage->GetBufferedRegion() );
displacementField->Allocate();
displacementField->FillBuffer( zeroVector );
typename DisplacementFieldType::Pointer inverseDisplacementField = DisplacementFieldType::New();
inverseDisplacementField->CopyInformation( fixedImage );
inverseDisplacementField->SetRegions( fixedImage->GetBufferedRegion() );
inverseDisplacementField->Allocate();
inverseDisplacementField->FillBuffer( zeroVector );
typedef typename TheItkFilterType::OutputTransformType OutputTransformType;
typename OutputTransformType::Pointer outputTransform = OutputTransformType::New();
outputTransform->SetDisplacementField( displacementField );
outputTransform->SetInverseDisplacementField( inverseDisplacementField );