diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkBatchedRegistration.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkBatchedRegistration.cpp index 878a0b9a1b..76b2fa7d55 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkBatchedRegistration.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkBatchedRegistration.cpp @@ -1,183 +1,178 @@ #include "mitkBatchedRegistration.h" #include "mitkPyramidImageRegistrationMethod.h" #include "mitkDiffusionImage.h" #include #include #include "itkB0ImageExtractionImageFilter.h" #include #include // VTK #include // DEBUG #include mitk::BatchedRegistration::BatchedRegistration() : m_RegisteredImagesValid(false) { } void mitk::BatchedRegistration::SetFixedImage(mitk::Image::Pointer& fixedImage) { m_FixedImage = fixedImage; } void mitk::BatchedRegistration::SetMovingReferenceImage(Image::Pointer &movingImage) { m_MovingReference = movingImage; m_RegisteredImagesValid = false; } void mitk::BatchedRegistration::SetBatch(std::vector imageBatch) { m_ImageBatch.clear(); m_ImageBatch = imageBatch; } std::vector mitk::BatchedRegistration::GetRegisteredImages() { if (!m_RegisteredImagesValid) { m_RegisteredImages.clear(); // First transform moving reference image - TransformType transf = GetTransformation(m_FixedImage, m_MovingReference); + RidgidTransformType transf = new double(6); + GetTransformation(m_FixedImage, m_MovingReference,transf); // store it as first element in vector ApplyTransformationToImage(m_MovingReference,transf); m_RegisteredImages.push_back(m_MovingReference); // apply transformation to whole batch std::vector::const_iterator itEnd = m_ImageBatch.end(); for (std::vector::iterator it = m_ImageBatch.begin(); it != itEnd; ++it) { ApplyTransformationToImage(*it,transf); m_RegisteredImages.push_back(*it); } } return m_RegisteredImages; } -void mitk::BatchedRegistration::ApplyTransformationToImage(mitk::Image::Pointer &img, const mitk::BatchedRegistration::TransformType &transformation) const +void mitk::BatchedRegistration::ApplyTransformationToImage(mitk::Image::Pointer &img, const mitk::BatchedRegistration::RidgidTransformType &transformation, mitk::Image::Pointer resampleReference, bool binary ) const { typedef mitk::DiffusionImage DiffusionImageType; - vtkMatrix4x4* transformationMatrix = vtkMatrix4x4::New(); - vnl_matrix_fixed vnlRotationMatrix; - for (int i = 0; i < 4;++i) - for (int j = 0; j < 4;++j) - { - transformationMatrix->Element[i][j] = transformation.get(i,j); - if (i < 3 && j < 3) - vnlRotationMatrix.set(i,j,transformation.get(i,j)); - } + mitk::Image::Pointer ref; + mitk::PyramidImageRegistrationMethod::Pointer registrationMethod = mitk::PyramidImageRegistrationMethod::New(); + registrationMethod->SetTransformToRigid(); + if (binary) + registrationMethod->SetUseNearestNeighborInterpolation(true); + + if (resampleReference.IsNotNull()) + { + registrationMethod->SetFixedImage( resampleReference ); + } + else + { + // clone image, to prevent recursive access on resampling .. + ref = img->Clone(); + registrationMethod->SetFixedImage( ref ); + } + + img = registrationMethod->GetResampledMovingImage(img, transformation); + - img->GetGeometry()->Compose( transformationMatrix); if (dynamic_cast (img.GetPointer()) != NULL) { DiffusionImageType::Pointer diffImages = dynamic_cast(img.GetPointer()); // Changes to default geometry does not affect the presentation // of Diffusion Images. (No idea why) // Workaround: // Extracting a B0-image // Casting B0 to an MITK image, so we have an image with // the same dimensions and space // Apply the transformation to this image // Cast back to ITK image // Copy Transformation informations from B0-ITK-Image to // diffusion vector image. itk::B0ImageExtractionImageFilter::Pointer b0Extraction = itk::B0ImageExtractionImageFilter::New(); b0Extraction->SetInput(diffImages->GetVectorImage()); b0Extraction->SetDirections(diffImages->GetDirections()); b0Extraction->Update(); mitk::Image::Pointer tmp = mitk::Image::New(); tmp->InitializeByItk(b0Extraction->GetOutput()); tmp->SetVolume(b0Extraction->GetOutput()->GetBufferPointer()); - tmp->GetGeometry()->Compose(transformationMatrix); + // TODO FIXME!!! + //tmp->GetGeometry()->Compose(transformationMatrix); itk::Image::Pointer itkTmp; mitk::CastToItkImage >(tmp, itkTmp); diffImages->GetVectorImage()->SetOrigin(itkTmp->GetOrigin()); diffImages->GetVectorImage()->SetDirection(itkTmp->GetDirection()); diffImages->GetVectorImage()->SetSpacing(itkTmp->GetSpacing()); mitk::DiffusionImageCorrectionFilter::Pointer correctionFilter = mitk::DiffusionImageCorrectionFilter::New(); // For Diff. Images: Need to rotate the gradients correctionFilter->SetImage(diffImages); - correctionFilter->CorrectDirections(vnlRotationMatrix); + //correctionFilter->CorrectDirections(vnlRotationMatrix); } } -mitk::BatchedRegistration::TransformType mitk::BatchedRegistration::GetTransformation(mitk::Image::Pointer fixedImage, mitk::Image::Pointer movingImage, mitk::Image::Pointer mask) +void mitk::BatchedRegistration::GetTransformation(mitk::Image::Pointer fixedImage, mitk::Image::Pointer movingImage, RidgidTransformType transformation, mitk::Image::Pointer mask) { // Handle case that fixed or moving image is a DWI image mitk::DiffusionImage* fixedDwi = dynamic_cast*> (fixedImage.GetPointer()); mitk::DiffusionImage* movingDwi = dynamic_cast*> (movingImage.GetPointer()); itk::B0ImageExtractionImageFilter::Pointer b0Extraction = itk::B0ImageExtractionImageFilter::New(); if (fixedDwi != NULL) { b0Extraction->SetInput(fixedDwi->GetVectorImage()); b0Extraction->SetDirections(fixedDwi->GetDirections()); b0Extraction->Update(); mitk::Image::Pointer tmp = mitk::Image::New(); tmp->InitializeByItk(b0Extraction->GetOutput()); tmp->SetVolume(b0Extraction->GetOutput()->GetBufferPointer()); fixedImage = tmp; } if (movingDwi != NULL) { b0Extraction->SetInput(movingDwi->GetVectorImage()); b0Extraction->SetDirections(movingDwi->GetDirections()); b0Extraction->Update(); mitk::Image::Pointer tmp = mitk::Image::New(); tmp->InitializeByItk(b0Extraction->GetOutput()); tmp->SetVolume(b0Extraction->GetOutput()->GetBufferPointer()); movingImage = tmp; } // Start registration mitk::PyramidImageRegistrationMethod::Pointer registrationMethod = mitk::PyramidImageRegistrationMethod::New(); registrationMethod->SetFixedImage( fixedImage ); if (mask.IsNotNull()) { registrationMethod->SetFixedImageMask(mask); registrationMethod->SetUseFixedImageMask(true); } else { registrationMethod->SetUseFixedImageMask(false); } registrationMethod->SetTransformToRigid(); registrationMethod->SetCrossModalityOn(); registrationMethod->SetMovingImage(movingImage); registrationMethod->Update(); - TransformType transformation; - mitk::PyramidImageRegistrationMethod::TransformMatrixType rotationMatrix; - rotationMatrix = registrationMethod->GetLastRotationMatrix().transpose(); - double param[6]; - registrationMethod->GetParameters(param); // first three: euler angles, last three translation - for (unsigned int i = 0; i < 3; i++) - { - for (unsigned int j = 0; j < 3; ++j) - { - double value = rotationMatrix.get(i,j); - transformation.set(i,j, value); - } - } - transformation.set(0,3,-param[3]); - transformation.set(1,3,-param[4]); - transformation.set(2,3,-param[5]); - transformation.set(3,3,1); - return transformation; + registrationMethod->GetParameters(transformation); // first three: euler angles, last three translation + } diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkBatchedRegistration.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkBatchedRegistration.h index 8d92ea1d9a..7beefe1305 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkBatchedRegistration.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkBatchedRegistration.h @@ -1,79 +1,80 @@ #ifndef MITKBATCHEDREGISTRATION_H #define MITKBATCHEDREGISTRATION_H // ITK #include // MITK #include #include "mitkCommon.h" #include "mitkImage.h" namespace mitk { /** * @brief The BatchedRegistration class Wrapper to calculate and apply a reference transformation to several images. * * Use if several pictures with the same world geometry are to be registered * to one reference image, the registration is only computed once (for the moving image) and the geometry transformed for the complete * image batch accordingly. Can handle image types that are usually not supported by registrations filters, e.g. fiber bundles and segmentations: * these can be registered if a "registerable" image such as B0/T2 from which they are derived is supplied, since the transformation can be calculated * on those and applied to the derived objects. */ class DiffusionCore_EXPORT BatchedRegistration : public itk::LightObject { public: - typedef vnl_matrix_fixed TransformType; + typedef double* RidgidTransformType; mitkClassMacro(BatchedRegistration, itk::LightObject) itkNewMacro(Self) void SetFixedImage(Image::Pointer &fixedImage); void SetMovingReferenceImage(mitk::Image::Pointer& movingImage); void SetBatch(std::vector imageBatch); /** * @brief GetRegisteredImages returns registered images , * * at position 0 the registered moving reference image is supplied followed all registered images from the batch. */ std::vector GetRegisteredImages(); - void ApplyTransformationToImage(mitk::Image::Pointer& img, const TransformType& transformation) const; + void ApplyTransformationToImage(mitk::Image::Pointer& img, const RidgidTransformType& transformation, mitk::Image::Pointer resampleReference = NULL , bool binary = false) const; + + void GetTransformation(mitk::Image::Pointer fixedImage , mitk::Image::Pointer movingImage, RidgidTransformType transformation, mitk::Image::Pointer mask = NULL); - TransformType GetTransformation(mitk::Image::Pointer fixedImage, mitk::Image::Pointer movingImage, mitk::Image::Pointer mask = NULL); protected: BatchedRegistration(); ~BatchedRegistration(){}; private: BatchedRegistration(const Self &); //purposely not implemented void operator=(const Self &); //purposely not implemented bool m_RegisteredImagesValid; mitk::Image::Pointer m_FixedImage; mitk::Image::Pointer m_MovingReference; /** * @brief m_ImageBatch List of images on which that the reference transformation is applied * */ std::vector m_ImageBatch; /** * @brief m_RegisteredImages List of references to registered images. * */ std::vector m_RegisteredImages; }; } #endif // MITKBATCHEDREGISTRATION_H diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidImageRegistrationMethod.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidImageRegistrationMethod.cpp index fe8fa3cf12..418a5a1bd7 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidImageRegistrationMethod.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidImageRegistrationMethod.cpp @@ -1,146 +1,180 @@ /*=================================================================== 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_UseWindowedSincInterpolator(false), + m_UseNearestNeighborInterpolator(false), m_UseMask(false), m_EstimatedParameters(NULL), m_Verbose(false), m_InitializeByGeometry(false) { } 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::SetFixedImageMask(mitk::Image::Pointer mask) { m_FixedImageMask = mask; } void mitk::PyramidImageRegistrationMethod::Update() { if( m_MovingImage.IsNull() ) { mitkThrow() << " Moving image is null"; } if( m_FixedImage.IsNull() ) { mitkThrow() << " Moving image is null"; } unsigned int allowedDimension = 3; if( m_FixedImage->GetDimension() != allowedDimension || m_MovingImage->GetDimension() != allowedDimension ) { mitkThrow() << " Only 3D Images supported."; } // // 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); } mitk::PyramidImageRegistrationMethod::TransformMatrixType mitk::PyramidImageRegistrationMethod ::GetLastRotationMatrix() { TransformMatrixType output; if( m_EstimatedParameters == NULL ) { output.set_identity(); return output; } 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; } 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; } return base_transform->GetMatrix().GetVnlMatrix(); } mitk::Image::Pointer mitk::PyramidImageRegistrationMethod ::GetResampledMovingImage() { mitk::Image::Pointer output = mitk::Image::New(); //output->Initialize( this->m_FixedImage ); AccessFixedDimensionByItk_1( this->m_MovingImage, ResampleMitkImage, 3, output ); return output; } + +mitk::Image::Pointer mitk::PyramidImageRegistrationMethod::GetResampledMovingImage(mitk::Image::Pointer movingImage, double* transform) +{ + mitk::Image::Pointer output = mitk::Image::New(); + + + unsigned int dim = 12; + if( !m_UseAffineTransform ) + dim = 6; + + if (m_EstimatedParameters == NULL) + m_EstimatedParameters = new double[dim]; + + double tmpParams[12]; + // save and set temporal transformation values + for( unsigned int i=0; i #include #include "itkImageMaskSpatialObject.h" #include "itkNotImageFilter.h" #include #include #include "mitkImage.h" #include "mitkBaseProcess.h" #include "mitkPyramidRegistrationMethodHelper.h" #include #include "mitkImageToItk.h" #include "mitkImageCast.h" #include "mitkImageCaster.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) /** Typedef for the transformation matrix, corresponds to the InternalMatrixType from ITK transforms */ typedef vnl_matrix_fixed< double, 3, 3> TransformMatrixType; + typedef itk::AffineTransform< double > AffineTransformType; + typedef itk::Euler3DTransform< double > RigidTransformType; + /** 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(); } /** Controll the verbosity of the registration output * * for true, each iteration step of the optimizer is watched and passed to the std::cout */ void SetVerbose(bool flag) { if( flag ) this->SetVerboseOn(); else this->SetVerboseOff(); } /** Turn verbosity on, \sa SetVerbose */ void SetVerboseOn() { m_Verbose = true; } /** Turn verbosity off, \sa SetVerbose */ void SetVerboseOff() { m_Verbose = false; } /** 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 ); - /** Fixed image mask, excludes the masked voxels from the registration metric*/ void SetFixedImageMask( mitk::Image::Pointer mask); void SetInitializeByGeometry(bool flag) { m_InitializeByGeometry = flag; } 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 * * It returns identity if the internal parameters are not-yet allocated ( i.e. the filter did not run yet ) * * @return \sa TransformMatrixType */ TransformMatrixType GetLastRotationMatrix(); protected: PyramidImageRegistrationMethod(); ~PyramidImageRegistrationMethod(); /** Fixed image, used as reference for registration */ mitk::Image::Pointer m_FixedImage; /** Moving image, will be transformed */ mitk::Image::Pointer m_MovingImage; mitk::Image::Pointer m_FixedImageMask; bool m_CrossModalityRegistration; bool m_UseAffineTransform; bool m_UseWindowedSincInterpolator; + bool m_UseNearestNeighborInterpolator; + bool m_UseMask; double* m_EstimatedParameters; /** Control the verbosity of the regsitistration output */ bool m_Verbose; bool m_InitializeByGeometry; template 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 = 6; float optminstep = 0.05f; unsigned int iterations = 100; if( min_value / max_schedule_val < 12 ) { max_schedule_val /= 2; optmaxstep *= 0.25f; optminstep *= 0.1f; } 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 ); // Initializing by Geometry if (!m_UseAffineTransform && m_InitializeByGeometry) { typedef itk::CenteredVersorTransformInitializer< FixedImageType, MovingImageType> TransformInitializerType; typename TransformInitializerType::TransformType::Pointer rigidTransform = TransformInitializerType::TransformType::New() ; MITK_INFO << "Initializer starting at : " << rigidTransform->GetParameters(); typename TransformInitializerType::Pointer initializer = TransformInitializerType::New(); initializer->SetTransform( rigidTransform); initializer->SetFixedImage( referenceImage.GetPointer() ); initializer->SetMovingImage( movingImage.GetPointer() ); initializer->MomentsOn(); initializer->InitializeTransform(); MITK_INFO << "Initialized Rigid position : " << rigidTransform->GetParameters(); initialParams[3]=rigidTransform->GetParameters()[3]; initialParams[4]=rigidTransform->GetParameters()[4]; initialParams[5]=rigidTransform->GetParameters()[5]; } // add observer tag if verbose */ unsigned long vopt_tag = 0; if(m_Verbose) { MITK_INFO << " :: Starting at " << initialParams; MITK_INFO << " :: Scales = " << optScales; OptimizerIterationCommand< OptimizerType >::Pointer iterationObserver = OptimizerIterationCommand::New(); vopt_tag = optimizer->AddObserver( itk::IterationEvent(), iterationObserver ); } // INPUT IMAGES registration->SetFixedImage( referenceImage ); registration->SetFixedImageRegion( maskedRegion ); registration->SetMovingImage( movingImage ); registration->SetSchedules( fixedSchedule, fixedSchedule); // SET MASKED AREA typedef itk::Image BinaryImageType; BinaryImageType::Pointer itkFixedImageMask = BinaryImageType::New(); itk::ImageMaskSpatialObject< 3 >::Pointer fixedMaskSpatialObject = itk::ImageMaskSpatialObject< 3 >::New(); if (m_UseMask) { CastToItkImage(m_FixedImageMask,itkFixedImageMask); itk::NotImageFilter::Pointer notFilter = itk::NotImageFilter::New(); notFilter->SetInput(itkFixedImageMask); notFilter->Update(); fixedMaskSpatialObject->SetImage( notFilter->GetOutput() ); metric->SetFixedImageMask( fixedMaskSpatialObject ); } // 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; iRemoveObserver( vopt_tag ); } } template< typename TPixel, unsigned int VDimension> 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 linear_interpolator = InterpolatorType::New(); + + typedef itk::NearestNeighborInterpolateImageFunction< ImageType, double > NearestNeighborInterpolatorType; + typename NearestNeighborInterpolatorType::Pointer nn_interpolator = NearestNeighborInterpolatorType::New(); + + typedef itk::WindowedSincInterpolateImageFunction< ImageType, 7> WindowedSincInterpolatorType; typename WindowedSincInterpolatorType::Pointer sinc_interpolator = WindowedSincInterpolatorType::New(); typename ImageType::Pointer reference_image = ImageType::New(); CastToItkImage(m_FixedImage,reference_image); 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(); + typedef itk::AffineTransform< double > TransformType; + TransformType::Pointer transform = TransformType::New(); - TransformType::ParametersType affine_params( TransformType::ParametersDimension ); - this->GetParameters( &affine_params[0] ); + TransformType::ParametersType affine_params( TransformType::ParametersDimension ); + this->GetParameters( &affine_params[0] ); - transform->SetParameters( affine_params ); + transform->SetParameters( affine_params ); - base_transform = transform; + base_transform = transform; } // Rigid else { - typedef itk::Euler3DTransform< double > RigidTransformType; - RigidTransformType::Pointer rtransform = RigidTransformType::New(); + typedef itk::Euler3DTransform< double > RigidTransformType; + RigidTransformType::Pointer rtransform = RigidTransformType::New(); - RigidTransformType::ParametersType rigid_params( RigidTransformType::ParametersDimension ); - this->GetParameters( &rigid_params[0] ); + RigidTransformType::ParametersType rigid_params( RigidTransformType::ParametersDimension ); + this->GetParameters( &rigid_params[0] ); - rtransform->SetParameters( rigid_params ); + rtransform->SetParameters( rigid_params ); - base_transform = rtransform; + base_transform = rtransform; } typename ResampleImageFilterType::Pointer resampler = ResampleImageFilterType::New(); + resampler->SetInterpolator( linear_interpolator ); if( m_UseWindowedSincInterpolator ) resampler->SetInterpolator( sinc_interpolator ); + if ( m_UseNearestNeighborInterpolator) + resampler->SetInterpolator ( nn_interpolator ); resampler->SetInput( itkImage ); resampler->SetTransform( base_transform ); resampler->SetReferenceImage( reference_image ); resampler->UseReferenceImageOn(); resampler->Update(); - mitk::GrabItkImageMemory( resampler->GetOutput(), outputImage); } }; } // end namespace #endif // MITKPYRAMIDIMAGEREGISTRATION_H