diff --git a/Modules/DiffusionImaging/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.cpp b/Modules/DiffusionImaging/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.cpp index beeba3d3e5..b78aaf57e8 100644 --- a/Modules/DiffusionImaging/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.cpp +++ b/Modules/DiffusionImaging/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.cpp @@ -1,775 +1,782 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __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 "itkRegularizedIVIMReconstructionFilter.h" #include #define IVIM_FOO -100000 namespace itk { - template< class TIn, class TOut> - DiffusionIntravoxelIncoherentMotionReconstructionImageFilter - ::DiffusionIntravoxelIncoherentMotionReconstructionImageFilter() : - m_GradientDirectionContainer(NULL), - m_Method(IVIM_DSTAR_FIX), - m_FitDStar(true), - m_Verbose(false) - { +template< class TIn, class TOut> +DiffusionIntravoxelIncoherentMotionReconstructionImageFilter +::DiffusionIntravoxelIncoherentMotionReconstructionImageFilter() : + m_GradientDirectionContainer(NULL), + 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() - { +template< class TIn, class TOut> +void DiffusionIntravoxelIncoherentMotionReconstructionImageFilter +::BeforeThreadedGenerateData() +{ // Input must be an itk::VectorImage. std::string gradientImageClassName( - this->ProcessObject::GetInput(0)->GetNameOfClass()); + 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 ); + itkExceptionMacro( << + "There is only one Gradient image. I expect that to be a VectorImage. " + << "But its of type: " << gradientImageClassName ); } - typename InputImageType::Pointer img = static_cast< InputImageType * >( - this->ProcessObject::GetInput(0) ); - // Compute the indicies of the baseline images and gradient images - GradientDirectionContainerType::ConstIterator gdcit = - this->m_GradientDirectionContainer->Begin(); + // 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() ) { - if(gdcit.Value().one_norm() <= 0.0) - { - 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; + 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; im_BThres) + 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) - { +template< class TIn, class TOut> +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, - int ) - { +template< class TIn, class TOut> +void DiffusionIntravoxelIncoherentMotionReconstructionImageFilter +::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, + int ) +{ typename OutputImageType::Pointer outputImage = - static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); + static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread); oit.GoToBegin(); typename OutputImageType::Pointer dImage = - static_cast< OutputImageType * >(this->ProcessObject::GetOutput(1)); + 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)); + 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 = NULL; // 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) ); + 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->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); + m_InternalVectorImage->SetVectorLength(m_Snap.Nhigh); + 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(); + InputVectorType measvec = iit.Get(); - typename NumericTraits::AccumulateType b0 = NumericTraits::Zero; + typename NumericTraits::AccumulateType b0 = NumericTraits::Zero; - m_Snap.meas.set_size(m_Snap.N); - m_Snap.allmeas.set_size(m_Snap.N); - if(!m_Snap.iterated_sequence) - { - // Average the baseline image pixels - for(unsigned int i = 0; i < m_Snap.baselineind.size(); ++i) + m_Snap.meas.set_size(m_Snap.N); + m_Snap.allmeas.set_size(m_Snap.N); + if(!m_Snap.iterated_sequence) { - b0 += measvec[m_Snap.baselineind[i]]; - } - if(m_Snap.baselineind.size()) - b0 /= m_Snap.baselineind.size(); + // 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) - { - 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); - } - else - { - m_Snap.meas[i] = IVIM_FOO; - } + // measurement vector + for(int i = 0; i < m_Snap.N; ++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); + } + else + { + m_Snap.meas[i] = IVIM_FOO; + } + } } - } - else - { - // measurement vector - for(int i = 0; i < m_Snap.N; ++i) + else { - 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); - } - else - { - m_Snap.meas[i] = IVIM_FOO; - } + // measurement vector + for(int i = 0; i < m_Snap.N; ++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); + } + else + { + m_Snap.meas[i] = IVIM_FOO; + } + } } - } - m_Snap.currentF = 0; - m_Snap.currentD = 0; - m_Snap.currentDStar = 0; + m_Snap.currentF = 0; + m_Snap.currentD = 0; + m_Snap.currentDStar = 0; - switch(m_Method) - { - - case IVIM_D_THEN_DSTAR: + switch(m_Method) { - 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); - 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_vector< double > 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]; - // 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]; - } + if(m_FitDStar) + { + MeasAndBvals input2 = ApplyS0Threshold(m_Snap.meas, 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; + break; } - case IVIM_DSTAR_FIX: + case IVIM_DSTAR_FIX: { - MeasAndBvals input = ApplyS0Threshold(m_Snap.meas, m_Snap.bvalues); - m_Snap.bvals1 = input.bvals; - m_Snap.meas1 = input.meas; - if (input.N < 2) break; + MeasAndBvals input = ApplyS0Threshold(m_Snap.meas, 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); - f_fixdstar.set_bvalues(input.bvals); - f_fixdstar.set_measurements(input.meas); + 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_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); + 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; + m_Snap.currentF = x[0]; + m_Snap.currentD = x[1]; + m_Snap.currentDStar = m_DStar; - break; + break; } - case IVIM_FIT_ALL: + case IVIM_FIT_ALL: { - MeasAndBvals input = ApplyS0Threshold(m_Snap.meas, m_Snap.bvalues); - m_Snap.bvals1 = input.bvals; - m_Snap.meas1 = input.meas; - if (input.N < 3) break; + MeasAndBvals input = ApplyS0Threshold(m_Snap.meas, 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); + 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_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); + 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]; + m_Snap.currentF = x[0]; + m_Snap.currentD = x[1]; + m_Snap.currentDStar = x[2]; - break; + break; } - case IVIM_LINEAR_D_THEN_F: + 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); - 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 + 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); + 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: + 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; + vnl_vector< double > 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() - { +template< class TIn, class TOut> +void DiffusionIntravoxelIncoherentMotionReconstructionImageFilter +::AfterThreadedGenerateData() +{ if(m_Method == IVIM_REGULARIZED) { - typename OutputImageType::Pointer outputImage = - static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); - 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) ); - - if(!m_Verbose) + typename OutputImageType::Pointer outputImage = + static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); + 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() ) { - // report the middle voxel - if( iit.GetIndex()[0] == m_CrossPosition[0] - && iit.GetIndex()[1] == m_CrossPosition[1] - && iit.GetIndex()[2] == m_CrossPosition[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]; - } - } + double f = iit.Get()[0]; + IVIM_CEIL( f, 0.0, 1.0 ); - ++oit0; - ++oit1; - ++oit2; - ++iit; - } + oit0.Set( myround(f * 100.0) ); + oit1.Set( myround(iit.Get()[1] * 10000.0) ); + oit2.Set( myround(iit.Get()[2] * 1000.0) ); + + 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] ) + { + 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) - { +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 ) - { +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 - { +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; + os << indent << "GradientDirectionContainer: " + << m_GradientDirectionContainer << std::endl; } else { - os << indent << - "GradientDirectionContainer: (Gradient directions not set)" << std::endl; + os << indent << + "GradientDirectionContainer: (Gradient directions not set)" << std::endl; } - } +} } #endif // __itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter_cpp diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/QmitkIVIMWidget.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/QmitkIVIMWidget.cpp index 2d4190365d..2eb8514478 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/QmitkIVIMWidget.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/QmitkIVIMWidget.cpp @@ -1,165 +1,167 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "QmitkIVIMWidget.h" #include "mitkHistogramGenerator.h" #include #include #include //#include #include #include #include #include #include #include #include #include #include #include #include #include #include #include //#include QmitkIVIMWidget::QmitkIVIMWidget( QWidget * parent ) : QmitkPlotWidget(parent) { // this->SetAxisTitle( QwtPlot::xBottom, "Grayvalue" ); // this->SetAxisTitle( QwtPlot::yLeft, "Probability" ); // this->Replot(); m_Plot->setCanvasLineWidth(0); m_Plot->setMargin(0); QwtLog10ScaleEngine* logScale = new QwtLog10ScaleEngine(); m_Plot->setAxisScaleEngine(0, logScale); m_Plot->setAxisScale( 0, 0.15, 1.0 ); } QmitkIVIMWidget::~QmitkIVIMWidget() { } void QmitkIVIMWidget::DrawGauss() { } void QmitkIVIMWidget::ClearItemModel() { } std::vector QmitkIVIMWidget::vec(vnl_vector vector) { std::vector retval(vector.size()); for(unsigned int i=0; iClear(); + if (snap.bvalues.empty()) + return; QString s("f=%1, D=%2, D*=%3"); s = s.arg(snap.currentF,4); s = s.arg(snap.currentD,4); s = s.arg(snap.currentDStar,4); int curveId = this->InsertCurve( s.toAscii() ); this->SetCurvePen( curveId, QPen( Qt::NoPen ) ); curveId = this->InsertCurve( "ignored measurement points" ); this->SetCurveData( curveId, vec(snap.bvalues), vec(snap.allmeas) ); this->SetCurvePen( curveId, QPen( Qt::NoPen ) ); QwtSymbol whiteSymbol(QwtSymbol::Diamond, QColor(Qt::white), QColor(Qt::black), QSize(10,10)); this->SetCurveSymbol(curveId, &whiteSymbol); if(snap.currentDStar != 0) { curveId = this->InsertCurve( "additional points second fit" ); this->SetCurveData( curveId, vec(snap.bvals2), vec(snap.meas2) ); this->SetCurvePen( curveId, QPen( Qt::NoPen ) ); QwtSymbol blackSymbol(QwtSymbol::Diamond, QColor(Qt::black), QColor(Qt::black), QSize(10,10)); this->SetCurveSymbol(curveId, &blackSymbol); } curveId = this->InsertCurve( "points first fit" ); this->SetCurveData( curveId, vec(snap.bvals1), vec(snap.meas1) ); this->SetCurvePen( curveId, QPen( Qt::NoPen ) ); QwtSymbol redSymbol(QwtSymbol::Diamond, QColor(Qt::red), QColor(Qt::red), QSize(10,10)); this->SetCurveSymbol(curveId, &redSymbol); QPen pen; pen.setColor( QColor(Qt::red) ); pen.setWidth(2); double maxb = snap.bvalues.max_value(); vnl_vector xvals(2); vnl_vector yvals(2); xvals[0] = 0; xvals[1] = maxb; yvals[0] = 1-snap.currentFunceiled; yvals[1] = yvals[0]*exp(-maxb * snap.currentD); curveId = this->InsertCurve( "contribution of D to the signal" ); this->SetCurveData( curveId, vec(xvals), vec(yvals) ); this->SetCurvePen( curveId, pen ); if(snap.currentDStar != 0) { pen.setColor(Qt::black); int nsampling = 50; xvals.set_size(nsampling); yvals.set_size(nsampling); double f = 1-snap.currentFunceiled; for(int i=0; iInsertCurve( "resulting fit of the model" ); this->SetCurveData( curveId, vec(xvals), vec(yvals) ); this->SetCurvePen( curveId, pen ); } // QMargins margins; // margins.setBottom(0); // margins.setLeft(0); // margins.setRight(0); // margins.setTop(0); QwtLegend* legend = new QwtLegend(); // legend->setContentsMargins(margins); m_Plot->insertLegend(legend, QwtPlot::BottomLegend); this->Replot(); } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkIVIMView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkIVIMView.cpp index 1ac4a8a560..dedf07d73f 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkIVIMView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkIVIMView.cpp @@ -1,816 +1,810 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ // Blueberry #include #include // Qmitk #include "QmitkIVIMView.h" #include "QmitkStdMultiWidget.h" // qt #include "qmessagebox.h" #include "qclipboard.h" // mitk #include "mitkDiffusionImage.h" #include "mitkImageCast.h" // itk #include "itkScalarImageToHistogramGenerator.h" #include "itkRegionOfInterestImageFilter.h" #include "itkImageRegionConstIteratorWithIndex.h" // itk/mitk #include "itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h" #include "itkRegularizedIVIMReconstructionFilter.h" #include "mitkImageCast.h" const std::string QmitkIVIMView::VIEW_ID = "org.mitk.views.ivim"; QmitkIVIMView::QmitkIVIMView() : QmitkFunctionality() , m_Controls( 0 ) , m_MultiWidget( NULL ) , m_Active(false) , m_SliceObserverTag1(0), m_SliceObserverTag2(0), m_SliceObserverTag3(0) , m_DiffusionImageNode(NULL) , m_MaskImageNode(NULL) { } QmitkIVIMView::~QmitkIVIMView() { // QmitkStdMultiWidget* MultiWidget = this->GetActiveStdMultiWidget(false); // if(MultiWidget) // { // //unregister observers when view is destroyed // if( MultiWidget->mitkWidget1 != NULL && m_SliceObserverTag1 != 0) // { // mitk::SliceNavigationController* slicer = MultiWidget->mitkWidget1->GetSliceNavigationController(); // slicer->RemoveObserver( m_SliceObserverTag1 ); // } // if( MultiWidget->mitkWidget2 != NULL && m_SliceObserverTag2 != 0) // { // mitk::SliceNavigationController* slicer = MultiWidget->mitkWidget2->GetSliceNavigationController(); // slicer->RemoveObserver( m_SliceObserverTag2 ); // } // if( MultiWidget->mitkWidget3!= NULL && m_SliceObserverTag3 != 0) // { // mitk::SliceNavigationController* slicer = MultiWidget->mitkWidget3->GetSliceNavigationController(); // slicer->RemoveObserver( m_SliceObserverTag3 ); // } // } } void QmitkIVIMView::CreateQtPartControl( QWidget *parent ) { // build up qt view, unless already done if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::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(MethodCombo(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_ChooseMethod, SIGNAL(clicked()), this, SLOT(ChooseMethod()) ); connect( m_Controls->m_CurveClipboard, SIGNAL(clicked()), this, SLOT(ClipboardCurveButtonClicked()) ); connect( m_Controls->m_ValuesClipboard, SIGNAL(clicked()), this, SLOT(ClipboardStatisticsButtonClicked()) ); } 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_MethodCombo->setVisible(m_Controls->m_ChooseMethod->isChecked()); m_Controls->m_Warning->setVisible(false); MethodCombo(m_Controls->m_MethodCombo->currentIndex()); } void QmitkIVIMView::Checkbox() { itk::StartEvent dummy; OnSliceChanged(dummy); } void QmitkIVIMView::MethodCombo(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; } itk::StartEvent dummy; OnSliceChanged(dummy); } void QmitkIVIMView::DStarSlider (int val) { QString sval = QString::number(val/1000.0); m_Controls->m_DStarLabel->setText(sval); itk::StartEvent dummy; OnSliceChanged(dummy); } void QmitkIVIMView::BThreshSlider (int val) { QString sval = QString::number(val*5.0); m_Controls->m_BThreshLabel->setText(sval); itk::StartEvent dummy; OnSliceChanged(dummy); } void QmitkIVIMView::S0ThreshSlider (int val) { QString sval = QString::number(val*0.5); m_Controls->m_S0ThreshLabel->setText(sval); itk::StartEvent dummy; OnSliceChanged(dummy); } void QmitkIVIMView::NumItsSlider (int val) { QString sval = QString::number(val); m_Controls->m_NumItsLabel->setText(sval); itk::StartEvent dummy; OnSliceChanged(dummy); } void QmitkIVIMView::LambdaSlider (int val) { QString sval = QString::number(val*.00001); m_Controls->m_LambdaLabel->setText(sval); itk::StartEvent dummy; OnSliceChanged(dummy); } void QmitkIVIMView::StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget) { m_MultiWidget = &stdMultiWidget; { mitk::SliceNavigationController* slicer = m_MultiWidget->mitkWidget1->GetSliceNavigationController(); itk::ReceptorMemberCommand::Pointer command = itk::ReceptorMemberCommand::New(); command->SetCallbackFunction( this, &QmitkIVIMView::OnSliceChanged ); m_SliceObserverTag1 = slicer->AddObserver( mitk::SliceNavigationController::GeometrySliceEvent(NULL, 0), command ); } { mitk::SliceNavigationController* slicer = m_MultiWidget->mitkWidget2->GetSliceNavigationController(); itk::ReceptorMemberCommand::Pointer command = itk::ReceptorMemberCommand::New(); command->SetCallbackFunction( this, &QmitkIVIMView::OnSliceChanged ); m_SliceObserverTag2 = slicer->AddObserver( mitk::SliceNavigationController::GeometrySliceEvent(NULL, 0), command ); } { mitk::SliceNavigationController* slicer = m_MultiWidget->mitkWidget3->GetSliceNavigationController(); itk::ReceptorMemberCommand::Pointer command = itk::ReceptorMemberCommand::New(); command->SetCallbackFunction( this, &QmitkIVIMView::OnSliceChanged ); m_SliceObserverTag3 = slicer->AddObserver( mitk::SliceNavigationController::GeometrySliceEvent(NULL, 0), command ); } } void QmitkIVIMView::StdMultiWidgetNotAvailable() { { mitk::SliceNavigationController* slicer = m_MultiWidget->mitkWidget1->GetSliceNavigationController(); slicer->RemoveObserver( m_SliceObserverTag1 ); } { mitk::SliceNavigationController* slicer = m_MultiWidget->mitkWidget2->GetSliceNavigationController(); slicer->RemoveObserver( m_SliceObserverTag2 ); } { mitk::SliceNavigationController* slicer = m_MultiWidget->mitkWidget3->GetSliceNavigationController(); slicer->RemoveObserver( m_SliceObserverTag3 ); } m_MultiWidget = NULL; } void QmitkIVIMView::OnSelectionChanged( std::vector nodes ) { bool foundOneDiffusionImage = false; m_Controls->m_InputData->setTitle("Please Select Input Data"); m_Controls->m_DiffusionImageLabel->setText("mandatory"); m_Controls->m_MaskImageLabel->setText("optional"); m_MaskImageNode = NULL; m_DiffusionImageNode = NULL; // iterate all selected objects, adjust warning visibility for( std::vector::iterator it = nodes.begin(); it != nodes.end(); ++it ) { mitk::DataNode::Pointer node = *it; if( node.IsNotNull() && dynamic_cast(node->GetData()) ) { if( dynamic_cast*>(node->GetData()) ) { m_DiffusionImageNode = node; foundOneDiffusionImage = true; m_Controls->m_DiffusionImageLabel->setText(node->GetName().c_str()); } else { bool isBinary = false; node->GetPropertyValue("binary", isBinary); if (isBinary) { m_MaskImageNode = node; m_Controls->m_MaskImageLabel->setText(node->GetName().c_str()); } } } } if (m_DiffusionImageNode.IsNotNull()) { m_Controls->m_VisualizeResultsWidget->setVisible(true); m_Controls->m_InputData->setTitle("Input Data"); } else m_Controls->m_VisualizeResultsWidget->setVisible(false); m_Controls->m_ButtonStart->setEnabled( foundOneDiffusionImage ); m_Controls->m_ButtonAutoThres->setEnabled( foundOneDiffusionImage ); m_Controls->m_ControlsFrame->setEnabled( foundOneDiffusionImage ); m_Controls->m_BottomControlsFrame->setEnabled( foundOneDiffusionImage ); itk::StartEvent dummy; OnSliceChanged(dummy); } void QmitkIVIMView::AutoThreshold() { std::vector nodes = this->GetDataManagerSelection(); if (nodes.empty()) return; if (!nodes.front()) { // Nothing selected. Inform the user and return QMessageBox::information( NULL, "Template", "Please load and select a diffusion image before starting image processing."); return; } typedef mitk::DiffusionImage DiffImgType; DiffImgType* dimg = dynamic_cast(nodes.front()->GetData()); if (!dimg) { // Nothing selected. Inform the user and return QMessageBox::information( NULL, "Template", "No valid diffusion image was found."); return; } // find bzero index int index = -1; DiffImgType::GradientDirectionContainerType::Pointer directions = dimg->GetDirections(); for(DiffImgType::GradientDirectionContainerType::ConstIterator it = directions->Begin(); it != directions->End(); ++it) { index++; DiffImgType::GradientDirectionType g = it.Value(); if(g[0] == 0 && g[1] == 0 && g[2] == 0 ) break; } typedef itk::VectorImage VecImgType; VecImgType::Pointer vecimg = dimg->GetVectorImage(); 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 = itw.Begin(); itk::ImageRegionConstIterator itr (vecimg, vecimg->GetLargestPossibleRegion() ); itr = itr.Begin(); 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->GetScalarValueMin() ); histogramGenerator->SetHistogramMax( dimg->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() { std::vector nodes = this->GetDataManagerSelection(); if (nodes.empty()) return; if (!nodes.front()) { // Nothing selected. Inform the user and return QMessageBox::information( NULL, "Template", "Please load and select a diffusion image before starting image processing."); return; } mitk::DiffusionImage* img = - dynamic_cast*>( + dynamic_cast*>( nodes.front()->GetData()); if (!img) { // Nothing selected. Inform the user and return QMessageBox::information( NULL, "Template", "No valid diffusion image was found."); return; } typedef itk::VectorImage VecImgType; VecImgType::Pointer vecimg = img->GetVectorImage(); OutImgType::IndexType dummy; + FittIVIM(vecimg, img->GetDirections(), img->GetB_Value(), true, dummy); OutputToDatastorage(nodes); } void QmitkIVIMView::OnSliceChanged(const itk::EventObject& /*e*/) { m_Controls->m_Warning->setVisible(false); if(!m_Controls || m_DiffusionImageNode.IsNull()) return; m_Controls->m_VisualizeResultsWidget->setVisible(false); mitk::DiffusionImage::Pointer diffusionImg = dynamic_cast*>(m_DiffusionImageNode->GetData()); mitk::Image::Pointer maskImg = NULL; if (m_MaskImageNode.IsNotNull()) maskImg = dynamic_cast(m_MaskImageNode->GetData()); - IVIMFilterType::GradientDirectionContainerType::ConstIterator gdcit = - diffusionImg->GetDirections()->Begin(); - bool foundB0 = false; - while( gdcit != diffusionImg->GetDirections()->End() ) - { - if(gdcit.Value().one_norm() <= 0.0) - foundB0 = true; - ++gdcit; - } - if(!foundB0) - { - m_Controls->m_Warning->setText(QString("No baseline (non diffusion-weighted) image found.. aborting:(")); - m_Controls->m_Warning->setVisible(true); - } - else - { - m_Controls->m_Warning->setVisible(false); - } - if (!m_MultiWidget) return; - m_Controls->m_VisualizeResultsWidget->setVisible(true); - typedef itk::VectorImage VecImgType; VecImgType::Pointer vecimg = (VecImgType*)diffusionImg->GetVectorImage().GetPointer(); VecImgType::Pointer roiImage = VecImgType::New(); + + bool success = false; if(maskImg.IsNull()) { int roisize = 0; if(m_Controls->m_MethodCombo->currentIndex() == 4) roisize = 5; mitk::Point3D pos = m_MultiWidget->GetCrossPosition(); VecImgType::IndexType crosspos; diffusionImg->GetTimeSlicedGeometry()->WorldToIndex(pos, crosspos); 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)); - FittIVIM(roiImage, diffusionImg->GetDirections(), diffusionImg->GetB_Value(), false, crosspos); + success = FittIVIM(roiImage, diffusionImg->GetDirections(), diffusionImg->GetB_Value(), 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); - FittIVIM(roiImage, diffusionImg->GetDirections(), diffusionImg->GetB_Value(), false, index); + success = FittIVIM(roiImage, diffusionImg->GetDirections(), diffusionImg->GetB_Value(), false, index); } vecimg->SetRegions( vecimg->GetLargestPossibleRegion() ); - m_Controls->m_VisualizeResultsWidget->SetParameters(m_Snap); - + if (success) + { + m_Controls->m_VisualizeResultsWidget->setVisible(true); + m_Controls->m_VisualizeResultsWidget->SetParameters(m_Snap); + } } -void QmitkIVIMView::FittIVIM(itk::VectorImage* vecimg, DirContainerType* dirs, float bval, bool multivoxel, OutImgType::IndexType &crosspos) +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); + m_Controls->m_Warning->setVisible(false); 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); - filter->Update(); - - m_Snap = filter->GetSnapshot(); - m_DStarMap = filter->GetOutput(2); - m_DMap = filter->GetOutput(1); - m_fMap = filter->GetOutput(0); + try{ + filter->Update(); + m_Snap = filter->GetSnapshot(); + m_DStarMap = filter->GetOutput(2); + m_DMap = filter->GetOutput(1); + m_fMap = filter->GetOutput(0); + } + 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(std::vector nodes) { // Outputs to Datastorage QString basename(nodes.front()->GetName().c_str()); 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 = basename; newname2 = newname2.append("_DStarMap_%1").arg(m_Controls->m_MethodCombo->currentText()); mitk::DataNode::Pointer node2=mitk::DataNode::New(); node2->SetData( dstarimage ); node2->SetName(newname2.toAscii()); GetDefaultDataStorage()->Add(node2); } 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 = basename; newname1 = newname1.append("_DMap_%1").arg(m_Controls->m_MethodCombo->currentText()); mitk::DataNode::Pointer node1=mitk::DataNode::New(); node1->SetData( dimage ); node1->SetName(newname1.toAscii()); GetDefaultDataStorage()->Add(node1); } 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 = basename; newname0 = newname0.append("_fMap_%1").arg(m_Controls->m_MethodCombo->currentText()); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( image ); node->SetName(newname0.toAscii()); GetDefaultDataStorage()->Add(node); } m_MultiWidget->RequestUpdate(); } void QmitkIVIMView::ChooseMethod() { m_Controls->m_MethodCombo->setVisible(m_Controls->m_ChooseMethod->isChecked()); } void QmitkIVIMView::ClipboardCurveButtonClicked() { if(true) { QString clipboard("Measurement Points\n"); for ( int i=0; isetText( clipboard, QClipboard::Clipboard ); } else { QApplication::clipboard()->clear(); } } void QmitkIVIMView::ClipboardStatisticsButtonClicked() { if ( true ) { QString clipboard( "f \t D \t D* \n" ); clipboard = clipboard.append( "%L1 \t %L2 \t %L3" ) .arg( m_Snap.currentF, 0, 'f', 10 ) .arg( m_Snap.currentD, 0, 'f', 10 ) .arg( m_Snap.currentDStar, 0, 'f', 10 ) ; QApplication::clipboard()->setText( clipboard, QClipboard::Clipboard ); } else { QApplication::clipboard()->clear(); } } void QmitkIVIMView::Activated() { m_Active = true; } void QmitkIVIMView::Deactivated() { m_Active = false; } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkIVIMView.h b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkIVIMView.h index 92bd30dc0e..7d06e23f65 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkIVIMView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkIVIMView.h @@ -1,115 +1,115 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef _QMITKIVIMVIEW_H_INCLUDED #define _QMITKIVIMVIEW_H_INCLUDED #include #include #include "ui_QmitkIVIMViewControls.h" #include "itkVectorImage.h" #include "itkImage.h" #include "mitkDiffusionImage.h" #include "itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h" /*! \brief QmitkIVIMView \warning This application module is not yet documented. Use "svn blame/praise/annotate" and ask the author to provide basic documentation. \sa QmitkFunctionality \ingroup Functionalities */ class QmitkIVIMView : public QmitkFunctionality { // this is needed for all Qt objects that should have a Qt meta-object // (everything that derives from QObject and wants to have signal/slots) Q_OBJECT public: static const std::string VIEW_ID; QmitkIVIMView(); virtual ~QmitkIVIMView(); typedef mitk::DiffusionImage::GradientDirectionContainerType DirContainerType; typedef itk::DiffusionIntravoxelIncoherentMotionReconstructionImageFilter IVIMFilterType; typedef itk::Image OutImgType; virtual void CreateQtPartControl(QWidget *parent); virtual void StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget); virtual void StdMultiWidgetNotAvailable(); void OnSliceChanged(const itk::EventObject& e); void OutputToDatastorage(std::vector nodes); - void FittIVIM(itk::VectorImage* vecimg, DirContainerType* dirs, float bval, bool multivoxel, OutImgType::IndexType &crosspos); + bool FittIVIM(itk::VectorImage* vecimg, DirContainerType* dirs, float bval, bool multivoxel, OutImgType::IndexType &crosspos); void Activated(); void Deactivated(); protected slots: /// \brief Called when the user clicks the GUI button void FittIVIMStart(); void AutoThreshold(); void MethodCombo(int val); void Checkbox(); void DStarSlider(int val); void BThreshSlider(int val); void S0ThreshSlider(int val); void NumItsSlider(int val); void LambdaSlider(int val); void ChooseMethod(); void ClipboardStatisticsButtonClicked(); void ClipboardCurveButtonClicked(); protected: /// \brief called by QmitkFunctionality when DataManager's selection has changed virtual void OnSelectionChanged( std::vector nodes ); Ui::QmitkIVIMViewControls* m_Controls; QmitkStdMultiWidget* m_MultiWidget; int m_SliceObserverTag1; int m_SliceObserverTag2; int m_SliceObserverTag3; OutImgType::Pointer m_DStarMap; OutImgType::Pointer m_DMap; OutImgType::Pointer m_fMap; IVIMFilterType::IVIMSnapshot m_Snap; mitk::DataNode::Pointer m_DiffusionImageNode; mitk::DataNode::Pointer m_MaskImageNode; bool m_Active; }; #endif // _QMITKIVIMVIEW_H_INCLUDED diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkIVIMViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkIVIMViewControls.ui index f8f6d9faac..2e17470388 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkIVIMViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkIVIMViewControls.ui @@ -1,816 +1,819 @@ QmitkIVIMViewControls 0 0 360 928 0 0 QmitkTemplate 9 Intra Voxel Incoherent Motion Estimation 0 9 0 0 Please Select Input Data DWI to analyze Raw DWI: DWI to analyze <html><head/><body><p><span style=" color:#ff0000;">mandatory</span></p></body></html> Optional ROI image ROI: Optional ROI image <html><head/><body><p><span style=" color:#969696;">optional</span></p></body></html> Parameters 9 QFrame::NoFrame QFrame::Raised 0 0 80 16777215 D* 100 60 Qt::Horizontal 51 16777215 200 Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter QFrame::NoFrame QFrame::Raised 0 0 80 16777215 neglect b< 250 34 Qt::Horizontal 51 16777215 46.5 Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter QFrame::NoFrame QFrame::Plain 0 0 80 16777215 neglect Si< 100 0 Qt::Horizontal 0 0 QFrame::NoFrame QFrame::Raised 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 80 16777215 #iterations 100 10 Qt::Horizontal 30 16777215 TextLabel Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter QFrame::NoFrame QFrame::Raised 0 0 80 16777215 lambda 1000 10 Qt::Horizontal 0 0 QFrame::NoFrame QFrame::Raised 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 QFrame::NoFrame QFrame::Raised 0 QFrame::NoFrame QFrame::Raised 0 QFrame::NoFrame QFrame::Raised 0 QFrame::NoFrame QFrame::Raised 0 80 0 Output Images f true D false D* false Generate Output Images warning display Qt::RichText + + true + QFrame::StyledPanel QFrame::Raised 0 0 true 0 0 0 400 0 QFrame::NoFrame QFrame::Raised 0 QFrame::NoFrame QFrame::Raised 0 QFrame::NoFrame QFrame::Raised 0 Datapoints to Clipboard Parameters to Clipboard Choose Method 2 3 Param. Fit Fit D & f with fixed D* value Fit D & f (high b), then fit D* Linearly fit D & f (high b), then fit D* Regularized Qt::Vertical QSizePolicy::Expanding 20 220 QmitkIVIMWidget QWidget
QmitkIVIMWidget.h
1