diff --git a/Modules/DiffusionImaging/DiffusionCore/src/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp b/Modules/DiffusionImaging/DiffusionCore/src/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp index e011bdb11c..6da075a982 100644 --- a/Modules/DiffusionImaging/DiffusionCore/src/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/src/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp @@ -1,317 +1,317 @@ /*=================================================================== 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_CPP #define MITKDWIHEADMOTIONCORRECTIONFILTER_CPP #include "mitkDWIHeadMotionCorrectionFilter.h" #include "itkSplitDWImageFilter.h" #include "itkB0ImageExtractionToSeparateImageFilter.h" #include "mitkImageTimeSelector.h" #include "mitkPyramidImageRegistrationMethod.h" //#include "mitkRegistrationMethodITK4.h" #include #include "mitkDiffusionImageCorrectionFilter.h" #include #include #include #include "mitkIOUtil.h" #include mitk::DWIHeadMotionCorrectionFilter::DWIHeadMotionCorrectionFilter() : m_CurrentStep(0), m_Steps(100), m_IsInValidState(true), m_AbortRegistration(false), m_AverageUnweighted(true) { } void mitk::DWIHeadMotionCorrectionFilter::GenerateData() { typedef itk::SplitDWImageFilter< DiffusionPixelType, DiffusionPixelType> SplitFilterType; InputImageType* input = const_cast(this->GetInput(0)); ITKDiffusionImageType::Pointer itkVectorImagePointer = ITKDiffusionImageType::New(); mitk::CastToItkImage(input, itkVectorImagePointer); typedef mitk::PyramidImageRegistrationMethod RegistrationMethod; // typedef mitk::RegistrationMethodITKV4 RegistrationMethod unsigned int numberOfSteps = itkVectorImagePointer->GetNumberOfComponentsPerPixel () ; m_Steps = numberOfSteps; // // (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; B0ExtractorType::Pointer b0_extractor = B0ExtractorType::New(); b0_extractor->SetInput( itkVectorImagePointer ); b0_extractor->SetDirections( static_cast( input->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer() ); 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(); const unsigned int numberOfb0Images = b0Image->GetTimeSteps(); // register b-zeros only if the flag to average is set, use the first one otherwise if( m_AverageUnweighted && numberOfb0Images > 1) { RegistrationMethod::Pointer registrationMethod = RegistrationMethod::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(); 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) { m_IsInValidState = false; mitkThrow() << "Failed to register the b0 images, the PyramidRegistration threw an exception: \n" << e.what(); } // import volume to the inter-results mitk::ImageWriteAccessor imac(registrationMethod->GetResampledMovingImage()); registeredB0Image->SetImportVolume( imac.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 // SplitFilterType::Pointer split_filter = SplitFilterType::New(); split_filter->SetInput (itkVectorImagePointer ); split_filter->SetExtractAllAboveThreshold(8, static_cast(input->GetProperty(mitk::DiffusionPropertyHelper::BVALUEMAPPROPERTYNAME.c_str()).GetPointer() )->GetBValueMap() ); try { split_filter->Update(); } catch( const itk::ExceptionObject &e) { m_IsInValidState = false; 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 // RegistrationMethod::Pointer weightedRegistrationMethod = RegistrationMethod::New(); weightedRegistrationMethod->SetTransformToAffine(); weightedRegistrationMethod->SetCrossModalityOn(); // // - (3.1) Set the reference image // - a single b0 image // - average over the registered b0 images if multiple present // weightedRegistrationMethod->SetFixedImage( b0referenceImage ); // use the advanced (windowed sinc) interpolation weightedRegistrationMethod->SetUseAdvancedInterpolation(true); - weightedRegistrationMethod->SetVerboseOn(); + weightedRegistrationMethod->SetVerboseOff(); // // - (3.2) Register all timesteps in the splitted image onto the first reference // unsigned int maxImageIdx = splittedImage->GetTimeSteps(); mitk::TimeGeometry* tsg = splittedImage->GetTimeGeometry(); mitk::ProportionalTimeGeometry* ptg = dynamic_cast(tsg); ptg->Expand(maxImageIdx+1); ptg->SetTimeStepGeometry( ptg->GetGeometryForTimeStep(0), maxImageIdx ); mitk::Image::Pointer registeredWeighted = mitk::Image::New(); registeredWeighted->Initialize( splittedImage->GetPixelType(0), *tsg ); // insert the first unweighted reference as the first volume // in own scope to release the accessor asap after copy { mitk::ImageWriteAccessor imac(b0referenceImage); registeredWeighted->SetImportVolume( imac.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 RegistrationMethod::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) { m_IsInValidState = false; mitkThrow() << "Failed to register the b0 images, the PyramidRegistration threw an exception: \n" << e.what(); } // allow expansion mitk::ImageWriteAccessor imac(weightedRegistrationMethod->GetResampledMovingImage()); registeredWeighted->SetImportVolume( imac.GetData(), i+1, 0, mitk::Image::CopyMemory); estimated_transforms.push_back( weightedRegistrationMethod->GetLastRotationMatrix() ); } // // (4) Cast the resulting image back to an diffusion weighted image // mitk::DiffusionPropertyHelper::GradientDirectionsContainerType *gradients = static_cast( input->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer(); mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer gradients_new = mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::New(); mitk::DiffusionPropertyHelper::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 ); SplitFilterType::IndexListType index_list = split_filter->GetIndexList(); SplitFilterType::IndexListType::const_iterator lIter = index_list.begin(); while( lIter != index_list.end() ) { gradients_new->push_back( gradients->at( *lIter ) ); ++lIter; } registeredWeighted->GetPropertyList()->ReplaceProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( gradients_new.GetPointer() )); registeredWeighted->GetPropertyList()->ReplaceProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( static_cast(input->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() )->GetValue() ) ); // // (5) Adapt the gradient directions according to the estimated transforms // typedef mitk::DiffusionImageCorrectionFilter CorrectionFilterType; CorrectionFilterType::Pointer corrector = CorrectionFilterType::New(); mitk::Image::Pointer output = registeredWeighted; corrector->SetImage( output ); corrector->CorrectDirections( estimated_transforms ); // // (6) Pass the corrected image to the filters output port // m_CurrentStep += 1; this->GetOutput()->Initialize( output ); mitk::DiffusionPropertyHelper::CopyProperties(output, this->GetOutput(), true); this->GetOutput()->GetPropertyList()->ReplaceProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( static_cast( output->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer() ) ); this->GetOutput()->GetPropertyList()->ReplaceProperty( mitk::DiffusionPropertyHelper::MEASUREMENTFRAMEPROPERTYNAME.c_str(), mitk::MeasurementFrameProperty::New( static_cast( output->GetProperty(mitk::DiffusionPropertyHelper::MEASUREMENTFRAMEPROPERTYNAME.c_str()).GetPointer() )->GetMeasurementFrame() ) ); this->GetOutput()->GetPropertyList()->ReplaceProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( static_cast(output->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() )->GetValue() ) ); mitk::DiffusionPropertyHelper propertyHelper( this->GetOutput() ); propertyHelper.InitializeImage(); this->GetOutput()->Modified(); } #endif // MITKDWIHEADMOTIONCORRECTIONFILTER_CPP diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.denoising/src/internal/QmitkDenoisingView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.denoising/src/internal/QmitkDenoisingView.cpp index 9d7e079b18..e3ed587f8d 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.denoising/src/internal/QmitkDenoisingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.denoising/src/internal/QmitkDenoisingView.cpp @@ -1,260 +1,260 @@ /*=================================================================== 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. ===================================================================*/ // Blueberry #include #include // Qmitk #include "QmitkDenoisingView.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include const std::string QmitkDenoisingView::VIEW_ID = "org.mitk.views.denoisingview"; QmitkDenoisingView::QmitkDenoisingView() : QmitkAbstractView() , m_Controls( 0 ) , m_ImageNode(nullptr) , m_SelectedFilter(TV) { } QmitkDenoisingView::~QmitkDenoisingView() { } void QmitkDenoisingView::CreateQtPartControl( QWidget *parent ) { // build up qt view, unless already done if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkDenoisingViewControls; m_Controls->setupUi( parent ); m_Controls->m_SelectFilterComboBox->clear(); m_Controls->m_SelectFilterComboBox->insertItem(TV, "Total-variation filter"); m_Controls->m_SelectFilterComboBox->insertItem(GAUSS, "Discrete gaussian filter"); m_Controls->m_SelectFilterComboBox->insertItem(NLM, "Non-local means filter"); // m_Controls->m_SelectFilterComboBox->insertItem(NLM_MORITZ, "Non-local means filter"); connect( (QObject*)(m_Controls->m_ApplyButton), SIGNAL(clicked()), this, SLOT(StartDenoising())); connect( (QObject*)(m_Controls->m_SelectFilterComboBox), SIGNAL(activated(int)), this, SLOT(UpdateGui(int))); connect( (QObject*)(m_Controls->m_InputImageBox), SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui(int))); m_Controls->m_InputImageBox->SetDataStorage(this->GetDataStorage()); mitk::TNodePredicateDataType::Pointer isImagePredicate = mitk::TNodePredicateDataType::New(); m_Controls->m_InputImageBox->SetPredicate( isImagePredicate ); UpdateGui(0); } } void QmitkDenoisingView::SetFocus() { m_Controls->m_SelectFilterComboBox->setFocus(); } void QmitkDenoisingView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList& ) { } void QmitkDenoisingView::StartDenoising() { typedef itk::NonLocalMeansDenoisingFilter< DiffusionPixelType > NlmRicianFilterType; typedef itk::DiscreteGaussianImageFilter < DwiVolumeType, DwiVolumeType > GaussianFilterType; typedef itk::PatchBasedDenoisingImageFilter < DwiVolumeType, DwiVolumeType > NlmFilterType; typedef itk::VectorImageToImageFilter < DiffusionPixelType > ExtractFilterType; typedef itk::ComposeImageFilter < itk::Image > ComposeFilterType; m_ImageNode = m_Controls->m_InputImageBox->GetSelectedNode(); mitk::Image::Pointer input_image = dynamic_cast(m_ImageNode->GetData()); if (!mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(input_image)) return; DwiImageType::Pointer itkVectorImagePointer = mitk::DiffusionPropertyHelper::GetItkVectorImage(input_image); mitk::Image::Pointer denoised_image = nullptr; mitk::DataNode::Pointer denoised_image_node = mitk::DataNode::New(); switch (m_SelectedFilter) { case TV: { ComposeFilterType::Pointer composer = ComposeFilterType::New(); for (unsigned int i=0; iGetVectorLength(); ++i) { ExtractFilterType::Pointer extractor = ExtractFilterType::New(); extractor->SetInput( itkVectorImagePointer ); extractor->SetIndex( i ); extractor->Update(); DwiVolumeType::Pointer gradient_volume = extractor->GetOutput(); itk::TotalVariationDenoisingImageFilter< DwiVolumeType, DwiVolumeType >::Pointer filter = itk::TotalVariationDenoisingImageFilter< DwiVolumeType, DwiVolumeType >::New(); filter->SetInput(gradient_volume); filter->SetLambda(m_Controls->m_TvLambdaBox->value()); filter->SetNumberIterations(m_Controls->m_TvIterationsBox->value()); filter->Update(); composer->SetInput(i, filter->GetOutput()); } composer->Update(); denoised_image = mitk::GrabItkImageMemory(composer->GetOutput()); denoised_image_node->SetName( "TotalVariation" ); break; } case GAUSS: { ExtractFilterType::Pointer extractor = ExtractFilterType::New(); extractor->SetInput( itkVectorImagePointer ); ComposeFilterType::Pointer composer = ComposeFilterType::New(); for (unsigned int i = 0; i < itkVectorImagePointer->GetVectorLength(); ++i) { extractor->SetIndex(i); extractor->Update(); GaussianFilterType::Pointer filter = GaussianFilterType::New(); filter->SetVariance(m_Controls->m_GaussVarianceBox->value()); filter->SetInput(extractor->GetOutput()); filter->Update(); composer->SetInput(i, filter->GetOutput()); } composer->Update(); denoised_image = mitk::GrabItkImageMemory(composer->GetOutput()); denoised_image_node->SetName( "Gauss" ); break; } case NLM: { - typedef itk::Statistics::GaussianRandomSpatialNeighborSubsampler< typename NlmFilterType::PatchSampleType, typename DwiVolumeType::RegionType > SamplerType; + typedef itk::Statistics::GaussianRandomSpatialNeighborSubsampler< NlmFilterType::PatchSampleType, DwiVolumeType::RegionType > SamplerType; // sampling the image to find similar patches - typename SamplerType::Pointer sampler = SamplerType::New(); + SamplerType::Pointer sampler = SamplerType::New(); sampler->SetRadius( m_Controls->m_NlmSearchRadiusBox->value() ); sampler->SetVariance( m_Controls->m_NlmSearchRadiusBox->value()*m_Controls->m_NlmSearchRadiusBox->value() ); sampler->SetNumberOfResultsRequested( m_Controls->m_NlmNumPatchesBox->value() ); MITK_INFO << "Starting NLM denoising"; ExtractFilterType::Pointer extractor = ExtractFilterType::New(); extractor->SetInput( itkVectorImagePointer ); ComposeFilterType::Pointer composer = ComposeFilterType::New(); for (unsigned int i = 0; i < itkVectorImagePointer->GetVectorLength(); ++i) { extractor->SetIndex(i); extractor->Update(); NlmFilterType::Pointer filter = NlmFilterType::New(); filter->SetInput(extractor->GetOutput()); filter->SetPatchRadius(m_Controls->m_NlmPatchSizeBox->value()); filter->SetNoiseModel(NlmFilterType::RICIAN); filter->UseSmoothDiscPatchWeightsOn(); filter->UseFastTensorComputationsOn(); filter->SetNumberOfIterations(m_Controls->m_NlmIterationsBox->value()); filter->SetSmoothingWeight( 1 ); filter->SetKernelBandwidthEstimation(true); filter->SetSampler( sampler ); filter->Update(); composer->SetInput(i, filter->GetOutput()); MITK_INFO << "Gradient " << i << " finished"; } composer->Update(); denoised_image = mitk::GrabItkImageMemory(composer->GetOutput()); denoised_image_node->SetName( "NLM" ); break; } case NLM_MORITZ: { NlmRicianFilterType::Pointer nlmFilter = NlmRicianFilterType::New(); nlmFilter->SetInputImage( itkVectorImagePointer ); // nlmFilter->SetUseRicianAdaption(m_Controls->m_RicianCheckbox->isChecked()); // nlmFilter->SetUseJointInformation(m_Controls->m_JointInformationCheckbox->isChecked()); nlmFilter->SetSearchRadius(m_Controls->m_NlmSearchRadiusBox->value()); nlmFilter->SetComparisonRadius(m_Controls->m_NlmPatchSizeBox->value()); // nlmFilter->SetVariance(m_Controls->m_NlmVarianceBox->value()); nlmFilter->Update(); denoised_image = mitk::GrabItkImageMemory(nlmFilter->GetOutput()); denoised_image_node->SetName( "NLM-Rician" ); break; } } denoised_image->GetPropertyList()->ReplaceProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( static_cast( input_image->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer() ) ); denoised_image->SetProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( static_cast(input_image->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() )->GetValue() ) ); mitk::DiffusionPropertyHelper propertyHelper( denoised_image ); propertyHelper.InitializeImage(); denoised_image_node->SetData( denoised_image ); GetDataStorage()->Add(denoised_image_node, m_ImageNode); } void QmitkDenoisingView::UpdateGui(int filter) { m_Controls->m_ApplyButton->setEnabled(false); m_Controls->m_TvFrame->setVisible(false); m_Controls->m_NlmFrame->setVisible(false); m_Controls->m_GaussFrame->setVisible(false); switch (filter) { case 0: { m_SelectedFilter = TV; m_Controls->m_TvFrame->setVisible(true); break; } case 1: { m_SelectedFilter = GAUSS; m_Controls->m_GaussFrame->setVisible(true); break; } case 2: { m_SelectedFilter = NLM; m_Controls->m_NlmFrame->setVisible(true); break; } case 3: { m_SelectedFilter = NLM_MORITZ; m_Controls->m_NlmFrame->setVisible(true); break; } } if (m_Controls->m_InputImageBox->GetSelectedNode().IsNotNull()) m_Controls->m_ApplyButton->setEnabled(true); }