diff --git a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp index d3eb88f170..5efbcdb46a 100644 --- a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp @@ -1,655 +1,655 @@ /*=================================================================== 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 __itkAnalyticalDiffusionQballReconstructionImageFilter_cpp #define __itkAnalyticalDiffusionQballReconstructionImageFilter_cpp #include #include #include #include #include #include #include #include #include #include #include "itkPointShell.h" namespace itk { template< class T, class TG, class TO, int ShOrder, int NrOdfDirections> AnalyticalDiffusionQballReconstructionImageFilter::AnalyticalDiffusionQballReconstructionImageFilter() : m_GradientDirectionContainer(nullptr), m_NumberOfGradientDirections(0), m_NumberOfBaselineImages(1), m_Threshold(NumericTraits< ReferencePixelType >::NonpositiveMin()), m_BValue(-1), m_Lambda(0.0), m_DirectionsDuplicated(false), m_Delta1(0.001), m_Delta2(0.001), m_UseMrtrixBasis(false) { // At least 1 inputs is necessary for a vector image. // For images added one at a time we need at least six this->SetNumberOfRequiredInputs( 1 ); } template< class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int ShOrder, int NrOdfDirections> typename itk::AnalyticalDiffusionQballReconstructionImageFilter< TReferenceImagePixelType,TGradientImagePixelType,TOdfPixelType, ShOrder,NrOdfDirections>::OdfPixelType itk::AnalyticalDiffusionQballReconstructionImageFilter::Normalize( OdfPixelType odf, typename NumericTraits::AccumulateType b0 ) { switch( m_NormalizationMethod ) { case QBAR_STANDARD: { TOdfPixelType sum = 0; for(int i=0; i0) odf /= sum; return odf; break; } case QBAR_B_ZERO_B_VALUE: { for(int i=0; i vnl_vector itk::AnalyticalDiffusionQballReconstructionImageFilter::PreNormalize( vnl_vector vec, typename NumericTraits::AccumulateType b0 ) { switch( m_NormalizationMethod ) { case QBAR_STANDARD: { int n = vec.size(); double b0f = (double)b0; for(int i=0; i=1) vec[i] = 1-m_Delta2/2; else if (vec[i]>=1-m_Delta2) vec[i] = 1-m_Delta2/2-(1-vec[i])*(1-vec[i])/(2*m_Delta2); vec[i] = log(-log(vec[i])); } return vec; break; } } return vec; } template< class T, class TG, class TO, int ShOrder, int NrOdfDirections> void AnalyticalDiffusionQballReconstructionImageFilter::BeforeThreadedGenerateData() { // If we have more than 2 inputs, then each input, except the first is a // gradient image. The number of gradient images must match the number of // gradient directions. //const unsigned int numberOfInputs = this->GetNumberOfInputs(); // There need to be at least 6 gradient directions to be able to compute the // tensor basis if( m_NumberOfGradientDirections < (ShOrder*ShOrder + ShOrder + 2)/2 + ShOrder ) { itkExceptionMacro( << "Not enough gradient directions supplied (" << m_NumberOfGradientDirections << "). At least " << (ShOrder*ShOrder + ShOrder + 2)/2 + ShOrder << " needed for SH-order " << ShOrder); } // Input must be an itk::VectorImage. 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 ); } this->ComputeReconstructionMatrix(); typename GradientImagesType::Pointer img = static_cast< GradientImagesType * >( this->ProcessObject::GetInput(0) ); m_BZeroImage = BZeroImageType::New(); m_BZeroImage->SetSpacing( img->GetSpacing() ); // Set the image spacing m_BZeroImage->SetOrigin( img->GetOrigin() ); // Set the image origin m_BZeroImage->SetDirection( img->GetDirection() ); // Set the image direction m_BZeroImage->SetLargestPossibleRegion( img->GetLargestPossibleRegion()); m_BZeroImage->SetBufferedRegion( img->GetLargestPossibleRegion() ); m_BZeroImage->Allocate(); m_ODFSumImage = BZeroImageType::New(); m_ODFSumImage->SetSpacing( img->GetSpacing() ); // Set the image spacing m_ODFSumImage->SetOrigin( img->GetOrigin() ); // Set the image origin m_ODFSumImage->SetDirection( img->GetDirection() ); // Set the image direction m_ODFSumImage->SetLargestPossibleRegion( img->GetLargestPossibleRegion()); m_ODFSumImage->SetBufferedRegion( img->GetLargestPossibleRegion() ); m_ODFSumImage->Allocate(); m_CoefficientImage = CoefficientImageType::New(); m_CoefficientImage->SetSpacing( img->GetSpacing() ); // Set the image spacing m_CoefficientImage->SetOrigin( img->GetOrigin() ); // Set the image origin m_CoefficientImage->SetDirection( img->GetDirection() ); // Set the image direction m_CoefficientImage->SetLargestPossibleRegion( img->GetLargestPossibleRegion()); m_CoefficientImage->SetBufferedRegion( img->GetLargestPossibleRegion() ); m_CoefficientImage->Allocate(); if(m_NormalizationMethod == QBAR_SOLID_ANGLE || m_NormalizationMethod == QBAR_NONNEG_SOLID_ANGLE) m_Lambda = 0.0; } template< class T, class TG, class TO, int ShOrder, int NrOdfDirections> void AnalyticalDiffusionQballReconstructionImageFilter ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, ThreadIdType ) { typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetPrimaryOutput()); ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread); oit.GoToBegin(); ImageRegionIterator< BZeroImageType > oit2(m_BZeroImage, outputRegionForThread); oit2.GoToBegin(); ImageRegionIterator< FloatImageType > oit3(m_ODFSumImage, outputRegionForThread); oit3.GoToBegin(); ImageRegionIterator< CoefficientImageType > oit4(m_CoefficientImage, outputRegionForThread); oit4.GoToBegin(); typedef ImageRegionConstIterator< GradientImagesType > GradientIteratorType; typedef typename GradientImagesType::PixelType GradientVectorType; typename GradientImagesType::Pointer gradientImagePointer = nullptr; // Would have liked a dynamic_cast here, but seems SGI doesn't like it // The enum will ensure that an inappropriate cast is not done gradientImagePointer = static_cast< GradientImagesType * >( this->ProcessObject::GetInput(0) ); GradientIteratorType git(gradientImagePointer, outputRegionForThread ); git.GoToBegin(); // Compute the indicies of the baseline images and gradient images std::vector baselineind; // contains the indicies of // the baseline images std::vector gradientind; // contains the indicies of // the gradient images for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin(); gdcit != this->m_GradientDirectionContainer->End(); ++gdcit) { float bval = gdcit.Value().two_norm(); bval = bval*bval*m_BValue; if(bval < 100) baselineind.push_back(gdcit.Index()); else gradientind.push_back(gdcit.Index()); } if( m_DirectionsDuplicated ) { int gradIndSize = gradientind.size(); for(int i=0; i::AccumulateType b0 = NumericTraits::Zero; // Average the baseline image pixels for(unsigned int i = 0; i < baselineind.size(); ++i) { b0 += b[baselineind[i]]; } b0 /= this->m_NumberOfBaselineImages; OdfPixelType odf(0.0); typename CoefficientImageType::PixelType coeffPixel(0.0); vnl_vector B(m_NumberOfGradientDirections); if( (b0 != 0) && (b0 >= m_Threshold) ) { for( unsigned int i = 0; i< m_NumberOfGradientDirections; i++ ) { B[i] = static_cast(b[gradientind[i]]); } B = PreNormalize(B, b0); if(m_NormalizationMethod == QBAR_SOLID_ANGLE) { vnl_vector coeffs(m_NumberCoefficients); coeffs = ( m_CoeffReconstructionMatrix * B ); coeffs[0] += 1.0/(2.0*sqrt(itk::Math::pi)); odf = ( m_SphericalHarmonicBasisMatrix * coeffs ).data_block(); coeffPixel = coeffs.data_block(); } else if(m_NormalizationMethod == QBAR_NONNEG_SOLID_ANGLE) { /** this would be the place to implement a non-negative * solver for quadratic programming problem: * min .5*|| Bc-s ||^2 subject to -CLPc <= 4*pi*ones * (refer to MICCAI 2009 Goh et al. "Estimating ODFs with PDF constraints") * .5*|| Bc-s ||^2 == .5*c'B'Bc - x'B's + .5*s's */ itkExceptionMacro( << "Nonnegative Solid Angle not yet implemented"); } else { vnl_vector coeffs(m_NumberCoefficients); coeffs = ( m_CoeffReconstructionMatrix * B ); coeffs[0] += 1.0/(2.0*sqrt(itk::Math::pi)); coeffPixel = coeffs.data_block(); odf = ( m_ReconstructionMatrix * B ).data_block(); } odf = Normalize(odf, b0); } oit.Set( odf ); oit2.Set( b0 ); float sum = 0; for (unsigned int k=0; k void AnalyticalDiffusionQballReconstructionImageFilter ::tofile2(vnl_matrix *pA, std::string fname) { vnl_matrix A = (*pA); std::ofstream myfile; std::locale C("C"); std::locale originalLocale = myfile.getloc(); myfile.imbue(C); myfile.open (fname.c_str()); myfile << "A1=["; for(unsigned int i=0; i void AnalyticalDiffusionQballReconstructionImageFilter::ComputeReconstructionMatrix() { m_NumberCoefficients = (ShOrder*ShOrder + ShOrder + 2)/2 + ShOrder; if( m_NumberOfGradientDirections < m_NumberCoefficients ) { itkExceptionMacro( << "Not enough gradient directions supplied (" << m_NumberOfGradientDirections << "). At least " << (ShOrder*ShOrder + ShOrder + 2)/2 + ShOrder << " needed for SH-order " << ShOrder); } // Gradient preprocessing { // check for duplicate diffusion gradients bool warning = false; for(GradientDirectionContainerType::ConstIterator gdcit1 = this->m_GradientDirectionContainer->Begin(); gdcit1 != this->m_GradientDirectionContainer->End(); ++gdcit1) { for(GradientDirectionContainerType::ConstIterator gdcit2 = this->m_GradientDirectionContainer->Begin(); gdcit2 != this->m_GradientDirectionContainer->End(); ++gdcit2) { if(gdcit1.Value() == gdcit2.Value() && gdcit1.Index() != gdcit2.Index()) { itkWarningMacro( << "Some of the Diffusion Gradients equal each other. Corresponding image data should be averaged before calling this filter." ); warning = true; break; } } if (warning) break; } // handle acquisition schemes where only half of the spherical // shell is sampled by the gradient directions. In this case, // each gradient direction is duplicated in negative direction. vnl_vector centerMass(3); centerMass.fill(0.0); int count = 0; for(GradientDirectionContainerType::ConstIterator gdcit1 = this->m_GradientDirectionContainer->Begin(); gdcit1 != this->m_GradientDirectionContainer->End(); ++gdcit1) { float bval = gdcit1.Value().two_norm(); bval = bval*bval*m_BValue; if(bval > 100) { centerMass += gdcit1.Value(); count ++; } } centerMass /= count; if(centerMass.two_norm() > 0.1) { m_DirectionsDuplicated = true; m_NumberOfGradientDirections *= 2; } } // Create 3xM matrix Q that contains the gradient vectors in polar coordinates vnl_matrix Q(3, m_NumberOfGradientDirections); { int i = 0; for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin(); gdcit != this->m_GradientDirectionContainer->End(); ++gdcit) { float bval = gdcit.Value().two_norm(); bval = bval*bval*m_BValue; if(bval > 100) { double x = gdcit.Value().get(0); double y = gdcit.Value().get(1); double z = gdcit.Value().get(2); double cart[3]; mitk::sh::Cart2Sph(x,y,z,cart); Q(0,i) = cart[0]; Q(1,i) = cart[1]; Q(2,i++) = cart[2]; } } if(m_DirectionsDuplicated) { for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin(); gdcit != this->m_GradientDirectionContainer->End(); ++gdcit) { float bval = gdcit.Value().two_norm(); bval = bval*bval*m_BValue; if(bval > 100) { double x = gdcit.Value().get(0); double y = gdcit.Value().get(1); double z = gdcit.Value().get(2); double cart[3]; mitk::sh::Cart2Sph(x,y,z,cart); Q(0,i) = cart[0]; Q(1,i) = cart[1]; Q(2,i++) = cart[2]; } } } } // Calcualte SH basis B vnl_matrix L(m_NumberCoefficients,m_NumberCoefficients, 0); vnl_vector lj(m_NumberCoefficients); vnl_matrix B(m_NumberOfGradientDirections,m_NumberCoefficients); for(unsigned int i=0; i P(m_NumberCoefficients,m_NumberCoefficients, 0); for(unsigned int i=0; i B_transpose(B.transpose()); vnl_matrix_inverse* pseudoInverse = new vnl_matrix_inverse( B_transpose * B + L*L*m_Lambda ); m_CoeffReconstructionMatrix = pseudoInverse->pinverse() * B_transpose; switch(m_NormalizationMethod) { case QBAR_ADC_ONLY: case QBAR_RAW_SIGNAL: break; case QBAR_STANDARD: case QBAR_B_ZERO_B_VALUE: case QBAR_B_ZERO: case QBAR_NONE: { m_CoeffReconstructionMatrix = P * m_CoeffReconstructionMatrix; break; } case QBAR_SOLID_ANGLE: { m_CoeffReconstructionMatrix = (float)(1.0/(8.0*itk::Math::pi)) * P * L * m_CoeffReconstructionMatrix; break; } case QBAR_NONNEG_SOLID_ANGLE: break; } // needed to calculate the ODF values from the SH coefficients vnl_matrix_fixed* U = itk::PointShell >::DistributePointShell(); m_SphericalHarmonicBasisMatrix = mitk::sh::CalcShBasisForDirections(ShOrder, U->as_matrix()); m_ReconstructionMatrix = m_SphericalHarmonicBasisMatrix * m_CoeffReconstructionMatrix; } template< class T, class TG, class TO, int ShOrder, int NrOdfDirections> void AnalyticalDiffusionQballReconstructionImageFilter ::SetGradientImage(const GradientDirectionContainerType *gradientDirection, const GradientImagesType *gradientImage ) { // Copy Gradient Direction Container this->m_GradientDirectionContainer = GradientDirectionContainerType::New(); for(GradientDirectionContainerType::ConstIterator it = gradientDirection->Begin(); it != gradientDirection->End(); it++) { this->m_GradientDirectionContainer->push_back(it.Value()); } if (m_BValue<0) - itkExceptionMacro("B-value needs to best before gradient image!"); + mitkThrow() << "B-value needs to be set before gradient image! " << m_BValue; unsigned int numImages = gradientDirection->Size(); this->m_NumberOfBaselineImages = 0; for(GradientDirectionContainerType::Iterator it = this->m_GradientDirectionContainer->Begin(); it != this->m_GradientDirectionContainer->End(); it++) { float bval = it.Value().two_norm(); bval = bval*bval*m_BValue; if(bval < 100) { this->m_NumberOfBaselineImages++; } else // Normalize non-zero gradient directions { it.Value() = it.Value() / it.Value().two_norm(); } } if (this->m_NumberOfBaselineImages==0) itkExceptionMacro("No baseline image detected (no b-zero image)"); this->m_NumberOfGradientDirections = numImages - this->m_NumberOfBaselineImages; // ensure that the gradient image we received has as many components as // the number of gradient directions if( gradientImage->GetVectorLength() != this->m_NumberOfBaselineImages + m_NumberOfGradientDirections ) { itkExceptionMacro( << m_NumberOfGradientDirections << " gradients + " << this->m_NumberOfBaselineImages << "baselines = " << m_NumberOfGradientDirections + this->m_NumberOfBaselineImages << " directions specified but image has " << gradientImage->GetVectorLength() << " components."); } this->ProcessObject::SetNthInput( 0, const_cast< GradientImagesType* >(gradientImage) ); } template< class T, class TG, class TO, int ShOrder, int NrOdfDirections> void AnalyticalDiffusionQballReconstructionImageFilter ::PrintSelf(std::ostream& os, Indent indent) const { std::locale C("C"); std::locale originalLocale = os.getloc(); os.imbue(C); Superclass::PrintSelf(os,indent); os << indent << "OdfReconstructionMatrix: " << m_ReconstructionMatrix << std::endl; if ( m_GradientDirectionContainer ) os << indent << "GradientDirectionContainer: " << m_GradientDirectionContainer << std::endl; else os << indent << "GradientDirectionContainer: (Gradient directions not set)" << std::endl; os << indent << "NumberOfGradientDirections: " << m_NumberOfGradientDirections << std::endl; os << indent << "NumberOfBaselineImages: " << m_NumberOfBaselineImages << std::endl; os << indent << "Threshold for reference B0 image: " << m_Threshold << std::endl; os << indent << "BValue: " << m_BValue << std::endl; os.imbue( originalLocale ); } } #endif // __itkAnalyticalDiffusionQballReconstructionImageFilter_cpp diff --git a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.h b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.h index d07decd32c..766ebbd08f 100644 --- a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.h @@ -1,269 +1,266 @@ /*=================================================================== 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 __itkAnalyticalDiffusionQballReconstructionImageFilter_h_ #define __itkAnalyticalDiffusionQballReconstructionImageFilter_h_ #include "itkImageToImageFilter.h" #include "vnl/vnl_vector_fixed.h" #include "vnl/vnl_matrix.h" #include "vnl/algo/vnl_svd.h" #include "itkVectorContainer.h" #include "itkVectorImage.h" namespace itk{ /** \class AnalyticalDiffusionQballReconstructionImageFilter * \brief This class takes as input one or more reference image (acquired in the * absence of diffusion sensitizing gradients) and 'n' diffusion * weighted images and their gradient directions and computes an image of * orientation distribution function coefficients in a spherical harmonic basis. * * \par Inputs and Usage * \par * When you have the 'n' gradient and one or more reference images in a single * multi-component image (VectorImage), you can specify the images as * \code * filter->SetGradientImage( directionsContainer, vectorImage ); * \endcode * Note that this method is used to specify both the reference and gradient images. * This is convenient when the DWI images are read in using the * NRRD * format. Like the Nrrd format, the reference images are those components of the * vectorImage whose gradient direction is (0,0,0). If more than one reference image * is present, they are averaged prior to the reconstruction. * * \par Outputs * The output image is an image of vectors that must be understood as ODFs: * \code * Image< Vector< TPixelType, OdfNrDirections >, 3 > * \endcode * * \par Parameters * \li Threshold - Threshold on the reference image data. The output ODF will * be a null pdf for pixels in the reference image that have a value less * than this. * \li BValue - See the documentation of SetBValue(). * \li At least 6 gradient images must be specified for the filter to be able * to run. If the input gradient directions g_i are majorly sampled on one half * of the sqhere, then each input image I_i will be duplicated and assign -g_i * in order to guarantee stability of the algorithm. * \li OdfDirections - directions of resulting orientation distribution function * \li EquatorNrSamplingPoints - number of sampling points on equator when * performing Funk Radeon Transform (FRT) * \li BasisFunctionCenters - the centers of the basis functions are used for * the sRBF (spherical radial basis functions interpolation). If not set, they * will be defaulted to equal m_EquatorNrSamplingPoints * * \par Template parameters * The class is templated over * \li the pixel type of the reference and gradient images * (expected to be scalar data types) * \li the internal representation of the ODF pixels (double, float etc). * \li the number of OdfDirections * \li the number of basis function centers for the sRBF * * \par References: * \li[1] * Tuch DS, * "Q-ball imaging", Magn Reson Med. 2004 Dec;52(6):1358-72. * */ template< class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int ShOrder, int NrOdfDirections> class AnalyticalDiffusionQballReconstructionImageFilter : public ImageToImageFilter< Image< TReferenceImagePixelType, 3 >, Image< Vector< TOdfPixelType, NrOdfDirections >, 3 > > { public: enum Normalization { QBAR_STANDARD, QBAR_B_ZERO_B_VALUE, QBAR_B_ZERO, QBAR_NONE, QBAR_ADC_ONLY, QBAR_RAW_SIGNAL, QBAR_SOLID_ANGLE, QBAR_NONNEG_SOLID_ANGLE }; typedef AnalyticalDiffusionQballReconstructionImageFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ImageToImageFilter< Image< TReferenceImagePixelType, 3>, Image< Vector< TOdfPixelType, NrOdfDirections >, 3 > > Superclass; /** Method for creation through the object factory. */ itkFactorylessNewMacro(Self) itkCloneMacro(Self) /** Runtime information support. */ itkTypeMacro(AnalyticalDiffusionQballReconstructionImageFilter, ImageToImageFilter) typedef TReferenceImagePixelType ReferencePixelType; typedef TGradientImagePixelType GradientPixelType; typedef Vector< TOdfPixelType, NrOdfDirections > OdfPixelType; typedef TOdfPixelType BZeroPixelType; /** Reference image data, This image is aquired in the absence * of a diffusion sensitizing field gradient */ typedef typename Superclass::InputImageType ReferenceImageType; typedef Image< OdfPixelType, 3 > OdfImageType; typedef OdfImageType OutputImageType; typedef Image< Vector< TOdfPixelType, (ShOrder*ShOrder + ShOrder + 2)/2 + ShOrder >, 3 > CoefficientImageType; typedef Image< BZeroPixelType, 3 > BZeroImageType; typedef typename Superclass::OutputImageRegionType OutputImageRegionType; /** Typedef defining one (of the many) gradient images. */ typedef Image< GradientPixelType, 3 > GradientImageType; /** An alternative typedef defining one (of the many) gradient images. * It will be assumed that the vectorImage has the same dimension as the * Reference image and a vector length parameter of \c n (number of * gradient directions)*/ typedef VectorImage< GradientPixelType, 3 > GradientImagesType; /** 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; /** 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 SetGradientImage( const GradientDirectionContainerType *, const GradientImagesType *image); /** Get reference image */ virtual ReferenceImageType * GetReferenceImage() { return ( static_cast< ReferenceImageType *>(this->ProcessObject::GetInput(0)) ); } /** 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 ); } static void tofile2(vnl_matrix *A, std::string fname); OdfPixelType Normalize(OdfPixelType odf, typename NumericTraits::AccumulateType b0 ); vnl_vector PreNormalize( vnl_vector vec, typename NumericTraits::AccumulateType b0 ); /** Threshold on the reference image data. The output ODF will be a null * pdf for pixels in the reference image that have a value less than this * threshold. */ itkSetMacro( Threshold, ReferencePixelType ) itkGetMacro( Threshold, ReferencePixelType ) itkSetMacro( NormalizationMethod, Normalization) itkGetMacro( NormalizationMethod, Normalization ) typedef Image FloatImageType; itkGetMacro( BZeroImage, typename BZeroImageType::Pointer) itkGetMacro( ODFSumImage, typename FloatImageType::Pointer) itkGetMacro( CoefficientImage, typename CoefficientImageType::Pointer) - itkSetMacro( BValue, TOdfPixelType) -#ifdef GetBValue -#undef GetBValue -#endif - itkGetConstReferenceMacro( BValue, TOdfPixelType) + itkSetMacro( BValue, float) + itkGetMacro( BValue, float) itkSetMacro( Lambda, double ) itkGetMacro( Lambda, double ) itkSetMacro( UseMrtrixBasis, bool ) #ifdef ITK_USE_CONCEPT_CHECKING /** Begin concept checking */ itkConceptMacro(ReferenceEqualityComparableCheck, (Concept::EqualityComparable)); itkConceptMacro(TensorEqualityComparableCheck, (Concept::EqualityComparable)); itkConceptMacro(GradientConvertibleToDoubleCheck, (Concept::Convertible)); itkConceptMacro(DoubleConvertibleToTensorCheck, (Concept::Convertible)); itkConceptMacro(GradientReferenceAdditiveOperatorsCheck, (Concept::AdditiveOperators)); itkConceptMacro(ReferenceOStreamWritableCheck, (Concept::OStreamWritable)); itkConceptMacro(TensorOStreamWritableCheck, (Concept::OStreamWritable)); /** End concept checking */ #endif protected: AnalyticalDiffusionQballReconstructionImageFilter(); ~AnalyticalDiffusionQballReconstructionImageFilter() override {}; void PrintSelf(std::ostream& os, Indent indent) const override; void ComputeReconstructionMatrix(); void BeforeThreadedGenerateData() override; void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType) override; private: vnl_matrix< float > m_ReconstructionMatrix; vnl_matrix< float > m_CoeffReconstructionMatrix; vnl_matrix< float > m_SphericalHarmonicBasisMatrix; /** container to hold gradient directions */ GradientDirectionContainerType::Pointer m_GradientDirectionContainer; /** Number of gradient measurements */ unsigned int m_NumberOfGradientDirections; /** Number of baseline images */ unsigned int m_NumberOfBaselineImages; /** Threshold on the reference image data */ ReferencePixelType m_Threshold; /** LeBihan's b-value for normalizing tensors */ - TOdfPixelType m_BValue; + float m_BValue; typename BZeroImageType::Pointer m_BZeroImage; double m_Lambda; bool m_DirectionsDuplicated; Normalization m_NormalizationMethod; unsigned int m_NumberCoefficients; FloatImageType::Pointer m_ODFSumImage; typename CoefficientImageType::Pointer m_CoefficientImage; TOdfPixelType m_Delta1; TOdfPixelType m_Delta2; bool m_UseMrtrixBasis; }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkAnalyticalDiffusionQballReconstructionImageFilter.cpp" #endif #endif //__itkAnalyticalDiffusionQballReconstructionImageFilter_h_ diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkQBallReconstructionView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkQBallReconstructionView.cpp index 5600269744..95f0309c19 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkQBallReconstructionView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkQBallReconstructionView.cpp @@ -1,874 +1,874 @@ /*=================================================================== 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. ===================================================================*/ //#define MBILOG_ENABLE_DEBUG #include "QmitkQBallReconstructionView.h" // qt includes #include // itk includes #include "itkTimeProbe.h" // mitk includes #include "mitkProgressBar.h" #include "mitkStatusBar.h" #include "mitkNodePredicateDataType.h" #include "QmitkDataStorageComboBox.h" #include "itkDiffusionQballReconstructionImageFilter.h" #include "itkAnalyticalDiffusionQballReconstructionImageFilter.h" #include "itkDiffusionMultiShellQballReconstructionImageFilter.h" #include "itkVectorContainer.h" #include "itkB0ImageExtractionImageFilter.h" #include #include "mitkOdfImage.h" #include "mitkProperties.h" #include "mitkVtkResliceInterpolationProperty.h" #include "mitkLookupTable.h" #include "mitkLookupTableProperty.h" #include "mitkTransferFunction.h" #include "mitkTransferFunctionProperty.h" #include "mitkDataNodeObject.h" #include "mitkOdfNormalizationMethodProperty.h" #include "mitkOdfScaleByProperty.h" #include #include "mitkDiffusionImagingConfigure.h" #include #include "berryIStructuredSelection.h" #include "berryIWorkbenchWindow.h" #include "berryISelectionService.h" #include #include #include #include const std::string QmitkQBallReconstructionView::VIEW_ID = "org.mitk.views.qballreconstruction"; typedef float TTensorPixelType; const int QmitkQBallReconstructionView::nrconvkernels = 252; struct QbrShellSelection { typedef mitk::DiffusionPropertyHelper::GradientDirectionType GradientDirectionType; typedef mitk::DiffusionPropertyHelper::GradientDirectionsContainerType GradientDirectionContainerType; typedef mitk::DiffusionPropertyHelper::BValueMapType BValueMapType; typedef itk::VectorImage< DiffusionPixelType, 3 > ITKDiffusionImageType; QmitkQBallReconstructionView* m_View; mitk::DataNode * m_Node; std::string m_NodeName; std::vector m_CheckBoxes; QLabel * m_Label; mitk::Image * m_Image; QbrShellSelection(QmitkQBallReconstructionView* view, mitk::DataNode * node) : m_View(view), m_Node(node), m_NodeName(node->GetName()) { m_Image = dynamic_cast (node->GetData()); if(!m_Image) { MITK_ERROR << "QmitkQBallReconstructionView::QbrShellSelection : no image selected"; return; } bool isDiffusionImage( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(m_Node->GetData())) ); if( !isDiffusionImage ) { MITK_ERROR << "QmitkQBallReconstructionView::QbrShellSelection : selected image contains no diffusion information"; return; } GenerateCheckboxes(); } void GenerateCheckboxes() { BValueMapType origMap = mitk::DiffusionPropertyHelper::GetBValueMap(m_Image); BValueMapType::iterator itStart = origMap.begin(); itStart++; BValueMapType::iterator itEnd = origMap.end(); m_Label = new QLabel(m_NodeName.c_str()); m_Label->setVisible(true); m_View->m_Controls->m_QBallSelectionBox->layout()->addWidget(m_Label); for(BValueMapType::iterator it = itStart ; it!= itEnd; it++) { QCheckBox * box = new QCheckBox(QString::number(it->first)); m_View->m_Controls->m_QBallSelectionBox->layout()->addWidget(box); box->setChecked(true); box->setCheckable(true); // box->setVisible(true); m_CheckBoxes.push_back(box); } } void SetVisible(bool vis) { foreach(QCheckBox * box, m_CheckBoxes) { box->setVisible(vis); } } BValueMapType GetBValueSelctionMap() { BValueMapType inputMap = mitk::DiffusionPropertyHelper::GetBValueMap(m_Image); BValueMapType outputMap; unsigned int val = 0; if(inputMap.find(0) == inputMap.end()){ MITK_INFO << "QbrShellSelection: return empty BValueMap from GUI Selection"; return outputMap; }else{ outputMap[val] = inputMap[val]; MITK_INFO << val; } foreach(QCheckBox * box, m_CheckBoxes) { if(box->isChecked()){ val = box->text().toDouble(); outputMap[val] = inputMap[val]; MITK_INFO << val; } } return outputMap; } ~QbrShellSelection() { m_View->m_Controls->m_QBallSelectionBox->layout()->removeWidget(m_Label); delete m_Label; for(std::vector::iterator it = m_CheckBoxes.begin() ; it!= m_CheckBoxes.end(); it++) { m_View->m_Controls->m_QBallSelectionBox->layout()->removeWidget((*it)); delete (*it); } m_CheckBoxes.clear(); } }; QmitkQBallReconstructionView::QmitkQBallReconstructionView() : QmitkAbstractView(), m_Controls(nullptr) { } QmitkQBallReconstructionView::~QmitkQBallReconstructionView() { } void QmitkQBallReconstructionView::CreateQtPartControl(QWidget *parent) { if (!m_Controls) { // create GUI widgets m_Controls = new Ui::QmitkQBallReconstructionViewControls; m_Controls->setupUi(parent); this->CreateConnections(); QStringList items; items << "2" << "4" << "6" << "8" << "10" << "12"; m_Controls->m_QBallReconstructionMaxLLevelComboBox->addItems(items); m_Controls->m_QBallReconstructionMaxLLevelComboBox->setCurrentIndex(1); MethodChoosen(m_Controls->m_QBallReconstructionMethodComboBox->currentIndex()); #ifndef DIFFUSION_IMAGING_EXTENDED m_Controls->m_QBallReconstructionMethodComboBox->removeItem(3); #endif } } void QmitkQBallReconstructionView::SetFocus() { } void QmitkQBallReconstructionView::CreateConnections() { if ( m_Controls ) { connect( (QObject*)(m_Controls->m_ButtonStandard), SIGNAL(clicked()), this, SLOT(ReconstructStandard()) ); connect( (QObject*)(m_Controls->m_QBallReconstructionMethodComboBox), SIGNAL(currentIndexChanged(int)), this, SLOT(MethodChoosen(int)) ); connect( (QObject*)(m_Controls->m_QBallReconstructionThreasholdEdit), SIGNAL(valueChanged(int)), this, SLOT(PreviewThreshold(int)) ); connect( (QObject*)(m_Controls->m_ConvertButton), SIGNAL(clicked()), this, SLOT(ConvertShImage()) ); m_Controls->m_ImageBox->SetDataStorage(this->GetDataStorage()); mitk::NodePredicateIsDWI::Pointer isDwi = mitk::NodePredicateIsDWI::New(); m_Controls->m_ImageBox->SetPredicate( isDwi ); m_Controls->m_ShImageBox->SetDataStorage(this->GetDataStorage()); mitk::TNodePredicateDataType::Pointer isSh = mitk::TNodePredicateDataType::New(); m_Controls->m_ShImageBox->SetPredicate( isSh ); connect( (QObject*)(m_Controls->m_ImageBox), SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); connect( (QObject*)(m_Controls->m_ShImageBox), SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); UpdateGui(); } } template void QmitkQBallReconstructionView::TemplatedConvertShImage(mitk::ShImage::Pointer mitkImage) { typedef itk::ShToOdfImageFilter< float, ShOrder > ShConverterType; typename ShConverterType::InputImageType::Pointer itkvol = ShConverterType::InputImageType::New(); mitk::CastToItkImage(mitkImage, itkvol); typename ShConverterType::Pointer converter = ShConverterType::New(); converter->SetInput(itkvol); converter->Update(); mitk::OdfImage::Pointer image = mitk::OdfImage::New(); image->InitializeByItk( converter->GetOutput() ); image->SetVolume( converter->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( image ); node->SetName(m_Controls->m_ShImageBox->GetSelectedNode()->GetName()); GetDataStorage()->Add(node, m_Controls->m_ShImageBox->GetSelectedNode()); } void QmitkQBallReconstructionView::ConvertShImage() { if (m_Controls->m_ShImageBox->GetSelectedNode().IsNotNull()) { mitk::ShImage::Pointer mitkImg = dynamic_cast(m_Controls->m_ShImageBox->GetSelectedNode()->GetData()); switch (mitkImg->ShOrder()) { case 2: TemplatedConvertShImage<2>(mitkImg); break; case 4: TemplatedConvertShImage<4>(mitkImg); break; case 6: TemplatedConvertShImage<6>(mitkImg); break; case 8: TemplatedConvertShImage<8>(mitkImg); break; case 10: TemplatedConvertShImage<10>(mitkImg); break; case 12: TemplatedConvertShImage<12>(mitkImg); break; default : QMessageBox::warning(nullptr, "Error", "Only spherical harmonics orders 2-12 are supported.", QMessageBox::Ok); } } } void QmitkQBallReconstructionView::UpdateGui() { m_Controls->m_ButtonStandard->setEnabled(false); if (m_Controls->m_ImageBox->GetSelectedNode().IsNotNull()) { m_Controls->m_ButtonStandard->setEnabled(true); GenerateShellSelectionUI(m_Controls->m_ImageBox->GetSelectedNode()); } m_Controls->m_ConvertButton->setEnabled(m_Controls->m_ShImageBox->GetSelectedNode().IsNotNull()); } void QmitkQBallReconstructionView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList& /*nodes*/) { UpdateGui(); } void QmitkQBallReconstructionView::Activated() { } void QmitkQBallReconstructionView::Deactivated() { } void QmitkQBallReconstructionView::Visible() { } void QmitkQBallReconstructionView::Hidden() { } void QmitkQBallReconstructionView::ReconstructStandard() { int index = m_Controls->m_QBallReconstructionMethodComboBox->currentIndex(); #ifndef DIFFUSION_IMAGING_EXTENDED if(index>=3) { index = index + 1; } #endif switch(index) { case 0: { // Numerical Reconstruct(0,0); break; } case 1: { // Standard Reconstruct(1,0); break; } case 2: { // Solid Angle Reconstruct(1,6); break; } case 3: { // Constrained Solid Angle Reconstruct(1,7); break; } case 4: { // ADC Reconstruct(1,4); break; } case 5: { // Raw Signal Reconstruct(1,5); break; } case 6: { // Q-Ball reconstruction Reconstruct(2,0); break; } } } void QmitkQBallReconstructionView::MethodChoosen(int method) { #ifndef DIFFUSION_IMAGING_EXTENDED if(method>=3) { method = method + 1; } #endif m_Controls->m_QBallSelectionBox->setHidden(true); m_Controls->m_OutputCoeffsImage->setHidden(true); if (method==0) m_Controls->m_ShFrame->setVisible(false); else m_Controls->m_ShFrame->setVisible(true); switch(method) { case 0: m_Controls->m_Description->setText("Numerical recon. (Tuch 2004)"); break; case 1: m_Controls->m_Description->setText("Spherical harmonics recon. (Descoteaux 2007)"); m_Controls->m_OutputCoeffsImage->setHidden(false); break; case 2: m_Controls->m_Description->setText("SH recon. with solid angle consideration (Aganj 2009)"); m_Controls->m_OutputCoeffsImage->setHidden(false); break; case 3: m_Controls->m_Description->setText("SH solid angle with non-neg. constraint (Goh 2009)"); m_Controls->m_OutputCoeffsImage->setHidden(false); break; case 4: m_Controls->m_Description->setText("SH recon. of the plain ADC-profiles"); m_Controls->m_OutputCoeffsImage->setHidden(false); break; case 5: m_Controls->m_Description->setText("SH recon. of the raw diffusion signal"); m_Controls->m_OutputCoeffsImage->setHidden(false); break; case 6: m_Controls->m_Description->setText("SH recon. of the multi shell diffusion signal (Aganj 2010)"); m_Controls->m_QBallSelectionBox->setHidden(false); m_Controls->m_OutputCoeffsImage->setHidden(false); break; } } void QmitkQBallReconstructionView::Reconstruct(int method, int normalization) { if (m_Controls->m_ImageBox->GetSelectedNode().IsNotNull()) { if(method == 0) { NumericalQBallReconstruction(m_Controls->m_ImageBox->GetSelectedNode(), normalization); } else if(method == 1) { AnalyticalQBallReconstruction(m_Controls->m_ImageBox->GetSelectedNode(), normalization); } else if(method == 2) { MultiQBallReconstruction(m_Controls->m_ImageBox->GetSelectedNode()); } } } void QmitkQBallReconstructionView::NumericalQBallReconstruction(mitk::DataNode::Pointer node, int normalization) { try { mitk::Image* vols = static_cast(node->GetData()); std::string nodename = node->GetName(); typedef itk::DiffusionQballReconstructionImageFilter QballReconstructionImageFilterType; ITKDiffusionImageType::Pointer itkVectorImagePointer = ITKDiffusionImageType::New(); mitk::CastToItkImage(vols, itkVectorImagePointer); QballReconstructionImageFilterType::Pointer filter = QballReconstructionImageFilterType::New(); - filter->SetGradientImage(mitk::DiffusionPropertyHelper::GetGradientContainer(vols), itkVectorImagePointer); filter->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(vols)); + filter->SetGradientImage(mitk::DiffusionPropertyHelper::GetGradientContainer(vols), itkVectorImagePointer); filter->SetThreshold( m_Controls->m_QBallReconstructionThreasholdEdit->value() ); std::string nodePostfix; switch(normalization) { case 0: { filter->SetNormalizationMethod(QballReconstructionImageFilterType::QBR_STANDARD); nodePostfix = "_Numerical_Qball"; break; } case 1: { filter->SetNormalizationMethod(QballReconstructionImageFilterType::QBR_B_ZERO_B_VALUE); nodePostfix = "_Numerical_ZeroBvalueNormalization_Qball"; break; } case 2: { filter->SetNormalizationMethod(QballReconstructionImageFilterType::QBR_B_ZERO); nodePostfix = "_NumericalQball_ZeroNormalization_Qball"; break; } case 3: { filter->SetNormalizationMethod(QballReconstructionImageFilterType::QBR_NONE); nodePostfix = "_NumericalQball_NoNormalization_Qball"; break; } default: { filter->SetNormalizationMethod(QballReconstructionImageFilterType::QBR_STANDARD); nodePostfix = "_NumericalQball_Qball"; } } filter->Update(); // ODFs TO DATATREE mitk::OdfImage::Pointer image = mitk::OdfImage::New(); image->InitializeByItk( filter->GetOutput() ); image->SetVolume( filter->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer new_node = mitk::DataNode::New(); new_node->SetData( image ); new_node->SetName(nodename+nodePostfix); mitk::ProgressBar::GetInstance()->Progress(); GetDataStorage()->Add(new_node, node); this->GetRenderWindowPart()->RequestUpdate(); } catch (itk::ExceptionObject &ex) { MITK_INFO << ex ; QMessageBox::information(0, "Reconstruction not possible:", ex.GetDescription()); return ; } } void QmitkQBallReconstructionView::AnalyticalQBallReconstruction( mitk::DataNode::Pointer node, int normalization) { try { float lambda = m_Controls->m_QBallReconstructionLambdaLineEdit->value(); switch(m_Controls->m_QBallReconstructionMaxLLevelComboBox->currentIndex()) { case 0: { TemplatedAnalyticalQBallReconstruction<2>(node, lambda, normalization); break; } case 1: { TemplatedAnalyticalQBallReconstruction<4>(node, lambda, normalization); break; } case 2: { TemplatedAnalyticalQBallReconstruction<6>(node, lambda, normalization); break; } case 3: { TemplatedAnalyticalQBallReconstruction<8>(node, lambda, normalization); break; } case 4: { TemplatedAnalyticalQBallReconstruction<10>(node, lambda, normalization); break; } case 5: { TemplatedAnalyticalQBallReconstruction<12>(node, lambda, normalization); break; } } this->GetRenderWindowPart()->RequestUpdate(); } catch (itk::ExceptionObject &ex) { MITK_INFO << ex; QMessageBox::information(0, "Reconstruction not possible:", ex.GetDescription()); return; } } template void QmitkQBallReconstructionView::TemplatedAnalyticalQBallReconstruction(mitk::DataNode* dataNodePointer, float lambda, int normalization) { typedef itk::AnalyticalDiffusionQballReconstructionImageFilter FilterType; typename FilterType::Pointer filter = FilterType::New(); mitk::Image* vols = dynamic_cast(dataNodePointer->GetData()); ITKDiffusionImageType::Pointer itkVectorImagePointer = ITKDiffusionImageType::New(); mitk::CastToItkImage(vols, itkVectorImagePointer); - filter->SetGradientImage(mitk::DiffusionPropertyHelper::GetGradientContainer(vols), itkVectorImagePointer); filter->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(vols)); + filter->SetGradientImage(mitk::DiffusionPropertyHelper::GetGradientContainer(vols), itkVectorImagePointer); filter->SetThreshold( m_Controls->m_QBallReconstructionThreasholdEdit->value() ); filter->SetLambda(lambda); std::string nodePostfix; switch(normalization) { case 0: { filter->SetNormalizationMethod(FilterType::QBAR_STANDARD); nodePostfix = "_SH_Qball"; break; } case 1: { filter->SetNormalizationMethod(FilterType::QBAR_B_ZERO_B_VALUE); nodePostfix = "_SH_1_Qball"; break; } case 2: { filter->SetNormalizationMethod(FilterType::QBAR_B_ZERO); nodePostfix = "_SH_2_Qball"; break; } case 3: { filter->SetNormalizationMethod(FilterType::QBAR_NONE); nodePostfix = "_SH_3_Qball"; break; } case 4: { filter->SetNormalizationMethod(FilterType::QBAR_ADC_ONLY); nodePostfix = "_AdcProfile"; break; } case 5: { filter->SetNormalizationMethod(FilterType::QBAR_RAW_SIGNAL); nodePostfix = "_RawSignal"; break; } case 6: { filter->SetNormalizationMethod(FilterType::QBAR_SOLID_ANGLE); nodePostfix = "_SH_CSA_Qball"; break; } case 7: { filter->SetNormalizationMethod(FilterType::QBAR_NONNEG_SOLID_ANGLE); nodePostfix = "_SH_NonNegCSA_Qball"; break; } default: { filter->SetNormalizationMethod(FilterType::QBAR_STANDARD); } } filter->Update(); // ODFs TO DATATREE mitk::OdfImage::Pointer image = mitk::OdfImage::New(); image->InitializeByItk( filter->GetOutput() ); image->SetVolume( filter->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( image ); node->SetName(dataNodePointer->GetName()+nodePostfix); GetDataStorage()->Add(node, dataNodePointer); if(m_Controls->m_OutputCoeffsImage->isChecked()) { mitk::Image::Pointer coeffsImage = dynamic_cast(mitk::ShImage::New().GetPointer()); coeffsImage->InitializeByItk( filter->GetCoefficientImage().GetPointer() ); coeffsImage->SetVolume( filter->GetCoefficientImage()->GetBufferPointer() ); mitk::DataNode::Pointer coeffsNode=mitk::DataNode::New(); coeffsNode->SetData( coeffsImage ); coeffsNode->SetProperty( "name", mitk::StringProperty::New(dataNodePointer->GetName()+"_SH-Coeffs") ); GetDataStorage()->Add(coeffsNode, node); } } void QmitkQBallReconstructionView::MultiQBallReconstruction(mitk::DataNode::Pointer node) { try { float lambda = m_Controls->m_QBallReconstructionLambdaLineEdit->value(); switch(m_Controls->m_QBallReconstructionMaxLLevelComboBox->currentIndex()) { case 0: { TemplatedMultiQBallReconstruction<2>(lambda, node); break; } case 1: { TemplatedMultiQBallReconstruction<4>(lambda, node); break; } case 2: { TemplatedMultiQBallReconstruction<6>(lambda, node); break; } case 3: { TemplatedMultiQBallReconstruction<8>(lambda, node); break; } case 4: { TemplatedMultiQBallReconstruction<10>(lambda, node); break; } case 5: { TemplatedMultiQBallReconstruction<12>(lambda, node); break; } } } catch (itk::ExceptionObject &ex) { MITK_INFO << ex ; QMessageBox::information(0, "Reconstruction not possible:", ex.GetDescription()); return ; } } template void QmitkQBallReconstructionView::TemplatedMultiQBallReconstruction(float lambda, mitk::DataNode* dataNodePointer) { typedef itk::DiffusionMultiShellQballReconstructionImageFilter FilterType; typename FilterType::Pointer filter = FilterType::New(); std::string nodename; dataNodePointer->GetStringProperty("name",nodename); mitk::Image* dwi = dynamic_cast(dataNodePointer->GetData()); BValueMapType currSelectionMap = m_ShellSelectorMap[dataNodePointer]->GetBValueSelctionMap(); if(currSelectionMap.size() != 4)// || currSelectionMap.find(0) != currSelectionMap.end()) { QMessageBox::information(0, "Reconstruction not possible:" ,QString("Only three equidistant shells are supported. (ImageName: " + QString(nodename.c_str()) + ")")); return; } BValueMapType::reverse_iterator it1 = currSelectionMap.rbegin(); BValueMapType::reverse_iterator it2 = currSelectionMap.rbegin(); ++it2; // Get average distance int avdistance = 0; for(; it2 != currSelectionMap.rend(); ++it1,++it2) avdistance += (int)it1->first - (int)it2->first; avdistance /= currSelectionMap.size()-1; // Check if all shells are using the same averae distance it1 = currSelectionMap.rbegin(); it2 = currSelectionMap.rbegin(); ++it2; for(; it2 != currSelectionMap.rend(); ++it1,++it2) { if(avdistance != (int)it1->first - (int)it2->first) { QMessageBox::information(0, "Reconstruction not possible:" ,QString("Selected Shells are not in a equidistant configuration. (ImageName: " + QString(nodename.c_str()) + ")")); return; } } ITKDiffusionImageType::Pointer itkVectorImagePointer = ITKDiffusionImageType::New(); mitk::CastToItkImage(dwi, itkVectorImagePointer); filter->SetBValueMap(m_ShellSelectorMap[dataNodePointer]->GetBValueSelctionMap()); filter->SetGradientImage(mitk::DiffusionPropertyHelper::GetGradientContainer(dwi), itkVectorImagePointer, mitk::DiffusionPropertyHelper::GetReferenceBValue(dwi)); filter->SetThreshold( m_Controls->m_QBallReconstructionThreasholdEdit->value() ); filter->SetLambda(lambda); filter->Update(); if(m_Controls->m_OutputCoeffsImage->isChecked()) { mitk::Image::Pointer coeffsImage = mitk::Image::New(); coeffsImage->InitializeByItk( filter->GetCoefficientImage().GetPointer() ); coeffsImage->SetVolume( filter->GetCoefficientImage()->GetBufferPointer() ); mitk::DataNode::Pointer coeffsNode=mitk::DataNode::New(); coeffsNode->SetData( coeffsImage ); coeffsNode->SetProperty( "name", mitk::StringProperty::New( QString(nodename.c_str()).append("_SH-Coefficients").toStdString()) ); GetDataStorage()->Add(coeffsNode, dataNodePointer); } // ODFs TO DATATREE mitk::OdfImage::Pointer image = mitk::OdfImage::New(); image->InitializeByItk( filter->GetOutput() ); image->SetVolume( filter->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( image ); node->SetName(nodename+"_SH_MultiShell_Qball"); GetDataStorage()->Add(node, dataNodePointer); } void QmitkQBallReconstructionView::GenerateShellSelectionUI(mitk::DataNode::Pointer node) { std::map tempMap; if(m_ShellSelectorMap.find( node.GetPointer() ) != m_ShellSelectorMap.end()) { tempMap[node.GetPointer()] = m_ShellSelectorMap[node.GetPointer()]; m_ShellSelectorMap.erase(node.GetPointer()); } else { tempMap[node.GetPointer()] = new QbrShellSelection(this, node ); tempMap[node.GetPointer()]->SetVisible(true); } for(std::map::iterator it = m_ShellSelectorMap.begin(); it != m_ShellSelectorMap.end();it ++) { delete it->second; } m_ShellSelectorMap.clear(); m_ShellSelectorMap = tempMap; } void QmitkQBallReconstructionView::PreviewThreshold(int threshold) { if (m_Controls->m_ImageBox->GetSelectedNode().IsNotNull()) { mitk::Image* vols = static_cast(m_Controls->m_ImageBox->GetSelectedNode()->GetData()); // Extract b0 image ITKDiffusionImageType::Pointer itkVectorImagePointer = ITKDiffusionImageType::New(); mitk::CastToItkImage(vols, itkVectorImagePointer); typedef itk::B0ImageExtractionImageFilter FilterType; FilterType::Pointer filterB0 = FilterType::New(); filterB0->SetInput( itkVectorImagePointer ); filterB0->SetDirections(mitk::DiffusionPropertyHelper::GetGradientContainer(vols)); filterB0->Update(); mitk::Image::Pointer mitkImage = mitk::Image::New(); typedef itk::Image ImageType; typedef itk::Image SegmentationType; typedef itk::BinaryThresholdImageFilter ThresholdFilterType; // apply threshold ThresholdFilterType::Pointer filterThreshold = ThresholdFilterType::New(); filterThreshold->SetInput(filterB0->GetOutput()); filterThreshold->SetLowerThreshold(threshold); filterThreshold->SetInsideValue(0); filterThreshold->SetOutsideValue(1); // mark cut off values red filterThreshold->Update(); mitkImage->InitializeByItk( filterThreshold->GetOutput() ); mitkImage->SetVolume( filterThreshold->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer node; if (this->GetDataStorage()->GetNamedDerivedNode("ThresholdOverlay", m_Controls->m_ImageBox->GetSelectedNode())) { node = this->GetDataStorage()->GetNamedDerivedNode("ThresholdOverlay", m_Controls->m_ImageBox->GetSelectedNode()); } else { // create a new node, to show thresholded values node = mitk::DataNode::New(); GetDataStorage()->Add( node, m_Controls->m_ImageBox->GetSelectedNode() ); node->SetProperty( "name", mitk::StringProperty::New("ThresholdOverlay")); node->SetBoolProperty("helper object", true); } node->SetData( mitkImage ); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.cpp index bb1f55848d..624c7d4342 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.cpp @@ -1,1028 +1,1042 @@ /*=================================================================== 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 #include // Qmitk #include "QmitkStreamlineTrackingView.h" #include "QmitkStdMultiWidget.h" // Qt #include // MITK #include #include #include #include #include #include #include #include #include #include #include #include +#include #include // VTK #include #include #include #include #include #include #include #include #include #include const std::string QmitkStreamlineTrackingView::VIEW_ID = "org.mitk.views.streamlinetracking"; const std::string id_DataManager = "org.mitk.views.datamanager"; using namespace berry; QmitkStreamlineTrackingWorker::QmitkStreamlineTrackingWorker(QmitkStreamlineTrackingView* view) : m_View(view) { } void QmitkStreamlineTrackingWorker::run() { m_View->m_Tracker->Update(); m_View->m_TrackingThread.quit(); } QmitkStreamlineTrackingView::QmitkStreamlineTrackingView() : m_TrackingWorker(this) , m_Controls(nullptr) , m_FirstTensorProbRun(true) , m_FirstInteractiveRun(true) , m_TrackingHandler(nullptr) , m_ThreadIsRunning(false) , m_DeleteTrackingHandler(false) , m_Visible(false) , m_LastPrior("") , m_TrackingPriorHandler(nullptr) { m_TrackingWorker.moveToThread(&m_TrackingThread); connect(&m_TrackingThread, SIGNAL(started()), this, SLOT(BeforeThread())); connect(&m_TrackingThread, SIGNAL(started()), &m_TrackingWorker, SLOT(run())); connect(&m_TrackingThread, SIGNAL(finished()), this, SLOT(AfterThread())); m_TrackingTimer = new QTimer(this); } // Destructor QmitkStreamlineTrackingView::~QmitkStreamlineTrackingView() { if (m_Tracker.IsNull()) return; m_Tracker->SetStopTracking(true); m_TrackingThread.wait(); } void QmitkStreamlineTrackingView::CreateQtPartControl( QWidget *parent ) { if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkStreamlineTrackingViewControls; m_Controls->setupUi( parent ); - m_Controls->m_FaImageBox->SetDataStorage(this->GetDataStorage()); - m_Controls->m_SeedImageBox->SetDataStorage(this->GetDataStorage()); - m_Controls->m_MaskImageBox->SetDataStorage(this->GetDataStorage()); - m_Controls->m_TargetImageBox->SetDataStorage(this->GetDataStorage()); - m_Controls->m_PriorImageBox->SetDataStorage(this->GetDataStorage()); - m_Controls->m_StopImageBox->SetDataStorage(this->GetDataStorage()); - m_Controls->m_ForestBox->SetDataStorage(this->GetDataStorage()); - m_Controls->m_ExclusionImageBox->SetDataStorage(this->GetDataStorage()); + m_Controls->m_FaImageSelectionWidget->SetDataStorage(this->GetDataStorage()); + m_Controls->m_SeedImageSelectionWidget->SetDataStorage(this->GetDataStorage()); + m_Controls->m_MaskImageSelectionWidget->SetDataStorage(this->GetDataStorage()); + m_Controls->m_TargetImageSelectionWidget->SetDataStorage(this->GetDataStorage()); + m_Controls->m_PriorImageSelectionWidget->SetDataStorage(this->GetDataStorage()); + m_Controls->m_StopImageSelectionWidget->SetDataStorage(this->GetDataStorage()); + m_Controls->m_ForestSelectionWidget->SetDataStorage(this->GetDataStorage()); + m_Controls->m_ExclusionImageSelectionWidget->SetDataStorage(this->GetDataStorage()); mitk::TNodePredicateDataType::Pointer isPeakImagePredicate = mitk::TNodePredicateDataType::New(); mitk::TNodePredicateDataType::Pointer isImagePredicate = mitk::TNodePredicateDataType::New(); mitk::TNodePredicateDataType::Pointer isTractographyForest = mitk::TNodePredicateDataType::New(); mitk::NodePredicateProperty::Pointer isBinaryPredicate = mitk::NodePredicateProperty::New("binary", mitk::BoolProperty::New(true)); mitk::NodePredicateNot::Pointer isNotBinaryPredicate = mitk::NodePredicateNot::New( isBinaryPredicate ); mitk::NodePredicateAnd::Pointer isNotABinaryImagePredicate = mitk::NodePredicateAnd::New( isImagePredicate, isNotBinaryPredicate ); mitk::NodePredicateDimension::Pointer dimensionPredicate = mitk::NodePredicateDimension::New(3); - m_Controls->m_ForestBox->SetPredicate(isTractographyForest); - m_Controls->m_FaImageBox->SetPredicate( mitk::NodePredicateAnd::New(isNotABinaryImagePredicate, dimensionPredicate) ); - m_Controls->m_FaImageBox->SetZeroEntryText("--"); - m_Controls->m_SeedImageBox->SetPredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); - m_Controls->m_SeedImageBox->SetZeroEntryText("--"); - m_Controls->m_MaskImageBox->SetPredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); - m_Controls->m_MaskImageBox->SetZeroEntryText("--"); - m_Controls->m_StopImageBox->SetPredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); - m_Controls->m_StopImageBox->SetZeroEntryText("--"); - m_Controls->m_TargetImageBox->SetPredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); - m_Controls->m_TargetImageBox->SetZeroEntryText("--"); - m_Controls->m_PriorImageBox->SetPredicate( isPeakImagePredicate ); - m_Controls->m_PriorImageBox->SetZeroEntryText("--"); - m_Controls->m_ExclusionImageBox->SetPredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); - m_Controls->m_ExclusionImageBox->SetZeroEntryText("--"); + m_Controls->m_ForestSelectionWidget->SetNodePredicate(isTractographyForest); + m_Controls->m_FaImageSelectionWidget->SetNodePredicate( mitk::NodePredicateAnd::New(isNotABinaryImagePredicate, dimensionPredicate) ); + m_Controls->m_FaImageSelectionWidget->SetEmptyInfo("--"); + m_Controls->m_FaImageSelectionWidget->SetSelectionIsOptional(true); + + m_Controls->m_SeedImageSelectionWidget->SetNodePredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); + m_Controls->m_SeedImageSelectionWidget->SetEmptyInfo("--"); + m_Controls->m_SeedImageSelectionWidget->SetSelectionIsOptional(true); + + m_Controls->m_MaskImageSelectionWidget->SetNodePredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); + m_Controls->m_MaskImageSelectionWidget->SetEmptyInfo("--"); + m_Controls->m_MaskImageSelectionWidget->SetSelectionIsOptional(true); + + m_Controls->m_StopImageSelectionWidget->SetNodePredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); + m_Controls->m_StopImageSelectionWidget->SetEmptyInfo("--"); + m_Controls->m_StopImageSelectionWidget->SetSelectionIsOptional(true); + + m_Controls->m_TargetImageSelectionWidget->SetNodePredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); + m_Controls->m_TargetImageSelectionWidget->SetEmptyInfo("--"); + m_Controls->m_TargetImageSelectionWidget->SetSelectionIsOptional(true); + + m_Controls->m_PriorImageSelectionWidget->SetNodePredicate( isPeakImagePredicate ); + m_Controls->m_PriorImageSelectionWidget->SetEmptyInfo("--"); + m_Controls->m_PriorImageSelectionWidget->SetSelectionIsOptional(true); + + m_Controls->m_ExclusionImageSelectionWidget->SetNodePredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); + m_Controls->m_ExclusionImageSelectionWidget->SetEmptyInfo("--"); + m_Controls->m_ExclusionImageSelectionWidget->SetSelectionIsOptional(true); connect( m_TrackingTimer, SIGNAL(timeout()), this, SLOT(TimerUpdate()) ); connect( m_Controls->commandLinkButton_2, SIGNAL(clicked()), this, SLOT(StopTractography()) ); connect( m_Controls->commandLinkButton, SIGNAL(clicked()), this, SLOT(DoFiberTracking()) ); connect( m_Controls->m_InteractiveBox, SIGNAL(stateChanged(int)), this, SLOT(ToggleInteractive()) ); connect( m_Controls->m_ModeBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui()) ); - connect( m_Controls->m_FaImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(DeleteTrackingHandler()) ); + connect( m_Controls->m_FaImageSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(DeleteTrackingHandler()) ); connect( m_Controls->m_ModeBox, SIGNAL(currentIndexChanged(int)), this, SLOT(DeleteTrackingHandler()) ); connect( m_Controls->m_OutputProbMap, SIGNAL(stateChanged(int)), this, SLOT(OutputStyleSwitched()) ); - connect( m_Controls->m_SeedImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); + connect( m_Controls->m_SeedImageSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_ModeBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); - connect( m_Controls->m_StopImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); - connect( m_Controls->m_TargetImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); - connect( m_Controls->m_PriorImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); - connect( m_Controls->m_ExclusionImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); - connect( m_Controls->m_MaskImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); - connect( m_Controls->m_FaImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); - connect( m_Controls->m_ForestBox, SIGNAL(currentIndexChanged(int)), this, SLOT(ForestSwitched()) ); - connect( m_Controls->m_ForestBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); + connect( m_Controls->m_StopImageSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(OnParameterChanged()) ); + connect( m_Controls->m_TargetImageSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(OnParameterChanged()) ); + connect( m_Controls->m_PriorImageSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(OnParameterChanged()) ); + connect( m_Controls->m_ExclusionImageSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(OnParameterChanged()) ); + connect( m_Controls->m_MaskImageSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(OnParameterChanged()) ); + connect( m_Controls->m_FaImageSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(OnParameterChanged()) ); + connect( m_Controls->m_ForestSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(ForestSwitched()) ); + connect( m_Controls->m_ForestSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SeedsPerVoxelBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_NumFibersBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_ScalarThresholdBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_OdfCutoffBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_StepSizeBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SamplingDistanceBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_AngularThresholdBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_MinTractLengthBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_fBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_gBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_NumSamplesBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SeedRadiusBox, SIGNAL(editingFinished()), this, SLOT(InteractiveSeedChanged()) ); connect( m_Controls->m_NumSeedsBox, SIGNAL(editingFinished()), this, SLOT(InteractiveSeedChanged()) ); connect( m_Controls->m_OutputProbMap, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SharpenOdfsBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_InterpolationBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_MaskInterpolationBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FlipXBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FlipYBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FlipZBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FrontalSamplesBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_StopVotesBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_LoopCheckBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_TrialsPerSeedBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_EpConstraintsBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); m_Controls->m_SeedsPerVoxelBox->editingFinished(); m_Controls->m_NumFibersBox->editingFinished(); m_Controls->m_ScalarThresholdBox->editingFinished(); m_Controls->m_OdfCutoffBox->editingFinished(); m_Controls->m_StepSizeBox->editingFinished(); m_Controls->m_SamplingDistanceBox->editingFinished(); m_Controls->m_AngularThresholdBox->editingFinished(); m_Controls->m_MinTractLengthBox->editingFinished(); m_Controls->m_fBox->editingFinished(); m_Controls->m_gBox->editingFinished(); m_Controls->m_NumSamplesBox->editingFinished(); m_Controls->m_SeedRadiusBox->editingFinished(); m_Controls->m_NumSeedsBox->editingFinished(); m_Controls->m_LoopCheckBox->editingFinished(); m_Controls->m_TrialsPerSeedBox->editingFinished(); StartStopTrackingGui(false); } UpdateGui(); } void QmitkStreamlineTrackingView::StopTractography() { if (m_Tracker.IsNull()) return; m_Tracker->SetStopTracking(true); } void QmitkStreamlineTrackingView::TimerUpdate() { if (m_Tracker.IsNull()) return; QString status_text(m_Tracker->GetStatusText().c_str()); m_Controls->m_StatusTextBox->setText(status_text); } void QmitkStreamlineTrackingView::BeforeThread() { m_TrackingTimer->start(1000); } void QmitkStreamlineTrackingView::AfterThread() { m_TrackingTimer->stop(); if (!m_Tracker->GetUseOutputProbabilityMap()) { vtkSmartPointer fiberBundle = m_Tracker->GetFiberPolyData(); if (!m_Controls->m_InteractiveBox->isChecked() && fiberBundle->GetNumberOfLines() == 0) { QMessageBox warnBox; warnBox.setWindowTitle("Warning"); warnBox.setText("No fiberbundle was generated!"); warnBox.setDetailedText("No fibers were generated using the chosen parameters. Typical reasons are:\n\n- Cutoff too high. Some images feature very low FA/GFA/peak size. Try to lower this parameter.\n- Angular threshold too strict. Try to increase this parameter.\n- A small step sizes also means many steps to go wrong. Especially in the case of probabilistic tractography. Try to adjust the angular threshold."); warnBox.setIcon(QMessageBox::Warning); warnBox.exec(); if (m_InteractivePointSetNode.IsNotNull()) m_InteractivePointSetNode->SetProperty("color", mitk::ColorProperty::New(1,1,1)); StartStopTrackingGui(false); if (m_DeleteTrackingHandler) DeleteTrackingHandler(); UpdateGui(); return; } mitk::FiberBundle::Pointer fib = mitk::FiberBundle::New(fiberBundle); fib->SetReferenceGeometry(dynamic_cast(m_ParentNode->GetData())->GetGeometry()); if (m_Controls->m_ResampleFibersBox->isChecked() && fiberBundle->GetNumberOfLines()>0) fib->Compress(m_Controls->m_FiberErrorBox->value()); fib->ColorFibersByOrientation(); m_Tracker->SetDicomProperties(fib); if (m_Controls->m_InteractiveBox->isChecked()) { if (m_InteractiveNode.IsNull()) { m_InteractiveNode = mitk::DataNode::New(); QString name("Interactive"); m_InteractiveNode->SetName(name.toStdString()); GetDataStorage()->Add(m_InteractiveNode); } m_InteractiveNode->SetData(fib); m_InteractiveNode->SetFloatProperty("Fiber2DSliceThickness", m_Tracker->GetMinVoxelSize()/2); if (auto renderWindowPart = this->GetRenderWindowPart()) renderWindowPart->RequestUpdate(); } else { mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(fib); QString name("FiberBundle_"); name += m_ParentNode->GetName().c_str(); name += "_Streamline"; node->SetName(name.toStdString()); node->SetFloatProperty("Fiber2DSliceThickness", m_Tracker->GetMinVoxelSize()/2); GetDataStorage()->Add(node, m_ParentNode); } } else { TrackerType::ItkDoubleImgType::Pointer outImg = m_Tracker->GetOutputProbabilityMap(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); if (m_Controls->m_InteractiveBox->isChecked()) { if (m_InteractiveNode.IsNull()) { m_InteractiveNode = mitk::DataNode::New(); QString name("Interactive"); m_InteractiveNode->SetName(name.toStdString()); GetDataStorage()->Add(m_InteractiveNode); } m_InteractiveNode->SetData(img); mitk::LookupTable::Pointer lut = mitk::LookupTable::New(); lut->SetType(mitk::LookupTable::JET_TRANSPARENT); mitk::LookupTableProperty::Pointer lut_prop = mitk::LookupTableProperty::New(); lut_prop->SetLookupTable(lut); m_InteractiveNode->SetProperty("LookupTable", lut_prop); m_InteractiveNode->SetProperty("opacity", mitk::FloatProperty::New(0.5)); m_InteractiveNode->SetFloatProperty("Fiber2DSliceThickness", m_Tracker->GetMinVoxelSize()/2); if (auto renderWindowPart = this->GetRenderWindowPart()) renderWindowPart->RequestUpdate(); } else { mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(img); QString name("ProbabilityMap_"); name += m_ParentNode->GetName().c_str(); node->SetName(name.toStdString()); mitk::LookupTable::Pointer lut = mitk::LookupTable::New(); lut->SetType(mitk::LookupTable::JET_TRANSPARENT); mitk::LookupTableProperty::Pointer lut_prop = mitk::LookupTableProperty::New(); lut_prop->SetLookupTable(lut); node->SetProperty("LookupTable", lut_prop); node->SetProperty("opacity", mitk::FloatProperty::New(0.5)); GetDataStorage()->Add(node, m_ParentNode); } } if (m_InteractivePointSetNode.IsNotNull()) m_InteractivePointSetNode->SetProperty("color", mitk::ColorProperty::New(1,1,1)); StartStopTrackingGui(false); if (m_DeleteTrackingHandler) DeleteTrackingHandler(); UpdateGui(); } void QmitkStreamlineTrackingView::InteractiveSeedChanged(bool posChanged) { if(!CheckAndStoreLastParams(sender()) && !posChanged) return; if (m_ThreadIsRunning || !m_Visible) return; if (!posChanged && (!m_Controls->m_InteractiveBox->isChecked() || !m_Controls->m_ParamUpdateBox->isChecked()) ) return; std::srand(std::time(0)); m_SeedPoints.clear(); itk::Point world_pos = this->GetRenderWindowPart()->GetSelectedPosition(); m_SeedPoints.push_back(world_pos); float radius = m_Controls->m_SeedRadiusBox->value(); int num = m_Controls->m_NumSeedsBox->value(); mitk::PointSet::Pointer pointset = mitk::PointSet::New(); pointset->InsertPoint(0, world_pos); m_InteractivePointSetNode->SetProperty("pointsize", mitk::FloatProperty::New(radius*2)); m_InteractivePointSetNode->SetProperty("point 2D size", mitk::FloatProperty::New(radius*2)); m_InteractivePointSetNode->SetData(pointset); for (int i=1; i p; p[0] = rand()%1000-500; p[1] = rand()%1000-500; p[2] = rand()%1000-500; p.Normalize(); p *= radius; m_SeedPoints.push_back(world_pos+p); } m_InteractivePointSetNode->SetProperty("color", mitk::ColorProperty::New(1,0,0)); DoFiberTracking(); } bool QmitkStreamlineTrackingView::CheckAndStoreLastParams(QObject* obj) { if (obj!=nullptr) { std::string new_val = ""; if(qobject_cast(obj)!=nullptr) new_val = boost::lexical_cast(qobject_cast(obj)->value()); else if (qobject_cast(obj)!=nullptr) new_val = boost::lexical_cast(qobject_cast(obj)->value()); else return true; if (m_LastTractoParams.find(obj->objectName())==m_LastTractoParams.end()) { m_LastTractoParams[obj->objectName()] = new_val; return false; } else if (m_LastTractoParams.at(obj->objectName()) != new_val) { m_LastTractoParams[obj->objectName()] = new_val; return true; } else if (m_LastTractoParams.at(obj->objectName()) == new_val) return false; } return true; } void QmitkStreamlineTrackingView::OnParameterChanged() { UpdateGui(); if(!CheckAndStoreLastParams(sender())) return; if (m_Controls->m_InteractiveBox->isChecked() && m_Controls->m_ParamUpdateBox->isChecked()) DoFiberTracking(); } void QmitkStreamlineTrackingView::ToggleInteractive() { UpdateGui(); m_Controls->m_SeedsPerVoxelBox->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->m_SeedsPerVoxelLabel->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); - m_Controls->m_SeedImageBox->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); + m_Controls->m_SeedImageSelectionWidget->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->label_6->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); if ( m_Controls->m_InteractiveBox->isChecked() ) { if (m_FirstInteractiveRun) { QMessageBox::information(nullptr, "Information", "Place and move a spherical seed region anywhere in the image by left-clicking and dragging. If the seed region is colored red, tracking is in progress. If the seed region is colored white, tracking is finished.\nPlacing the seed region for the first time in a newly selected dataset might cause a short delay, since the tracker needs to be initialized."); m_FirstInteractiveRun = false; } QApplication::setOverrideCursor(Qt::PointingHandCursor); QApplication::processEvents(); m_InteractivePointSetNode = mitk::DataNode::New(); m_InteractivePointSetNode->SetProperty("color", mitk::ColorProperty::New(1,1,1)); m_InteractivePointSetNode->SetName("InteractiveSeedRegion"); mitk::PointSetShapeProperty::Pointer shape_prop = mitk::PointSetShapeProperty::New(); shape_prop->SetValue(mitk::PointSetShapeProperty::PointSetShape::CIRCLE); m_InteractivePointSetNode->SetProperty("Pointset.2D.shape", shape_prop); GetDataStorage()->Add(m_InteractivePointSetNode); m_SliceChangeListener.RenderWindowPartActivated(this->GetRenderWindowPart()); connect(&m_SliceChangeListener, SIGNAL(SliceChanged()), this, SLOT(OnSliceChanged())); } else { QApplication::restoreOverrideCursor(); QApplication::processEvents(); m_InteractiveNode = nullptr; m_InteractivePointSetNode = nullptr; m_SliceChangeListener.RenderWindowPartActivated(this->GetRenderWindowPart()); disconnect(&m_SliceChangeListener, SIGNAL(SliceChanged()), this, SLOT(OnSliceChanged())); } } void QmitkStreamlineTrackingView::Activated() { } void QmitkStreamlineTrackingView::Deactivated() { } void QmitkStreamlineTrackingView::Visible() { m_Visible = true; } void QmitkStreamlineTrackingView::Hidden() { m_Visible = false; m_Controls->m_InteractiveBox->setChecked(false); ToggleInteractive(); } void QmitkStreamlineTrackingView::OnSliceChanged() { InteractiveSeedChanged(true); } void QmitkStreamlineTrackingView::SetFocus() { } void QmitkStreamlineTrackingView::DeleteTrackingHandler() { if (!m_ThreadIsRunning && m_TrackingHandler != nullptr) { delete m_TrackingHandler; m_TrackingHandler = nullptr; m_DeleteTrackingHandler = false; m_LastPrior = ""; if (m_TrackingPriorHandler != nullptr) delete m_TrackingPriorHandler; } else if (m_ThreadIsRunning) { m_DeleteTrackingHandler = true; } } void QmitkStreamlineTrackingView::ForestSwitched() { DeleteTrackingHandler(); } void QmitkStreamlineTrackingView::OutputStyleSwitched() { if (m_InteractiveNode.IsNotNull()) GetDataStorage()->Remove(m_InteractiveNode); m_InteractiveNode = nullptr; } void QmitkStreamlineTrackingView::OnSelectionChanged( berry::IWorkbenchPart::Pointer , const QList& nodes ) { std::vector< mitk::DataNode::Pointer > last_nodes = m_InputImageNodes; m_InputImageNodes.clear(); m_InputImages.clear(); m_AdditionalInputImages.clear(); bool retrack = false; for( auto node : nodes ) { if( node.IsNotNull() && dynamic_cast(node->GetData()) ) { if( dynamic_cast(node->GetData()) ) { m_InputImageNodes.push_back(node); m_InputImages.push_back(dynamic_cast(node->GetData())); retrack = true; } else if ( dynamic_cast(node->GetData()) ) { m_InputImageNodes.push_back(node); m_InputImages.push_back(dynamic_cast(node->GetData())); retrack = true; } else if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(node->GetData())) ) { m_InputImageNodes.push_back(node); m_InputImages.push_back(dynamic_cast(node->GetData())); retrack = true; } else { mitk::Image* img = dynamic_cast(node->GetData()); if (img!=nullptr) { int dim = img->GetDimension(); unsigned int* dimensions = img->GetDimensions(); if (dim==4 && dimensions[3]%3==0) { m_InputImageNodes.push_back(node); m_InputImages.push_back(dynamic_cast(node->GetData())); retrack = true; } else if (dim==3) { m_AdditionalInputImages.push_back(dynamic_cast(node->GetData())); } } } } } // sometimes the OnSelectionChanged event is sent twice and actually no selection has changed for the first event. We need to catch that. if (last_nodes.size() == m_InputImageNodes.size()) { bool same_nodes = true; for (unsigned int i=0; im_TensorImageLabel->setText("select in data-manager"); m_Controls->m_fBox->setEnabled(false); m_Controls->m_fLabel->setEnabled(false); m_Controls->m_gBox->setEnabled(false); m_Controls->m_gLabel->setEnabled(false); - m_Controls->m_FaImageBox->setEnabled(true); + m_Controls->m_FaImageSelectionWidget->setEnabled(true); m_Controls->mFaImageLabel->setEnabled(true); m_Controls->m_OdfCutoffBox->setEnabled(false); m_Controls->m_OdfCutoffLabel->setEnabled(false); m_Controls->m_SharpenOdfsBox->setEnabled(false); - m_Controls->m_ForestBox->setVisible(false); + m_Controls->m_ForestSelectionWidget->setVisible(false); m_Controls->m_ForestLabel->setVisible(false); m_Controls->commandLinkButton->setEnabled(false); m_Controls->m_TrialsPerSeedBox->setEnabled(false); m_Controls->m_TrialsPerSeedLabel->setEnabled(false); - m_Controls->m_TargetImageBox->setEnabled(false); + m_Controls->m_TargetImageSelectionWidget->setEnabled(false); m_Controls->m_TargetImageLabel->setEnabled(false); if (m_Controls->m_InteractiveBox->isChecked()) { m_Controls->m_InteractiveSeedingFrame->setVisible(true); m_Controls->m_StaticSeedingFrame->setVisible(false); m_Controls->commandLinkButton_2->setVisible(false); m_Controls->commandLinkButton->setVisible(false); } else { m_Controls->m_InteractiveSeedingFrame->setVisible(false); m_Controls->m_StaticSeedingFrame->setVisible(true); m_Controls->commandLinkButton_2->setVisible(m_ThreadIsRunning); m_Controls->commandLinkButton->setVisible(!m_ThreadIsRunning); } if (m_Controls->m_EpConstraintsBox->currentIndex()>0) { - m_Controls->m_TargetImageBox->setEnabled(true); + m_Controls->m_TargetImageSelectionWidget->setEnabled(true); m_Controls->m_TargetImageLabel->setEnabled(true); } // trials per seed are only important for probabilistic tractography if (m_Controls->m_ModeBox->currentIndex()==1) { m_Controls->m_TrialsPerSeedBox->setEnabled(true); m_Controls->m_TrialsPerSeedLabel->setEnabled(true); } if(!m_InputImageNodes.empty()) { if (m_InputImageNodes.size()>1) m_Controls->m_TensorImageLabel->setText( ( std::to_string(m_InputImageNodes.size()) + " images selected").c_str() ); else m_Controls->m_TensorImageLabel->setText(m_InputImageNodes.at(0)->GetName().c_str()); m_Controls->commandLinkButton->setEnabled(!m_Controls->m_InteractiveBox->isChecked() && !m_ThreadIsRunning); m_Controls->m_ScalarThresholdBox->setEnabled(true); m_Controls->m_FaThresholdLabel->setEnabled(true); if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { m_Controls->m_fBox->setEnabled(true); m_Controls->m_fLabel->setEnabled(true); m_Controls->m_gBox->setEnabled(true); m_Controls->m_gLabel->setEnabled(true); } else if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { m_Controls->m_OdfCutoffBox->setEnabled(true); m_Controls->m_OdfCutoffLabel->setEnabled(true); m_Controls->m_SharpenOdfsBox->setEnabled(true); } else if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(m_InputImageNodes.at(0)->GetData())) ) { - m_Controls->m_ForestBox->setVisible(true); + m_Controls->m_ForestSelectionWidget->setVisible(true); m_Controls->m_ForestLabel->setVisible(true); m_Controls->m_ScalarThresholdBox->setEnabled(false); m_Controls->m_FaThresholdLabel->setEnabled(false); } } } void QmitkStreamlineTrackingView::StartStopTrackingGui(bool start) { m_ThreadIsRunning = start; if (!m_Controls->m_InteractiveBox->isChecked()) { m_Controls->commandLinkButton_2->setVisible(start); m_Controls->commandLinkButton->setVisible(!start); m_Controls->m_InteractiveBox->setEnabled(!start); m_Controls->m_StatusTextBox->setVisible(start); } } void QmitkStreamlineTrackingView::DoFiberTracking() { if (m_InputImages.empty()) { QMessageBox::information(nullptr, "Information", "Please select an input image in the datamaneger (tensor, ODF, peak or dMRI image)!"); return; } if (m_ThreadIsRunning || !m_Visible) return; if (m_Controls->m_InteractiveBox->isChecked() && m_SeedPoints.empty()) return; StartStopTrackingGui(true); m_Tracker = TrackerType::New(); if( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { if (m_Controls->m_ModeBox->currentIndex()==1) { if (m_InputImages.size()>1) { QMessageBox::information(nullptr, "Information", "Probabilistic tensor tractography is only implemented for single-tensor mode!"); StartStopTrackingGui(false); return; } if (m_TrackingHandler==nullptr) { m_TrackingHandler = new mitk::TrackingHandlerOdf(); mitk::TensorImage::ItkTensorImageType::Pointer itkImg = mitk::TensorImage::ItkTensorImageType::New(); mitk::CastToItkImage(m_InputImages.at(0), itkImg); typedef itk::TensorImageToOdfImageFilter< float, float > FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetInput( itkImg ); filter->Update(); dynamic_cast(m_TrackingHandler)->SetOdfImage(filter->GetOutput()); - if (m_Controls->m_FaImageBox->GetSelectedNode().IsNotNull()) + if (m_Controls->m_FaImageSelectionWidget->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer itkImg = ItkFloatImageType::New(); - mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageBox->GetSelectedNode()->GetData()), itkImg); + mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageSelectionWidget->GetSelectedNode()->GetData()), itkImg); dynamic_cast(m_TrackingHandler)->SetGfaImage(itkImg); } } dynamic_cast(m_TrackingHandler)->SetGfaThreshold(m_Controls->m_ScalarThresholdBox->value()); dynamic_cast(m_TrackingHandler)->SetOdfThreshold(0); dynamic_cast(m_TrackingHandler)->SetSharpenOdfs(true); dynamic_cast(m_TrackingHandler)->SetIsOdfFromTensor(true); } else { if (m_TrackingHandler==nullptr) { m_TrackingHandler = new mitk::TrackingHandlerTensor(); for (int i=0; i<(int)m_InputImages.size(); i++) { typedef mitk::ImageToItk< mitk::TrackingHandlerTensor::ItkTensorImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(m_InputImages.at(i)); caster->Update(); mitk::TrackingHandlerTensor::ItkTensorImageType::ConstPointer itkImg = caster->GetOutput(); dynamic_cast(m_TrackingHandler)->AddTensorImage(itkImg); } - if (m_Controls->m_FaImageBox->GetSelectedNode().IsNotNull()) + if (m_Controls->m_FaImageSelectionWidget->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer itkImg = ItkFloatImageType::New(); - mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageBox->GetSelectedNode()->GetData()), itkImg); + mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageSelectionWidget->GetSelectedNode()->GetData()), itkImg); dynamic_cast(m_TrackingHandler)->SetFaImage(itkImg); } } dynamic_cast(m_TrackingHandler)->SetFaThreshold(m_Controls->m_ScalarThresholdBox->value()); dynamic_cast(m_TrackingHandler)->SetF((float)m_Controls->m_fBox->value()); dynamic_cast(m_TrackingHandler)->SetG((float)m_Controls->m_gBox->value()); } } else if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { if (m_TrackingHandler==nullptr) { m_TrackingHandler = new mitk::TrackingHandlerOdf(); mitk::TrackingHandlerOdf::ItkOdfImageType::Pointer itkImg = mitk::TrackingHandlerOdf::ItkOdfImageType::New(); mitk::CastToItkImage(m_InputImages.at(0), itkImg); dynamic_cast(m_TrackingHandler)->SetOdfImage(itkImg); - if (m_Controls->m_FaImageBox->GetSelectedNode().IsNotNull()) + if (m_Controls->m_FaImageSelectionWidget->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer itkImg = ItkFloatImageType::New(); - mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageBox->GetSelectedNode()->GetData()), itkImg); + mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageSelectionWidget->GetSelectedNode()->GetData()), itkImg); dynamic_cast(m_TrackingHandler)->SetGfaImage(itkImg); } } dynamic_cast(m_TrackingHandler)->SetGfaThreshold(m_Controls->m_ScalarThresholdBox->value()); dynamic_cast(m_TrackingHandler)->SetOdfThreshold(m_Controls->m_OdfCutoffBox->value()); dynamic_cast(m_TrackingHandler)->SetSharpenOdfs(m_Controls->m_SharpenOdfsBox->isChecked()); } else if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(m_InputImageNodes.at(0)->GetData())) ) { - if ( m_Controls->m_ForestBox->GetSelectedNode().IsNull() ) + if ( m_Controls->m_ForestSelectionWidget->GetSelectedNode().IsNull() ) { QMessageBox::information(nullptr, "Information", "Not random forest for machine learning based tractography (raw dMRI tractography) selected. Did you accidentally select the raw diffusion-weighted image in the datamanager?"); StartStopTrackingGui(false); return; } if (m_TrackingHandler==nullptr) { - mitk::TractographyForest::Pointer forest = dynamic_cast(m_Controls->m_ForestBox->GetSelectedNode()->GetData()); + mitk::TractographyForest::Pointer forest = dynamic_cast(m_Controls->m_ForestSelectionWidget->GetSelectedNode()->GetData()); mitk::Image::Pointer dwi = dynamic_cast(m_InputImageNodes.at(0)->GetData()); std::vector< std::vector< ItkFloatImageType::Pointer > > additionalFeatureImages; additionalFeatureImages.push_back(std::vector< ItkFloatImageType::Pointer >()); for (auto img : m_AdditionalInputImages) { ItkFloatImageType::Pointer itkimg = ItkFloatImageType::New(); mitk::CastToItkImage(img, itkimg); additionalFeatureImages.at(0).push_back(itkimg); } bool forest_valid = false; if (forest->GetNumFeatures()>=100) { int num_previous_directions = (forest->GetNumFeatures() - (100 + additionalFeatureImages.at(0).size()))/3; m_TrackingHandler = new mitk::TrackingHandlerRandomForest<6, 100>(); dynamic_cast*>(m_TrackingHandler)->AddDwi(dwi); dynamic_cast*>(m_TrackingHandler)->SetAdditionalFeatureImages(additionalFeatureImages); dynamic_cast*>(m_TrackingHandler)->SetForest(forest); dynamic_cast*>(m_TrackingHandler)->SetNumPreviousDirections(num_previous_directions); forest_valid = dynamic_cast*>(m_TrackingHandler)->IsForestValid(); } else { int num_previous_directions = (forest->GetNumFeatures() - (28 + additionalFeatureImages.at(0).size()))/3; m_TrackingHandler = new mitk::TrackingHandlerRandomForest<6, 28>(); dynamic_cast*>(m_TrackingHandler)->AddDwi(dwi); dynamic_cast*>(m_TrackingHandler)->SetAdditionalFeatureImages(additionalFeatureImages); dynamic_cast*>(m_TrackingHandler)->SetForest(forest); dynamic_cast*>(m_TrackingHandler)->SetNumPreviousDirections(num_previous_directions); forest_valid = dynamic_cast*>(m_TrackingHandler)->IsForestValid(); } if (!forest_valid) { QMessageBox::information(nullptr, "Information", "Random forest is invalid. The forest signatue does not match the parameters of TrackingHandlerRandomForest."); StartStopTrackingGui(false); return; } } } else { if (m_Controls->m_ModeBox->currentIndex()==1) { QMessageBox::information(nullptr, "Information", "Probabilstic tractography is not implemented for peak images."); StartStopTrackingGui(false); return; } try { if (m_TrackingHandler==nullptr) { typedef mitk::ImageToItk< mitk::TrackingHandlerPeaks::PeakImgType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(m_InputImages.at(0)); caster->SetCopyMemFlag(true); caster->Update(); mitk::TrackingHandlerPeaks::PeakImgType::Pointer itkImg = caster->GetOutput(); m_TrackingHandler = new mitk::TrackingHandlerPeaks(); dynamic_cast(m_TrackingHandler)->SetPeakImage(itkImg); } dynamic_cast(m_TrackingHandler)->SetPeakThreshold(m_Controls->m_ScalarThresholdBox->value()); } catch(...) { QMessageBox::information(nullptr, "Error", "Peak tracker could not be initialized. Is your input image in the correct format (4D float image, peaks in the 4th dimension)?"); StartStopTrackingGui(false); return; } } m_TrackingHandler->SetFlipX(m_Controls->m_FlipXBox->isChecked()); m_TrackingHandler->SetFlipY(m_Controls->m_FlipYBox->isChecked()); m_TrackingHandler->SetFlipZ(m_Controls->m_FlipZBox->isChecked()); m_TrackingHandler->SetInterpolate(m_Controls->m_InterpolationBox->isChecked()); switch (m_Controls->m_ModeBox->currentIndex()) { case 0: m_TrackingHandler->SetMode(mitk::TrackingDataHandler::MODE::DETERMINISTIC); break; case 1: m_TrackingHandler->SetMode(mitk::TrackingDataHandler::MODE::PROBABILISTIC); break; default: m_TrackingHandler->SetMode(mitk::TrackingDataHandler::MODE::DETERMINISTIC); } if (m_Controls->m_InteractiveBox->isChecked()) { m_Tracker->SetSeedPoints(m_SeedPoints); } - else if (m_Controls->m_SeedImageBox->GetSelectedNode().IsNotNull()) + else if (m_Controls->m_SeedImageSelectionWidget->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); - mitk::CastToItkImage(dynamic_cast(m_Controls->m_SeedImageBox->GetSelectedNode()->GetData()), mask); + mitk::CastToItkImage(dynamic_cast(m_Controls->m_SeedImageSelectionWidget->GetSelectedNode()->GetData()), mask); m_Tracker->SetSeedImage(mask); } - if (m_Controls->m_MaskImageBox->GetSelectedNode().IsNotNull()) + if (m_Controls->m_MaskImageSelectionWidget->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); - mitk::CastToItkImage(dynamic_cast(m_Controls->m_MaskImageBox->GetSelectedNode()->GetData()), mask); + mitk::CastToItkImage(dynamic_cast(m_Controls->m_MaskImageSelectionWidget->GetSelectedNode()->GetData()), mask); m_Tracker->SetMaskImage(mask); } - if (m_Controls->m_StopImageBox->GetSelectedNode().IsNotNull()) + if (m_Controls->m_StopImageSelectionWidget->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); - mitk::CastToItkImage(dynamic_cast(m_Controls->m_StopImageBox->GetSelectedNode()->GetData()), mask); + mitk::CastToItkImage(dynamic_cast(m_Controls->m_StopImageSelectionWidget->GetSelectedNode()->GetData()), mask); m_Tracker->SetStoppingRegions(mask); } - if (m_Controls->m_TargetImageBox->GetSelectedNode().IsNotNull()) + if (m_Controls->m_TargetImageSelectionWidget->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); - mitk::CastToItkImage(dynamic_cast(m_Controls->m_TargetImageBox->GetSelectedNode()->GetData()), mask); + mitk::CastToItkImage(dynamic_cast(m_Controls->m_TargetImageSelectionWidget->GetSelectedNode()->GetData()), mask); m_Tracker->SetTargetRegions(mask); } - if (m_Controls->m_PriorImageBox->GetSelectedNode().IsNotNull()) + if (m_Controls->m_PriorImageSelectionWidget->GetSelectedNode().IsNotNull()) { - if (m_LastPrior!=m_Controls->m_PriorImageBox->GetSelectedNode()->GetUID() || m_TrackingPriorHandler==nullptr) + if (m_LastPrior!=m_Controls->m_PriorImageSelectionWidget->GetSelectedNode()->GetUID() || m_TrackingPriorHandler==nullptr) { typedef mitk::ImageToItk< mitk::TrackingHandlerPeaks::PeakImgType > CasterType; CasterType::Pointer caster = CasterType::New(); - caster->SetInput(dynamic_cast(m_Controls->m_PriorImageBox->GetSelectedNode()->GetData())); + caster->SetInput(dynamic_cast(m_Controls->m_PriorImageSelectionWidget->GetSelectedNode()->GetData())); caster->SetCopyMemFlag(true); caster->Update(); mitk::TrackingHandlerPeaks::PeakImgType::Pointer itkImg = caster->GetOutput(); m_TrackingPriorHandler = new mitk::TrackingHandlerPeaks(); dynamic_cast(m_TrackingPriorHandler)->SetPeakImage(itkImg); dynamic_cast(m_TrackingPriorHandler)->SetPeakThreshold(0.0); - m_LastPrior = m_Controls->m_PriorImageBox->GetSelectedNode()->GetUID(); + m_LastPrior = m_Controls->m_PriorImageSelectionWidget->GetSelectedNode()->GetUID(); } m_TrackingPriorHandler->SetInterpolate(m_Controls->m_InterpolationBox->isChecked()); m_TrackingPriorHandler->SetMode(mitk::TrackingDataHandler::MODE::DETERMINISTIC); m_Tracker->SetTrackingPriorHandler(m_TrackingPriorHandler); m_Tracker->SetTrackingPriorWeight(m_Controls->m_PriorWeightBox->value()); m_Tracker->SetTrackingPriorAsMask(m_Controls->m_PriorAsMaskBox->isChecked()); m_Tracker->SetIntroduceDirectionsFromPrior(m_Controls->m_NewDirectionsFromPriorBox->isChecked()); } - else if (m_Controls->m_PriorImageBox->GetSelectedNode().IsNull()) + else if (m_Controls->m_PriorImageSelectionWidget->GetSelectedNode().IsNull()) m_Tracker->SetTrackingPriorHandler(nullptr); - if (m_Controls->m_ExclusionImageBox->GetSelectedNode().IsNotNull()) + if (m_Controls->m_ExclusionImageSelectionWidget->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); - mitk::CastToItkImage(dynamic_cast(m_Controls->m_ExclusionImageBox->GetSelectedNode()->GetData()), mask); + mitk::CastToItkImage(dynamic_cast(m_Controls->m_ExclusionImageSelectionWidget->GetSelectedNode()->GetData()), mask); m_Tracker->SetExclusionRegions(mask); } // Endpoint constraints switch (m_Controls->m_EpConstraintsBox->currentIndex()) { case 0: m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::NONE); m_Tracker->SetTargetRegions(nullptr); break; case 1: m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_TARGET); break; case 2: m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_TARGET_LABELDIFF); break; case 3: m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_SEED_AND_TARGET); break; case 4: m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::MIN_ONE_EP_IN_TARGET); break; case 5: m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::ONE_EP_IN_TARGET); break; case 6: m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::NO_EP_IN_TARGET); break; } - if (m_Tracker->GetEndpointConstraint()!=itk::StreamlineTrackingFilter::EndpointConstraints::NONE && m_Controls->m_TargetImageBox->GetSelectedNode().IsNull()) + if (m_Tracker->GetEndpointConstraint()!=itk::StreamlineTrackingFilter::EndpointConstraints::NONE && m_Controls->m_TargetImageSelectionWidget->GetSelectedNode().IsNull()) { QMessageBox::information(nullptr, "Error", "Endpoint constraints are used but no target image is set!"); StartStopTrackingGui(false); return; } else if (m_Tracker->GetEndpointConstraint()==itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_SEED_AND_TARGET - && (m_Controls->m_SeedImageBox->GetSelectedNode().IsNull()|| m_Controls->m_TargetImageBox->GetSelectedNode().IsNull()) ) + && (m_Controls->m_SeedImageSelectionWidget->GetSelectedNode().IsNull()|| m_Controls->m_TargetImageSelectionWidget->GetSelectedNode().IsNull()) ) { QMessageBox::information(nullptr, "Error", "Endpoint constraint EPS_IN_SEED_AND_TARGET is used but no target or no seed image is set!"); StartStopTrackingGui(false); return; } m_Tracker->SetInterpolateMasks(m_Controls->m_MaskInterpolationBox->isChecked()); m_Tracker->SetVerbose(!m_Controls->m_InteractiveBox->isChecked()); m_Tracker->SetSeedsPerVoxel(m_Controls->m_SeedsPerVoxelBox->value()); m_Tracker->SetStepSize(m_Controls->m_StepSizeBox->value()); m_Tracker->SetSamplingDistance(m_Controls->m_SamplingDistanceBox->value()); m_Tracker->SetUseStopVotes(m_Controls->m_StopVotesBox->isChecked()); m_Tracker->SetOnlyForwardSamples(m_Controls->m_FrontalSamplesBox->isChecked()); m_Tracker->SetTrialsPerSeed(m_Controls->m_TrialsPerSeedBox->value()); m_Tracker->SetMaxNumTracts(m_Controls->m_NumFibersBox->value()); m_Tracker->SetNumberOfSamples(m_Controls->m_NumSamplesBox->value()); m_Tracker->SetTrackingHandler(m_TrackingHandler); m_Tracker->SetLoopCheck(m_Controls->m_LoopCheckBox->value()); m_Tracker->SetAngularThreshold(m_Controls->m_AngularThresholdBox->value()); m_Tracker->SetMinTractLength(m_Controls->m_MinTractLengthBox->value()); m_Tracker->SetUseOutputProbabilityMap(m_Controls->m_OutputProbMap->isChecked()); m_ParentNode = m_InputImageNodes.at(0); m_TrackingThread.start(QThread::LowestPriority); } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingViewControls.ui index 0ea54e5f85..0c3d0f1200 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingViewControls.ui @@ -1,1575 +1,1491 @@ QmitkStreamlineTrackingViewControls 0 0 453 859 0 0 QmitkTemplate QCommandLinkButton:disabled { border: none; } QGroupBox { background-color: transparent; } 3 3 0 40 QFrame::NoFrame QFrame::Raised 0 15 0 0 6 15 true 0 0 true QFrame::NoFrame QFrame::Raised 0 0 0 0 Input Image. ODF, tensor and peak images are currently supported. Input Image: Input Image. ODF, tensor, peak, and, in case of ML tractography, raw diffusion-weighted images are currently supported. <html><head/><body><p><span style=" color:#ff0000;">select image in data-manager</span></p></body></html> true Tractography Forest: - - - Random forest for machine learning based tractography. - - - QComboBox::AdjustToMinimumContentsLength + + + - - - - - - true Stop tractography and return all fibers reconstructed until now. Stop Tractography false Start Tractography 0 0 0 0 75 true 0 0 0 421 - 262 + 254 Seeding Specify how, where and how many tractography seed points are placed. QFrame::NoFrame QFrame::Raised 0 0 0 0 QFrame::NoFrame QFrame::Raised 0 0 0 0 - Number of seed points equally distributed around selected position. + Number of seed points equally distributed around selected position. 1 9999999 50 Radius: Seedpoints are equally distributed within a sphere centered at the selected position with the specified radius (in mm). 2 50.000000000000000 0.100000000000000 2.000000000000000 Num. Seeds: true When checked, parameter changes cause instant retracking while in interactive mode. Update on Parameter Change true QFrame::NoFrame QFrame::Raised 0 0 0 0 Try each seed N times until a valid streamline is obtained (only for probabilistic tractography). Minimum fiber length (in mm) 1 999 10 Trials Per Seed: Max. Num. Fibers: Tractography is stopped after the desired number of fibers is reached, even before all seed points are processed (-1 means no limit). -1 999999999 -1 QFrame::NoFrame QFrame::Raised 0 0 0 0 - - + + - Number of seed points placed in each voxel. - - - 1 + - - 9999999 + + Seeds per Voxel: - - + + - Seeds per Voxel: + Seed Image: - - + + - Seed points are only placed inside the regions defined in the seed image. If no seed image is selected, the whole image is seeded. + Number of seed points placed in each voxel. - - QComboBox::AdjustToMinimumContentsLength + + 1 + + + 9999999 - - - - - - - - - - - - - Seed Image: + + + + true Dynamically pick a seed location by click into image. Enable Interactive Tractography Qt::Vertical 20 40 0 0 435 - 224 + 232 ROI Constraints Specify various ROI and mask images to constrain the tractography process. Mask Image: - - - - Fibers that enter a region defined in this image will stop immediately. - - - QComboBox::AdjustToMinimumContentsLength - - - - - - - - - Select which fibers should be accepted or rejected based on the location of their endpoints. QComboBox::AdjustToMinimumContentsLength No Constraints on EP locations Both EPs in Target Image Both EPs in Target Image But Different Label One EP in Seed Image and One EP in Target Image At Least One EP in Target Image Exactly One EP in Target Image No EP in Target Image Endpoint Constraints: Stop ROI Image: - - - - The target image is used for the endpoint constraint strategy defined above. - - - QComboBox::AdjustToMinimumContentsLength - - - - - - - - - Exclusion ROI Image: - - - - Fibers that leave the regions defined in this image will stop immediately. - - - QComboBox::AdjustToMinimumContentsLength - - - - - - - - - Target ROI Image: - - - - Fibers that enter a region defined in this image will be discarded. - - - QComboBox::AdjustToMinimumContentsLength - - - - - - - - - Qt::Vertical 20 40 + + + + + + + + + + + + 0 0 421 - 351 + 343 Tractography Parameters Specify the behavior of the tractography at each streamline integration step (step size, deterministic/probabilistic, ...). Qt::Vertical 20 40 f=1 + g=0 means FACT (depending on the chosen interpolation). f=0 and g=1 means TEND (disable interpolation for this mode!). 2 1.000000000000000 0.100000000000000 0.000000000000000 Toggle between deterministic and probabilistic tractography. Some modes might not be available for all types of tractography. Deterministic Probabilistic Cutoff: FA/GFA Image: Mode: Angular Threshold: Step size (in voxels) 2 0.010000000000000 10.000000000000000 0.100000000000000 0.500000000000000 Maximum allowed angular SDTEV over 4 voxel lengths. Default: no loop check. -1 180 -1 - - - - If an image is selected, the stopping criterion is not calculated from the input image but instead the selected image is used. - - - QComboBox::AdjustToMinimumContentsLength - - - - - - - - - Step Size: Additional threshold on the ODF magnitude. This is useful in case of CSD fODF tractography. For fODFs a good default value is 0.1, for normalized dODFs, e.g. Q-ball ODFs, this threshold should be very low (0.00025) or 0. 5 1.000000000000000 0.100000000000000 0.000250000000000 f=1 + g=0 means FACT (depending on the chosen interpolation). f=0 and g=1 means TEND (disable interpolation for this mode!). 2 1.000000000000000 0.100000000000000 1.000000000000000 If you are using dODF images as input, it is advisable to sharpen the ODFs (min-max normalize and raise to the power of 4). This is not necessary for CSD fODFs, since they are naturally much sharper. f parameter of tensor tractography. f=1 + g=0 means FACT (depending on the chosen interpolation). f=0 and g=1 means TEND (disable interpolation for this mode!). f: Min. Tract Length: Threshold on peak magnitude, FA, GFA, ... 5 1.000000000000000 0.100000000000000 0.100000000000000 ODF Cutoff: Minimum tract length in mm. Shorter fibers are discarded. Minimum fiber length (in mm) 1 999.000000000000000 1.000000000000000 20.000000000000000 Loop Check: Angular threshold between two steps (in degree). Default: 90° * step_size -1 90 1 -1 g: Sharpen ODFs: + + + 0 0 435 - 224 + 232 Tractography Prior Restrict to Prior: Weight: Peak Image: Weighting factor between prior and data. 1.000000000000000 0.100000000000000 0.500000000000000 Restrict tractography to regions where the prior is valid. true - - - - - - - Qt::Vertical 20 40 New Directions from Prior: If unchecked, the prior cannot create directions where there are none in the data. true + + + 0 0 435 - 224 + 232 Neighborhood Sampling Specify if and how information about the current streamline neighborhood should be used. Only neighborhood samples in front of the current streamline position are considered. Use Only Frontal Samples false If checked, the majority of sampling points has to place a stop-vote for the streamline to terminate. If not checked, all sampling positions have to vote for a streamline termination. Use Stop-Votes false QFrame::NoFrame QFrame::Raised 0 0 0 0 Num. Samples: Number of neighborhood samples that are used to determine the next fiber progression direction. 50 Sampling Distance: Sampling distance (in voxels) 2 10.000000000000000 0.100000000000000 0.250000000000000 Qt::Vertical 20 40 0 0 435 - 224 + 232 Data Handling Specify interpolation and direction flips. QFrame::NoFrame QFrame::Raised 0 0 0 0 Trilinearly interpolate the input image used for tractography. Interpolate Tractography Data true Trilinearly interpolate the ROI images used to constrain the tractography. Interpolate ROI Images true QFrame::NoFrame QFrame::Raised 0 0 0 0 QFrame::NoFrame QFrame::Raised 0 0 0 0 Internally flips progression directions. This might be necessary depending on the input data. x Internally flips progression directions. This might be necessary depending on the input data. y Internally flips progression directions. This might be necessary depending on the input data. z Flip directions: Qt::Vertical 20 40 0 0 435 - 224 + 232 Output and Postprocessing Specify the tractography output (streamlines or probability maps) and postprocessing steps. QFrame::NoFrame QFrame::Raised 0 0 0 0 Compress fibers using the specified error constraint. Compress Fibers true Qt::StrongFocus Lossy fiber compression. Recommended for large tractograms. Maximum error in mm. 3 10.000000000000000 0.010000000000000 0.100000000000000 Output map with voxel-wise visitation counts instead of streamlines. Output Probability Map false Qt::Vertical 20 40 - QmitkDataStorageComboBox - QComboBox -
QmitkDataStorageComboBox.h
- - - QmitkDataStorageComboBoxWithSelectNone - QComboBox -
QmitkDataStorageComboBoxWithSelectNone.h
+ QmitkSingleNodeSelectionWidget + QWidget +
QmitkSingleNodeSelectionWidget.h
+ 1