diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp index 6959af5467..e82b5c81fd 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp @@ -1,253 +1,275 @@ /*=================================================================== 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 MITKDIFFUSIONIMAGETODIFFUSIONIMAGEFILTER_CPP #define MITKDIFFUSIONIMAGETODIFFUSIONIMAGEFILTER_CPP #include "mitkDWIHeadMotionCorrectionFilter.h" #include "itkSplitDWImageFilter.h" #include "itkB0ImageExtractionToSeparateImageFilter.h" #include "mitkImageTimeSelector.h" #include "mitkPyramidImageRegistrationMethod.h" #include "mitkImageToDiffusionImageSource.h" +#include "mitkDiffusionImageCorrectionFilter.h" + +#include + #include "mitkIOUtil.h" template< typename DiffusionPixelType> mitk::DWIHeadMotionCorrectionFilter ::DWIHeadMotionCorrectionFilter() { } template< typename DiffusionPixelType> void mitk::DWIHeadMotionCorrectionFilter ::GenerateData() { typedef itk::SplitDWImageFilter< DiffusionPixelType, DiffusionPixelType> SplitFilterType; DiffusionImageType* input = const_cast(this->GetInput(0)); // // (1) Extract the b-zero images to a 3d+t image, register them by NCorr metric and // rigid registration : they will then be used are reference image for registering // the gradient images // typedef itk::B0ImageExtractionToSeparateImageFilter< DiffusionPixelType, DiffusionPixelType> B0ExtractorType; typename B0ExtractorType::Pointer b0_extractor = B0ExtractorType::New(); b0_extractor->SetInput( input->GetVectorImage() ); b0_extractor->SetDirections( input->GetDirections() ); b0_extractor->Update(); mitk::Image::Pointer b0Image = mitk::Image::New(); b0Image->InitializeByItk( b0_extractor->GetOutput() ); b0Image->SetImportChannel( b0_extractor->GetOutput()->GetBufferPointer(), mitk::Image::CopyMemory ); // (2.1) Use the extractor to access the extracted b0 volumes mitk::ImageTimeSelector::Pointer t_selector = mitk::ImageTimeSelector::New(); t_selector->SetInput( b0Image ); t_selector->SetTimeNr(0); t_selector->Update(); // first unweighted image as reference space for the registration mitk::Image::Pointer b0referenceImage = t_selector->GetOutput(); mitk::PyramidImageRegistrationMethod::Pointer registrationMethod = mitk::PyramidImageRegistrationMethod::New(); registrationMethod->SetFixedImage( b0referenceImage ); registrationMethod->SetTransformToRigid(); // the unweighted images are of same modality registrationMethod->SetCrossModalityOff(); // Initialize the temporary output image mitk::Image::Pointer registeredB0Image = b0Image->Clone(); const unsigned int numberOfb0Images = b0Image->GetTimeSteps(); mitk::ImageTimeSelector::Pointer t_selector2 = mitk::ImageTimeSelector::New(); t_selector2->SetInput( b0Image ); for( unsigned int i=1; iSetTimeNr(i); t_selector2->Update(); registrationMethod->SetMovingImage( t_selector2->GetOutput() ); try { MITK_INFO << " === (" << i <<"/"<< numberOfb0Images-1 << ") :: Starting registration"; registrationMethod->Update(); } catch( const itk::ExceptionObject& e) { mitkThrow() << "Failed to register the b0 images, the PyramidRegistration threw an exception: \n" << e.what(); } // import volume to the inter-results registeredB0Image->SetImportVolume( registrationMethod->GetResampledMovingImage()->GetData(), i, 0, mitk::Image::ReferenceMemory ); } // // (2) Split the diffusion image into a 3d+t regular image, extract only the weighted images // typename SplitFilterType::Pointer split_filter = SplitFilterType::New(); split_filter->SetInput (input->GetVectorImage() ); split_filter->SetExtractAllAboveThreshold(20, input->GetB_ValueMap() ); try { split_filter->Update(); } catch( const itk::ExceptionObject &e) { mitkThrow() << " Caught exception from SplitImageFilter : " << e.what(); } mitk::Image::Pointer splittedImage = mitk::Image::New(); splittedImage->InitializeByItk( split_filter->GetOutput() ); splittedImage->SetImportChannel( split_filter->GetOutput()->GetBufferPointer(), mitk::Image::CopyMemory ); // // (3) Use again the time-selector to access the components separately in order // to perform the registration of Image -> unweighted reference // mitk::PyramidImageRegistrationMethod::Pointer weightedRegistrationMethod = mitk::PyramidImageRegistrationMethod::New(); weightedRegistrationMethod->SetTransformToAffine(); weightedRegistrationMethod->SetCrossModalityOn(); // // - (3.1) Create a reference image by averaging the aligned b0 images // // !!!FIXME: For rapid prototyping using the first one // weightedRegistrationMethod->SetFixedImage( b0referenceImage ); // // - (3.2) Register all timesteps in the splitted image onto the first reference // unsigned int maxImageIdx = splittedImage->GetTimeSteps(); mitk::TimeSlicedGeometry* tsg = splittedImage->GetTimeSlicedGeometry(); tsg->ExpandToNumberOfTimeSteps( maxImageIdx+1 ); mitk::Image::Pointer registeredWeighted = mitk::Image::New(); registeredWeighted->Initialize( splittedImage->GetPixelType(0), *tsg ); // insert the first unweighted reference as the first volume registeredWeighted->SetImportVolume( b0referenceImage->GetData(), 0,0, mitk::Image::CopyMemory ); // mitk::Image::Pointer registeredWeighted = splittedImage->Clone(); // this time start at 0, we have only gradient images in the 3d+t file // the reference image comes form an other image mitk::ImageTimeSelector::Pointer t_selector_w = mitk::ImageTimeSelector::New(); t_selector_w->SetInput( splittedImage ); + // store the rotation parts of the transformations in a vector + typedef mitk::PyramidImageRegistrationMethod::TransformMatrixType MatrixType; + std::vector< MatrixType > estimated_transforms; + for( unsigned int i=0; iSetTimeNr(i); t_selector_w->Update(); weightedRegistrationMethod->SetMovingImage( t_selector_w->GetOutput() ); try { MITK_INFO << " === (" << i+1 <<"/"<< maxImageIdx << ") :: Starting registration"; weightedRegistrationMethod->Update(); } catch( const itk::ExceptionObject& e) { mitkThrow() << "Failed to register the b0 images, the PyramidRegistration threw an exception: \n" << e.what(); } // allow expansion registeredWeighted->SetImportVolume( weightedRegistrationMethod->GetResampledMovingImage()->GetData(), i+1, 0, mitk::Image::CopyMemory); + + estimated_transforms.push_back( weightedRegistrationMethod->GetLastRotationMatrix() ); } // // (4) Cast the resulting image back to an diffusion weighted image // typename DiffusionImageType::GradientDirectionContainerType *gradients = input->GetDirections(); typename DiffusionImageType::GradientDirectionContainerType::Pointer gradients_new = DiffusionImageType::GradientDirectionContainerType::New(); typename DiffusionImageType::GradientDirectionType bzero_vector; bzero_vector.fill(0); // compose the direction vector // - no direction for the first image // - correct ordering of the directions based on the index list gradients_new->push_back( bzero_vector ); typename SplitFilterType::IndexListType index_list = split_filter->GetIndexList(); typename SplitFilterType::IndexListType::const_iterator lIter = index_list.begin(); while( lIter != index_list.end() ) { gradients_new->push_back( gradients->at( *lIter ) ); ++lIter; } typename mitk::ImageToDiffusionImageSource< DiffusionPixelType >::Pointer caster = mitk::ImageToDiffusionImageSource< DiffusionPixelType >::New(); caster->SetImage( registeredWeighted ); caster->SetBValue( input->GetB_Value() ); caster->SetGradientDirections( gradients_new.GetPointer() ); try { caster->Update(); } catch( const itk::ExceptionObject& e) { MITK_ERROR << "Casting back to diffusion image failed: "; mitkThrow() << "Subprocess failed with exception: " << e.what(); } + // + // (5) Adapt the gradient directions according to the estimated transforms + // + typedef mitk::DiffusionImageCorrectionFilter< DiffusionPixelType > CorrectionFilterType; + typename CorrectionFilterType::Pointer corrector = CorrectionFilterType::New(); + OutputImagePointerType output = caster->GetOutput(); + corrector->SetImage( output ); + corrector->CorrectDirections( estimated_transforms ); + + // + // (6) Pass the corrected image to the filters output port + // this->GetOutput()->SetVectorImage(output->GetVectorImage()); this->GetOutput()->SetB_Value(output->GetB_Value()); this->GetOutput()->SetDirections(output->GetDirections()); this->GetOutput()->SetMeasurementFrame(output->GetMeasurementFrame()); this->GetOutput()->InitializeFromVectorImage(); this->GetOutput()->Modified(); } #endif // MITKDIFFUSIONIMAGETODIFFUSIONIMAGEFILTER_CPP diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidImageRegistrationMethod.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidImageRegistrationMethod.cpp index 1c7332ead3..d6b00b34ea 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidImageRegistrationMethod.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidImageRegistrationMethod.cpp @@ -1,97 +1,139 @@ /*=================================================================== 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_EstimatedParameters(NULL), m_Verbose(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::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; } diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidImageRegistrationMethod.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidImageRegistrationMethod.h index 2cdcaecb0b..ad36bcabe1 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidImageRegistrationMethod.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidImageRegistrationMethod.h @@ -1,412 +1,430 @@ /*=================================================================== 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 #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) + /** Typedef for the transformation matrix, corresponds to the InternalMatrixType from ITK transforms */ + typedef vnl_matrix_fixed< double, 3, 3> TransformMatrixType; + /** 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 ); 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; bool m_CrossModalityRegistration; bool m_UseAffineTransform; double* m_EstimatedParameters; /** Control the verbosity of the regsitistration output */ bool m_Verbose; 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 ); // 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); // 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 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/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.cpp b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.cpp new file mode 100644 index 0000000000..42112e7b0f --- /dev/null +++ b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.cpp @@ -0,0 +1,101 @@ +/*=================================================================== + +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 MITKDIFFUSIONIMAGECORRECTIONFILTER_CPP +#define MITKDIFFUSIONIMAGECORRECTIONFILTER_CPP + +#include "mitkDiffusionImageCorrectionFilter.h" + +#include +#include + +template< typename TPixelType > +mitk::DiffusionImageCorrectionFilter::DiffusionImageCorrectionFilter() +{ + +} + +template< typename TPixelType > +typename mitk::DiffusionImageCorrectionFilter::TransformMatrixType +mitk::DiffusionImageCorrectionFilter +::GetRotationComponent(const TransformMatrixType &A) +{ + TransformMatrixType B; + + B = A * A.transpose(); + + // get the eigenvalues and eigenvectors + typedef double MType; + vnl_vector< MType > eigvals; + vnl_matrix< MType > eigvecs; + vnl_symmetric_eigensystem_compute< MType > ( B, eigvecs, eigvals ); + + vnl_matrix_fixed< MType, 3, 3 > eigvecs_fixed; + eigvecs_fixed.set_columns(0, eigvecs ); + + TransformMatrixType C; + C.set_identity(); + + vnl_vector_fixed< MType, 3 > eigvals_sqrt; + for(unsigned int i=0; i<3; i++) + { + C(i,i) = std::sqrt( eigvals[i] ); + } + + TransformMatrixType S = vnl_inverse( eigvecs_fixed * C * vnl_inverse( eigvecs_fixed )) * A; + + return S; +} + +template< typename TPixelType > +void mitk::DiffusionImageCorrectionFilter +::CorrectDirections( const TransformsVectorType& transformations) +{ + if( m_SourceImage.IsNull() ) + { + mitkThrow() << " No diffusion image given! "; + } + + GradientDirectionContainerPointerType directions = m_SourceImage->GetDirections(); + GradientDirectionContainerPointerType corrected_directions = + GradientDirectionContainerType::New(); + + unsigned int transformed = 0; + for(size_t i=0; i< directions->Size(); i++ ) + { + + // skip b-zero images + if( directions->ElementAt(i).one_norm() <= 0.0 ) + { + corrected_directions->push_back( directions->ElementAt(i) ); + continue; + } + + GradientDirectionType corrected = GetRotationComponent( + transformations.at(transformed)) + * directions->ElementAt(i); + // store the corrected direction + corrected_directions->push_back( corrected ); + transformed++; + } + + // replace the old directions with the corrected ones + m_SourceImage->SetDirections( corrected_directions ); + +} + +#endif diff --git a/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.h b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.h new file mode 100644 index 0000000000..159cd5321b --- /dev/null +++ b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.h @@ -0,0 +1,90 @@ +/*=================================================================== + +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 MITKDIFFUSIONIMAGECORRECTIONFILTER_H +#define MITKDIFFUSIONIMAGECORRECTIONFILTER_H + +#include "mitkDiffusionImageSource.h" + +namespace mitk +{ +/** + * @class DiffusionImageCorrectionFilter + */ +template< typename TPixelType > +class DiffusionImageCorrectionFilter + : public DiffusionImageSource< TPixelType > +{ +public: + /** class macros */ + mitkClassMacro( DiffusionImageCorrectionFilter, + DiffusionImageSource ) + + itkSimpleNewMacro(Self) + + typedef vnl_vector_fixed< double, 3 > GradientDirectionType; + typedef vnl_matrix_fixed< double, 3, 3 > TransformMatrixType; + typedef itk::VectorContainer< unsigned int, GradientDirectionType > + GradientDirectionContainerType; + typedef GradientDirectionContainerType::Pointer GradientDirectionContainerPointerType; + + typedef std::vector< TransformMatrixType > TransformsVectorType; + + typedef typename Superclass::OutputType DiffusionImageType; + typedef typename DiffusionImageType::Pointer DiffusionImageTypePointer; + typedef itk::VectorImage ImageType; + + /** + * @brief Set the mitk image ( a 3d+t image ) which is to be reinterpreted as dw image + * @param mitkImage + */ + void SetImage( DiffusionImageTypePointer input ) + { + m_SourceImage = input; + } + + /** + * @brief Correct each gradient direction according to the given transform + * + * The size of the input is expected to correspond to the count of gradient images in the image. + */ + void CorrectDirections( const TransformsVectorType& ); + + +protected: + DiffusionImageCorrectionFilter(); + virtual ~DiffusionImageCorrectionFilter() {} + + /** + * @brief Get the rotation component following the Finite Strain + * + * For a given transformation \f$A\f$ its rotation component is defined as \f$ (AA^{T})^{-1/2}\f$. + * + * The computation first computes \f$ B = AA^T \f$ and then estimates the square root. Square root of + * diagonal matrices is defined as + * \f$ S = Q * \sqrt{C} * Q^{-1} \f$ with \f$ C \f$ having the eigenvalues on the diagonal. + * + */ + TransformMatrixType GetRotationComponent( const TransformMatrixType& ); + + DiffusionImageTypePointer m_SourceImage; +}; + +} + +#include "mitkDiffusionImageCorrectionFilter.cpp" + +#endif // MITKDIFFUSIONIMAGECORRECTIONFILTER_H diff --git a/Modules/DiffusionImaging/DiffusionCore/files.cmake b/Modules/DiffusionImaging/DiffusionCore/files.cmake index 1740783a51..063f891c42 100644 --- a/Modules/DiffusionImaging/DiffusionCore/files.cmake +++ b/Modules/DiffusionImaging/DiffusionCore/files.cmake @@ -1,134 +1,135 @@ set(CPP_FILES # DicomImport DicomImport/mitkDicomDiffusionImageReader.cpp # DicomImport/mitkGroupDiffusionHeadersFilter.cpp DicomImport/mitkDicomDiffusionImageHeaderReader.cpp DicomImport/mitkGEDicomDiffusionImageHeaderReader.cpp DicomImport/mitkPhilipsDicomDiffusionImageHeaderReader.cpp DicomImport/mitkSiemensDicomDiffusionImageHeaderReader.cpp DicomImport/mitkSiemensMosaicDicomDiffusionImageHeaderReader.cpp # DataStructures IODataStructures/mitkDiffusionCoreObjectFactory.cpp # DataStructures -> DWI IODataStructures/DiffusionWeightedImages/mitkDiffusionImageHeaderInformation.cpp IODataStructures/DiffusionWeightedImages/mitkDiffusionImageSource.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageReader.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageWriter.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageIOFactory.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageWriterFactory.cpp IODataStructures/DiffusionWeightedImages/mitkDiffusionImageSerializer.cpp IODataStructures/DiffusionWeightedImages/mitkImageToDiffusionImageSource.cpp + IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.cpp # DataStructures -> QBall IODataStructures/QBallImages/mitkQBallImageSource.cpp IODataStructures/QBallImages/mitkNrrdQBallImageReader.cpp IODataStructures/QBallImages/mitkNrrdQBallImageWriter.cpp IODataStructures/QBallImages/mitkNrrdQBallImageIOFactory.cpp IODataStructures/QBallImages/mitkNrrdQBallImageWriterFactory.cpp IODataStructures/QBallImages/mitkQBallImage.cpp IODataStructures/QBallImages/mitkQBallImageSerializer.cpp # DataStructures -> Tensor IODataStructures/TensorImages/mitkTensorImageSource.cpp IODataStructures/TensorImages/mitkNrrdTensorImageReader.cpp IODataStructures/TensorImages/mitkNrrdTensorImageWriter.cpp IODataStructures/TensorImages/mitkNrrdTensorImageIOFactory.cpp IODataStructures/TensorImages/mitkNrrdTensorImageWriterFactory.cpp IODataStructures/TensorImages/mitkTensorImage.cpp IODataStructures/TensorImages/mitkTensorImageSerializer.cpp #IODataStructures/mitkRegistrationObject.cpp # Rendering Rendering/vtkMaskedProgrammableGlyphFilter.cpp Rendering/mitkCompositeMapper.cpp Rendering/mitkVectorImageVtkGlyphMapper3D.cpp Rendering/vtkOdfSource.cxx Rendering/vtkThickPlane.cxx Rendering/mitkOdfNormalizationMethodProperty.cpp Rendering/mitkOdfScaleByProperty.cpp Rendering/mitkPlanarFigureMapper3D.cpp # Algorithms Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.cpp Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.cpp Algorithms/itkDwiGradientLengthCorrectionFilter.cpp # Registration Algorithms & Co. Algorithms/Registration/mitkPyramidImageRegistrationMethod.cpp # Function Collection mitkDiffusionFunctionCollection.cpp ) set(H_FILES # function Collection mitkDiffusionFunctionCollection.h # Rendering Rendering/mitkDiffusionImageMapper.h Rendering/mitkOdfVtkMapper2D.h Rendering/mitkPlanarFigureMapper3D.h # Reconstruction Algorithms/Reconstruction/itkDiffusionQballReconstructionImageFilter.h Algorithms/Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.h Algorithms/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.h Algorithms/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h Algorithms/Reconstruction/itkPointShell.h Algorithms/Reconstruction/itkOrientationDistributionFunction.h Algorithms/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.h Algorithms/Reconstruction/itkMultiShellRadialAdcKurtosisImageFilter.h # IO Datastructures IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h # Algorithms Algorithms/itkDiffusionQballGeneralizedFaImageFilter.h Algorithms/itkDiffusionQballPrepareVisualizationImageFilter.h Algorithms/itkTensorDerivedMeasurementsFilter.h Algorithms/itkBrainMaskExtractionImageFilter.h Algorithms/itkB0ImageExtractionImageFilter.h Algorithms/itkB0ImageExtractionToSeparateImageFilter.h Algorithms/itkTensorImageToDiffusionImageFilter.h Algorithms/itkTensorToL2NormImageFilter.h Algorithms/itkGaussianInterpolateImageFunction.h Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.h Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.h Algorithms/itkDiffusionTensorPrincipalDirectionImageFilter.h Algorithms/itkCartesianToPolarVectorImageFilter.h Algorithms/itkPolarToCartesianVectorImageFilter.h Algorithms/itkDistanceMapFilter.h Algorithms/itkProjectionFilter.h Algorithms/itkResidualImageFilter.h Algorithms/itkExtractChannelFromRgbaImageFilter.h Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.h Algorithms/itkMergeDiffusionImagesFilter.h Algorithms/itkDwiPhantomGenerationFilter.h Algorithms/itkFiniteDiffOdfMaximaExtractionFilter.h Algorithms/itkMrtrixPeakImageConverter.h Algorithms/itkFslPeakImageConverter.h Algorithms/itkShCoefficientImageImporter.h Algorithms/itkOdfMaximaExtractionFilter.h Algorithms/itkResampleDwiImageFilter.h Algorithms/itkDwiGradientLengthCorrectionFilter.h Algorithms/itkAdcImageFilter.h Algorithms/itkSplitDWImageFilter.h Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.h Algorithms/mitkDiffusionImageToDiffusionImageFilter.h ) set( TOOL_FILES )