diff --git a/Modules/DiffusionCore/Algorithms/Reconstruction/itkDiffusionKurtosisReconstructionImageFilter.cxx b/Modules/DiffusionCore/Algorithms/Reconstruction/itkDiffusionKurtosisReconstructionImageFilter.cxx index 37c66c7..8a26054 100644 --- a/Modules/DiffusionCore/Algorithms/Reconstruction/itkDiffusionKurtosisReconstructionImageFilter.cxx +++ b/Modules/DiffusionCore/Algorithms/Reconstruction/itkDiffusionKurtosisReconstructionImageFilter.cxx @@ -1,484 +1,493 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. 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 DIFFUSIONKURTOSISRECONSTRUCTIONIMAGEFILTER_CXX #define DIFFUSIONKURTOSISRECONSTRUCTIONIMAGEFILTER_CXX #include "itkDiffusionKurtosisReconstructionImageFilter.h" #include #include #include #include #include #include #include template< class TInputPixelType> static void FitSingleVoxel( const itk::VariableLengthVector< TInputPixelType > &input, const vnl_vector& bvalues, vnl_vector& result, itk::KurtosisFitConfiguration kf_config) { // check for length assert( input.Size() == bvalues.size() ); // assembly data vectors for fitting auto bvalueIter = bvalues.begin(); unsigned int unused_values = 0; while( bvalueIter != bvalues.end() ) { if( *bvalueIter < vnl_math::eps && kf_config.omit_bzero ) { ++unused_values; } ++bvalueIter; } // initialize data vectors with the estimated size (after filtering) vnl_vector fit_measurements( input.Size() - unused_values, 0 ); vnl_vector fit_bvalues( input.Size() - unused_values, 0 ); bvalueIter = bvalues.begin(); unsigned int running_index = 0; unsigned int skip_count = 0; while( bvalueIter != bvalues.end() ) { // skip b=0 if the corresponding flag is set if( ( kf_config.omit_bzero && *bvalueIter < vnl_math::eps ) // skip bvalues higher than the limit provided, if the flag is activated || ( kf_config.exclude_high_b && *bvalueIter > kf_config.b_upper_threshold ) ) { ++skip_count; } else { fit_measurements[ running_index - skip_count ] = input.GetElement(running_index); fit_bvalues[ running_index - skip_count] = *bvalueIter; } ++running_index; ++bvalueIter; } MITK_DEBUG("KurtosisFilter.FitSingleVoxel.Meas") << fit_measurements; MITK_DEBUG("KurtosisFilter.FitSingleVoxel.Bval") << fit_bvalues; // perform fit on data vectors if( kf_config.omit_bzero ) { itk::kurtosis_fit_omit_unweighted kurtosis_cost_fn( fit_measurements.size() ); // configuration kurtosis_cost_fn.set_fit_logscale( static_cast(kf_config.fit_scale) ); kurtosis_cost_fn.initialize( fit_measurements, fit_bvalues ); if( kf_config.use_K_limits) { kurtosis_cost_fn.set_K_bounds( kf_config.K_limits ); } vnl_levenberg_marquardt nonlinear_fit( kurtosis_cost_fn ); nonlinear_fit.minimize( result ); } else { itk::kurtosis_fit_lsq_function kurtosis_cost_fn( fit_measurements.size() ); // configuration kurtosis_cost_fn.set_fit_logscale( static_cast(kf_config.fit_scale) ); kurtosis_cost_fn.initialize( fit_measurements, fit_bvalues ); if( kf_config.use_K_limits) { kurtosis_cost_fn.set_K_bounds( kf_config.K_limits ); } vnl_levenberg_marquardt nonlinear_fit( kurtosis_cost_fn ); nonlinear_fit.minimize(result); } MITK_DEBUG("KurtosisFilter.FitSingleVoxel.Rslt") << result; } template< class TInputPixelType, class TOutputPixelType> itk::DiffusionKurtosisReconstructionImageFilter ::DiffusionKurtosisReconstructionImageFilter() : m_ReferenceBValue(-1), m_OmitBZero(false), m_ApplyPriorSmoothing(false), m_SmoothingSigma(1.5), m_UseKBounds( false ), m_MaxFitBValue( 3000 ), m_ScaleForFitting( STRAIGHT ) { this->m_InitialPosition = vnl_vector(3, 0); this->m_InitialPosition[2] = 1000.0; // S_0 this->m_InitialPosition[0] = 0.001; // D this->m_InitialPosition[1] = 1; // K this->m_KurtosisBounds.fill(0); // single input image this->SetNumberOfRequiredInputs(1); // two output images (D, K) this->SetNumberOfRequiredOutputs(2); typename OutputImageType::Pointer outputPtr1 = OutputImageType::New(); typename OutputImageType::Pointer outputPtr2 = OutputImageType::New(); this->SetNthOutput(0, outputPtr1.GetPointer() ); this->SetNthOutput(1, outputPtr2.GetPointer() ); } template< class TInputPixelType, class TOutputPixelType> void itk::DiffusionKurtosisReconstructionImageFilter ::GenerateOutputInformation() { // first call superclass Superclass::GenerateOutputInformation(); } template< class TInputPixelType, class TOutputPixelType> void itk::DiffusionKurtosisReconstructionImageFilter ::SetImageMask(MaskImageType::Pointer mask) { this->m_MaskImage = mask; } template< class TInputPixelType, class TOutputPixelType> void itk::DiffusionKurtosisReconstructionImageFilter ::BeforeThreadedGenerateData() { // if we have a region set, convert it to a mask image, which is to be used as default for telling // we need to set the image anyway, so by default the mask is overall 1 if( m_MaskImage.IsNull() ) { m_MaskImage = MaskImageType::New(); m_MaskImage->SetRegions( this->GetInput()->GetLargestPossibleRegion() ); m_MaskImage->CopyInformation( this->GetInput() ); m_MaskImage->Allocate(); if( this->m_MapOutputRegion.GetNumberOfPixels() > 0 ) { m_MaskImage->FillBuffer(0); typedef itk::ImageRegionIteratorWithIndex< MaskImageType > MaskIteratorType; MaskIteratorType maskIter( this->m_MaskImage, this->m_MapOutputRegion ); maskIter.GoToBegin(); while( !maskIter.IsAtEnd() ) { maskIter.Set( 1 ); ++maskIter; } } else { m_MaskImage->FillBuffer(1); } } // apply smoothing to the input image if( this->m_ApplyPriorSmoothing ) { // filter typedefs typedef itk::DiscreteGaussianImageFilter< itk::Image, itk::Image > GaussianFilterType; typedef itk::VectorIndexSelectionCastImageFilter< InputImageType, itk::Image > IndexSelectionType; typedef itk::ComposeImageFilter< itk::Image, InputImageType > ComposeFilterType; auto vectorImage = this->GetInput(); typename IndexSelectionType::Pointer indexSelectionFilter = IndexSelectionType::New(); indexSelectionFilter->SetInput( vectorImage ); typename ComposeFilterType::Pointer vec_composer = ComposeFilterType::New(); for( unsigned int i=0; iGetVectorLength(); ++i) { typename GaussianFilterType::Pointer gaussian_filter = GaussianFilterType::New(); indexSelectionFilter->SetIndex( i ); gaussian_filter->SetInput( indexSelectionFilter->GetOutput() ); gaussian_filter->SetVariance( m_SmoothingSigma ); vec_composer->SetInput(i, gaussian_filter->GetOutput() ); gaussian_filter->Update(); } try { vec_composer->Update(); } catch(const itk::ExceptionObject &e) { mitkThrow() << "[VectorImage.GaussianSmoothing] !! Failed with ITK Exception: " << e.what(); } this->m_ProcessedInputImage = vec_composer->GetOutput(); } else { this->m_ProcessedInputImage = const_cast( this->GetInput() ); } } template< class TInputPixelType, class TOutputPixelType> void itk::DiffusionKurtosisReconstructionImageFilter ::AfterThreadedGenerateData() { /* // initialize buffer to zero overall, but don't forget the requested region pointer for( unsigned int i=0; iGetNumberOfOutputs(); ++i) { typename OutputImageType::Pointer output = this->GetOutput(i); // after generating, set the buffered region to max, we have taken care about filling it with valid data in // UpdateOutputInformation, if not set, a writer (or any filter using the LargestPossibleRegion() during update may // complain about not getting the requested region output->SetBufferedRegion( output->GetLargestPossibleRegion() ); std::cout << "[DiffusionKurtosisReconstructionImageFilter.After]" << output->GetLargestPossibleRegion() << std::endl; }*/ } template< class TInputPixelType, class TOutputPixelType> typename itk::DiffusionKurtosisReconstructionImageFilter::KurtosisSnapshot itk::DiffusionKurtosisReconstructionImageFilter ::GetSnapshot(const itk::VariableLengthVector &input, GradientDirectionContainerType::Pointer gradients, float bvalue, KurtosisFitConfiguration kf_conf) { // initialize bvalues from reference value and the gradients provided on input this->SetReferenceBValue(bvalue); this->SetGradientDirections( gradients ); // call the other method return this->GetSnapshot( input, this->m_BValues, kf_conf ); } template< class TInputPixelType, class TOutputPixelType> typename itk::DiffusionKurtosisReconstructionImageFilter::KurtosisSnapshot itk::DiffusionKurtosisReconstructionImageFilter ::GetSnapshot(const itk::VariableLengthVector &input, vnl_vector bvalues, KurtosisFitConfiguration kf_conf) { // initialize vnl_vector initial_position; bool omit_bzero = kf_conf.omit_bzero; if( !omit_bzero ) { initial_position.set_size(2); initial_position[0] = this->m_InitialPosition[0]; initial_position[1] = this->m_InitialPosition[1]; } else { initial_position.set_size(3); initial_position = this->m_InitialPosition; } // fit FitSingleVoxel( input, bvalues, initial_position, kf_conf ); // assembly result snapshot KurtosisSnapshot result; result.m_D = initial_position[0]; result.m_K = initial_position[1]; if( omit_bzero ) { result.m_f = 1 - initial_position[2] / std::fmax(0.01, input.GetElement(0)); result.m_BzeroFit = initial_position[2]; } else result.m_f = 1; // assembly data vectors for fitting auto bvalueIter = bvalues.begin(); unsigned int unused_values = 0; while( bvalueIter != bvalues.end() ) { if( *bvalueIter < vnl_math::eps && omit_bzero ) { ++unused_values; } ++bvalueIter; } // initialize data vectors with the estimated size (after filtering) vnl_vector fit_measurements( input.Size() - unused_values, 0 ); vnl_vector fit_bvalues( input.Size() - unused_values, 0 ); // original data vnl_vector orig_measurements( input.Size(), 0 ); bvalueIter = bvalues.begin(); unsigned int running_index = 0; unsigned int skip_count = 0; + unsigned int bzero_count = 0; + result.m_Bzero = 0; while( bvalueIter != bvalues.end() ) { if( *bvalueIter < vnl_math::eps && omit_bzero ) { ++skip_count; } else { fit_measurements[ running_index - skip_count ] = input.GetElement(running_index); fit_bvalues[ running_index - skip_count ] = *bvalueIter; } + if( *bvalueIter < vnl_math::eps ) + { + result.m_Bzero += input.GetElement(running_index); + ++bzero_count; + } + orig_measurements[ running_index ] = input.GetElement(running_index); ++running_index; ++bvalueIter; } + result.m_Bzero /= bzero_count; result.fit_bvalues = fit_bvalues; result.fit_measurements = fit_measurements; result.bvalues = bvalues; result.measurements = orig_measurements; result.m_fittedBZero = omit_bzero; return result; } template< class TInputPixelType, class TOutputPixelType> void itk::DiffusionKurtosisReconstructionImageFilter ::SetGradientDirections(GradientDirectionContainerType::Pointer gradients) { if( this->m_ReferenceBValue < 0) { itkExceptionMacro( << "Specify reference b-value prior to setting the gradient directions." ); } if( gradients->Size() == 0 ) { itkExceptionMacro( << "Empty gradient directions container retrieved" ); } vnl_vector vnl_bvalues( gradients->Size(), 0 ); // initialize the b-values auto grIter = gradients->Begin(); unsigned int index = 0; while( grIter != gradients->End() ) { const double twonorm = grIter.Value().two_norm(); vnl_bvalues( index++ ) = this->m_ReferenceBValue * twonorm * twonorm; ++grIter; } this->m_BValues = vnl_bvalues; } template< class TInputPixelType, class TOutputPixelType> void itk::DiffusionKurtosisReconstructionImageFilter ::SetInitialSolution(const vnl_vector& x0 ) { assert( x0.size() == (2 + static_cast( this->m_OmitBZero )) ); this->m_InitialPosition = x0; } template< class TInputPixelType, class TOutputPixelType> void itk::DiffusionKurtosisReconstructionImageFilter ::ThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, ThreadIdType /*threadId*/) { typename OutputImageType::Pointer dImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); itk::ImageRegionIteratorWithIndex< OutputImageType > dImageIt(dImage, outputRegionForThread); dImageIt.GoToBegin(); typename OutputImageType::Pointer kImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(1)); itk::ImageRegionIteratorWithIndex< OutputImageType > kImageIt(kImage, outputRegionForThread); kImageIt.GoToBegin(); typedef itk::ImageRegionConstIteratorWithIndex< InputImageType > InputIteratorType; InputIteratorType inputIter( m_ProcessedInputImage, outputRegionForThread ); inputIter.GoToBegin(); typedef itk::ImageRegionConstIteratorWithIndex< MaskImageType > MaskIteratorType; MaskIteratorType maskIter( this->m_MaskImage, outputRegionForThread ); maskIter.GoToBegin(); KurtosisFitConfiguration fit_config; fit_config.omit_bzero = this->m_OmitBZero; fit_config.fit_scale = this->m_ScaleForFitting; fit_config.use_K_limits = this->m_UseKBounds; fit_config.K_limits = this->m_KurtosisBounds; vnl_vector initial_position; if( !this->m_OmitBZero ) { initial_position.set_size(2); initial_position[0] = this->m_InitialPosition[0]; initial_position[1] = this->m_InitialPosition[1]; } else { initial_position.set_size(3); initial_position = this->m_InitialPosition; } while( !inputIter.IsAtEnd() ) { // set (reset) each iteration vnl_vector result = initial_position; // fit single voxel (if inside mask ) if( maskIter.Get() > 0 ) { FitSingleVoxel( inputIter.Get(), this->m_BValues, result, fit_config ); } else { result.fill(0); } // regardless the fit type, the parameters are always in the first two position // of the results vector dImageIt.Set( result[0] ); kImageIt.Set( result[1] ); //std::cout << "[Kurtosis.Fit]" << inputIter.GetIndex() << " --> " << dImageIt.GetIndex() << " result: " << result << std::endl; ++maskIter; ++inputIter; ++dImageIt; ++kImageIt; } } #endif // guards diff --git a/Modules/DiffusionCore/Algorithms/Reconstruction/itkDiffusionKurtosisReconstructionImageFilter.h b/Modules/DiffusionCore/Algorithms/Reconstruction/itkDiffusionKurtosisReconstructionImageFilter.h index 734c040..8ea52d8 100644 --- a/Modules/DiffusionCore/Algorithms/Reconstruction/itkDiffusionKurtosisReconstructionImageFilter.h +++ b/Modules/DiffusionCore/Algorithms/Reconstruction/itkDiffusionKurtosisReconstructionImageFilter.h @@ -1,442 +1,443 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. 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 DIFFUSIONKURTOSISRECONSTRUCTIONIMAGEFILTER_H #define DIFFUSIONKURTOSISRECONSTRUCTIONIMAGEFILTER_H #include "itkImageToImageFilter.h" #include "itkVectorImage.h" #include "mitkDiffusionPropertyHelper.h" // vnl includes #include #include namespace itk { // Fitting routines /** @struct Kurtosis_fit_lsq_function @brief A base least-squares function for the diffusion kurtosis fit (non-IVIM) The basic function fits the signal to S_0 = S * exp [ -b * D + -b^2 * D^2 * K^2 ] */ struct kurtosis_fit_lsq_function : public vnl_least_squares_function { public: /** full lsq_function constructor */ kurtosis_fit_lsq_function( unsigned int num_params, unsigned int num_measurements, UseGradient g=no_gradient) : vnl_least_squares_function( num_params, num_measurements, g), m_use_bounds(false), m_use_logscale(false), m_skip_fit(false) {} /** simplified constructor for the 2-parameters fit */ kurtosis_fit_lsq_function( unsigned int number_measurements) : kurtosis_fit_lsq_function( 2, number_measurements, no_gradient ) {} /** Initialize the function by setting measurements and the corresponding b-values */ void initialize( vnl_vector< double > const& _meas, vnl_vector< double> const& _bvals ) { meas = _meas; if( m_use_logscale ) { for( unsigned int i=0; i< meas.size(); ++i) { // would produce NaN values, skip the fit // using the virtual function from the superclass (sets a boolean flag) if( meas[i] < vnl_math::eps ) { m_skip_fit = true; throw_failure(); continue; } meas[i] = log( meas[i] ); } } bvalues = _bvals; } /** use penalty terms on fitting to force the parameters stay within the default bounds */ void use_bounds() { m_use_bounds = true; // initialize bounds double upper_bounds[2] = {4e-3, 4 }; kurtosis_upper_bounds = vnl_vector(2, 2, upper_bounds); kurtosis_lower_bounds = vnl_vector(2, 0); } void set_fit_logscale( bool flag ) { this->m_use_logscale = flag; } void set_K_bounds( const vnl_vector_fixed k_bounds ) { // init this->use_bounds(); // override K bounds kurtosis_lower_bounds[1] = k_bounds[0]; kurtosis_upper_bounds[1] = k_bounds[1]; } virtual void f(const vnl_vector &x, vnl_vector &fx) { for ( unsigned int s=0; s < fx.size(); s++ ) { const double factor = ( meas[s] - M(x, s) ); fx[s] = factor * factor + penalty_term(x); } MITK_DEBUG("Fit.x_and_f") << x << " | " << fx; } protected: /** Formula for diffusion term, use for internal computations */ double Diff( double x1, double x2, double b) { const double quotient = -1. * b * x1 + b*b * x1 * x1 * x2 / 6; if( m_use_logscale ) return quotient; else return exp(quotient); } /** The fitting measurement function, has to be reimplemented in the classes */ virtual double M( vnl_vector const& x, unsigned int idx ) { const double bvalue = bvalues[idx]; double result = Diff( x[0], x[1], bvalue); if( m_use_logscale ) return meas[0] + result; else return meas[0] * result ; } /** Penalty term on D and K during fitting, make sure the vector that is passed in contains (D, K) in this ordering */ virtual double penalty_term( vnl_vector const& x) { double penalty = 0; // skip when turned off if( !m_use_bounds ) return penalty; // we have bounds for D and K only (the first two params ) for( unsigned int i=0; i< 2; i++) { // 5% penalty boundary // use exponential function to scale the penalty (max when x[i] == bounds ) double penalty_boundary = 0.02 * (kurtosis_upper_bounds[i] - kurtosis_lower_bounds[i]); if( x[i] < kurtosis_lower_bounds[i] + penalty_boundary ) { penalty += 1e6 * exp( -1 * ( x[i] - kurtosis_lower_bounds[i]) / penalty_boundary ); } else if ( x[i] > kurtosis_upper_bounds[i] - penalty_boundary ) { penalty += 1e6 * exp( -1 * ( kurtosis_upper_bounds[i] - x[i]) / penalty_boundary ); } } MITK_DEBUG("Fit.Orig.Penalty") << x << " || penalty: " << penalty; return penalty; } bool m_use_bounds; bool m_use_logscale; bool m_skip_fit; vnl_vector kurtosis_upper_bounds; vnl_vector kurtosis_lower_bounds; vnl_vector meas; vnl_vector bvalues; }; /** @struct kurtosis_fit_omit_unweighted @brief A fitting function handling the unweighted signal b_0 as a fitted parameter */ struct kurtosis_fit_omit_unweighted : public kurtosis_fit_lsq_function { public: /** simplified constructor for the 3-parameters fit */ kurtosis_fit_omit_unweighted( unsigned int number_measurements) : kurtosis_fit_lsq_function( 3, number_measurements, no_gradient ) {} protected: virtual double M(const vnl_vector &x, unsigned int idx) override { const double bvalue = bvalues[idx]; double result = Diff( x[0], x[1], bvalue); if( m_use_logscale ) return log( x[2] ) + result; else return x[2] * result ; } }; enum FitScale { STRAIGHT = 0, LOGARITHMIC }; struct KurtosisFitConfiguration { KurtosisFitConfiguration() : omit_bzero(false) , use_K_limits(false) , exclude_high_b(false) , b_upper_threshold(10e9) {} bool omit_bzero; FitScale fit_scale; bool use_K_limits; vnl_vector_fixed K_limits; bool exclude_high_b; double b_upper_threshold; }; /** @class DiffusionKurtosisReconstructionImageFilter @brief This filter provides the fit of the kurtosis (non-IVIM) signal to the data It has two main modes of operation, either as an image filter to compute the D and K maps, i.e. fitting the values to each voxel or a computation on a single voxel or a voxel group (with mask) can be triggered by @sa GetSnapshot, GetCurrentSnapshot methods. */ template< class TInputPixelType, class TOutputPixelType > class DiffusionKurtosisReconstructionImageFilter : public ImageToImageFilter< VectorImage< TInputPixelType, 3>, Image > { public: /** @struct KurtosisSnapshot @brief Struct describing a result (and the data) of a Kurtosis model fit */ struct KurtosisSnapshot { KurtosisSnapshot() : m_f(1), m_BzeroFit(1), m_D(0.001), m_K(0) {} // input data structures //vnl_vector filtered_measurements; vnl_vector bvalues; vnl_vector measurements; vnl_vector fit_bvalues; vnl_vector fit_measurements; vnl_vector weighted_image_indices; bool m_fittedBZero; // variables holding the fitted values double m_f; double m_BzeroFit; double m_D; double m_K; + double m_Bzero; // mean baseline signal }; //-- class typedefs typedef DiffusionKurtosisReconstructionImageFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ImageToImageFilter< VectorImage< TInputPixelType, 3>, Image< TOutputPixelType,3 > > Superclass; /** Method for creation through the object factory. */ itkFactorylessNewMacro(Self) itkCloneMacro(Self) /** Runtime information support. */ itkTypeMacro(DiffusionKurtosisReconstructionImageFilter, ImageToImageFilter) typedef TOutputPixelType OutputPixelType; typedef TInputPixelType InputPixelType; typedef typename Superclass::InputImageType InputImageType; typedef Image< OutputPixelType, 3 > OutputImageType; typedef itk::Image< short , 3> MaskImageType; typedef typename Superclass::OutputImageRegionType OutputImageRegionType; /** Holds each magnetic field gradient used to acquire one DWImage */ typedef mitk::DiffusionPropertyHelper::GradientDirectionsContainerType GradientDirectionContainerType; // vector image typedefs for regularized fit typedef itk::VectorImage VectorImageType; typedef itk::Image, 3> InitialFitImageType; //-- input (Set) methods /** Set the initial solution for fitting, make sure the length and the values correspond to the parameters * x0 = ( S_0, ADC_0, AKC_0 ) when also the S_0 is estimated * x0 = ( ADC_0, AKC_0 ) when the S_0 is used in fitting */ void SetInitialSolution(const vnl_vector& x0 ); /** Set whether the S_0 value is fitted or used in fitting */ void SetOmitUnweightedValue( bool flag) { this->m_OmitBZero = flag; } /** Trigger a single computation of the Kurtosis values from the given input vector and the bvalues and returns the result as a KurtosisSnapshot object */ KurtosisSnapshot GetSnapshot( const itk::VariableLengthVector< TInputPixelType > &input, vnl_vector bvalues, KurtosisFitConfiguration kf_conf); /** Trigger a single computation of the kurtosis values, first the bvalues vector is computed internally but then also stored into the returend snapshot */ KurtosisSnapshot GetSnapshot( const itk::VariableLengthVector< TInputPixelType > &input, GradientDirectionContainerType::Pointer, float bvalue, KurtosisFitConfiguration kf_conf); /** * Returns the value of the current data presented to the filter. If a mask is set, the voxels are first averaged before passed to the fitting procedure */ KurtosisSnapshot GetCurrentSnapshot(bool omit_bzero); /** Set the reference bvalue of the input DW image */ void SetReferenceBValue( double bvalue ) { this->m_ReferenceBValue = bvalue; } /** Set the gradient directions */ void SetGradientDirections( GradientDirectionContainerType::Pointer gradients ); /** Restrict map generation to an image region */ void SetMapOutputRegion( OutputImageRegionType region ) { m_MapOutputRegion = region; this->m_ApplyPriorSmoothing = true; } void SetImageMask( MaskImageType::Pointer mask ); /** Set smoothing sigma (default = 1.5 ), automatically enables smoothing prior to fitting */ void SetSmoothingSigma( double sigma ) { this->m_SmoothingSigma = sigma; this->m_ApplyPriorSmoothing = true; } /** Activate/Deactivate the gaussian smoothing applied to the input prior to fitting ( default = off ) */ void SetUseSmoothingPriorToFitting( bool flag) { this->m_ApplyPriorSmoothing = flag; } /** Set boundaries enforced by penalty terms in the fitting procedure */ void SetBoundariesForKurtosis( double lower, double upper ) { m_UseKBounds = true; m_KurtosisBounds[0] = lower; m_KurtosisBounds[1] = upper; } /** Exclude measurements associated with b-values higher than max_bvalue from fitting */ void SetMaximalBValueUsedForFitting( double max_bvalue ) { m_MaxFitBValue = max_bvalue; } /** Select the method used in fitting of the data STRAIHT - fit the exponential signal equation S / S_0 = exp [ ... ] LOGARITHMIC - fit the logarithmic signal equation ln( S / S_0 ) = [] */ void SetFittingScale( FitScale scale ) { m_ScaleForFitting = scale; } protected: DiffusionKurtosisReconstructionImageFilter(); virtual ~DiffusionKurtosisReconstructionImageFilter() {} void GenerateOutputInformation() override; void AfterThreadedGenerateData() override; void BeforeThreadedGenerateData() override; void ThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, ThreadIdType threadId) override; double m_ReferenceBValue; vnl_vector m_BValues; vnl_vector m_InitialPosition; bool m_OmitBZero; OutputImageRegionType m_MapOutputRegion; MaskImageType::Pointer m_MaskImage; typename InputImageType::Pointer m_ProcessedInputImage; bool m_ApplyPriorSmoothing; double m_SmoothingSigma; bool m_UseKBounds; vnl_vector_fixed m_KurtosisBounds; double m_MaxFitBValue; FitScale m_ScaleForFitting; private: }; } //end namespace itk #ifndef ITK_MANUAL_INSTANTIATION #include "itkDiffusionKurtosisReconstructionImageFilter.cxx" #endif #endif // DIFFUSIONKURTOSISRECONSTRUCTIONIMAGEFILTER_H diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.ivim/src/internal/QmitkIVIMView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.ivim/src/internal/QmitkIVIMView.cpp index ebf9220..4d39745 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.ivim/src/internal/QmitkIVIMView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.ivim/src/internal/QmitkIVIMView.cpp @@ -1,1098 +1,1063 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. 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 "QmitkIVIMView.h" // qt #include "qmessagebox.h" #include "qclipboard.h" // mitk #include "mitkImage.h" #include "mitkImageCast.h" #include "mitkLookupTable.h" #include "mitkLookupTableProperty.h" #include #include // itk #include "itkScalarImageToHistogramGenerator.h" #include "itkRegionOfInterestImageFilter.h" #include "itkImageRegionConstIteratorWithIndex.h" // itk/mitk #include "itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h" #include "itkRegularizedIVIMReconstructionFilter.h" #include "mitkImageCast.h" #include #include #include #include #include #include const std::string QmitkIVIMView::VIEW_ID = "org.mitk.views.ivim"; QmitkIVIMView::QmitkIVIMView() : QmitkAbstractView() , m_Controls( 0 ) , m_Active(false) , m_Visible(false) , m_HoldUpdate(false) { } QmitkIVIMView::~QmitkIVIMView() { } void QmitkIVIMView::CreateQtPartControl( QWidget *parent ) { // hold update untill all elements are set this->m_HoldUpdate = true; // build up qt view, unless already done if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkIVIMViewControls; m_Controls->setupUi( parent ); connect( m_Controls->m_ButtonStart, SIGNAL(clicked()), this, SLOT(FittIVIMStart()) ); connect( m_Controls->m_ButtonAutoThres, SIGNAL(clicked()), this, SLOT(AutoThreshold()) ); connect( m_Controls->m_MethodCombo, SIGNAL(currentIndexChanged(int)), this, SLOT(OnIvimFitChanged(int)) ); connect( m_Controls->m_DStarSlider, SIGNAL(valueChanged(int)), this, SLOT(DStarSlider(int)) ); connect( m_Controls->m_BThreshSlider, SIGNAL(valueChanged(int)), this, SLOT(BThreshSlider(int)) ); connect( m_Controls->m_S0ThreshSlider, SIGNAL(valueChanged(int)), this, SLOT(S0ThreshSlider(int)) ); connect( m_Controls->m_NumItSlider, SIGNAL(valueChanged(int)), this, SLOT(NumItsSlider(int)) ); connect( m_Controls->m_LambdaSlider, SIGNAL(valueChanged(int)), this, SLOT(LambdaSlider(int)) ); connect( m_Controls->m_CheckDStar, SIGNAL(clicked()), this, SLOT(Checkbox()) ); connect( m_Controls->m_CheckD, SIGNAL(clicked()), this, SLOT(Checkbox()) ); connect( m_Controls->m_Checkf, SIGNAL(clicked()), this, SLOT(Checkbox()) ); connect( m_Controls->m_CurveClipboard, SIGNAL(clicked()), this, SLOT(ClipboardCurveButtonClicked()) ); connect( m_Controls->m_ValuesClipboard, SIGNAL(clicked()), this, SLOT(ClipboardStatisticsButtonClicked()) ); connect( m_Controls->m_SavePlot, SIGNAL(clicked()), this, SLOT(SavePlotButtonClicked()) ); // connect all kurtosis actions to a recompute connect( m_Controls->m_KurtosisRangeWidget, SIGNAL( rangeChanged(double, double)), this, SLOT(OnKurtosisParamsChanged() ) ); connect( m_Controls->m_OmitBZeroCB, SIGNAL( stateChanged(int) ), this, SLOT( OnKurtosisParamsChanged() ) ); connect( m_Controls->m_KurtosisFitScale, SIGNAL( currentIndexChanged(int)), this, SLOT( OnKurtosisParamsChanged() ) ); connect( m_Controls->m_UseKurtosisBoundsCB, SIGNAL(clicked() ), this, SLOT( OnKurtosisParamsChanged() ) ); m_Controls->m_DwiBox->SetDataStorage(this->GetDataStorage()); mitk::NodePredicateIsDWI::Pointer isDwi = mitk::NodePredicateIsDWI::New(); m_Controls->m_DwiBox->SetPredicate( isDwi ); connect( (QObject*)(m_Controls->m_DwiBox), SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); m_Controls->m_MaskBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_MaskBox->SetZeroEntryText("--"); mitk::TNodePredicateDataType::Pointer isImagePredicate = mitk::TNodePredicateDataType::New(); mitk::NodePredicateProperty::Pointer isBinaryPredicate = mitk::NodePredicateProperty::New("binary", mitk::BoolProperty::New(true)); mitk::NodePredicateDimension::Pointer is3D = mitk::NodePredicateDimension::New(3); m_Controls->m_MaskBox->SetPredicate( mitk::NodePredicateAnd::New(isBinaryPredicate, mitk::NodePredicateAnd::New(isImagePredicate, is3D)) ); connect( (QObject*)(m_Controls->m_MaskBox), SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); connect( (QObject*)(m_Controls->m_ModelTabSelectionWidget), SIGNAL(currentChanged(int)), this, SLOT(OnModelTabChanged(int))); } QString dstar = QString::number(m_Controls->m_DStarSlider->value()/1000.0); m_Controls->m_DStarLabel->setText(dstar); QString bthresh = QString::number(m_Controls->m_BThreshSlider->value()*5.0); m_Controls->m_BThreshLabel->setText(bthresh); QString s0thresh = QString::number(m_Controls->m_S0ThreshSlider->value()*0.5); m_Controls->m_S0ThreshLabel->setText(s0thresh); QString numits = QString::number(m_Controls->m_NumItSlider->value()); m_Controls->m_NumItsLabel->setText(numits); QString lambda = QString::number(m_Controls->m_LambdaSlider->value()*.00001); m_Controls->m_LambdaLabel->setText(lambda); m_Controls->m_Warning->setVisible(false); OnIvimFitChanged(m_Controls->m_MethodCombo->currentIndex()); m_Controls->m_KurtosisRangeWidget->setSingleStep(0.1); m_Controls->m_KurtosisRangeWidget->setRange( 0.0, 10.0 ); m_Controls->m_KurtosisRangeWidget->setMaximumValue( 5.0 ); // LogScale not working yet, have to fix that first // m_Controls->m_KurtosisFitScale->setEnabled(false); //m_Controls->m_MaximalBValueWidget->setVisible( false ); // release update block after the UI-elements were all set this->m_HoldUpdate = false; QmitkIVIMView::InitChartIvim(); m_ListenerActive = false; if (this->GetRenderWindowPart()) { m_SliceChangeListener.RenderWindowPartActivated(this->GetRenderWindowPart()); connect(&m_SliceChangeListener, SIGNAL(SliceChanged()), this, SLOT(OnSliceChanged())); m_ListenerActive = true; } } void QmitkIVIMView::OnModelTabChanged(int tab) { if (tab==0) InitChartIvim(); else if (tab==1) InitChartKurtosis(); UpdateGui(); } void QmitkIVIMView::AddSecondFitPlot() { if (m_Controls->m_ChartWidget->GetDataElementByLabel("signal values (used for second fit)") == nullptr) { - std::map< double, double > init_data; + std::vector< std::pair > init_data; m_Controls->m_ChartWidget->AddData2D(init_data, "signal values (used for second fit)", QmitkChartWidget::ChartType::scatter); m_Controls->m_ChartWidget->SetColor("signal values (used for second fit)", "white"); m_Controls->m_ChartWidget->SetMarkerSymbol("signal values (used for second fit)", QmitkChartWidget::MarkerSymbol::x_thin); m_Controls->m_ChartWidget->Show(true); m_Controls->m_ChartWidget->SetShowDataPoints(false); m_Controls->m_ChartWidget->SetShowSubchart(false); } } void QmitkIVIMView::RemoveSecondFitPlot() { if (m_Controls->m_ChartWidget->GetDataElementByLabel("signal values (used for second fit)") != nullptr) { m_Controls->m_ChartWidget->RemoveData("signal values (used for second fit)"); m_Controls->m_ChartWidget->Show(true); m_Controls->m_ChartWidget->SetShowDataPoints(false); m_Controls->m_ChartWidget->SetShowSubchart(false); } } void QmitkIVIMView::InitChartIvim() { m_Controls->m_ChartWidget->Clear(); - std::map< double, double > init_data; + std::vector< std::pair > init_data; m_Controls->m_ChartWidget->AddData2D(init_data, "D-part of fitted model", QmitkChartWidget::ChartType::line); m_Controls->m_ChartWidget->AddData2D(init_data, "fitted model", QmitkChartWidget::ChartType::line); m_Controls->m_ChartWidget->AddData2D(init_data, "signal values", QmitkChartWidget::ChartType::scatter); m_Controls->m_ChartWidget->SetColor("fitted model", "red"); m_Controls->m_ChartWidget->SetColor("signal values", "white"); m_Controls->m_ChartWidget->SetColor("D-part of fitted model", "cyan"); m_Controls->m_ChartWidget->SetYAxisScale(QmitkChartWidget::AxisScale::log); m_Controls->m_ChartWidget->SetYAxisLabel("S/S0"); m_Controls->m_ChartWidget->SetXAxisLabel("b-value"); m_Controls->m_ChartWidget->SetLineStyle("fitted model", QmitkChartWidget::LineStyle::solid); m_Controls->m_ChartWidget->SetLineStyle("D-part of fitted model", QmitkChartWidget::LineStyle::dashed); m_Controls->m_ChartWidget->SetMarkerSymbol("signal values", QmitkChartWidget::MarkerSymbol::diamond); m_Controls->m_ChartWidget->Show(true); m_Controls->m_ChartWidget->SetShowDataPoints(false); m_Controls->m_ChartWidget->SetShowSubchart(false); } void QmitkIVIMView::InitChartKurtosis() { m_Controls->m_ChartWidget->Clear(); - std::map< double, double > init_data; + std::vector< std::pair > init_data; m_Controls->m_ChartWidget->AddData2D(init_data, "D-part of fitted model", QmitkChartWidget::ChartType::line); m_Controls->m_ChartWidget->AddData2D(init_data, "fitted model", QmitkChartWidget::ChartType::line); m_Controls->m_ChartWidget->AddData2D(init_data, "signal values", QmitkChartWidget::ChartType::scatter); m_Controls->m_ChartWidget->SetColor("fitted model", "red"); m_Controls->m_ChartWidget->SetColor("signal values", "white"); m_Controls->m_ChartWidget->SetColor("D-part of fitted model", "cyan"); m_Controls->m_ChartWidget->SetYAxisScale(QmitkChartWidget::AxisScale::log); m_Controls->m_ChartWidget->SetYAxisLabel("S"); m_Controls->m_ChartWidget->SetXAxisLabel("b-value"); m_Controls->m_ChartWidget->SetLineStyle("fitted model", QmitkChartWidget::LineStyle::solid); m_Controls->m_ChartWidget->SetLineStyle("D-part of fitted model", QmitkChartWidget::LineStyle::dashed); m_Controls->m_ChartWidget->SetMarkerSymbol("signal values", QmitkChartWidget::MarkerSymbol::diamond); m_Controls->m_ChartWidget->Show(true); m_Controls->m_ChartWidget->SetShowDataPoints(false); m_Controls->m_ChartWidget->SetShowSubchart(false); } void QmitkIVIMView::SetFocus() { m_Controls->m_ButtonAutoThres->setFocus(); } void QmitkIVIMView::Checkbox() { OnSliceChanged(); } void QmitkIVIMView::OnIvimFitChanged(int val) { switch(val) { case 0: m_Controls->m_DstarFrame->setVisible(false); m_Controls->m_NeglSiFrame->setVisible(true); m_Controls->m_NeglBframe->setVisible(false); m_Controls->m_IterationsFrame->setVisible(false); m_Controls->m_LambdaFrame->setVisible(false); break; case 1: m_Controls->m_DstarFrame->setVisible(true); m_Controls->m_NeglSiFrame->setVisible(true); m_Controls->m_NeglBframe->setVisible(false); m_Controls->m_IterationsFrame->setVisible(false); m_Controls->m_LambdaFrame->setVisible(false); break; case 2: m_Controls->m_DstarFrame->setVisible(false); m_Controls->m_NeglSiFrame->setVisible(true); m_Controls->m_NeglBframe->setVisible(true); m_Controls->m_IterationsFrame->setVisible(false); m_Controls->m_LambdaFrame->setVisible(false); break; case 3: m_Controls->m_DstarFrame->setVisible(false); m_Controls->m_NeglSiFrame->setVisible(true); m_Controls->m_NeglBframe->setVisible(true); m_Controls->m_IterationsFrame->setVisible(false); m_Controls->m_LambdaFrame->setVisible(false); break; case 4: m_Controls->m_DstarFrame->setVisible(false); m_Controls->m_NeglSiFrame->setVisible(false); m_Controls->m_NeglBframe->setVisible(false); m_Controls->m_IterationsFrame->setVisible(false); m_Controls->m_LambdaFrame->setVisible(false); break; } OnSliceChanged(); } void QmitkIVIMView::DStarSlider (int val) { QString sval = QString::number(val/1000.0); m_Controls->m_DStarLabel->setText(sval); OnSliceChanged(); } void QmitkIVIMView::BThreshSlider (int val) { QString sval = QString::number(val*5.0); m_Controls->m_BThreshLabel->setText(sval); OnSliceChanged(); } void QmitkIVIMView::S0ThreshSlider (int val) { QString sval = QString::number(val*0.5); m_Controls->m_S0ThreshLabel->setText(sval); OnSliceChanged(); } void QmitkIVIMView::NumItsSlider (int val) { QString sval = QString::number(val); m_Controls->m_NumItsLabel->setText(sval); OnSliceChanged(); } void QmitkIVIMView::LambdaSlider (int val) { QString sval = QString::number(val*.00001); m_Controls->m_LambdaLabel->setText(sval); OnSliceChanged(); } void QmitkIVIMView::UpdateGui() { m_Controls->m_FittedParamsLabel->setText(""); if (m_Controls->m_DwiBox->GetSelectedNode().IsNotNull()) { m_Controls->m_ChartWidget->setVisible(true); m_HoldUpdate = false; } else { m_Controls->m_ChartWidget->setVisible(false); } m_Controls->m_ButtonStart->setEnabled( m_Controls->m_DwiBox->GetSelectedNode().IsNotNull() ); m_Controls->m_ButtonAutoThres->setEnabled( m_Controls->m_DwiBox->GetSelectedNode().IsNotNull() ); m_Controls->m_ControlsFrame->setEnabled( m_Controls->m_DwiBox->GetSelectedNode().IsNotNull() ); m_Controls->m_BottomControlsFrame->setEnabled( m_Controls->m_DwiBox->GetSelectedNode().IsNotNull() ); OnSliceChanged(); } void QmitkIVIMView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList& ) { // UpdateGui(); } void QmitkIVIMView::AutoThreshold() { if (m_Controls->m_DwiBox->GetSelectedNode().IsNull()) { // Nothing selected. Inform the user and return QMessageBox::information( nullptr, "Template", "Please load and select a diffusion image before starting image processing."); return; } mitk::Image* dimg = dynamic_cast(m_Controls->m_DwiBox->GetSelectedNode()->GetData()); if (!dimg) { // Nothing selected. Inform the user and return QMessageBox::information( nullptr, "Template", "No valid diffusion image was found."); return; } // find bzero index int index = -1; auto directions = mitk::DiffusionPropertyHelper::GetGradientContainer(dimg); for(DirContainerType::ConstIterator it = directions->Begin(); it != directions->End(); ++it) { index++; GradientDirectionType g = it.Value(); if(g[0] == 0 && g[1] == 0 && g[2] == 0 ) break; } VecImgType::Pointer vecimg = VecImgType::New(); mitk::CastToItkImage(dimg, vecimg); int vecLength = vecimg->GetVectorLength(); index = index > vecLength-1 ? vecLength-1 : index; MITK_INFO << "Performing Histogram Analysis on Channel" << index; typedef itk::Image ImgType; ImgType::Pointer img = ImgType::New(); mitk::CastToItkImage(dimg, img); itk::ImageRegionIterator itw (img, img->GetLargestPossibleRegion() ); itw.GoToBegin(); itk::ImageRegionConstIterator itr (vecimg, vecimg->GetLargestPossibleRegion() ); itr.GoToBegin(); while(!itr.IsAtEnd()) { itw.Set(itr.Get().GetElement(index)); ++itr; ++itw; } typedef itk::Statistics::ScalarImageToHistogramGenerator< ImgType > HistogramGeneratorType; typedef HistogramGeneratorType::HistogramType HistogramType; HistogramGeneratorType::Pointer histogramGenerator = HistogramGeneratorType::New(); histogramGenerator->SetInput( img ); histogramGenerator->SetMarginalScale( 10 ); // Defines y-margin width of histogram histogramGenerator->SetNumberOfBins( 100 ); // CT range [-1024, +2048] --> bin size 4 values histogramGenerator->SetHistogramMin( dimg->GetStatistics()->GetScalarValueMin() ); histogramGenerator->SetHistogramMax( dimg->GetStatistics()->GetScalarValueMax() * .5 ); histogramGenerator->Compute(); HistogramType::ConstIterator iter = histogramGenerator->GetOutput()->Begin(); float maxFreq = 0; float maxValue = 0; while ( iter != histogramGenerator->GetOutput()->End() ) { if(iter.GetFrequency() > maxFreq) { maxFreq = iter.GetFrequency(); maxValue = iter.GetMeasurementVector()[0]; } ++iter; } maxValue *= 2; int sliderPos = maxValue * 2; m_Controls->m_S0ThreshSlider->setValue(sliderPos); S0ThreshSlider(sliderPos); } void QmitkIVIMView::FittIVIMStart() { if (m_Controls->m_DwiBox->GetSelectedNode().IsNull()) { QMessageBox::information( nullptr, "Template", "No valid diffusion-weighted image selected."); return; } mitk::Image* img = dynamic_cast(m_Controls->m_DwiBox->GetSelectedNode()->GetData()); VecImgType::Pointer vecimg = VecImgType::New(); mitk::CastToItkImage(img, vecimg); OutImgType::IndexType dummy; if( m_Controls->m_ModelTabSelectionWidget->currentIndex() ) { // KURTOSIS KurtosisFilterType::Pointer filter = KurtosisFilterType::New(); filter->SetInput(vecimg); filter->SetReferenceBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(img)); filter->SetGradientDirections(mitk::DiffusionPropertyHelper::GetGradientContainer(img)); filter->SetSmoothingSigma( m_Controls->m_SigmaSpinBox->value() ); if( m_Controls->m_UseKurtosisBoundsCB->isChecked() ) filter->SetBoundariesForKurtosis( m_Controls->m_KurtosisRangeWidget->minimumValue(), m_Controls->m_KurtosisRangeWidget->maximumValue() ); filter->SetFittingScale( static_cast(m_Controls->m_KurtosisFitScale->currentIndex() ) ); if( m_Controls->m_MaskBox->GetSelectedNode().IsNotNull() ) { mitk::Image::Pointer maskImg = dynamic_cast(m_Controls->m_MaskBox->GetSelectedNode()->GetData()); typedef itk::Image MaskImgType; MaskImgType::Pointer maskItk; CastToItkImage( maskImg, maskItk ); filter->SetImageMask( maskItk ); } filter->Update(); mitk::LookupTable::Pointer kurt_map_lut = mitk::LookupTable::New(); kurt_map_lut->SetType( mitk::LookupTable::JET ); mitk::LookupTableProperty::Pointer kurt_lut_prop = mitk::LookupTableProperty::New(); kurt_lut_prop->SetLookupTable( kurt_map_lut ); mitk::Image::Pointer dimage = mitk::Image::New(); dimage->InitializeByItk( filter->GetOutput(0) ); dimage->SetVolume( filter->GetOutput(0)->GetBufferPointer()); mitk::Image::Pointer kimage = mitk::Image::New(); kimage->InitializeByItk( filter->GetOutput(1) ); kimage->SetVolume( filter->GetOutput(1)->GetBufferPointer()); QString new_dname = "Kurtosis_DMap"; new_dname.append("_Method-"+m_Controls->m_KurtosisFitScale->currentText()); QString new_kname = "Kurtosis_KMap"; new_kname.append("_Method-"+m_Controls->m_KurtosisFitScale->currentText()); if( m_Controls->m_CheckKurtD->isChecked() ) { mitk::DataNode::Pointer dnode = mitk::DataNode::New(); dnode->SetData( dimage ); dnode->SetName(new_dname.toLatin1()); dnode->SetProperty("LookupTable", kurt_lut_prop ); GetDataStorage()->Add(dnode, m_Controls->m_DwiBox->GetSelectedNode()); } if( m_Controls->m_CheckKurtK->isChecked() ) { mitk::DataNode::Pointer knode = mitk::DataNode::New(); knode->SetData( kimage ); knode->SetName(new_kname.toLatin1()); knode->SetProperty("LookupTable", kurt_lut_prop ); GetDataStorage()->Add(knode, m_Controls->m_DwiBox->GetSelectedNode()); } } else { FittIVIM(vecimg, mitk::DiffusionPropertyHelper::GetGradientContainer(img), mitk::DiffusionPropertyHelper::GetReferenceBValue(img), true, dummy); OutputToDatastorage(m_Controls->m_DwiBox->GetSelectedNode()); } } void QmitkIVIMView::OnKurtosisParamsChanged() { OnSliceChanged(); } void QmitkIVIMView::OnSliceChanged() { if(m_HoldUpdate || !m_Visible) return; m_Controls->m_Warning->setVisible(false); if(m_Controls->m_DwiBox->GetSelectedNode().IsNull()) return; mitk::Image::Pointer diffusionImg = dynamic_cast(m_Controls->m_DwiBox->GetSelectedNode()->GetData()); mitk::Image::Pointer maskImg = nullptr; if (m_Controls->m_MaskBox->GetSelectedNode().IsNotNull()) maskImg = dynamic_cast(m_Controls->m_MaskBox->GetSelectedNode()->GetData()); if (!this->GetRenderWindowPart()) return; if (!m_ListenerActive) { m_SliceChangeListener.RenderWindowPartActivated(this->GetRenderWindowPart()); connect(&m_SliceChangeListener, SIGNAL(SliceChanged()), this, SLOT(OnSliceChanged())); m_ListenerActive = true; } VecImgType::Pointer vecimg = VecImgType::New(); mitk::CastToItkImage(diffusionImg, vecimg); VecImgType::Pointer roiImage = VecImgType::New(); if(maskImg.IsNull()) { int roisize = 0; if(m_Controls->m_MethodCombo->currentIndex() == 4) roisize = 5; mitk::Point3D pos = this->GetRenderWindowPart()->GetSelectedPosition(); VecImgType::IndexType crosspos; diffusionImg->GetGeometry()->WorldToIndex(pos, crosspos); if (!vecimg->GetLargestPossibleRegion().IsInside(crosspos)) { m_Controls->m_Warning->setText(QString("Crosshair position not inside of selected diffusion weighted image. Reinit needed!")); m_Controls->m_Warning->setVisible(true); return; } else m_Controls->m_Warning->setVisible(false); VecImgType::IndexType index; index[0] = crosspos[0] - roisize; index[0] = index[0] < 0 ? 0 : index[0]; index[1] = crosspos[1] - roisize; index[1] = index[1] < 0 ? 0 : index[1]; index[2] = crosspos[2] - roisize; index[2] = index[2] < 0 ? 0 : index[2]; VecImgType::SizeType size; size[0] = roisize*2+1; size[1] = roisize*2+1; size[2] = roisize*2+1; VecImgType::SizeType maxSize = vecimg->GetLargestPossibleRegion().GetSize(); size[0] = index[0]+size[0] > maxSize[0] ? maxSize[0]-index[0] : size[0]; size[1] = index[1]+size[1] > maxSize[1] ? maxSize[1]-index[1] : size[1]; size[2] = index[2]+size[2] > maxSize[2] ? maxSize[2]-index[2] : size[2]; VecImgType::RegionType region; region.SetSize( size ); region.SetIndex( index ); vecimg->SetRequestedRegion( region ); VecImgType::IndexType newstart; newstart.Fill(0); VecImgType::RegionType newregion; newregion.SetSize( size ); newregion.SetIndex( newstart ); roiImage->CopyInformation( vecimg ); roiImage->SetRegions( newregion ); roiImage->SetOrigin( pos ); roiImage->Allocate(); roiImage->SetPixel(newstart, vecimg->GetPixel(index)); if( m_Controls->m_ModelTabSelectionWidget->currentIndex() ) { FitKurtosis(roiImage, mitk::DiffusionPropertyHelper::GetGradientContainer(diffusionImg), mitk::DiffusionPropertyHelper::GetReferenceBValue(diffusionImg), newstart); } else { FittIVIM(roiImage, mitk::DiffusionPropertyHelper::GetGradientContainer(diffusionImg), mitk::DiffusionPropertyHelper::GetReferenceBValue(diffusionImg), false, crosspos); } } else { typedef itk::Image MaskImgType; MaskImgType::Pointer maskItk; CastToItkImage( maskImg, maskItk ); mitk::Point3D pos; pos[0] = 0; pos[1] = 0; pos[2] = 0; VecImgType::IndexType index; index[0] = 0; index[1] = 0; index[2] = 0; VecImgType::SizeType size; size[0] = 1; size[1] = 1; size[2] = 1; VecImgType::RegionType region; region.SetSize( size ); region.SetIndex( index ); vecimg->SetRequestedRegion( region ); // iterators over output and input itk::ImageRegionConstIteratorWithIndex vecit(vecimg, vecimg->GetLargestPossibleRegion()); itk::VariableLengthVector avg(vecimg->GetVectorLength()); avg.Fill(0); float numPixels = 0; while ( ! vecit.IsAtEnd() ) { VecImgType::PointType point; vecimg->TransformIndexToPhysicalPoint(vecit.GetIndex(), point); MaskImgType::IndexType index; maskItk->TransformPhysicalPointToIndex(point, index); if(maskItk->GetPixel(index) != 0) { avg += vecit.Get(); numPixels += 1.0; } // update iterators ++vecit; } avg /= numPixels; m_Controls->m_Warning->setText(QString("Averaging ")+QString::number((int)numPixels)+QString(" voxels!")); m_Controls->m_Warning->setVisible(true); roiImage->CopyInformation( vecimg ); roiImage->SetRegions( region ); roiImage->SetOrigin( pos ); roiImage->Allocate(); roiImage->SetPixel(index, avg); if( m_Controls->m_ModelTabSelectionWidget->currentIndex() ) { FitKurtosis(roiImage, mitk::DiffusionPropertyHelper::GetGradientContainer(diffusionImg), mitk::DiffusionPropertyHelper::GetReferenceBValue(diffusionImg), index); } else { FittIVIM(roiImage, mitk::DiffusionPropertyHelper::GetGradientContainer(diffusionImg), mitk::DiffusionPropertyHelper::GetReferenceBValue(diffusionImg), false, index); } // do not update until selection changed, the values will remain the same as long as the mask is selected! m_HoldUpdate = true; } vecimg->SetRegions( vecimg->GetLargestPossibleRegion() ); } bool QmitkIVIMView::FitKurtosis( itk::VectorImage *vecimg, DirContainerType *dirs, float bval, OutImgType::IndexType &crosspos ) { KurtosisFilterType::Pointer filter = KurtosisFilterType::New(); itk::KurtosisFitConfiguration fit_config; fit_config.omit_bzero = m_Controls->m_OmitBZeroCB->isChecked(); if( m_Controls->m_UseKurtosisBoundsCB->isChecked() ) { fit_config.use_K_limits = true; vnl_vector_fixed k_limits; k_limits[0] = m_Controls->m_KurtosisRangeWidget->minimumValue(); k_limits[1] = m_Controls->m_KurtosisRangeWidget->maximumValue(); fit_config.K_limits = k_limits; } fit_config.fit_scale = static_cast(m_Controls->m_KurtosisFitScale->currentIndex() ); m_KurtosisSnap = filter->GetSnapshot( vecimg->GetPixel( crosspos ), dirs, bval, fit_config ); QString param_label_text("K=%2, D=%1"); param_label_text = param_label_text.arg( m_KurtosisSnap.m_K, 4); param_label_text = param_label_text.arg( m_KurtosisSnap.m_D, 4); m_Controls->m_FittedParamsLabel->setText(param_label_text); const double maxb = m_KurtosisSnap.bvalues.max_value(); - double S0 = m_KurtosisSnap.measurements[0]; + double S0 = m_KurtosisSnap.m_Bzero; + if (m_KurtosisSnap.m_fittedBZero) + S0 = m_KurtosisSnap.m_BzeroFit; - std::map d_line; - d_line[0] = S0; - d_line[maxb] = d_line[0]*exp(-maxb * m_KurtosisSnap.m_D); + std::vector< std::pair > d_line; + d_line.push_back(std::pair(0, S0)); + d_line.push_back(std::pair(maxb, S0*exp(-maxb * m_KurtosisSnap.m_D))); m_Controls->m_ChartWidget->UpdateData2D(d_line, "D-part of fitted model"); const unsigned int num_samples = 50; - std::map y; + std::vector< std::pair > y; for( unsigned int i=0; i<=num_samples; ++i) { double b = (((1.0)*i)/(1.0*num_samples))*maxb; - y[b] = S0 * exp( -b * m_KurtosisSnap.m_D + b*b * m_KurtosisSnap.m_D * m_KurtosisSnap.m_D * m_KurtosisSnap.m_K / 6.0 ); + y.push_back(std::pair(b, S0 * exp( -b * m_KurtosisSnap.m_D + b*b * m_KurtosisSnap.m_D * m_KurtosisSnap.m_D * m_KurtosisSnap.m_K / 6.0 ))); } m_Controls->m_ChartWidget->UpdateData2D(y, "fitted model"); - std::map y_meas; + std::vector< std::pair > y_meas; for (unsigned int i=0; im_OmitBZeroCB->isChecked() || m_KurtosisSnap.bvalues[i] > 0.01) - y_meas[m_KurtosisSnap.bvalues[i]] = m_KurtosisSnap.measurements[i]; + y_meas.push_back(std::pair(m_KurtosisSnap.bvalues[i], m_KurtosisSnap.measurements[i])); + } m_Controls->m_ChartWidget->UpdateData2D(y_meas, "signal values"); return true; } bool QmitkIVIMView::FittIVIM(itk::VectorImage* vecimg, DirContainerType* dirs, float bval, bool multivoxel, OutImgType::IndexType &crosspos) { IVIMFilterType::Pointer filter = IVIMFilterType::New(); filter->SetInput(vecimg); filter->SetGradientDirections(dirs); filter->SetBValue(bval); switch(m_Controls->m_MethodCombo->currentIndex()) { case 0: filter->SetMethod(IVIMFilterType::IVIM_FIT_ALL); filter->SetS0Thres(m_Controls->m_S0ThreshLabel->text().toDouble()); break; case 1: filter->SetMethod(IVIMFilterType::IVIM_DSTAR_FIX); filter->SetDStar(m_Controls->m_DStarLabel->text().toDouble()); filter->SetS0Thres(m_Controls->m_S0ThreshLabel->text().toDouble()); break; case 2: filter->SetMethod(IVIMFilterType::IVIM_D_THEN_DSTAR); filter->SetBThres(m_Controls->m_BThreshLabel->text().toDouble()); filter->SetS0Thres(m_Controls->m_S0ThreshLabel->text().toDouble()); filter->SetFitDStar(m_Controls->m_CheckDStar->isChecked()); break; case 3: filter->SetMethod(IVIMFilterType::IVIM_LINEAR_D_THEN_F); filter->SetBThres(m_Controls->m_BThreshLabel->text().toDouble()); filter->SetS0Thres(m_Controls->m_S0ThreshLabel->text().toDouble()); filter->SetFitDStar(m_Controls->m_CheckDStar->isChecked()); break; case 4: filter->SetMethod(IVIMFilterType::IVIM_REGULARIZED); filter->SetBThres(m_Controls->m_BThreshLabel->text().toDouble()); filter->SetS0Thres(m_Controls->m_S0ThreshLabel->text().toDouble()); filter->SetNumberIterations(m_Controls->m_NumItsLabel->text().toInt()); filter->SetLambda(m_Controls->m_LambdaLabel->text().toDouble()); filter->SetFitDStar(m_Controls->m_CheckDStar->isChecked()); break; } if(!multivoxel) { filter->SetFitDStar(true); } filter->SetNumberOfThreads(1); filter->SetVerbose(false); filter->SetCrossPosition(crosspos); try{ filter->Update(); m_IvimSnap = filter->GetSnapshot(); m_DStarMap = filter->GetOutput(2); m_DMap = filter->GetOutput(1); m_fMap = filter->GetOutput(); QString param_label_text("f=%1, D=%2, D*=%3"); param_label_text = param_label_text.arg(m_IvimSnap.currentF,4); param_label_text = param_label_text.arg(m_IvimSnap.currentD,4); param_label_text = param_label_text.arg(m_IvimSnap.currentDStar,4); m_Controls->m_FittedParamsLabel->setText(param_label_text); double maxb = m_IvimSnap.bvalues.max_value(); - std::map d_line; - d_line[0] = 1-m_IvimSnap.currentFunceiled; - d_line[maxb] = d_line[0]*exp(-maxb * m_IvimSnap.currentD); + std::vector< std::pair > d_line; + d_line.push_back(std::pair(0, 1-m_IvimSnap.currentFunceiled)); + d_line.push_back(std::pair(maxb, d_line[0].second*exp(-maxb * m_IvimSnap.currentD))); m_Controls->m_ChartWidget->UpdateData2D(d_line, "D-part of fitted model"); if(m_IvimSnap.currentDStar != 0) { - std::map y; + std::vector< std::pair > y; int nsampling = 50; double f = 1-m_IvimSnap.currentFunceiled; for(int i=0; i<=nsampling; i++) { double x = (((1.0)*i)/(1.0*nsampling))*maxb; - y[x] = f*exp(- x * m_IvimSnap.currentD) + (1-f)*exp(- x * (m_IvimSnap.currentD+m_IvimSnap.currentDStar)); + y.push_back(std::pair(x, f*exp(- x * m_IvimSnap.currentD) + (1-f)*exp(- x * (m_IvimSnap.currentD+m_IvimSnap.currentDStar)))); } m_Controls->m_ChartWidget->UpdateData2D(y, "fitted model"); - std::map y_meas; + std::vector< std::pair > y_meas; for (unsigned int i=0; i(m_IvimSnap.bvals1[i], m_IvimSnap.meas1[i])); m_Controls->m_ChartWidget->UpdateData2D(y_meas, "signal values"); if(!m_IvimSnap.high_indices.empty()) { - MITK_INFO << "m_Snap.high_indices found"; AddSecondFitPlot(); - std::map additonal_meas; + std::vector< std::pair > additonal_meas; for (int i=0; i(m_IvimSnap.bvals2[i], m_IvimSnap.meas2[i])); m_Controls->m_ChartWidget->UpdateData2D(additonal_meas, "signal values (used for second fit)"); } else { - MITK_INFO << "NO m_Snap.high_indices found"; RemoveSecondFitPlot(); } } } catch (itk::ExceptionObject &ex) { MITK_INFO << ex ; m_Controls->m_Warning->setText(QString("IVIM fit not possible: ")+ex.GetDescription()); m_Controls->m_Warning->setVisible(true); return false; } return true; } void QmitkIVIMView::OutputToDatastorage(mitk::DataNode::Pointer node) { mitk::LookupTable::Pointer lut = mitk::LookupTable::New(); lut->SetType( mitk::LookupTable::JET ); mitk::LookupTableProperty::Pointer lut_prop = mitk::LookupTableProperty::New(); lut_prop->SetLookupTable( lut ); if(m_Controls->m_CheckDStar->isChecked()) { mitk::Image::Pointer dstarimage = mitk::Image::New(); dstarimage->InitializeByItk(m_DStarMap.GetPointer()); dstarimage->SetVolume(m_DStarMap->GetBufferPointer()); QString newname2 = ""; newname2 = newname2.append("IVIM_DStarMap_Method-%1").arg(m_Controls->m_MethodCombo->currentText()); mitk::DataNode::Pointer node2=mitk::DataNode::New(); node2->SetData( dstarimage ); node2->SetName(newname2.toLatin1()); node2->SetProperty("LookupTable", lut_prop ); GetDataStorage()->Add(node2, node); } if(m_Controls->m_CheckD->isChecked()) { mitk::Image::Pointer dimage = mitk::Image::New(); dimage->InitializeByItk(m_DMap.GetPointer()); dimage->SetVolume(m_DMap->GetBufferPointer()); QString newname1 = ""; newname1 = newname1.append("IVIM_DMap_Method-%1").arg(m_Controls->m_MethodCombo->currentText()); mitk::DataNode::Pointer node1=mitk::DataNode::New(); node1->SetData( dimage ); node1->SetName(newname1.toLatin1()); node1->SetProperty("LookupTable", lut_prop ); GetDataStorage()->Add(node1, node); } if(m_Controls->m_Checkf->isChecked()) { mitk::Image::Pointer image = mitk::Image::New(); image->InitializeByItk(m_fMap.GetPointer()); image->SetVolume(m_fMap->GetBufferPointer()); QString newname0 = ""; newname0 = newname0.append("IVIM_fMap_Method-%1").arg(m_Controls->m_MethodCombo->currentText()); mitk::DataNode::Pointer node3=mitk::DataNode::New(); node3->SetData( image ); node3->SetName(newname0.toLatin1()); node3->SetProperty("LookupTable", lut_prop ); GetDataStorage()->Add(node3, node); } this->GetRenderWindowPart()->RequestUpdate(); } void QmitkIVIMView::ClipboardCurveButtonClicked() { // Kurtosis if ( m_Controls->m_ModelTabSelectionWidget->currentIndex() ) { std::stringstream ss; QString clipboard("Measurement Points\n"); ss << m_KurtosisSnap.bvalues << "\n" << m_KurtosisSnap.measurements << "\n\n"; ss << "Fitted Values ( D K [b_0] ) \n" << m_KurtosisSnap.m_D << " " << m_KurtosisSnap.m_K; + if( m_KurtosisSnap.m_fittedBZero ) ss << " " << m_KurtosisSnap.m_BzeroFit; ss << "\n\n"; clipboard.append( QString( ss.str().c_str() )); ss.str( std::string() ); ss.clear(); - QApplication::clipboard()->setText( - clipboard, QClipboard::Clipboard ); + QApplication::clipboard()->setText(clipboard, QClipboard::Clipboard ); } else { - - QString clipboard("Measurement Points\n"); - for ( unsigned int i=0; isetText( - clipboard, QClipboard::Clipboard ); + QApplication::clipboard()->setText(clipboard, QClipboard::Clipboard ); } } void QmitkIVIMView::SavePlotButtonClicked() { m_Controls->m_ChartWidget->SavePlotAsImage(); } void QmitkIVIMView::ClipboardStatisticsButtonClicked() { // Kurtosis if ( m_Controls->m_ModelTabSelectionWidget->currentIndex() ) { QString clipboard( "D \t K \n" ); clipboard = clipboard.append( "%L1 \t %L2" ) .arg( m_KurtosisSnap.m_D, 0, 'f', 10 ) .arg( m_KurtosisSnap.m_K, 0, 'f', 10 ) ; - QApplication::clipboard()->setText( - clipboard, QClipboard::Clipboard ); + QApplication::clipboard()->setText(clipboard, QClipboard::Clipboard ); } else { QString clipboard( "f \t D \t D* \n" ); clipboard = clipboard.append( "%L1 \t %L2 \t %L3" ) .arg( m_IvimSnap.currentF, 0, 'f', 10 ) .arg( m_IvimSnap.currentD, 0, 'f', 10 ) .arg( m_IvimSnap.currentDStar, 0, 'f', 10 ) ; - QApplication::clipboard()->setText( - clipboard, QClipboard::Clipboard ); + QApplication::clipboard()->setText(clipboard, QClipboard::Clipboard ); } } void QmitkIVIMView::Activated() { m_Active = true; } void QmitkIVIMView::Deactivated() { m_Active = false; } void QmitkIVIMView::Visible() { m_Visible = true; if (this->GetRenderWindowPart() and !m_ListenerActive) { m_SliceChangeListener.RenderWindowPartActivated(this->GetRenderWindowPart()); connect(&m_SliceChangeListener, SIGNAL(SliceChanged()), this, SLOT(OnSliceChanged())); m_ListenerActive = true; } } void QmitkIVIMView::Hidden() { m_Visible = false; } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.ivim/src/internal/QmitkIVIMViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging.ivim/src/internal/QmitkIVIMViewControls.ui index 4c9b10e..0ec0c50 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.ivim/src/internal/QmitkIVIMViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.ivim/src/internal/QmitkIVIMViewControls.ui @@ -1,1218 +1,1218 @@ QmitkIVIMViewControls 0 0 - 423 - 648 + 424 + 689 0 0 QmitkTemplate QCommandLinkButton:disabled { border: none; } QGroupBox { background-color: transparent; } 9 9 9 9 9 QFrame::NoFrame QFrame::Raised 0 0 0 0 warning display Qt::RichText true 0 0 16 16 QFrame::NoFrame QFrame::Raised 0 0 0 0 0 0 16 16 QFrame::NoFrame QFrame::Raised 0 0 0 0 0 0 16 16 QFrame::NoFrame QFrame::Raised 0 0 0 0 0 0 0 0 16777215 16777215 Generate Parameter Maps QFrame::StyledPanel QFrame::Raised 0 0 0 0 0 0 0 IVIM Parameters 9 QFrame::NoFrame QFrame::Raised 0 0 0 0 0 51 16777215 200 Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter 80 16777215 D* 100 60 Qt::Horizontal QFrame::NoFrame QFrame::Plain 0 0 0 0 0 130 0 130 16777215 Signal threshold: 100 0 Qt::Horizontal 0 0 QFrame::NoFrame QFrame::Raised 0 0 0 0 0 0 0 30 16777215 TextLabel Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter 15 16777215 Calculate threshold from histogram * QFrame::NoFrame QFrame::Raised 0 0 0 0 0 130 0 130 16777215 Ignore b< (first fit) 250 34 Qt::Horizontal 51 16777215 46.5 Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter QFrame::NoFrame QFrame::Raised 0 0 0 0 0 80 16777215 #iterations: 100 10 Qt::Horizontal 30 16777215 TextLabel Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter QFrame::NoFrame QFrame::Raised 0 0 0 0 0 80 16777215 lambda: 1000 10 Qt::Horizontal 0 0 QFrame::NoFrame QFrame::Raised 0 0 0 0 0 0 0 30 16777215 TextLabel Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter 15 16777215 Calculate threshold from histogram * QFrame::NoFrame QFrame::Raised 0 0 0 0 Fit method: 2 3 Param. Fit Fit D & f with fixed D* value Fit D & f (high b), then fit D* Linearly fit D & f (high b), then fit D* Regularized fit QFrame::NoFrame QFrame::Raised 0 0 0 0 80 0 Output Maps f true D false D* false Kurtosis QFrame::StyledPanel QFrame::Raised 2 2 2 2 2 Smoothing sigma Select Fit Type Omit b=0 Measurements - true + false 80 0 Output Maps Force the fitting of K to remain within the given boundaries Boundaries for K Select if the data is fitted directly (straight) or the logarithmic equation is used Straight Fit Logarithmic Fit 2 QLayout::SetMaximumSize D false K true Signa for gaussian smoothing applied prior to map computation 0.000000000000000 5.000000000000000 0.100000000000000 On Input Data 6 6 6 6 Optional ROI image ROI: DWI to analyze Raw DWI: Qt::Vertical 20 40 Copy the signal values in the current voxel to the clipboard. Save plot as image QFrame::NoFrame QFrame::Raised 0 0 0 0 QFrame::NoFrame QFrame::Raised 0 0 0 0 Copy the signal values in the current voxel to the clipboard. Signal values to clipboard Copy the parameters of the curve fit in the current voxel to the clipboard. Fit-parameters to clipboard Intra Voxel Incoherent Motion Estimation 6 9 6 6 Qt::AlignCenter QmitkDataStorageComboBox QComboBox
QmitkDataStorageComboBox.h
 QmitkDataStorageComboBoxWithSelectNone QComboBox
QmitkDataStorageComboBoxWithSelectNone.h
 ctkRangeWidget QWidget
ctkRangeWidget.h
1 QmitkChartWidget QWidget
QmitkChartWidget.h
1