diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/mitkPyramidImageRegistrationMethod.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/mitkPyramidImageRegistrationMethod.cpp index 0be0933ca4..f3633a0050 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/mitkPyramidImageRegistrationMethod.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/mitkPyramidImageRegistrationMethod.cpp @@ -1,72 +1,83 @@ /*=================================================================== 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 "mitkPyramidImageRegistrationMethod.h" #include "mitkException.h" #include "mitkImageAccessByItk.h" mitk::PyramidImageRegistrationMethod::PyramidImageRegistrationMethod() : m_FixedImage(NULL), m_MovingImage(NULL), m_CrossModalityRegistration(true), - m_UseAffineTransform(true) + m_UseAffineTransform(true), + m_EstimatedParameters(NULL) { } mitk::PyramidImageRegistrationMethod::~PyramidImageRegistrationMethod() { - + if( m_EstimatedParameters != NULL) + { + delete [] m_EstimatedParameters; + } } void mitk::PyramidImageRegistrationMethod::SetFixedImage(mitk::Image::Pointer fixed) { if( fixed.IsNotNull() ) { m_FixedImage = fixed; } } void mitk::PyramidImageRegistrationMethod::SetMovingImage(mitk::Image::Pointer moving) { if( moving.IsNotNull() ) { m_MovingImage = moving; } } void mitk::PyramidImageRegistrationMethod::Update() { if( m_MovingImage.IsNull() ) { mitkThrow() << " Moving image is null"; } if( m_FixedImage.IsNull() ) { mitkThrow() << " Moving image is null"; } - if( m_FixedImage->GetDimension() != 3 || - m_MovingImage->GetDimension() != 3 ) + unsigned int allowedDimension = 3; + + if( m_FixedImage->GetDimension() != allowedDimension || + m_MovingImage->GetDimension() != allowedDimension ) { mitkThrow() << " Only 3D Images supported."; } - AccessTwoImagesFixedDimensionByItk( m_FixedImage, m_MovingImage, RegisterTwoImages, 3); + // + // One possibility: use the FixedDimesnionByItk, but this instantiates ALL possible + // pixel type combinations! + // AccessTwoImagesFixedDimensionByItk( m_FixedImage, m_MovingImage, RegisterTwoImages, 3); + // in helper: TypeSubset : short, float + AccessTwoImagesFixedDimensionTypeSubsetByItk( m_FixedImage, m_MovingImage, RegisterTwoImages, 3); } diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/mitkPyramidImageRegistrationMethod.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/mitkPyramidImageRegistrationMethod.h index 6408431ab0..e52a70de27 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/mitkPyramidImageRegistrationMethod.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/mitkPyramidImageRegistrationMethod.h @@ -1,244 +1,328 @@ /*=================================================================== 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" 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(); + +/* + * moment initialization not default settings, to be used as soon as initialization settings available + typedef typename itk::ImageMomentsCalculator< FixedImageType > FImageMoments; + typedef typename itk::ImageMomentsCalculator< MovingImageType > MImageMoments; + + typename FImageMoments::Pointer f_moments = FImageMoments::New(); + f_moments->SetImage( itkImage1 ); + f_moments->Compute(); + + typename MImageMoments::Pointer m_moments = MImageMoments::New(); + m_moments->SetImage( itkImage2 ); + m_moments->Compute(); + + itk::Vector< double, 3> f_cog = f_moments->GetCenterOfGravity(); + itk::Vector< double, 3> m_cog = m_moments->GetCenterOfGravity(); +*/ + 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::ParametersType initialParams(paramDim); - initialParams.Fill(0); - 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; + } + + /* + * moment initialization not default settings, to be used as soon as initialization settings available + initialParams[paramDim-3] = m_cog[2] - c_cog[2]; + initialParams[paramDim-2] = m_cog[1] - c_cog[1]; + initialParams[paramDim-1] = m_cog[0] - c_cog[0]; +*/ 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 = 4; unsigned int max_schedule_val = 8; 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 = 16; + float optmaxstep = 12; float optminstep = 0.1f; if( min_value / max_schedule_val < 8 ) { //max_pyramid_lvl--; max_schedule_val /= 2; optmaxstep *= 0.25f; optminstep *= 0.1f; //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( 100 ); optimizer->SetGradientMagnitudeTolerance( 1e-4 ); optimizer->SetRelaxationFactor( 0.7 ); optimizer->SetMaximumStepLength( optmaxstep ); optimizer->SetMinimumStepLength( optminstep ); // INPUT IMAGES registration->SetFixedImage( referenceImage ); - //registration->SetFixedImageRegion( maskedRegion ); + 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: "; - MITK_ERROR << e; + 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