diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp index 2eb5e274cc..d55b2ca14a 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp @@ -1,284 +1,290 @@ /*=================================================================== 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" #include 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(); // use the advanced (windowed sinc) interpolation registrationMethod->SetUseAdvancedInterpolation(true); // 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; i 1) { - t_selector2->SetTimeNr(i); - t_selector2->Update(); + mitk::ImageTimeSelector::Pointer t_selector2 = + mitk::ImageTimeSelector::New(); - registrationMethod->SetMovingImage( t_selector2->GetOutput() ); + t_selector2->SetInput( b0Image ); - try + for( unsigned int i=1; iUpdate(); - } - catch( const itk::ExceptionObject& e) - { - mitkThrow() << "Failed to register the b0 images, the PyramidRegistration threw an exception: \n" << e.what(); + + t_selector2->SetTimeNr(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 ); + } - // import volume to the inter-results - registeredB0Image->SetImportVolume( registrationMethod->GetResampledMovingImage()->GetData(), - i, 0, mitk::Image::ReferenceMemory ); + // use the accumulateImageFilter as provided by the ItkAccumulateFilter method in the header file + AccessFixedDimensionByItk_1(registeredB0Image, ItkAccumulateFilter, (4), b0referenceImage ); } // // (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 + // - (3.1) Set the reference image + // - a single b0 image + // - average over the registered b0 images if multiple present // - - // use the accumulateImageFilter as provided by the ItkAccumulateFilter method in the header file - AccessFixedDimensionByItk_1(registeredB0Image, ItkAccumulateFilter, (4), b0referenceImage ); - weightedRegistrationMethod->SetFixedImage( b0referenceImage ); // use the advanced (windowed sinc) interpolation weightedRegistrationMethod->SetUseAdvancedInterpolation(true); // // - (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/mitkDWIHeadMotionCorrectionFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.h index 92f2b08587..196383c65f 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.h @@ -1,105 +1,119 @@ /*=================================================================== 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 MITKDWIHEADMOTIONCORRECTIONFILTER_H #define MITKDWIHEADMOTIONCORRECTIONFILTER_H #include "mitkDiffusionImageToDiffusionImageFilter.h" #include +#include + #include "mitkITKImageImport.h" namespace mitk { /** * @class DWIHeadMotionCorrectionFilter * * @brief Performs standard head-motion correction by using affine registration of the gradient images. * * (Head) motion correction is a essential pre-processing step before performing any further analysis of a diffusion-weighted * images since all model fits ( tensor, QBI ) rely on an aligned diffusion-weighted dataset. The correction is done in two steps. First the * unweighted images ( if multiple present ) are separately registered on the first one by means of rigid registration and normalized correlation * as error metric. Second, the weighted gradient images are registered to the unweighted reference ( computed as average from the aligned images from first step ) * by an affine transformation using the MattesMutualInformation metric as optimizer guidance. * */ template< typename DiffusionPixelType> class DWIHeadMotionCorrectionFilter : public DiffusionImageToDiffusionImageFilter< DiffusionPixelType > { public: // class macros mitkClassMacro( DWIHeadMotionCorrectionFilter, DiffusionImageToDiffusionImageFilter ) itkNewMacro(Self) // public typedefs typedef typename Superclass::InputImageType DiffusionImageType; typedef typename Superclass::InputImagePointerType DiffusionImagePointerType; typedef typename Superclass::OutputImageType OutputImageType; typedef typename Superclass::OutputImagePointerType OutputImagePointerType; protected: DWIHeadMotionCorrectionFilter(); virtual ~DWIHeadMotionCorrectionFilter() {} virtual void GenerateData(); /** * @brief Averages an 3d+t image along the time axis. * * The method uses the AccumulateImageFilter as provided by ITK and collapses the given 3d+t image * to an 3d image while computing the average over the time axis for each of the spatial voxels. */ template< typename TPixel, unsigned int VDimensions> static void ItkAccumulateFilter( const itk::Image< TPixel, VDimensions>* image, mitk::Image::Pointer& output) { // input 3d+t --> output 3d typedef itk::Image< TPixel, 4> InputItkType; typedef itk::Image< TPixel, 3> OutputItkType; - typedef typename itk::AccumulateImageFilter< InputItkType, OutputItkType > FilterType; + // the accumulate filter requires the same dimension in output and input image + typedef typename itk::AccumulateImageFilter< InputItkType, InputItkType > FilterType; typename FilterType::Pointer filter = FilterType::New(); filter->SetInput( image ); filter->SetAccumulateDimension( 3 ); filter->SetAverage( true ); + // we need to extract the volume to reduce the size from 4 to 3 for further processing + typedef typename itk::ExtractImageFilter< InputItkType, OutputItkType > ExtractFilterType; + typename ExtractFilterType::Pointer extractor = ExtractFilterType::New(); + extractor->SetInput( filter->GetOutput() ); + extractor->SetDirectionCollapseToIdentity(); + + typename InputItkType::RegionType extractRegion = image->GetLargestPossibleRegion(); + // crop out the time axis + extractRegion.SetSize( 3, 0); + extractor->SetExtractionRegion( extractRegion ); + try { - filter->Update(); + extractor->Update(); } catch( const itk::ExceptionObject& e) { mitkThrow() << " Exception while averaging: " << e.what(); } - output = mitk::GrabItkImageMemory( filter->GetOutput() ); + output = mitk::GrabItkImageMemory( extractor->GetOutput() ); } }; } //end namespace mitk #include "mitkDWIHeadMotionCorrectionFilter.cpp" #endif // MITKDWIHEADMOTIONCORRECTIONFILTER_H