diff --git a/Modules/DiffusionCore/Algorithms/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.cpp b/Modules/DiffusionCore/Algorithms/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.cpp index 027d3f1..658a943 100644 --- a/Modules/DiffusionCore/Algorithms/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.cpp +++ b/Modules/DiffusionCore/Algorithms/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.cpp @@ -1,781 +1,712 @@ /*=================================================================== 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 __itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter_cpp #define __itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter_cpp #include "itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h" #include "itkImageRegionConstIterator.h" #include "itkImageRegionConstIteratorWithIndex.h" #include "itkImageRegionIterator.h" #include "vnl/vnl_matrix.h" #include "vnl/algo/vnl_symmetric_eigensystem.h" +#include #include "itkRegularizedIVIMReconstructionFilter.h" #include #define IVIM_FOO -100000 namespace itk { template< class TIn, class TOut> DiffusionIntravoxelIncoherentMotionReconstructionImageFilter ::DiffusionIntravoxelIncoherentMotionReconstructionImageFilter() : m_GradientDirectionContainer(nullptr), m_Method(IVIM_DSTAR_FIX), m_FitDStar(true), m_Verbose(false) { this->SetNumberOfRequiredInputs( 1 ); this->SetNumberOfRequiredOutputs( 3 ); typename OutputImageType::Pointer outputPtr1 = OutputImageType::New(); this->SetNthOutput(0, outputPtr1.GetPointer()); typename OutputImageType::Pointer outputPtr2 = OutputImageType::New(); this->SetNthOutput(1, outputPtr2.GetPointer()); typename OutputImageType::Pointer outputPtr3 = OutputImageType::New(); this->SetNthOutput(2, outputPtr3.GetPointer()); } template< class TIn, class TOut> void DiffusionIntravoxelIncoherentMotionReconstructionImageFilter ::BeforeThreadedGenerateData() { // Input must be an itk::VectorImage. - std::string gradientImageClassName( - this->ProcessObject::GetInput(0)->GetNameOfClass()); + std::string gradientImageClassName(this->ProcessObject::GetInput(0)->GetNameOfClass()); if ( strcmp(gradientImageClassName.c_str(),"VectorImage") != 0 ) { itkExceptionMacro( << "There is only one Gradient image. I expect that to be a VectorImage. " << "But its of type: " << gradientImageClassName ); } // Compute the indicies of the baseline images and gradient images // If no b=0 mm/s² gradients ar found, the next lowest b-value is used. GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin(); double minNorm = itk::NumericTraits::max(); while( gdcit != this->m_GradientDirectionContainer->End() ) { double norm = gdcit.Value().one_norm(); if (normm_GradientDirectionContainer->Begin(); while( gdcit != this->m_GradientDirectionContainer->End() ) { if(gdcit.Value().one_norm() <= minNorm) - { m_Snap.baselineind.push_back(gdcit.Index()); - } else - { m_Snap.gradientind.push_back(gdcit.Index()); - double twonorm = gdcit.Value().two_norm(); - m_Snap.bvals.push_back( m_BValue*twonorm*twonorm ); - } ++gdcit; } if (m_Snap.gradientind.size()==0) itkExceptionMacro("Only one b-value supplied. At least two needed for IVIM fit."); // check sind die grad und base gleichlang? alle grad gerade und base ungerade? dann iterierende aufnahme!! m_Snap.iterated_sequence = false; if(m_Snap.baselineind.size() == m_Snap.gradientind.size()) { int size = m_Snap.baselineind.size(); int sum_b = 0, sum_g = 0; for(int i=0; iGetElement(m_Snap.gradientind[i]).two_norm(); + m_Snap.bvalues[i] = m_BValue*twonorm*twonorm; } if(m_Verbose) { std::cout << "ref bval: " << m_BValue << "; b-values: "; - for(int i=0; im_BThres) - { m_Snap.high_indices.push_back(i); - } - } } - m_Snap.Nhigh = m_Snap.high_indices.size(); - m_Snap.high_bvalues.set_size(m_Snap.Nhigh); - m_Snap.high_meas.set_size(m_Snap.Nhigh); - for(int i=0; i MeasAndBvals DiffusionIntravoxelIncoherentMotionReconstructionImageFilter ::ApplyS0Threshold(vnl_vector &meas, vnl_vector &bvals) { std::vector newmeas; std::vector newbvals; int N = meas.size(); for(int i=0; i void DiffusionIntravoxelIncoherentMotionReconstructionImageFilter -::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, - ThreadIdType ) +::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, ThreadIdType ) { - typename OutputImageType::Pointer outputImage = - static_cast< OutputImageType * >(this->ProcessObject::GetPrimaryOutput()); + typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetPrimaryOutput()); ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread); oit.GoToBegin(); - typename OutputImageType::Pointer dImage = - static_cast< OutputImageType * >(this->ProcessObject::GetOutput(1)); + typename OutputImageType::Pointer dImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(1)); ImageRegionIterator< OutputImageType > oit1(dImage, outputRegionForThread); oit1.GoToBegin(); - typename OutputImageType::Pointer dstarImage = - static_cast< OutputImageType * >(this->ProcessObject::GetOutput(2)); + typename OutputImageType::Pointer dstarImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(2)); ImageRegionIterator< OutputImageType > oit2(dstarImage, outputRegionForThread); oit2.GoToBegin(); typedef ImageRegionConstIterator< InputImageType > InputIteratorType; typedef typename InputImageType::PixelType InputVectorType; typename InputImageType::Pointer inputImagePointer = nullptr; // Would have liked a dynamic_cast here, but seems SGI doesn't like it // The enum will DiffusionIntravoxelIncoherentMotionReconstructionImageFilterensure that an inappropriate cast is not done inputImagePointer = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) ); InputIteratorType iit(inputImagePointer, outputRegionForThread ); iit.GoToBegin(); // init internal vector image for regularized fit m_InternalVectorImage = VectorImageType::New(); m_InternalVectorImage->SetSpacing( inputImagePointer->GetSpacing() ); // Set the image spacing m_InternalVectorImage->SetOrigin( inputImagePointer->GetOrigin() ); // Set the image origin m_InternalVectorImage->SetDirection( inputImagePointer->GetDirection() ); // Set the image direction m_InternalVectorImage->SetRegions( inputImagePointer->GetLargestPossibleRegion() ); m_InitialFitImage = InitialFitImageType::New(); m_InitialFitImage->SetSpacing( inputImagePointer->GetSpacing() ); // Set the image spacing m_InitialFitImage->SetOrigin( inputImagePointer->GetOrigin() ); // Set the image origin m_InitialFitImage->SetDirection( inputImagePointer->GetDirection() ); // Set the image direction m_InitialFitImage->SetRegions( inputImagePointer->GetLargestPossibleRegion() ); if(m_Method == IVIM_REGULARIZED) { - m_InternalVectorImage->SetVectorLength(m_Snap.Nhigh); + m_InternalVectorImage->SetVectorLength(m_Snap.num_high); m_InternalVectorImage->Allocate(); - VectorImageType::PixelType varvec(m_Snap.Nhigh); - for(int i=0; iFillBuffer(varvec); m_InitialFitImage->Allocate(); InitialFitImageType::PixelType vec; vec[0] = 0.5; vec[1] = 0.01; vec[2]=0.001; m_InitialFitImage->FillBuffer(vec); } typedef itk::ImageRegionIterator VectorIteratorType; VectorIteratorType vecit(m_InternalVectorImage, outputRegionForThread ); vecit.GoToBegin(); typedef itk::ImageRegionIterator InitIteratorType; InitIteratorType initit(m_InitialFitImage, outputRegionForThread ); initit.GoToBegin(); while( !iit.IsAtEnd() ) { InputVectorType measvec = iit.Get(); typename NumericTraits::AccumulateType b0 = NumericTraits::Zero; - m_Snap.meas.set_size(m_Snap.N); - m_Snap.allmeas.set_size(m_Snap.N); + m_Snap.meas_for_threshold.set_size(m_Snap.num_weighted); + m_Snap.allmeas.set_size(m_Snap.num_weighted); if(!m_Snap.iterated_sequence) { // Average the baseline image pixels for(unsigned int i = 0; i < m_Snap.baselineind.size(); ++i) - { b0 += measvec[m_Snap.baselineind[i]]; - } + if(m_Snap.baselineind.size()) b0 /= m_Snap.baselineind.size(); // measurement vector - for(int i = 0; i < m_Snap.N; ++i) + for(int i = 0; i < m_Snap.num_weighted; ++i) { m_Snap.allmeas[i] = measvec[m_Snap.gradientind[i]] / (b0+.0001); if(measvec[m_Snap.gradientind[i]] > m_S0Thres) - { - m_Snap.meas[i] = measvec[m_Snap.gradientind[i]] / (b0+.0001); - } + m_Snap.meas_for_threshold[i] = measvec[m_Snap.gradientind[i]] / (b0+.0001); else - { - m_Snap.meas[i] = IVIM_FOO; - } + m_Snap.meas_for_threshold[i] = IVIM_FOO; } } else { // measurement vector - for(int i = 0; i < m_Snap.N; ++i) + for(int i = 0; i < m_Snap.num_weighted; ++i) { b0 = measvec[m_Snap.baselineind[i]]; - m_Snap.allmeas[i] = measvec[m_Snap.gradientind[i]] / (b0+.0001); if(measvec[m_Snap.gradientind[i]] > m_S0Thres) - { - m_Snap.meas[i] = measvec[m_Snap.gradientind[i]] / (b0+.0001); - } + m_Snap.meas_for_threshold[i] = measvec[m_Snap.gradientind[i]] / (b0+.0001); else - { - m_Snap.meas[i] = IVIM_FOO; - } + m_Snap.meas_for_threshold[i] = IVIM_FOO; } } m_Snap.currentF = 0; m_Snap.currentD = 0; m_Snap.currentDStar = 0; switch(m_Method) { case IVIM_D_THEN_DSTAR: { - - for(int i=0; i x_donly(2); x_donly[0] = 0.001; x_donly[1] = 0.1; // f 0.1 Dstar 0.01 D 0.001 vnl_levenberg_marquardt lm_donly(f_donly); lm_donly.set_f_tolerance(0.0001); lm_donly.minimize(x_donly); m_Snap.currentD = x_donly[0]; m_Snap.currentF = x_donly[1]; - if(m_FitDStar) { - MeasAndBvals input2 = ApplyS0Threshold(m_Snap.meas, m_Snap.bvalues); + MeasAndBvals input2 = ApplyS0Threshold(m_Snap.meas_for_threshold, m_Snap.bvalues); m_Snap.bvals2 = input2.bvals; m_Snap.meas2 = input2.meas; if (input2.N < 2) break; IVIM_dstar_only f_dstar_only(input2.N,m_Snap.currentD,m_Snap.currentF); f_dstar_only.set_bvalues(input2.bvals); f_dstar_only.set_measurements(input2.meas); vnl_vector< double > x_dstar_only(1); vnl_vector< double > fx_dstar_only(input2.N); double opt = 1111111111111111.0; int opt_idx = -1; int num_its = 100; double min_val = .001; double max_val = .15; for(int i=0; i x_fixd(2); - // x_fixd[0] = 0.1; - // x_fixd[1] = 0.01; - // // f 0.1 Dstar 0.01 D 0.001 - - // vnl_levenberg_marquardt lm_fixd(f_fixd); - // lm_fixd.set_f_tolerance(0.0001); - // lm_fixd.minimize(x_fixd); - - // m_Snap.currentF = x_fixd[0]; - // m_Snap.currentDStar = x_fixd[1]; } break; } case IVIM_DSTAR_FIX: { - MeasAndBvals input = ApplyS0Threshold(m_Snap.meas, m_Snap.bvalues); + MeasAndBvals input = ApplyS0Threshold(m_Snap.meas_for_threshold, m_Snap.bvalues); m_Snap.bvals1 = input.bvals; m_Snap.meas1 = input.meas; if (input.N < 2) break; - IVIM_fixdstar f_fixdstar(input.N,m_DStar); + IVIM_fixdstar f_fixdstar(input.N, m_DStar); f_fixdstar.set_bvalues(input.bvals); f_fixdstar.set_measurements(input.meas); vnl_vector< double > x(2); x[0] = 0.1; x[1] = 0.001; // f 0.1 Dstar 0.01 D 0.001 vnl_levenberg_marquardt lm(f_fixdstar); lm.set_f_tolerance(0.0001); lm.minimize(x); m_Snap.currentF = x[0]; m_Snap.currentD = x[1]; m_Snap.currentDStar = m_DStar; break; } case IVIM_FIT_ALL: { - MeasAndBvals input = ApplyS0Threshold(m_Snap.meas, m_Snap.bvalues); + MeasAndBvals input = ApplyS0Threshold(m_Snap.meas_for_threshold, m_Snap.bvalues); m_Snap.bvals1 = input.bvals; m_Snap.meas1 = input.meas; if (input.N < 3) break; IVIM_3param f_3param(input.N); f_3param.set_bvalues(input.bvals); f_3param.set_measurements(input.meas); vnl_vector< double > x(3); x[0] = 0.1; x[1] = 0.001; x[2] = 0.01; // f 0.1 Dstar 0.01 D 0.001 vnl_levenberg_marquardt lm(f_3param); lm.set_f_tolerance(0.0001); lm.minimize(x); m_Snap.currentF = x[0]; m_Snap.currentD = x[1]; m_Snap.currentDStar = x[2]; break; } case IVIM_LINEAR_D_THEN_F: { - - // // neglect zero-measurements - // bool zero = false; - // for(int i=0; i X(input.N,2); for(int i=0; i XX = X.transpose() * X; vnl_symmetric_eigensystem eigs(XX); vnl_vector eig; if(eigs.get_eigenvalue(0) > eigs.get_eigenvalue(1)) eig = eigs.get_eigenvector(0); else eig = eigs.get_eigenvector(1); m_Snap.currentF = 1 - exp( meas_m - bval_m*(eig(1)/eig(0)) ); m_Snap.currentD = -eig(1)/eig(0); if(m_FitDStar) { - MeasAndBvals input2 = ApplyS0Threshold(m_Snap.meas, m_Snap.bvalues); + MeasAndBvals input2 = ApplyS0Threshold(m_Snap.meas_for_threshold, m_Snap.bvalues); m_Snap.bvals2 = input2.bvals; m_Snap.meas2 = input2.meas; if (input2.N < 2) break; IVIM_dstar_only f_dstar_only(input2.N,m_Snap.currentD,m_Snap.currentF); f_dstar_only.set_bvalues(input2.bvals); f_dstar_only.set_measurements(input2.meas); vnl_vector< double > x_dstar_only(1); vnl_vector< double > fx_dstar_only(input2.N); double opt = 1111111111111111.0; int opt_idx = -1; int num_its = 100; double min_val = .001; double max_val = .15; for(int i=0; i " << DStar; // x_dstar_only[0] = 0.01; // // f 0.1 Dstar 0.01 D 0.001 // vnl_levenberg_marquardt lm_dstar_only(f_dstar_only); // lm_dstar_only.set_f_tolerance(0.0001); // lm_dstar_only.minimize(x_dstar_only); // DStar = x_dstar_only[0]; break; } case IVIM_REGULARIZED: { //m_Snap.high_meas, m_Snap.high_bvalues; - for(int i=0; i x_donly(2); x_donly[0] = 0.001; x_donly[1] = 0.1; if(input.N >= 2) { IVIM_d_and_f f_donly(input.N); f_donly.set_bvalues(input.bvals); f_donly.set_measurements(input.meas); //MITK_INFO << "initial fit N=" << input.N << ", min-b = " << input.bvals[0] << ", max-b = " << input.bvals[input.N-1]; vnl_levenberg_marquardt lm_donly(f_donly); lm_donly.set_f_tolerance(0.0001); lm_donly.minimize(x_donly); } typename InitialFitImageType::PixelType initvec; initvec[0] = x_donly[1]; initvec[1] = x_donly[0]; initit.Set(initvec); //MITK_INFO << "Init vox " << initit.GetIndex() << " with " << initvec[0] << "; " << initvec[1]; ++initit; int N = m_Snap.high_meas.size(); typename VectorImageType::PixelType vec(N); for(int i=0; i void DiffusionIntravoxelIncoherentMotionReconstructionImageFilter ::AfterThreadedGenerateData() { if(m_Method == IVIM_REGULARIZED) { typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetPrimaryOutput()); ImageRegionIterator< OutputImageType > oit0(outputImage, outputImage->GetLargestPossibleRegion()); oit0.GoToBegin(); typename OutputImageType::Pointer dImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(1)); ImageRegionIterator< OutputImageType > oit1(dImage, dImage->GetLargestPossibleRegion()); oit1.GoToBegin(); typename OutputImageType::Pointer dstarImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(2)); ImageRegionIterator< OutputImageType > oit2(dstarImage, dstarImage->GetLargestPossibleRegion()); oit2.GoToBegin(); typedef itk::RegularizedIVIMReconstructionFilter RegFitType; RegFitType::Pointer filter = RegFitType::New(); filter->SetInput(m_InitialFitImage); filter->SetReferenceImage(m_InternalVectorImage); filter->SetBValues(m_Snap.high_bvalues); filter->SetNumberIterations(m_NumberIterations); filter->SetNumberOfThreads(1); filter->SetLambda(m_Lambda); filter->Update(); typename RegFitType::OutputImageType::Pointer outimg = filter->GetOutput(); ImageRegionConstIterator< RegFitType::OutputImageType > iit(outimg, outimg->GetLargestPossibleRegion()); iit.GoToBegin(); while( !iit.IsAtEnd() ) { double f = iit.Get()[0]; IVIM_CEIL( f, 0.0, 1.0 ); - oit0.Set( myround(f * 100.0) ); - oit1.Set( myround(iit.Get()[1] * 10000.0) ); - oit2.Set( myround(iit.Get()[2] * 1000.0) ); + oit0.Set( myround(f) ); + oit1.Set( myround(iit.Get()[1]) ); + oit2.Set( myround(iit.Get()[2]) ); if(!m_Verbose) { // report the middle voxel - if( iit.GetIndex()[0] == m_CrossPosition[0] - && iit.GetIndex()[1] == m_CrossPosition[1] - && iit.GetIndex()[2] == m_CrossPosition[2] ) + if( iit.GetIndex()[0] == static_cast(outimg->GetLargestPossibleRegion().GetSize(0)-1)/2 + && iit.GetIndex()[1] == static_cast(outimg->GetLargestPossibleRegion().GetSize(2)-1)/2 + && iit.GetIndex()[2] == static_cast(outimg->GetLargestPossibleRegion().GetSize(1)-1)/2 ) { m_Snap.currentF = f; m_Snap.currentD = iit.Get()[1]; m_Snap.currentDStar = iit.Get()[2]; m_Snap.allmeas = m_tmp_allmeas; MITK_INFO << "setting " << f << ";" << iit.Get()[1] << ";" << iit.Get()[2]; } } ++oit0; ++oit1; ++oit2; ++iit; } } } template< class TIn, class TOut> double DiffusionIntravoxelIncoherentMotionReconstructionImageFilter ::myround(double number) { return number < 0.0 ? ceil(number - 0.5) : floor(number + 0.5); } template< class TIn, class TOut> void DiffusionIntravoxelIncoherentMotionReconstructionImageFilter ::SetGradientDirections( GradientDirectionContainerType *gradientDirection ) { this->m_GradientDirectionContainer = gradientDirection; this->m_NumberOfGradientDirections = gradientDirection->Size(); } template< class TIn, class TOut> void DiffusionIntravoxelIncoherentMotionReconstructionImageFilter ::PrintSelf(std::ostream& os, Indent indent) const { Superclass::PrintSelf(os,indent); if ( m_GradientDirectionContainer ) { os << indent << "GradientDirectionContainer: " << m_GradientDirectionContainer << std::endl; } else { os << indent << "GradientDirectionContainer: (Gradient directions not set)" << std::endl; } } } #endif // __itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter_cpp diff --git a/Modules/DiffusionCore/Algorithms/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h b/Modules/DiffusionCore/Algorithms/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h index d97f4fd..33e2165 100644 --- a/Modules/DiffusionCore/Algorithms/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h +++ b/Modules/DiffusionCore/Algorithms/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h @@ -1,381 +1,403 @@ /*=================================================================== 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 __itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter_h #define __itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter_h #include "itkImageToImageFilter.h" //#include "vnl/vnl_matrix.h" #include "vnl/vnl_vector_fixed.h" #include "vnl/vnl_matrix.h" #include "vnl/algo/vnl_svd.h" #include "itkVectorContainer.h" #include "itkVectorImage.h" //#include "QuadProg.h" #include "itkVectorImage.h" #include "vnl/vnl_least_squares_function.h" +#include "vnl/vnl_cost_function.h" #include "vnl/algo/vnl_levenberg_marquardt.h" #include "vnl/vnl_math.h" #define IVIM_CEIL(val,u,o) (val) = \ ( (val) < (u) ) ? ( (u) ) : ( ( (val)>(o) ) ? ( (o) ) : ( (val) ) ); namespace itk{ /** baseclass for IVIM fitting algorithms */ struct IVIM_base { void set_measurements(const vnl_vector& x) { measurements.set_size(x.size()); measurements.copy_in(x.data_block()); } void set_bvalues(const vnl_vector& x) { bvalues.set_size(x.size()); bvalues.copy_in(x.data_block()); } vnl_vector measurements; vnl_vector bvalues; int N; }; /** Fitt all three parameters */ struct IVIM_3param : public IVIM_base, vnl_least_squares_function { IVIM_3param(unsigned int number_of_measurements) : vnl_least_squares_function(3, number_of_measurements, no_gradient) { N = get_number_of_residuals(); } void f(const vnl_vector& x, vnl_vector& fx) override { double ef = x[0]; double D = x[1]; double Dstar = x[2]; - for(int s=0; s& x, vnl_vector& fx) override { double ef = x[0]; double D = x[1]; for(int s=0; s& x, vnl_vector& fx) override { double D = x[0]; double f = x[1]; for(int s=0; s& x, vnl_vector& fx) override { double ef = x[0]; double Dstar = x[1]; for(int s=0; s& x, vnl_vector& fx) override { double Dstar = x[0]; for(int s=0; s meas; vnl_vector bvals; int N; }; /** \class DiffusionIntravoxelIncoherentMotionReconstructionImageFilter */ template< class TInputPixelType, class TOutputPixelType> class DiffusionIntravoxelIncoherentMotionReconstructionImageFilter : public ImageToImageFilter< VectorImage< TInputPixelType, 3 >, Image< TOutputPixelType, 3 > > { public: struct IVIMSnapshot { double currentF; double currentFunceiled; double currentD; double currentDStar; bool iterated_sequence; // wether each measurement has its own b0-acqu. std::vector baselineind; // baseline image indicies - - int N; // total number of measurements std::vector gradientind; // gradient image indicies - std::vector bvals; // b-values != 0 - vnl_vector bvalues; // copy of bvalues != 0 - vnl_vector meas; // all measurements, thresholded blanked out + vnl_vector bvalues; // bvalues != 0 + + int num_weighted; // total number of measurements + + vnl_vector meas_for_threshold; // all measurements, thresholded blanked out vnl_vector allmeas; // all measurements - int Nhigh; // number of used measurements + int num_high; // number of used measurements std::vector high_indices; // indices of used measurements vnl_vector high_bvalues; // bvals of used measurements vnl_vector high_meas; // used measurements + // fit 1 vnl_vector meas1; vnl_vector bvals1; + // fit 2 vnl_vector meas2; vnl_vector bvals2; - }; enum IVIM_Method { IVIM_FIT_ALL, IVIM_DSTAR_FIX, IVIM_D_THEN_DSTAR, IVIM_LINEAR_D_THEN_F, IVIM_REGULARIZED }; typedef DiffusionIntravoxelIncoherentMotionReconstructionImageFilter 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(DiffusionIntravoxelIncoherentMotionReconstructionImageFilter, ImageToImageFilter); typedef TOutputPixelType OutputPixelType; typedef TInputPixelType InputPixelType; /** Reference image data, This image is aquired in the absence * of a diffusion sensitizing field gradient */ typedef typename Superclass::InputImageType InputImageType; typedef Image< OutputPixelType, 3 > OutputImageType; typedef typename Superclass::OutputImageRegionType OutputImageRegionType; /** Holds each magnetic field gradient used to acquire one DWImage */ typedef vnl_vector_fixed< double, 3 > GradientDirectionType; /** Container to hold gradient directions of the 'n' DW measurements */ typedef VectorContainer< unsigned int, GradientDirectionType > GradientDirectionContainerType; // vector image typedefs for regularized fit typedef itk::VectorImage VectorImageType; typedef itk::Image, 3> InitialFitImageType; /** set method to add gradient directions and its corresponding * image. The image here is a VectorImage. The user is expected to pass the * gradient directions in a container. The ith element of the container * corresponds to the gradient direction of the ith component image the * VectorImage. For the baseline image, a vector of all zeros * should be set.*/ void SetGradientDirections( GradientDirectionContainerType * ); void SetBValue(double bval){m_BValue = bval;} void SetBThres(double bval){m_BThres = bval;} void SetS0Thres(double val){m_S0Thres = val;} void SetDStar(double dstar){m_DStar = dstar;} void SetFitDStar(bool fit){m_FitDStar = fit;} void SetVerbose(bool verbose){m_Verbose = verbose;} void SetNumberIterations(int num){m_NumberIterations = num;} void SetLambda(double lambda){m_Lambda = lambda;} void SetCrossPosition(typename InputImageType::IndexType crosspos){this->m_CrossPosition = crosspos;} void SetMethod(IVIM_Method method){m_Method = method;} IVIMSnapshot GetSnapshot(){return m_Snap;} /** Return the gradient direction. idx is 0 based */ virtual GradientDirectionType GetGradientDirection( unsigned int idx) const { if( idx >= m_NumberOfGradientDirections ) { itkExceptionMacro( << "Gradient direction " << idx << "does not exist" ); } return m_GradientDirectionContainer->ElementAt( idx+1 ); } protected: DiffusionIntravoxelIncoherentMotionReconstructionImageFilter(); ~DiffusionIntravoxelIncoherentMotionReconstructionImageFilter() {}; void PrintSelf(std::ostream& os, Indent indent) const; void BeforeThreadedGenerateData(); void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType); void AfterThreadedGenerateData(); MeasAndBvals ApplyS0Threshold(vnl_vector &meas, vnl_vector &bvals); private: double myround(double number); /** container to hold gradient directions */ GradientDirectionContainerType::Pointer m_GradientDirectionContainer; /** Number of gradient measurements */ unsigned int m_NumberOfGradientDirections; double m_BValue; double m_BThres; double m_S0Thres; double m_DStar; IVIM_Method m_Method; typename OutputImageType::Pointer m_DMap; typename OutputImageType::Pointer m_DStarMap; bool m_FitDStar; IVIMSnapshot m_Snap; bool m_Verbose; typename VectorImageType::Pointer m_InternalVectorImage; typename InitialFitImageType::Pointer m_InitialFitImage; int m_NumberIterations; // for total variation vnl_vector m_tmp_allmeas; double m_Lambda; typename InputImageType::IndexType m_CrossPosition; }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.cpp" #endif #endif //__itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter_h diff --git a/Modules/MriSimulation/Algorithms/itkFitFibersToImageFilter.cpp b/Modules/MriSimulation/Algorithms/itkFitFibersToImageFilter.cpp index 05d6bd8..e0c67cd 100644 --- a/Modules/MriSimulation/Algorithms/itkFitFibersToImageFilter.cpp +++ b/Modules/MriSimulation/Algorithms/itkFitFibersToImageFilter.cpp @@ -1,993 +1,965 @@ #include "itkFitFibersToImageFilter.h" #include #include namespace itk{ FitFibersToImageFilter::FitFibersToImageFilter() : m_PeakImage(nullptr) , m_DiffImage(nullptr) , m_ScalarImage(nullptr) , m_MaskImage(nullptr) , m_FitIndividualFibers(true) , m_GradientTolerance(1e-5) , m_Lambda(1.0) , m_MaxIterations(20) , m_Coverage(0) , m_Overshoot(0) , m_RMSE(0.0) , m_FilterOutliers(false) , m_MeanWeight(1.0) , m_MedianWeight(1.0) , m_MinWeight(1.0) , m_MaxWeight(1.0) , m_Verbose(true) , m_NumUnknowns(0) , m_NumResiduals(0) , m_NumCoveredDirections(0) , m_SignalModel(nullptr) , sz_x(0) , sz_y(0) , sz_z(0) , m_MeanTractDensity(0) , m_MeanSignal(0) , fiber_count(0) , m_Regularization(VnlCostFunction::REGU::VOXEL_VARIANCE) { this->SetNumberOfRequiredOutputs(3); } FitFibersToImageFilter::~FitFibersToImageFilter() { } void FitFibersToImageFilter::CreateDiffSystem() { sz_x = m_DiffImage->GetLargestPossibleRegion().GetSize(0); sz_y = m_DiffImage->GetLargestPossibleRegion().GetSize(1); sz_z = m_DiffImage->GetLargestPossibleRegion().GetSize(2); dim_four_size = m_DiffImage->GetVectorLength(); int num_voxels = sz_x*sz_y*sz_z; if (m_MaskImage.IsNotNull() && (m_MaskImage->GetLargestPossibleRegion().GetSize()[0]!=sz_x || m_MaskImage->GetLargestPossibleRegion().GetSize()[1]!=sz_y || m_MaskImage->GetLargestPossibleRegion().GetSize()[2]!=sz_z)) mitkThrow() << "Mask and peak image must have same size!"; auto spacing = m_DiffImage->GetSpacing(); m_NumResiduals = num_voxels * dim_four_size; MITK_INFO << "Num. unknowns: " << m_NumUnknowns; MITK_INFO << "Num. residuals: " << m_NumResiduals; MITK_INFO << "Creating system ..."; A.set_size(m_NumResiduals, m_NumUnknowns); b.set_size(m_NumResiduals); b.fill(0.0); m_MeanTractDensity = 0; m_MeanSignal = 0; m_NumCoveredDirections = 0; fiber_count = 0; vnl_vector voxel_indicator; voxel_indicator.set_size(sz_x*sz_y*sz_z); voxel_indicator.fill(0); unsigned int numFibers = 0; for (unsigned int bundle=0; bundleGetNumFibers(); boost::progress_display disp(numFibers); m_GroupSizes.clear(); for (unsigned int bundle=0; bundle polydata = m_Tractograms.at(bundle)->GetFiberPolyData(); m_GroupSizes.push_back(m_Tractograms.at(bundle)->GetNumFibers()); for (unsigned int i=0; iGetNumFibers(); ++i) { ++disp; vtkCell* cell = polydata->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (numPoints<2) MITK_INFO << "FIBER WITH ONLY ONE POINT ENCOUNTERED!"; for (int j=0; jGetPoint(j)); itk::Index<3> startIndex; itk::ContinuousIndex startIndexCont; m_DiffImage->TransformPhysicalPointToIndex(startVertex, startIndex); m_DiffImage->TransformPhysicalPointToContinuousIndex(startVertex, startIndexCont); PointType3 endVertex = mitk::imv::GetItkPoint(points->GetPoint(j+1)); itk::Index<3> endIndex; itk::ContinuousIndex endIndexCont; m_DiffImage->TransformPhysicalPointToIndex(endVertex, endIndex); m_DiffImage->TransformPhysicalPointToContinuousIndex(endVertex, endIndexCont); mitk::DiffusionSignalModel<>::GradientType fiber_dir; fiber_dir[0] = endVertex[0]-startVertex[0]; fiber_dir[1] = endVertex[1]-startVertex[1]; fiber_dir[2] = endVertex[2]-startVertex[2]; fiber_dir.Normalize(); std::vector< std::pair< itk::Index<3>, double > > segments = mitk::imv::IntersectImage(spacing, startIndex, endIndex, startIndexCont, endIndexCont); for (std::pair< itk::Index<3>, double > seg : segments) { if (!m_DiffImage->GetLargestPossibleRegion().IsInside(seg.first) || (m_MaskImage.IsNotNull() && m_MaskImage->GetPixel(seg.first)==0)) continue; int x = seg.first[0]; int y = seg.first[1]; int z = seg.first[2]; mitk::DiffusionSignalModel<>::PixelType simulated_pixel = m_SignalModel->SimulateMeasurement(fiber_dir)*seg.second; VectorImgType::PixelType measured_pixel = m_DiffImage->GetPixel(seg.first); double simulated_mean = 0; double measured_mean = 0; int num_nonzero_g = 0; for (int g=0; gGetGradientDirection(g).GetNorm()GetLargestPossibleRegion().GetSize(0); sz_y = m_PeakImage->GetLargestPossibleRegion().GetSize(1); sz_z = m_PeakImage->GetLargestPossibleRegion().GetSize(2); dim_four_size = m_PeakImage->GetLargestPossibleRegion().GetSize(3)/3 + 1; // +1 for zero - peak int num_voxels = sz_x*sz_y*sz_z; itk::Vector< double, 3 > spacing; spacing[0] = m_PeakImage->GetSpacing()[0]; spacing[1] = m_PeakImage->GetSpacing()[1]; spacing[2] = m_PeakImage->GetSpacing()[2]; PointType3 origin; origin[0] = m_PeakImage->GetOrigin()[0]; origin[1] = m_PeakImage->GetOrigin()[1]; origin[2] = m_PeakImage->GetOrigin()[2]; UcharImgType::RegionType::SizeType size; size[0] = m_PeakImage->GetLargestPossibleRegion().GetSize()[0]; size[1] = m_PeakImage->GetLargestPossibleRegion().GetSize()[1]; size[2] = m_PeakImage->GetLargestPossibleRegion().GetSize()[2]; UcharImgType::RegionType imageRegion; imageRegion.SetSize(size); itk::Matrix direction; for (int r=0; r<3; r++) for (int c=0; c<3; c++) direction[r][c] = m_PeakImage->GetDirection()[r][c]; if (m_MaskImage.IsNull()) { m_MaskImage = UcharImgType::New(); m_MaskImage->SetSpacing( spacing ); m_MaskImage->SetOrigin( origin ); m_MaskImage->SetDirection( direction ); m_MaskImage->SetRegions( imageRegion ); m_MaskImage->Allocate(); m_MaskImage->FillBuffer(1); } else if (m_MaskImage->GetLargestPossibleRegion().GetSize()[0]!=sz_x || m_MaskImage->GetLargestPossibleRegion().GetSize()[1]!=sz_y || m_MaskImage->GetLargestPossibleRegion().GetSize()[2]!=sz_z) mitkThrow() << "Mask and peak image must have same size!"; m_NumResiduals = num_voxels * dim_four_size; MITK_INFO << "Num. unknowns: " << m_NumUnknowns; MITK_INFO << "Num. residuals: " << m_NumResiduals; MITK_INFO << "Creating system ..."; A.set_size(m_NumResiduals, m_NumUnknowns); b.set_size(m_NumResiduals); b.fill(0.0); m_MeanTractDensity = 0; m_MeanSignal = 0; m_NumCoveredDirections = 0; fiber_count = 0; m_GroupSizes.clear(); int numFibers = 0; for (unsigned int bundle=0; bundleGetNumFibers(); boost::progress_display disp(numFibers); for (unsigned int bundle=0; bundle polydata = m_Tractograms.at(bundle)->GetFiberPolyData(); m_GroupSizes.push_back(m_Tractograms.at(bundle)->GetNumFibers()); for (unsigned int i=0; iGetNumFibers(); ++i) { ++disp; vtkCell* cell = polydata->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (numPoints<2) MITK_INFO << "FIBER WITH ONLY ONE POINT ENCOUNTERED!"; for (int j=0; jGetPoint(j)); itk::Index<3> startIndex; itk::ContinuousIndex startIndexCont; m_MaskImage->TransformPhysicalPointToIndex(startVertex, startIndex); m_MaskImage->TransformPhysicalPointToContinuousIndex(startVertex, startIndexCont); PointType3 endVertex = mitk::imv::GetItkPoint(points->GetPoint(j+1)); itk::Index<3> endIndex; itk::ContinuousIndex endIndexCont; m_MaskImage->TransformPhysicalPointToIndex(endVertex, endIndex); m_MaskImage->TransformPhysicalPointToContinuousIndex(endVertex, endIndexCont); vnl_vector_fixed fiber_dir; fiber_dir[0] = endVertex[0]-startVertex[0]; fiber_dir[1] = endVertex[1]-startVertex[1]; fiber_dir[2] = endVertex[2]-startVertex[2]; fiber_dir.normalize(); std::vector< std::pair< itk::Index<3>, double > > segments = mitk::imv::IntersectImage(spacing, startIndex, endIndex, startIndexCont, endIndexCont); for (std::pair< itk::Index<3>, double > seg : segments) { if (!m_MaskImage->GetLargestPossibleRegion().IsInside(seg.first) || m_MaskImage->GetPixel(seg.first)==0) continue; itk::Index<4> idx4; idx4[0]=seg.first[0]; idx4[1]=seg.first[1]; idx4[2]=seg.first[2]; idx4[3]=0; double w = 1; int peak_id = dim_four_size-1; double peak_mag = 0; GetClosestPeak(idx4, m_PeakImage, fiber_dir, peak_id, w, peak_mag); w *= seg.second; int x = idx4[0]; int y = idx4[1]; int z = idx4[2]; unsigned int linear_index = x + sz_x*y + sz_x*sz_y*z + sz_x*sz_y*sz_z*peak_id; if (b[linear_index] == 0 && peak_idGetLargestPossibleRegion().GetSize(0); sz_y = m_ScalarImage->GetLargestPossibleRegion().GetSize(1); sz_z = m_ScalarImage->GetLargestPossibleRegion().GetSize(2); int num_voxels = sz_x*sz_y*sz_z; if (m_MaskImage.IsNotNull() && (m_MaskImage->GetLargestPossibleRegion().GetSize()[0]!=sz_x || m_MaskImage->GetLargestPossibleRegion().GetSize()[1]!=sz_y || m_MaskImage->GetLargestPossibleRegion().GetSize()[2]!=sz_z)) mitkThrow() << "Mask and peak image must have same size!"; auto spacing = m_ScalarImage->GetSpacing(); m_NumResiduals = num_voxels; MITK_INFO << "Num. unknowns: " << m_NumUnknowns; MITK_INFO << "Num. residuals: " << m_NumResiduals; MITK_INFO << "Creating system ..."; A.set_size(m_NumResiduals, m_NumUnknowns); b.set_size(m_NumResiduals); b.fill(0.0); m_MeanTractDensity = 0; m_MeanSignal = 0; int numCoveredVoxels = 0; fiber_count = 0; int numFibers = 0; for (unsigned int bundle=0; bundleGetNumFibers(); boost::progress_display disp(numFibers); m_GroupSizes.clear(); for (unsigned int bundle=0; bundle polydata = m_Tractograms.at(bundle)->GetFiberPolyData(); m_GroupSizes.push_back(m_Tractograms.at(bundle)->GetNumFibers()); for (unsigned int i=0; iGetNumFibers(); ++i) { ++disp; vtkCell* cell = polydata->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; jGetPoint(j)); itk::Index<3> startIndex; itk::ContinuousIndex startIndexCont; m_ScalarImage->TransformPhysicalPointToIndex(startVertex, startIndex); m_ScalarImage->TransformPhysicalPointToContinuousIndex(startVertex, startIndexCont); PointType3 endVertex = mitk::imv::GetItkPoint(points->GetPoint(j+1)); itk::Index<3> endIndex; itk::ContinuousIndex endIndexCont; m_ScalarImage->TransformPhysicalPointToIndex(endVertex, endIndex); m_ScalarImage->TransformPhysicalPointToContinuousIndex(endVertex, endIndexCont); std::vector< std::pair< itk::Index<3>, double > > segments = mitk::imv::IntersectImage(spacing, startIndex, endIndex, startIndexCont, endIndexCont); for (std::pair< itk::Index<3>, double > seg : segments) { if (!m_ScalarImage->GetLargestPossibleRegion().IsInside(seg.first) || (m_MaskImage.IsNotNull() && m_MaskImage->GetPixel(seg.first)==0)) continue; float image_value = m_ScalarImage->GetPixel(seg.first); int x = seg.first[0]; int y = seg.first[1]; int z = seg.first[2]; unsigned int linear_index = x + sz_x*y + sz_x*sz_y*z; if (b[linear_index] == 0) { numCoveredVoxels++; m_MeanSignal += image_value; } m_MeanTractDensity += seg.second; if (m_FitIndividualFibers) { b[linear_index] = image_value; A.put(linear_index, fiber_count, A.get(linear_index, fiber_count) + seg.second); } else { b[linear_index] = image_value; A.put(linear_index, bundle, A.get(linear_index, bundle) + seg.second); } } } ++fiber_count; } } if (numCoveredVoxels==0 || m_MeanSignal<0.0001) mitkThrow() << "No overlap between fibers and non-zero image!"; m_MeanTractDensity /= numCoveredVoxels; m_MeanSignal /= numCoveredVoxels; b /= m_MeanSignal; b *= m_MeanTractDensity; } void FitFibersToImageFilter::GenerateData() { m_NumUnknowns = m_Tractograms.size(); if (m_FitIndividualFibers) { m_NumUnknowns = 0; for (unsigned int bundle=0; bundleGetNumFibers(); } else m_FilterOutliers = false; if (m_NumUnknowns<1) { MITK_INFO << "No fibers in tractogram."; return; } fiber_count = 0; sz_x = 0; sz_y = 0; sz_z = 0; m_MeanTractDensity = 0; m_MeanSignal = 0; if (m_PeakImage.IsNotNull()) CreatePeakSystem(); else if (m_DiffImage.IsNotNull()) CreateDiffSystem(); else if (m_ScalarImage.IsNotNull()) CreateScalarSystem(); else mitkThrow() << "No input image set!"; MITK_INFO << "Initializing optimizer"; double init_lambda = fiber_count; // initialization for lambda estimation itk::TimeProbe clock; clock.Start(); cost = VnlCostFunction(m_NumUnknowns); cost.SetProblem(A, b, init_lambda, m_Regularization); cost.SetGroupSizes(m_GroupSizes); m_Weights.set_size(m_NumUnknowns); m_Weights.fill( 1.0/m_NumUnknowns ); vnl_lbfgsb minimizer(cost); vnl_vector l; l.set_size(m_NumUnknowns); l.fill(0); vnl_vector bound_selection; bound_selection.set_size(m_NumUnknowns); bound_selection.fill(1); minimizer.set_bound_selection(bound_selection); minimizer.set_lower_bound(l); minimizer.set_projected_gradient_tolerance(m_GradientTolerance); if (m_Regularization==VnlCostFunction::REGU::MSM) MITK_INFO << "Regularization type: MSM"; else if (m_Regularization==VnlCostFunction::REGU::VARIANCE) MITK_INFO << "Regularization type: VARIANCE"; else if (m_Regularization==VnlCostFunction::REGU::LASSO) MITK_INFO << "Regularization type: LASSO"; else if (m_Regularization==VnlCostFunction::REGU::VOXEL_VARIANCE) MITK_INFO << "Regularization type: VOXEL_VARIANCE"; else if (m_Regularization==VnlCostFunction::REGU::GROUP_LASSO) MITK_INFO << "Regularization type: GROUP_LASSO"; else if (m_Regularization==VnlCostFunction::REGU::GROUP_VARIANCE) MITK_INFO << "Regularization type: GROUP_VARIANCE"; else if (m_Regularization==VnlCostFunction::REGU::NONE) MITK_INFO << "Regularization type: NONE"; -// if (m_Regularization!=VnlCostFunction::REGU::NONE) // REMOVE FOR NEW FIT AND SET cost.m_Lambda = m_Lambda -// { -// MITK_INFO << "Estimating regularization"; -// minimizer.set_trace(false); -// minimizer.set_max_function_evals(2); -// minimizer.minimize(m_Weights); -// vnl_vector dx; dx.set_size(m_NumUnknowns); dx.fill(0.0); -// cost.calc_regularization_gradient(m_Weights, dx); - -// if (m_Weights.magnitude()==0) -// { -// MITK_INFO << "Regularization estimation failed. Using default value."; -// cost.m_Lambda = fiber_count*m_Lambda; -// } -// else -// { -// double r = dx.magnitude()/m_Weights.magnitude(); // wtf??? -// cost.m_Lambda *= m_Lambda*55.0/r; -// MITK_INFO << r << " - " << m_Lambda*55.0/r; -// if (cost.m_Lambda>10e7) -// { -// MITK_INFO << "Regularization estimation failed. Using default value."; -// cost.m_Lambda = fiber_count*m_Lambda; -// } -// } -// } -// else -// cost.m_Lambda = 0; cost.m_Lambda = m_Lambda; MITK_INFO << "Using regularization factor of " << cost.m_Lambda << " (λ: " << m_Lambda << ")"; MITK_INFO << "Fitting fibers"; minimizer.set_trace(m_Verbose); minimizer.set_max_function_evals(m_MaxIterations); minimizer.minimize(m_Weights); std::vector< double > weights; if (m_FilterOutliers) { for (auto w : m_Weights) weights.push_back(w); std::sort(weights.begin(), weights.end()); MITK_INFO << "Setting upper weight bound to " << weights.at(m_NumUnknowns*0.99); vnl_vector u; u.set_size(m_NumUnknowns); u.fill(weights.at(m_NumUnknowns*0.99)); minimizer.set_upper_bound(u); bound_selection.fill(2); minimizer.set_bound_selection(bound_selection); minimizer.minimize(m_Weights); weights.clear(); } for (auto w : m_Weights) weights.push_back(w); std::sort(weights.begin(), weights.end()); m_MeanWeight = m_Weights.mean(); m_MedianWeight = weights.at(m_NumUnknowns*0.5); m_MinWeight = weights.at(0); m_MaxWeight = weights.at(m_NumUnknowns-1); MITK_INFO << "*************************"; MITK_INFO << "Weight statistics"; MITK_INFO << "Sum: " << m_Weights.sum(); MITK_INFO << "Mean: " << m_MeanWeight; MITK_INFO << "1% quantile: " << weights.at(m_NumUnknowns*0.01); MITK_INFO << "5% quantile: " << weights.at(m_NumUnknowns*0.05); MITK_INFO << "25% quantile: " << weights.at(m_NumUnknowns*0.25); MITK_INFO << "Median: " << m_MedianWeight; MITK_INFO << "75% quantile: " << weights.at(m_NumUnknowns*0.75); MITK_INFO << "95% quantile: " << weights.at(m_NumUnknowns*0.95); MITK_INFO << "99% quantile: " << weights.at(m_NumUnknowns*0.99); MITK_INFO << "Min: " << m_MinWeight; MITK_INFO << "Max: " << m_MaxWeight; MITK_INFO << "*************************"; MITK_INFO << "NumEvals: " << minimizer.get_num_evaluations(); MITK_INFO << "NumIterations: " << minimizer.get_num_iterations(); MITK_INFO << "Residual cost: " << minimizer.get_end_error(); m_RMSE = cost.S->get_rms_error(m_Weights); MITK_INFO << "Final RMSE: " << m_RMSE; clock.Stop(); int h = clock.GetTotal()/3600; int m = ((int)clock.GetTotal()%3600)/60; int s = (int)clock.GetTotal()%60; MITK_INFO << "Optimization took " << h << "h, " << m << "m and " << s << "s"; MITK_INFO << "Weighting fibers"; m_RmsDiffPerBundle.set_size(m_Tractograms.size()); std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); if (m_FitIndividualFibers) { unsigned int fiber_count = 0; for (unsigned int bundle=0; bundle temp_weights; temp_weights.set_size(m_Weights.size()); temp_weights.copy_in(m_Weights.data_block()); for (unsigned int i=0; iGetNumFibers(); i++) { m_Tractograms.at(bundle)->SetFiberWeight(i, m_Weights[fiber_count]); temp_weights[fiber_count] = 0; ++fiber_count; } double d_rms = cost.S->get_rms_error(temp_weights) - m_RMSE; m_RmsDiffPerBundle[bundle] = d_rms; } } else { for (unsigned int i=0; i temp_weights; temp_weights.set_size(m_Weights.size()); temp_weights.copy_in(m_Weights.data_block()); temp_weights[i] = 0; double d_rms = cost.S->get_rms_error(temp_weights) - m_RMSE; m_RmsDiffPerBundle[i] = d_rms; m_Tractograms.at(i)->SetFiberWeights(m_Weights[i]); } } std::cout.rdbuf (old); // transform back A *= m_MeanSignal/100.0; b *= m_MeanSignal/100.0; MITK_INFO << "Generating output images ..."; if (m_PeakImage.IsNotNull()) GenerateOutputPeakImages(); else if (m_DiffImage.IsNotNull()) GenerateOutputDiffImages(); else if (m_ScalarImage.IsNotNull()) GenerateOutputScalarImages(); m_Coverage = m_Coverage/m_MeanSignal; m_Overshoot = m_Overshoot/m_MeanSignal; MITK_INFO << std::fixed << "Coverage: " << setprecision(2) << 100.0*m_Coverage << "%"; MITK_INFO << std::fixed << "Overshoot: " << setprecision(2) << 100.0*m_Overshoot << "%"; } void FitFibersToImageFilter::GenerateOutputDiffImages() { VectorImgType::PixelType pix; pix.SetSize(m_DiffImage->GetVectorLength()); pix.Fill(0); itk::ImageDuplicator< VectorImgType >::Pointer duplicator = itk::ImageDuplicator< VectorImgType >::New(); duplicator->SetInputImage(m_DiffImage); duplicator->Update(); m_UnderexplainedImageDiff = duplicator->GetOutput(); m_UnderexplainedImageDiff->FillBuffer(pix); duplicator->SetInputImage(m_UnderexplainedImageDiff); duplicator->Update(); m_OverexplainedImageDiff = duplicator->GetOutput(); m_OverexplainedImageDiff->FillBuffer(pix); duplicator->SetInputImage(m_OverexplainedImageDiff); duplicator->Update(); m_ResidualImageDiff = duplicator->GetOutput(); m_ResidualImageDiff->FillBuffer(pix); duplicator->SetInputImage(m_ResidualImageDiff); duplicator->Update(); m_FittedImageDiff = duplicator->GetOutput(); m_FittedImageDiff->FillBuffer(pix); vnl_vector fitted_b; fitted_b.set_size(b.size()); cost.S->multiply(m_Weights, fitted_b); itk::ImageRegionIterator it1 = itk::ImageRegionIterator(m_DiffImage, m_DiffImage->GetLargestPossibleRegion()); itk::ImageRegionIterator it2 = itk::ImageRegionIterator(m_FittedImageDiff, m_FittedImageDiff->GetLargestPossibleRegion()); itk::ImageRegionIterator it3 = itk::ImageRegionIterator(m_ResidualImageDiff, m_ResidualImageDiff->GetLargestPossibleRegion()); itk::ImageRegionIterator it4 = itk::ImageRegionIterator(m_UnderexplainedImageDiff, m_UnderexplainedImageDiff->GetLargestPossibleRegion()); itk::ImageRegionIterator it5 = itk::ImageRegionIterator(m_OverexplainedImageDiff, m_OverexplainedImageDiff->GetLargestPossibleRegion()); m_MeanSignal = 0; m_Coverage = 0; m_Overshoot = 0; while( !it2.IsAtEnd() ) { itk::Index<3> idx3 = it2.GetIndex(); VectorImgType::PixelType original_pix =it1.Get(); VectorImgType::PixelType fitted_pix =it2.Get(); VectorImgType::PixelType residual_pix =it3.Get(); VectorImgType::PixelType underexplained_pix =it4.Get(); VectorImgType::PixelType overexplained_pix =it5.Get(); int num_nonzero_g = 0; double original_mean = 0; for (int g=0; gGetGradientDirection(g).GetNorm()>=mitk::eps ) { original_mean += original_pix[g]; ++num_nonzero_g; } } original_mean /= num_nonzero_g; for (int g=0; g=0) { underexplained_pix[g] = residual_pix[g]; m_Coverage += fitted_b[linear_index] + original_mean; } m_MeanSignal += b[linear_index] + original_mean; } it2.Set(fitted_pix); it3.Set(residual_pix); it4.Set(underexplained_pix); it5.Set(overexplained_pix); ++it1; ++it2; ++it3; ++it4; ++it5; } } void FitFibersToImageFilter::GenerateOutputScalarImages() { itk::ImageDuplicator< DoubleImgType >::Pointer duplicator = itk::ImageDuplicator< DoubleImgType >::New(); duplicator->SetInputImage(m_ScalarImage); duplicator->Update(); m_UnderexplainedImageScalar = duplicator->GetOutput(); m_UnderexplainedImageScalar->FillBuffer(0); duplicator->SetInputImage(m_UnderexplainedImageScalar); duplicator->Update(); m_OverexplainedImageScalar = duplicator->GetOutput(); m_OverexplainedImageScalar->FillBuffer(0); duplicator->SetInputImage(m_OverexplainedImageScalar); duplicator->Update(); m_ResidualImageScalar = duplicator->GetOutput(); m_ResidualImageScalar->FillBuffer(0); duplicator->SetInputImage(m_ResidualImageScalar); duplicator->Update(); m_FittedImageScalar = duplicator->GetOutput(); m_FittedImageScalar->FillBuffer(0); vnl_vector fitted_b; fitted_b.set_size(b.size()); cost.S->multiply(m_Weights, fitted_b); itk::ImageRegionIterator it1 = itk::ImageRegionIterator(m_ScalarImage, m_ScalarImage->GetLargestPossibleRegion()); itk::ImageRegionIterator it2 = itk::ImageRegionIterator(m_FittedImageScalar, m_FittedImageScalar->GetLargestPossibleRegion()); itk::ImageRegionIterator it3 = itk::ImageRegionIterator(m_ResidualImageScalar, m_ResidualImageScalar->GetLargestPossibleRegion()); itk::ImageRegionIterator it4 = itk::ImageRegionIterator(m_UnderexplainedImageScalar, m_UnderexplainedImageScalar->GetLargestPossibleRegion()); itk::ImageRegionIterator it5 = itk::ImageRegionIterator(m_OverexplainedImageScalar, m_OverexplainedImageScalar->GetLargestPossibleRegion()); m_MeanSignal = 0; m_Coverage = 0; m_Overshoot = 0; while( !it2.IsAtEnd() ) { itk::Index<3> idx3 = it2.GetIndex(); DoubleImgType::PixelType original_pix =it1.Get(); DoubleImgType::PixelType fitted_pix =it2.Get(); DoubleImgType::PixelType residual_pix =it3.Get(); DoubleImgType::PixelType underexplained_pix =it4.Get(); DoubleImgType::PixelType overexplained_pix =it5.Get(); unsigned int linear_index = idx3[0] + sz_x*idx3[1] + sz_x*sz_y*idx3[2]; fitted_pix = fitted_b[linear_index]; residual_pix = original_pix - fitted_pix; if (residual_pix<0) { overexplained_pix = residual_pix; m_Coverage += b[linear_index]; m_Overshoot -= residual_pix; } else if (residual_pix>=0) { underexplained_pix = residual_pix; m_Coverage += fitted_b[linear_index]; } m_MeanSignal += b[linear_index]; it2.Set(fitted_pix); it3.Set(residual_pix); it4.Set(underexplained_pix); it5.Set(overexplained_pix); ++it1; ++it2; ++it3; ++it4; ++it5; } } VnlCostFunction::REGU FitFibersToImageFilter::GetRegularization() const { return m_Regularization; } void FitFibersToImageFilter::SetRegularization(const VnlCostFunction::REGU &Regularization) { m_Regularization = Regularization; } void FitFibersToImageFilter::GenerateOutputPeakImages() { itk::ImageDuplicator< PeakImgType >::Pointer duplicator = itk::ImageDuplicator< PeakImgType >::New(); duplicator->SetInputImage(m_PeakImage); duplicator->Update(); m_UnderexplainedImage = duplicator->GetOutput(); m_UnderexplainedImage->FillBuffer(0.0); duplicator->SetInputImage(m_UnderexplainedImage); duplicator->Update(); m_OverexplainedImage = duplicator->GetOutput(); m_OverexplainedImage->FillBuffer(0.0); duplicator->SetInputImage(m_OverexplainedImage); duplicator->Update(); m_ResidualImage = duplicator->GetOutput(); m_ResidualImage->FillBuffer(0.0); duplicator->SetInputImage(m_ResidualImage); duplicator->Update(); m_FittedImage = duplicator->GetOutput(); m_FittedImage->FillBuffer(0.0); vnl_vector fitted_b; fitted_b.set_size(b.size()); cost.S->multiply(m_Weights, fitted_b); for (unsigned int r=0; r idx4; unsigned int linear_index = r; idx4[0] = linear_index % sz_x; linear_index /= sz_x; idx4[1] = linear_index % sz_y; linear_index /= sz_y; idx4[2] = linear_index % sz_z; linear_index /= sz_z; int peak_id = linear_index % dim_four_size; if (peak_id peak_dir; idx4[3] = peak_id*3; peak_dir[0] = m_PeakImage->GetPixel(idx4); idx4[3] += 1; peak_dir[1] = m_PeakImage->GetPixel(idx4); idx4[3] += 1; peak_dir[2] = m_PeakImage->GetPixel(idx4); peak_dir.normalize(); peak_dir *= fitted_b[r]; idx4[3] = peak_id*3; m_FittedImage->SetPixel(idx4, peak_dir[0]); idx4[3] += 1; m_FittedImage->SetPixel(idx4, peak_dir[1]); idx4[3] += 1; m_FittedImage->SetPixel(idx4, peak_dir[2]); } } m_MeanSignal = 0; m_Coverage = 0; m_Overshoot = 0; itk::Index<4> idx4; for (idx4[0]=0; idx4[0] idx3; idx3[0] = idx4[0]; idx3[1] = idx4[1]; idx3[2] = idx4[2]; if (m_MaskImage.IsNotNull() && m_MaskImage->GetPixel(idx3)==0) continue; vnl_vector_fixed peak_dir; vnl_vector_fixed fitted_dir; vnl_vector_fixed overshoot_dir; for (idx4[3]=0; idx4[3]<(itk::IndexValueType)m_PeakImage->GetLargestPossibleRegion().GetSize(3); ++idx4[3]) { peak_dir[idx4[3]%3] = m_PeakImage->GetPixel(idx4); fitted_dir[idx4[3]%3] = m_FittedImage->GetPixel(idx4); m_ResidualImage->SetPixel(idx4, m_PeakImage->GetPixel(idx4) - m_FittedImage->GetPixel(idx4)); if (idx4[3]%3==2) { m_MeanSignal += peak_dir.magnitude(); itk::Index<4> tidx= idx4; if (peak_dir.magnitude()>fitted_dir.magnitude()) { m_Coverage += fitted_dir.magnitude(); m_UnderexplainedImage->SetPixel(tidx, peak_dir[2]-fitted_dir[2]); tidx[3]--; m_UnderexplainedImage->SetPixel(tidx, peak_dir[1]-fitted_dir[1]); tidx[3]--; m_UnderexplainedImage->SetPixel(tidx, peak_dir[0]-fitted_dir[0]); } else { overshoot_dir[0] = fitted_dir[0]-peak_dir[0]; overshoot_dir[1] = fitted_dir[1]-peak_dir[1]; overshoot_dir[2] = fitted_dir[2]-peak_dir[2]; m_Coverage += peak_dir.magnitude(); m_Overshoot += overshoot_dir.magnitude(); m_OverexplainedImage->SetPixel(tidx, overshoot_dir[2]); tidx[3]--; m_OverexplainedImage->SetPixel(tidx, overshoot_dir[1]); tidx[3]--; m_OverexplainedImage->SetPixel(tidx, overshoot_dir[0]); } } } } } void FitFibersToImageFilter::GetClosestPeak(itk::Index<4> idx, PeakImgType::Pointer peak_image , vnl_vector_fixed fiber_dir, int& id, double& w, double& peak_mag ) { int m_NumDirs = peak_image->GetLargestPossibleRegion().GetSize()[3]/3; vnl_vector_fixed out_dir; out_dir.fill(0); float angle = 0.9; for (int i=0; i dir; idx[3] = i*3; dir[0] = peak_image->GetPixel(idx); idx[3] += 1; dir[1] = peak_image->GetPixel(idx); idx[3] += 1; dir[2] = peak_image->GetPixel(idx); float mag = dir.magnitude(); if (magangle) { angle = a; w = angle; peak_mag = mag; id = i; } } } std::vector FitFibersToImageFilter::GetTractograms() const { return m_Tractograms; } void FitFibersToImageFilter::SetTractograms(const std::vector &tractograms) { m_Tractograms = tractograms; } void FitFibersToImageFilter::SetSignalModel(mitk::DiffusionSignalModel<> *SignalModel) { m_SignalModel = SignalModel; } } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.ivim/documentation/UserManual/QmitkUserIVIMViewManual.dox b/Plugins/org.mitk.gui.qt.diffusionimaging.ivim/documentation/UserManual/QmitkUserIVIMViewManual.dox index 6549fe1..79938ad 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.ivim/documentation/UserManual/QmitkUserIVIMViewManual.dox +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.ivim/documentation/UserManual/QmitkUserIVIMViewManual.dox @@ -1,41 +1,34 @@ /** \page org_mitk_views_ivim Intra-voxel incoherent motion estimation (IVIM) -The required input for the "Intra-voxel incoherent motion estimation" (IVIM) is a diffusion weighted image (.dwi or .hdwi) that was acquired with several different b-values. +This view enables the estimation of "Intra-voxel incoherent motion estimation" (IVIM) and diffusion kurtosis parameters from diffusion weighted images. The required input is a diffusion weighted image that was acquired with multiple different b-values (unidirectional). -\imageMacro{ivimview.png,"The IVIM View",13.91} +The view enables an interactive exploration of the selected image (click around in the image and watch the estimated parameters in the figure of the view) as well as generation of parameter maps. -Once an input image is selected in the datamanager, the IVIM view allows for interactive exploration of the dataset (click around in the image and watch the estimated parameters in the figure of the view) as well as generation of f-, D-, and D*-maps (activate the checkmarks and press "Generate Output Images"). -The "neglect b<" threshold allows you to ignore b-values smaller then a threshold for the initial fit of f and D. D* is then estimated using all measurements. +\section QmitkIVIMViewAdvanced IVIM fit methods -The exact values of the current fit are always given in the legend underneath the figure. +The IVIM model is a bi-exponential model defined by the parameters f, D and D*: \f$(1-f)*exp(-b*D)+f*exp(-b*(D+D*))\f$. All fit methods available in the view ultimatle use this model but they differ how the three parameters are estimated: -\section QmitkIVIMViewROIAnalysis Region of interest analysis +\li Jointly fit D, f and D*: All three parameters are estimated simultaneously. +\li Fit D & f with fixed D* value: Only f and D are estimated, D* is fixed (user defined). +\li Fit D & f (high b), then fit D*: First fit f and D (monoexponentially \f$(1-f)*exp(-b*D)\f$) and use these then fixed parameters in a second fit of D* with the complete bi-exponential model. +\li Linearly fit D & f (high b), then fit D*: First fit f and D linearly and use these then fixed parameters in a second fit of D* with the complete bi-exponential model. -Create region of interest: To create a new segmentatin, open the "quantification" perspective, select the tab "Segmentation", and create a segmentation of the structure of interest. Alternatively, of course, you may also load a binary image from file or generate your segmentation in any other possible way. +Negative values for D and D* are penalized during the fit. -IVIM in region of interset: Go back to the "IVIM" perspective and select both, the diffusion image and the segmentation (holding the CTRL key). A red message should appear "Averaging N voxels". - -\section QmitkIVIMViewExport Export - -All model parameters and corresponding curves can be exported to clipboard using the buttons underneath the figure. - -\section QmitkIVIMViewAdvanced Advanced Settings - -Advanced users, that know what they are doing, can change the method for the model-fit under "Advanced Settings" on the very bottom of the view. 3-param fit, linear fit of f/D, and fix D* are among the options. \section QmitkIVIMViewReferences Suggested Readings Toward an optimal distribution of b values for intravoxel incoherent motion imaging. Lemke A, Stieltjes B, Schad LR, Laun FB. Magn Reson Imaging. 2011 Jul;29(6):766-76. Epub 2011 May 5. PMID: 21549538 Differentiation of pancreas carcinoma from healthy pancreatic tissue using multiple b-values: comparison of apparent diffusion coefficient and intravoxel incoherent motion derived parameters. Lemke A, Laun FB, Klauss M, Re TJ, Simon D, Delorme S, Schad LR, Stieltjes B. Invest Radiol. 2009 Dec;44(12):769-75. PMID: 19838121 */ 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 6adcaf1..c6bfd96 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,1063 +1,1076 @@ /*=================================================================== 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::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::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::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; + roisize = 3; 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)); +// roiImage->SetPixel(newstart, vecimg->GetPixel(index)); + + typedef itk::ImageRegionIterator VectorIteratorType; + VectorIteratorType vecit(vecimg, vecimg->GetRequestedRegion() ); + vecit.GoToBegin(); + + typedef itk::ImageRegionIterator VectorIteratorType; + VectorIteratorType vecit2(roiImage, roiImage->GetLargestPossibleRegion() ); + vecit2.GoToBegin(); + + while( !vecit.IsAtEnd() ) + { + vecit2.Set(vecit.Get()); + + ++vecit; + ++vecit2; + } 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.m_Bzero; if (m_KurtosisSnap.m_fittedBZero) S0 = m_KurtosisSnap.m_BzeroFit; std::vector< std::pair > d_line; d_line.emplace_back(0, S0); d_line.emplace_back(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::vector< std::pair > y; for( unsigned int i=0; i<=num_samples; ++i) { double b = (((1.0)*i)/(1.0*num_samples))*maxb; y.emplace_back(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::vector< std::pair > y_meas; for (unsigned int i=0; im_OmitBZeroCB->isChecked() || m_KurtosisSnap.bvalues[i] > 0.01) y_meas.emplace_back(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::vector< std::pair > d_line; d_line.emplace_back(0, 1-m_IvimSnap.currentFunceiled); d_line.emplace_back(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::vector< std::pair > y; + int nsampling = 50; + double f = 1-m_IvimSnap.currentFunceiled; + for(int i=0; i<=nsampling; i++) { - 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.emplace_back(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"); + double x = (((1.0)*i)/(1.0*nsampling))*maxb; + y.emplace_back(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::vector< std::pair > y_meas; - for (unsigned int i=0; i > y_meas; + for (unsigned int i=0; im_ChartWidget->UpdateData2D(y_meas, "signal values"); + m_Controls->m_ChartWidget->UpdateData2D(y_meas, "signal values"); - if(!m_IvimSnap.high_indices.empty()) - { - AddSecondFitPlot(); - std::vector< std::pair > additonal_meas; - for (int i=0; i > additonal_meas; + for (int i=0; im_ChartWidget->UpdateData2D(additonal_meas, "signal values (used for second fit)"); - } - else - { - RemoveSecondFitPlot(); - } + m_Controls->m_ChartWidget->UpdateData2D(additonal_meas, "signal values (used for second fit)"); + } + else + { + 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 ); } else { std::stringstream ss; QString clipboard("Normalized Measurement Points\n"); ss << m_IvimSnap.bvalues << "\n" << m_IvimSnap.allmeas << "\n\n"; ss << "Fitted Values ( f D D* ) \n" << m_IvimSnap.currentF << " " << m_IvimSnap.currentD << " " << m_IvimSnap.currentDStar; ss << "\n\n"; clipboard.append( QString( ss.str().c_str() )); ss.str( std::string() ); ss.clear(); 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 ); } 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 ); } } void QmitkIVIMView::Activated() { m_Active = true; } void QmitkIVIMView::Deactivated() { m_Active = false; } void QmitkIVIMView::Visible() { m_Visible = true; if (this->GetRenderWindowPart() && !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 0ec0c50..9d36f29 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,1213 @@ QmitkIVIMViewControls 0 0 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 + 0 - 3 Param. Fit + Jointly fit D, f and D* 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 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