diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/mitkPyramidImageRegistrationMethod.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/mitkPyramidImageRegistrationMethod.h new file mode 100644 index 0000000000..fa01f85077 --- /dev/null +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/mitkPyramidImageRegistrationMethod.h @@ -0,0 +1,375 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#ifndef MITKPYRAMIDIMAGEREGISTRATION_H +#define MITKPYRAMIDIMAGEREGISTRATION_H + +#include + +#include + +#include "mitkImage.h" +#include "mitkBaseProcess.h" + +#include "mitkPyramidRegistrationMethodHelper.h" + +#include "mitkImageToItk.h" +#include "mitkITKImageImport.h" + + +namespace mitk +{ + +/** + * @brief The PyramidImageRegistration class implements a multi-scale registration method + * + * The PyramidImageRegistration class is suitable for aligning (f.e.) brain MR images. The method offers two + * transform types + * - Rigid: optimizing translation and rotation only and + * - Affine ( default ): with scaling in addition ( 12 DOF ) + * + * The error metric is internally chosen based on the selected task type ( @sa SetCrossModalityOn ) + * It uses + * - MattesMutualInformation for CrossModality=on ( default ) and + * - NormalizedCorrelation for CrossModality=off. + */ +class DiffusionCore_EXPORT PyramidImageRegistrationMethod : + public itk::Object +{ +public: + + /** Typedefs */ + mitkClassMacro(PyramidImageRegistrationMethod, itk::Object) + + /** Smart pointer support */ + itkNewMacro(Self) + + /** Registration is between modalities - will take MattesMutualInformation as error metric */ + void SetCrossModalityOn() + { + m_CrossModalityRegistration = true; + } + + /** Registration is between modalities - will take NormalizedCorrelation as error metric */ + void SetCrossModalityOff() + { + m_CrossModalityRegistration = false; + } + + /** Turn the cross-modality on/off */ + void SetCrossModality(bool flag) + { + if( flag ) + this->SetCrossModalityOn(); + else + this->SetCrossModalityOff(); + } + + /** A rigid ( 6dof : translation + rotation ) transform is optimized */ + void SetTransformToRigid() + { + m_UseAffineTransform = false; + } + + /** An affine ( 12dof : Rigid + scale + skew ) transform is optimized */ + void SetTransformToAffine() + { + m_UseAffineTransform = true; + } + + /** Input image, the reference one */ + void SetFixedImage( mitk::Image::Pointer ); + + /** Input image, the one to be transformed */ + void SetMovingImage( mitk::Image::Pointer ); + + void Update(); + + /** + * @brief Get the number of parameters optimized ( 12 or 6 ) + * + * @return number of paramters + */ + unsigned int GetNumberOfParameters() + { + unsigned int retValue = 12; + if(!m_UseAffineTransform) + retValue = 6; + + return retValue; + } + + /** + * @brief Copies the estimated parameters to the given array + * @param paramArray, target array for copy, make sure to allocate enough space + */ + void GetParameters( double* paramArray) + { + if( m_EstimatedParameters == NULL ) + { + mitkThrow() << "No parameters were estimated yet, call Update() first."; + } + + unsigned int dim = 12; + if( !m_UseAffineTransform ) + dim = 6; + + for( unsigned int i=0; i + void RegisterTwoImages(itk::Image* itkImage1, itk::Image* itkImage2) + { + typedef typename itk::Image< TPixel1, VImageDimension1> FixedImageType; + typedef typename itk::Image< TPixel2, VImageDimension2> MovingImageType; + typedef typename itk::MultiResolutionImageRegistrationMethod< FixedImageType, MovingImageType > RegistrationType; + typedef itk::RegularStepGradientDescentOptimizer OptimizerType; + + typedef itk::AffineTransform< double > AffineTransformType; + typedef itk::Euler3DTransform< double > RigidTransformType; + typedef itk::MatrixOffsetTransformBase< double, VImageDimension1, VImageDimension2 > BaseTransformType; + + typedef typename itk::MattesMutualInformationImageToImageMetric< FixedImageType, MovingImageType > MMIMetricType; + typedef typename itk::NormalizedCorrelationImageToImageMetric< FixedImageType, MovingImageType > NCMetricType; + typedef typename itk::ImageToImageMetric< FixedImageType, MovingImageType> BaseMetricType; + + typename itk::LinearInterpolateImageFunction::Pointer interpolator = + itk::LinearInterpolateImageFunction::New(); + + + typename BaseMetricType::Pointer metric; + + if( m_CrossModalityRegistration ) + { + metric = MMIMetricType::New(); + } + else + { + metric = NCMetricType::New(); + } + + + + // initial parameter dimension ( affine = default ) + unsigned int paramDim = 12; + typename BaseTransformType::Pointer transform; + if( m_UseAffineTransform ) + { + transform = AffineTransformType::New(); + } + else + { + transform = RigidTransformType::New(); + paramDim = 6; + } + + typename BaseTransformType::ParametersType initialParams(paramDim); + + initialParams.Fill(0); + if( m_UseAffineTransform ) + { + initialParams[0] = initialParams[4] = initialParams[8] = 1; + } + + + typename FixedImageType::Pointer referenceImage = itkImage1; + typename MovingImageType::Pointer movingImage = itkImage2; + typename FixedImageType::RegionType maskedRegion = referenceImage->GetLargestPossibleRegion(); + + typename RegistrationType::Pointer registration = RegistrationType::New(); + unsigned int max_pyramid_lvl = 3; + unsigned int max_schedule_val = 4; + typename FixedImageType::RegionType::SizeType image_size = referenceImage->GetLargestPossibleRegion().GetSize(); + unsigned int min_value = std::min( image_size[0], std::min( image_size[1], image_size[2])); + + // condition for the top level pyramid image + float optmaxstep = 12; + float optminstep = 0.1f; + unsigned int iterations = 40; + if( min_value / max_schedule_val < 16 ) + { + //max_pyramid_lvl--; + max_schedule_val /= 2; + optmaxstep *= 0.5f; + optminstep *= 0.2f; + iterations *= 1.5; + + //std::cout << "Changed default pyramid: lvl " << max_pyramid_lvl << std::endl; + } + + typename RegistrationType::ScheduleType fixedSchedule(max_pyramid_lvl,3); + fixedSchedule.Fill(1); + + fixedSchedule[0][0] = max_schedule_val; + fixedSchedule[0][1] = max_schedule_val; + fixedSchedule[0][2] = max_schedule_val; + + for( unsigned int i=1; iSetScales( optScales ); + optimizer->SetInitialPosition( initialParams ); + optimizer->SetNumberOfIterations( iterations ); + optimizer->SetGradientMagnitudeTolerance( 1e-4 ); + optimizer->SetRelaxationFactor( 0.7 ); + optimizer->SetMaximumStepLength( optmaxstep ); + optimizer->SetMinimumStepLength( optminstep ); + + OptimizerIterationCommand< OptimizerType >::Pointer iterationObserver = + OptimizerIterationCommand::New(); + + unsigned long vopt_tag = optimizer->AddObserver( itk::IterationEvent(), iterationObserver ); + + + // INPUT IMAGES + registration->SetFixedImage( referenceImage ); + registration->SetFixedImageRegion( maskedRegion ); + registration->SetMovingImage( movingImage ); + registration->SetSchedules( fixedSchedule, fixedSchedule); + // OTHER INPUTS + registration->SetMetric( metric ); + registration->SetOptimizer( optimizer ); + registration->SetTransform( transform.GetPointer() ); + registration->SetInterpolator( interpolator ); + registration->SetInitialTransformParameters( initialParams ); + + typename PyramidOptControlCommand::Pointer pyramid_observer = + PyramidOptControlCommand::New(); + + registration->AddObserver( itk::IterationEvent(), pyramid_observer ); + + try + { + registration->Update(); + } + catch (itk::ExceptionObject &e) + { + MITK_ERROR << "[Registration Update] Caught ITK exception: "; + mitkThrow() << "Registration failed with exception: " << e.what(); + } + + if( m_EstimatedParameters != NULL) + { + delete [] m_EstimatedParameters; + } + + m_EstimatedParameters = new double[paramDim]; + + typename BaseTransformType::ParametersType finalParameters = registration->GetLastTransformParameters(); + for( unsigned int i=0; i + void ResampleMitkImage( itk::Image* itkImage, + mitk::Image::Pointer& outputImage ) + { + typedef typename itk::Image< TPixel, VDimension> ImageType; + typedef typename itk::ResampleImageFilter< ImageType, ImageType, double> ResampleImageFilterType; + + typedef itk::LinearInterpolateImageFunction< ImageType, double > InterpolatorType; + typename InterpolatorType::Pointer interpolator = InterpolatorType::New(); + + typename mitk::ImageToItk< ImageType >::Pointer reference_image = mitk::ImageToItk< ImageType >::New(); + reference_image->SetInput( this->m_FixedImage ); + reference_image->Update(); + + typedef itk::MatrixOffsetTransformBase< double, 3, 3> BaseTransformType; + BaseTransformType::Pointer base_transform = BaseTransformType::New(); + + if( this->m_UseAffineTransform ) + { + typedef itk::AffineTransform< double > TransformType; + TransformType::Pointer transform = TransformType::New(); + + TransformType::ParametersType affine_params( TransformType::ParametersDimension ); + this->GetParameters( &affine_params[0] ); + + transform->SetParameters( affine_params ); + + base_transform = transform; + } + // Rigid + else + { + typedef itk::Euler3DTransform< double > RigidTransformType; + RigidTransformType::Pointer rtransform = RigidTransformType::New(); + + RigidTransformType::ParametersType rigid_params( RigidTransformType::ParametersDimension ); + this->GetParameters( &rigid_params[0] ); + + rtransform->SetParameters( rigid_params ); + + base_transform = rtransform; + } + + typename ResampleImageFilterType::Pointer resampler = ResampleImageFilterType::New(); + resampler->SetInput( itkImage ); + resampler->SetTransform( base_transform ); + resampler->SetReferenceImage( reference_image->GetOutput() ); + resampler->UseReferenceImageOn(); + + resampler->Update(); + + + mitk::GrabItkImageMemory( resampler->GetOutput(), outputImage); + + } + +}; + +} // end namespace + +#endif // MITKPYRAMIDIMAGEREGISTRATION_H diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/mitkPyramidRegistrationMethodHelper.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/mitkPyramidRegistrationMethodHelper.h new file mode 100644 index 0000000000..84c8374bcd --- /dev/null +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/mitkPyramidRegistrationMethodHelper.h @@ -0,0 +1,142 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#ifndef MITKPYRAMIDREGISTRATIONMETHODHELPER_H +#define MITKPYRAMIDREGISTRATIONMETHODHELPER_H + +#include + +#include + +#include +#include +#include + +#include +#include + +#include +#include + +#include "mitkImageAccessByItk.h" + +/** + * @brief Provides same functionality as \a AccessTwoImagesFixedDimensionByItk for a subset of types + * + * For now, the subset defined is only short and float. + */ +#define AccessTwoImagesFixedDimensionTypeSubsetByItk(mitkImage1, mitkImage2, itkImageTypeFunction, dimension) \ +{ \ + const mitk::PixelType& pixelType1 = mitkImage1->GetPixelType(); \ + const mitk::PixelType& pixelType2 = mitkImage2->GetPixelType(); \ + const mitk::Image* constImage1 = mitkImage1; \ + const mitk::Image* constImage2 = mitkImage2; \ + mitk::Image* nonConstImage1 = const_cast(constImage1); \ + mitk::Image* nonConstImage2 = const_cast(constImage2); \ + nonConstImage1->Update(); \ + nonConstImage2->Update(); \ + _checkSpecificDimension(mitkImage1, (dimension)); \ + _checkSpecificDimension(mitkImage2, (dimension)); \ + _accessTwoImagesByItkForEach(itkImageTypeFunction, ((short, dimension))((float, dimension)), ((short, dimension))((float, dimension)) ) \ + { \ + std::string msg("Pixel type "); \ + msg.append(pixelType1.GetComponentTypeAsString() ); \ + msg.append(" or pixel type "); \ + msg.append(pixelType2.GetComponentTypeAsString() ); \ + msg.append(" is not in " MITK_PP_STRINGIZE(MITK_ACCESSBYITK_TYPES_DIMN_SEQ(dimension))); \ + throw mitk::AccessByItkException(msg); \ + } \ +} + + +/** + * @brief The PyramidOptControlCommand class stears the step lenght of the + * gradient descent optimizer in multi-scale registration methods + */ +template +class PyramidOptControlCommand : public itk::Command +{ +public: + + typedef itk::RegularStepGradientDescentOptimizer OptimizerType; + + itkNewMacro( PyramidOptControlCommand ) + + void Execute(itk::Object *caller, const itk::EventObject & /*event*/) + { + RegistrationType* registration = dynamic_cast< RegistrationType* >( caller ); + + std::cout << "\t - Pyramid level " << registration->GetCurrentLevel() << std::endl; + if( registration->GetCurrentLevel() == 0 ) + return; + + + + OptimizerType* optimizer = dynamic_cast< OptimizerType* >(registration->GetOptimizer()); + + std::cout << optimizer->GetStopConditionDescription() << std::endl; + + optimizer->SetMaximumStepLength( optimizer->GetMaximumStepLength() * 0.5f ); + optimizer->SetMinimumStepLength( optimizer->GetMinimumStepLength() * 0.1f ); + optimizer->SetNumberOfIterations( optimizer->GetNumberOfIterations() * 1.5f ); + } + + void Execute(const itk::Object * /*object*/, const itk::EventObject & /*event*/){} +}; + + +template +class OptimizerIterationCommand : public itk::Command +{ +public: + itkNewMacro( OptimizerIterationCommand ) + + void Execute(itk::Object *caller, const itk::EventObject & event) + { + OptimizerType* optimizer = dynamic_cast< OptimizerType* >( caller ); + + unsigned int currentIter = optimizer->GetCurrentIteration(); + std::cout << "[" << currentIter << "] : " << optimizer->GetValue() << " : " << optimizer->GetCurrentPosition() << std::endl; + +/* typedef typename OptimizerType::CostFunctionType CostFunctionType; + const CostFunctionType *cfcn = optimizer->GetCostFunction(); + + const std::string match("MattesMutualInformationImageToImageMetric"); + const std::string name( cfcn->GetNameOfClass() ); + + double compareValue = -0.1; + + if ( name.compare( match ) == 0 ) + { + compareValue = -0.2; + } + + if (optimizer->GetValue() > compareValue ) + { + optimizer->StopOptimization(); + itkGenericExceptionMacro( << " Unrealistic metric value reached. Stopping optimization! "); + } +*/ + } + + void Execute(const itk::Object * object, const itk::EventObject & event) + { + + } +}; + + +#endif // MITKPYRAMIDREGISTRATIONMETHODHELPER_H diff --git a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkPyramidImageRegistrationMethodTest.cpp b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkPyramidImageRegistrationMethodTest.cpp new file mode 100644 index 0000000000..6aa9dab239 --- /dev/null +++ b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkPyramidImageRegistrationMethodTest.cpp @@ -0,0 +1,156 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + + +#include "mitkTestingMacros.h" +#include "mitkIOUtil.h" + +#include "mitkPyramidImageRegistrationMethod.h" + +#include + + +int mitkPyramidImageRegistrationMethodTest( int argc, char* argv[] ) +{ + if( argc < 4 ) + { + MITK_ERROR << "Not enough input \n Usage: fixed moving type [output_image [output_transform]]" + << "\n \t fixed : the path to the fixed image \n" + << " \t moving : path to the image to be registered" + << " \t type : Affine or Rigid defining the type of the transformation" + << " \t output_image : output file optional, (full) path, and optionally output_transform : also (full)path to file"; + return EXIT_FAILURE; + } + + MITK_TEST_BEGIN("PyramidImageRegistrationMethodTest"); + + mitk::Image::Pointer fixedImage = mitk::IOUtil::LoadImage( argv[1] ); + mitk::Image::Pointer movingImage = mitk::IOUtil::LoadImage( argv[2] ); + + std::string type_flag( argv[3] ); + + mitk::PyramidImageRegistrationMethod::Pointer registrationMethod = mitk::PyramidImageRegistrationMethod::New(); + registrationMethod->SetFixedImage( fixedImage ); + registrationMethod->SetMovingImage( movingImage ); + + if( type_flag == "Rigid" ) + { + registrationMethod->SetTransformToRigid(); + } + else if( type_flag == "Affine" ) + { + registrationMethod->SetTransformToAffine(); + } + else + { + MITK_WARN << " No type specified, using 'Affine' ."; + } + + registrationMethod->Update(); + + bool imageOutput = false; + bool transformOutput = false; + + std::string image_out_filename, transform_out_filename; + + std::string first_output( argv[4] ); + // check for txt, otherwise suppose it is an image + if( first_output.find(".txt") != std::string::npos ) + { + transformOutput = true; + transform_out_filename = first_output; + } + else + { + imageOutput = true; + image_out_filename = first_output; + } + + if( argc > 4 ) + { + std::string second_output( argv[5] ); + if( second_output.find(".txt") != std::string::npos ) + { + transformOutput = true; + transform_out_filename = second_output; + } + } + + MITK_INFO << " Selected output: " << transform_out_filename << " " << image_out_filename; + + try{ + + unsigned int paramCount = registrationMethod->GetNumberOfParameters(); + double* params = new double[ paramCount ]; + registrationMethod->GetParameters( ¶ms[0] ); + + std::cout << "Parameters: "; + for( unsigned int i=0; i< paramCount; i++) + { + std::cout << params[ i ] << " "; + } + std::cout << std::endl; + + if( imageOutput ) + { + mitk::IOUtil::SaveImage( registrationMethod->GetResampledMovingImage(), image_out_filename.c_str() ); + } + + + if( transformOutput ) + { + + itk::TransformFileWriter::Pointer writer = itk::TransformFileWriter::New(); + + // Get transform parameter for resampling / saving + // Affine + if( paramCount == 12 ) + { + typedef itk::AffineTransform< double > TransformType; + TransformType::Pointer transform = TransformType::New(); + + TransformType::ParametersType affine_params( paramCount ); + registrationMethod->GetParameters( &affine_params[0] ); + + transform->SetParameters( affine_params ); + writer->SetInput( transform ); + } + // Rigid + else + { + typedef itk::Euler3DTransform< double > RigidTransformType; + RigidTransformType::Pointer rtransform = RigidTransformType::New(); + + RigidTransformType::ParametersType rigid_params( paramCount ); + registrationMethod->GetParameters( &rigid_params[0] ); + + rtransform->SetParameters( rigid_params ); + writer->SetInput( rtransform ); + } + + writer->SetFileName( transform_out_filename ); + writer->Update(); + } + + } + catch( const std::exception &e) + { + MITK_ERROR << "Caught exception: " << e.what(); + } + + + MITK_TEST_END(); +}