diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStochasticTractographyFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStochasticTractographyFilter.h index 03ef745c80..c8770d5a56 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStochasticTractographyFilter.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStochasticTractographyFilter.h @@ -1,304 +1,305 @@ /*=================================================================== 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 __itkStochasticTractographyFilter_h__ #define __itkStochasticTractographyFilter_h__ #include "itkImageToImageFilter.h" #include "vnl/vnl_random.h" #include "vnl/vnl_vector_fixed.h" #include "vnl/vnl_matrix.h" #include "itkArray.h" #include "itkVectorContainer.h" #include "vnl/algo/vnl_qr.h" #include "itkVariableLengthVector.h" #include "StochasticTracking/itkSlowPolyLineParametricPath.h" #include "itkSimpleFastMutexLock.h" #include "itkRealTimeClock.h" #include "itkDiffusionTensor3D.h" #include namespace itk{ /** * \brief Performs probabilistic streamline tracking on the input dwi. */ /**Types for Probability Distribution **/ typedef Image< Array< double >, 3 > ProbabilityDistributionImageType; template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage > class ITK_EXPORT StochasticTractographyFilter : public ImageToImageFilter< TInputDWIImage, TOutputConnectivityImage >{ public: typedef StochasticTractographyFilter Self; typedef ImageToImageFilter< TInputDWIImage, TOutputConnectivityImage > Superclass; typedef SmartPointer< Self > Pointer; typedef SmartPointer< const Self > ConstPointer; - itkNewMacro(Self); + itkNewMacro(Self) itkTypeMacro( StochasticTractographyFilter, - ImageToImageFilter ); + ImageToImageFilter ) /** Types for the DWI Input Image **/ typedef TInputDWIImage InputDWIImageType; /** Types for the Connectivity Output Image**/ typedef TOutputConnectivityImage OutputConnectivityImageType; /** Types for the Mask Image **/ typedef TInputWhiteMatterProbabilityImage InputWhiteMatterProbabilityImageType; /** Tract Types **/ typedef SlowPolyLineParametricPath< 3 > TractType; /** Types for the TractContainer **/ typedef VectorContainer< unsigned int, typename TractType::Pointer > TractContainerType; /** Types for Tensor Output Image **/ typedef Image< DiffusionTensor3D< double >, 3 > OutputTensorImageType; /** Types for the Image-wide Magnetic Field Gradient Directions **/ typedef VectorContainer< unsigned int, vnl_vector_fixed< double, 3 > > GradientDirectionContainerType; /** Types for the Image-wide bValues **/ typedef double bValueType; typedef VectorContainer< unsigned int, bValueType > bValueContainerType; /** Types for the Measurement Frame of the Gradients **/ typedef vnl_matrix_fixed< double, 3, 3 > MeasurementFrameType; /** Type for the sample directions **/ typedef VectorContainer< unsigned int, vnl_vector_fixed< double, 3 > > TractOrientationContainerType; /** the number of Tracts to generate **/ - itkSetMacro( TotalTracts, unsigned int); - itkGetMacro( TotalTracts, unsigned int); + itkSetMacro( TotalTracts, unsigned int) + itkGetMacro( TotalTracts, unsigned int) /** the maximum length of Tract **/ - itkSetMacro( MaxTractLength, unsigned int ); - itkGetMacro( MaxTractLength, unsigned int ); + itkSetMacro( MaxTractLength, unsigned int ) + itkGetMacro( MaxTractLength, unsigned int ) /** Set/Get bvalues **/ - itkSetConstObjectMacro( bValues, bValueContainerType ); - itkGetConstObjectMacro( bValues, bValueContainerType ); + itkSetConstObjectMacro( bValues, bValueContainerType ) + itkGetConstObjectMacro( bValues, bValueContainerType ) /** Set/Get of gradient directions **/ - itkSetConstObjectMacro( Gradients, GradientDirectionContainerType ); - itkGetConstObjectMacro( Gradients, GradientDirectionContainerType ); + itkSetConstObjectMacro( Gradients, GradientDirectionContainerType ) + itkGetConstObjectMacro( Gradients, GradientDirectionContainerType ) /** Set/Get the White Matter Probability Input image **/ /* At each voxel specifies the probability of a mylinated fiber existing at that location. This probability is interpreted to be the probability that a fiber tract passes through that region. */ - itkSetInputMacro(WhiteMatterProbabilityImage, InputWhiteMatterProbabilityImageType); - itkGetInputMacro(WhiteMatterProbabilityImage, InputWhiteMatterProbabilityImageType); + itkSetInputMacro(WhiteMatterProbabilityImage, InputWhiteMatterProbabilityImageType) + itkGetInputMacro(WhiteMatterProbabilityImage, InputWhiteMatterProbabilityImageType) //overide the built in set input function //we need to create a new cache everytime we change the input image //but we need to preserve it when the input image is the same + using Superclass::SetInput; void SetInput( typename InputDWIImageType::Pointer dwiimagePtr ){ Superclass::SetInput( dwiimagePtr ); //update the likelihood cache this->m_LikelihoodCachePtr = ProbabilityDistributionImageType::New(); this->m_LikelihoodCachePtr->CopyInformation( this->GetInput() ); this->m_LikelihoodCachePtr->SetBufferedRegion( this->GetInput()->GetBufferedRegion() ); this->m_LikelihoodCachePtr->SetRequestedRegion( this->GetInput()->GetRequestedRegion() ); this->m_LikelihoodCachePtr->Allocate(); this->m_CurrentLikelihoodCacheElements = 0; //update the likelihoodcache mutex image this->m_LikelihoodCacheMutexImagePtr = LikelihoodCacheMutexImageType::New(); this->m_LikelihoodCacheMutexImagePtr->CopyInformation( this->GetInput() ); this->m_LikelihoodCacheMutexImagePtr->SetBufferedRegion( this->GetInput()->GetBufferedRegion() ); this->m_LikelihoodCacheMutexImagePtr->SetRequestedRegion( this->GetInput()->GetRequestedRegion() ); this->m_LikelihoodCacheMutexImagePtr->Allocate(); } /** Set/Get the seed index **/ - itkSetMacro( SeedIndex, typename InputDWIImageType::IndexType ); - itkGetMacro( SeedIndex, typename InputDWIImageType::IndexType ); + itkSetMacro( SeedIndex, typename InputDWIImageType::IndexType ) + itkGetMacro( SeedIndex, typename InputDWIImageType::IndexType ) /** Set/Get the list of directions to sample **/ - itkSetConstObjectMacro( SampleDirections, TractOrientationContainerType ); - itkGetConstObjectMacro( SampleDirections, TractOrientationContainerType ); + itkSetConstObjectMacro( SampleDirections, TractOrientationContainerType ) + itkGetConstObjectMacro( SampleDirections, TractOrientationContainerType ) /** Set/Get the Measurement Frame **/ - itkSetMacro( MeasurementFrame, MeasurementFrameType ); - itkGetMacro( MeasurementFrame, MeasurementFrameType ); + itkSetMacro( MeasurementFrame, MeasurementFrameType ) + itkGetMacro( MeasurementFrame, MeasurementFrameType ) /** Set/Get the Maximum Likelihood Cache Size, the max num. of cached voxels **/ - itkSetMacro( MaxLikelihoodCacheSize, unsigned int ); - itkGetMacro( MaxLikelihoodCacheSize, unsigned int ); + itkSetMacro( MaxLikelihoodCacheSize, unsigned int ) + itkGetMacro( MaxLikelihoodCacheSize, unsigned int ) /** Get the Tracts that are generated **/ - itkGetObjectMacro( OutputTractContainer, TractContainerType ); + itkGetObjectMacro( OutputTractContainer, TractContainerType ) /** Get TensorImage **/ - itkGetObjectMacro( OutputTensorImage, OutputTensorImageType ); + itkGetObjectMacro( OutputTensorImage, OutputTensorImageType ) void GenerateData(); void GenerateTractContainerOutput( void ); void GenerateTensorImageOutput( void ); protected: /** Convenience Types used only inside the filter **/ /**Types for the parameters of the Tensor Model **/ typedef vnl_vector_fixed< double, 7 > TensorModelParamType; /**Types for the parameters of the Constrained Model **/ typedef vnl_vector_fixed< double, 6 > ConstrainedModelParamType; /**Type to hold generated DWI values**/ typedef Image< VariableLengthVector< double >, 3 > DWIVectorImageType; /**Types for Probability Distribution **/ typedef Image< Array< double >, 3 > ProbabilityDistributionImageType; /** Types for the Image of Mutexes of the Likelihood distribution **/ typedef Image< SimpleFastMutexLock, 3 > LikelihoodCacheMutexImageType; StochasticTractographyFilter(); virtual ~StochasticTractographyFilter(); /** Load the default Sample Directions**/ void LoadDefaultSampleDirections( void ); /** Randomly chose a neighboring pixel weighted on distance **/ void ProbabilisticallyInterpolate( vnl_random& randomgenerator, const TractType::ContinuousIndexType& cindex, typename InputDWIImageType::IndexType& index); /** Functions and data related to fitting the tensor model at each pixel **/ void UpdateGradientDirections(void); void UpdateTensorModelFittingMatrices( void ); void CalculateTensorModelParameters( const DWIVectorImageType::PixelType& dwivalues, vnl_diag_matrix& W, TensorModelParamType& tensormodelparams); void CalculateConstrainedModelParameters( const TensorModelParamType& tensormodelparams, ConstrainedModelParamType& constrainedmodelparams); void CalculateNoiseFreeDWIFromConstrainedModel( const ConstrainedModelParamType& constrainedmodelparams, DWIVectorImageType::PixelType& noisefreedwi); void CalculateResidualVariance( const DWIVectorImageType::PixelType& noisydwi, const DWIVectorImageType::PixelType& noisefreedwi, const vnl_diag_matrix< double >& W, const unsigned int numberofparameters, double& residualvariance); void CalculateLikelihood( const DWIVectorImageType::PixelType &dwipixel, TractOrientationContainerType::ConstPointer orientations, ProbabilityDistributionImageType::PixelType& likelihood); void CalculatePrior( TractOrientationContainerType::Element v_prev, TractOrientationContainerType::ConstPointer orientations, ProbabilityDistributionImageType::PixelType& prior ); void CalculatePosterior( const ProbabilityDistributionImageType::PixelType& likelihood, const ProbabilityDistributionImageType::PixelType& prior, ProbabilityDistributionImageType::PixelType& posterior); void SampleTractOrientation( vnl_random& randomgenerator, const ProbabilityDistributionImageType::PixelType& posterior, TractOrientationContainerType::ConstPointer orientations, TractOrientationContainerType::Element& choosendirection ); void StochasticTractGeneration( typename InputDWIImageType::ConstPointer dwiimagePtr, typename InputWhiteMatterProbabilityImageType::ConstPointer maskimagePtr, typename InputDWIImageType::IndexType seedindex, unsigned long randomseed, TractType::Pointer tract ); /** Callback routine used by the threading library. This routine just calls the ThreadedGenerateData method after setting the correct region for this thread. **/ static ITK_THREAD_RETURN_TYPE StochasticTractGenerationCallback( void *arg ); struct StochasticTractGenerationCallbackStruct{ Pointer Filter; }; /** Thread Safe Function to check/update an entry in the likelihood cache **/ ProbabilityDistributionImageType::PixelType& AccessLikelihoodCache( typename InputDWIImageType::IndexType index ); /** Thread Safe Function to delegate a tract and obtain a randomseed to start tracking **/ bool DelegateTract(unsigned long& randomseed); /** Function to write a tract to the connectivity map **/ void TractContainerToConnectivityMap(TractContainerType::Pointer tractcontainer); /** Thread Safe Function to store a tract to a TractContainer **/ void StoreTract(TractType::Pointer tract); /** Randomly samples the existence of a fiber tract in the current voxel **/ bool FiberExistenceTest( vnl_random& randomgenerator, typename InputWhiteMatterProbabilityImageType::ConstPointer wmpimage, typename InputWhiteMatterProbabilityImageType::IndexType index ); MeasurementFrameType m_MeasurementFrame; LikelihoodCacheMutexImageType::Pointer m_LikelihoodCacheMutexImagePtr; unsigned int m_TotalTracts; unsigned int m_MaxTractLength; GradientDirectionContainerType::ConstPointer m_Gradients; GradientDirectionContainerType::Pointer m_TransformedGradients; bValueContainerType::ConstPointer m_bValues; typename InputDWIImageType::IndexType m_SeedIndex; TractOrientationContainerType::ConstPointer m_SampleDirections; //these will be the same for every pixel in the image so //go ahead and do a QR decomposition to optimize the //LS fitting process for estimating the weighing matrix W //in this case we solve instead: //R*Beta = Q'logPhi vnl_matrix< double >* m_A; vnl_qr< double >* m_Aqr; ProbabilityDistributionImageType::Pointer m_LikelihoodCachePtr; unsigned long m_MaxLikelihoodCacheSize; //in Megabytes unsigned long m_MaxLikelihoodCacheElements; //in Elements (Voxels) unsigned long m_CurrentLikelihoodCacheElements; SimpleFastMutexLock m_LikelihoodCacheMutex; RealTimeClock::Pointer m_ClockPtr; unsigned int m_TotalDelegatedTracts; SimpleFastMutexLock m_TotalDelegatedTractsMutex; //unsigned long m_RandomSeed; SimpleFastMutexLock m_OutputImageMutex; TractContainerType::Pointer m_OutputTractContainer; SimpleFastMutexLock m_OutputTractContainerMutex; OutputTensorImageType::Pointer m_OutputTensorImage; vnl_random m_RandomGenerator; }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkStochasticTractographyFilter.txx" #include "StochasticTracking/itkStochasticTractographyFilter_SD.txx" #endif #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStochasticTractographyFilter.txx b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStochasticTractographyFilter.txx index d01b5e2b04..fd2b9670a9 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStochasticTractographyFilter.txx +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStochasticTractographyFilter.txx @@ -1,700 +1,700 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "itkStochasticTractographyFilter.h" #include "vnl/vnl_math.h" #include "vnl/vnl_matrix_fixed.h" #include "vnl/vnl_vector_fixed.h" #include "vnl/vnl_matrix.h" #include "vnl/vnl_sym_matrix.h" #include "vnl/vnl_vector.h" #include "vnl/vnl_diag_matrix.h" #include "vnl/algo/vnl_qr.h" //#include "vnl/algo/vnl_svd.h" #include "vnl/algo/vnl_matrix_inverse.h" //#include "vnl/algo/vnl_symmetric_eigensystem.h" #include "itkSymmetricEigenAnalysis.h" #include "vnl/vnl_transpose.h" #include "itkVariableSizeMatrix.h" #include "itkPathIterator.h" #include "itkImageRegionIterator.h" #include "itkImageRegionConstIterator.h" #include "itkImageRegionConstIteratorWithIndex.h" namespace itk{ template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage > StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage > ::StochasticTractographyFilter(): m_TotalTracts(0),m_MaxTractLength(0),m_Gradients(NULL), m_TransformedGradients(NULL),m_bValues(NULL), m_SampleDirections(NULL), m_A(NULL), m_Aqr(NULL), m_LikelihoodCachePtr(NULL), m_MaxLikelihoodCacheSize(0), m_CurrentLikelihoodCacheElements(0), m_ClockPtr(NULL), m_TotalDelegatedTracts(0), m_OutputTractContainer(NULL){ this->m_SeedIndex[0]=0; this->m_SeedIndex[1]=0; this->m_SeedIndex[2]=0; this->m_MeasurementFrame.set_identity(); - this->SetNumberOfRequiredInputs(2); //Filter needs a DWI image and a Mask Image + this->SetNumberOfRequiredInputs(1); //Filter needs a DWI image and a Mask Image, but check will be outside itkProcessObject and inside GenerateData() m_ClockPtr = RealTimeClock::New(); this->m_RandomGenerator.reseed( ((unsigned long) this->m_ClockPtr->GetTimeInSeconds()) ); //load in default sample directions this->LoadDefaultSampleDirections(); } template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage > StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage > ::~StochasticTractographyFilter(){ delete this->m_A; delete this->m_Aqr; } template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage > void StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage > ::ProbabilisticallyInterpolate( vnl_random& randomgenerator, const TractType::ContinuousIndexType& cindex, typename InputDWIImageType::IndexType& index){ for(int i=0; i<3; i++){ if ((vcl_ceil(cindex[i]+vnl_math::eps)-cindex[i]) < randomgenerator.drand64()) index[i]=(int)vcl_ceil(cindex[i]); else index[i]=(int)vcl_floor(cindex[i]); } } template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage > void StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage > ::UpdateGradientDirections(void){ //the gradient direction is transformed into IJK space //by moving into the image space and then to IJK space this->m_TransformedGradients = GradientDirectionContainerType::New(); unsigned int N = this->m_Gradients->Size(); for(unsigned int i=0; im_MeasurementFrame * this->m_Gradients->GetElement(i); /** The correction to LPS space is not neccessary as of itk 3.2 **/ //g_i[0] = -g_i[0]; //g_i[1] = -g_i[1]; g_i = this->GetInput()->GetDirection().GetInverse() * g_i; this->m_TransformedGradients->InsertElement(i, g_i); } } template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage > void StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage > ::UpdateTensorModelFittingMatrices( void ){ //std::cout<<"UpdateTensorFittingMatrix\n"; //estimate the parameters using linear LS estimation //using convention specified by Salvador //solve for Beta in: logPhi=X*Beta //number of rows of the matrix depends on the number of inputs, //i.e. the number of measurements of the voxel (n) unsigned int N = this->m_TransformedGradients->Size(); if(this->m_A!=NULL) delete this->m_A; this->m_A = new vnl_matrix< double >(N, 7); //potential memory leak here vnl_matrix< double >& A = *(this->m_A); for(unsigned int j=0; j< N ; j++){ GradientDirectionContainerType::Element g = m_TransformedGradients->GetElement(j); const bValueType& b_i = m_bValues->GetElement(j); A(j,0)=1.0; A(j,1)=-1*b_i*(g[0]*g[0]); A(j,2)=-1*b_i*(g[1]*g[1]); A(j,3)=-1*b_i*(g[2]*g[2]); A(j,4)=-1*b_i*(2*g[0]*g[1]); A(j,5)=-1*b_i*(2*g[0]*g[2]); A(j,6)=-1*b_i*(2*g[1]*g[2]); } //Store a QR decomposition to quickly estimate //the weighing matrix for each voxel if(this->m_Aqr!=NULL) delete this->m_Aqr; this->m_Aqr = new vnl_qr< double >(A); //potential memory leak here } template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage > void StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage > ::CalculateTensorModelParameters( const DWIVectorImageType::PixelType& dwivalues, vnl_diag_matrix& W, TensorModelParamType& tensormodelparams){ unsigned int N = this->m_TransformedGradients->Size(); //setup const references for code clarity const vnl_matrix< double >& A = *(this->m_A); const vnl_qr< double >& Aqr = *(this->m_Aqr); //vnl_vector is used because the itk vector is limited in its methods and does not //contain an internal vnl class like VariableSizematrix //also itk_matrix has methods which are compatible with vnl_vectors vnl_vector< double > logPhi( N ); for(unsigned int j=0; j< N ; j++){ //fill up the logPhi vector using log(dwi) values logPhi.put(j, vcl_log(static_cast(dwivalues[j]) + vnl_math::eps)); } /** Find WLS estimate of the parameters of the Tensor model **/ // First estimate W by LS estimation of the intensities //vnl_matrix< double > Q = Aqr.Q(); //vnl_vector< double > QtB = Aqr.Q().transpose()*logPhi; //vnl_vector< double > QTB = Aqr.QtB(logPhi); //vnl_matrix< double > R = Aqr.R(); W = A* vnl_qr< double >(Aqr.R()).solve(Aqr.QtB(logPhi)); //W = A * Aqr.solve(logPhi); for(vnl_diag_matrix< double >::iterator i = W.begin();i!=W.end(); i++){ *i = vcl_exp( *i ); } // Now solve for parameters using the estimated weighing matrix tensormodelparams = vnl_qr< double >((W*A).transpose()*W*A).solve( (W*A).transpose()*W*logPhi); //int a; //tensormodelparams = vnl_qr< double >((W*A)).solve(W*logPhi); } template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage > void StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage > ::CalculateConstrainedModelParameters( const TensorModelParamType& tensormodelparams, ConstrainedModelParamType& constrainedmodelparams){ vnl_sym_matrix< double > D( 3, 0 ); double alpha =0; double beta=0; //set the tensor model parameters into a Diffusion tensor D(0,0) = tensormodelparams[1]; D(0,1) = tensormodelparams[4]; D(0,2) = tensormodelparams[5]; D(1,0) = tensormodelparams[4]; D(1,1) = tensormodelparams[2]; D(1,2) = tensormodelparams[6]; D(2,0) = tensormodelparams[5]; D(2,1) = tensormodelparams[6]; D(2,2) = tensormodelparams[3]; //pass through the no gradient intensity Z_0 and //calculate alpha, beta and v hat (the eigenvector //associated with the largest eigenvalue) vnl_matrix_fixed< double, 3, 3 > S(0.0); vnl_vector_fixed< double, 3 > Lambda(0.0); SymmetricEigenAnalysis< vnl_sym_matrix< double >, vnl_vector_fixed< double, 3 >, vnl_matrix_fixed< double, 3, 3 > > eigensystem( 3 ); eigensystem.ComputeEigenValuesAndVectors( D, Lambda, S ); //need to take abs to get rid of negative eigenvalues alpha = (vcl_abs(Lambda[0]) + vcl_abs(Lambda[1])) / 2; beta = vcl_abs(Lambda[2]) - alpha; constrainedmodelparams[0] = tensormodelparams[0]; constrainedmodelparams[1] = alpha; constrainedmodelparams[2] = beta; constrainedmodelparams[3] = S[2][0]; constrainedmodelparams[4] = S[2][1]; constrainedmodelparams[5] = S[2][2]; } template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage > void StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage > ::CalculateNoiseFreeDWIFromConstrainedModel( const ConstrainedModelParamType& constrainedmodelparams, DWIVectorImageType::PixelType& noisefreedwi){ unsigned int N = this->m_TransformedGradients->Size(); const double& z_0 = constrainedmodelparams[0]; const double& alpha = constrainedmodelparams[1]; const double& beta = constrainedmodelparams[2]; TractOrientationContainerType::Element v_hat( constrainedmodelparams[3], constrainedmodelparams[4], constrainedmodelparams[5]); for(unsigned int i=0; i < N ; i++ ){ const double& b_i = this->m_bValues->GetElement(i); const GradientDirectionContainerType::Element& g_i = this->m_TransformedGradients->GetElement(i); noisefreedwi.SetElement(i, vcl_exp(z_0-(alpha*b_i+beta*b_i*vnl_math_sqr(dot_product(g_i, v_hat))))); } } template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage > void StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage > ::CalculateResidualVariance( const DWIVectorImageType::PixelType& noisydwi, const DWIVectorImageType::PixelType& noisefreedwi, const vnl_diag_matrix< double >& W, const unsigned int numberofparameters, double& residualvariance){ unsigned int N = this->m_TransformedGradients->Size(); residualvariance=0; /** Not sure if we should be taking difference of log or nonlog intensities **/ /** residual variance is too low if we take the difference of log intensities **/ /** perhaps using WLS will correct this problem **/ for(unsigned int i=0; i void StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage > ::CalculateLikelihood( const DWIVectorImageType::PixelType &dwipixel, TractOrientationContainerType::ConstPointer orientations, ProbabilityDistributionImageType::PixelType& likelihood){ unsigned int N = this->m_TransformedGradients->Size(); TensorModelParamType tensorparams( 0.0 ); vnl_diag_matrix< double > W(N,0); ConstrainedModelParamType constrainedparams( 0.0 ); DWIVectorImageType::PixelType noisefreedwi(N); double residualvariance=0; double jointlikelihood=1; CalculateTensorModelParameters( dwipixel, W, tensorparams ); CalculateConstrainedModelParameters( tensorparams, constrainedparams ); CalculateNoiseFreeDWIFromConstrainedModel( constrainedparams, noisefreedwi ); CalculateResidualVariance( dwipixel, noisefreedwi, W, 6, residualvariance ); for(unsigned int i=0; i < orientations->Size(); i++){ /** Vary the entry corresponding to the estimated Tract orientation over the selected sample directions, while preserving the best estimate for the other parameters **/ TractOrientationContainerType::Element currentdir = orientations->GetElement(i); /** Incorporate the current sample direction into the secondary parameters **/ constrainedparams[3]=currentdir[0]; constrainedparams[4]=currentdir[1]; constrainedparams[5]=currentdir[2]; /** Obtain the estimated intensity for this choice of Tract direction **/ CalculateNoiseFreeDWIFromConstrainedModel(constrainedparams, noisefreedwi); jointlikelihood = 1.0; for(unsigned int j=0; j void StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage > ::CalculatePrior( TractOrientationContainerType::Element v_prev, TractOrientationContainerType::ConstPointer orientations, ProbabilityDistributionImageType::PixelType& prior ){ const double gamma = 1; for(unsigned int i=0; i < orientations->Size(); i++){ if(v_prev.squared_magnitude()==0){ prior[i]=1.0; } else{ prior[i] = dot_product(orientations->GetElement(i),v_prev);; if(prior[i]<0){ prior[i]=0; } else{ prior[i]=vcl_pow(prior[i],gamma); } } } } template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage > void StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage > ::CalculatePosterior( const ProbabilityDistributionImageType::PixelType& likelihood, const ProbabilityDistributionImageType::PixelType& prior, ProbabilityDistributionImageType::PixelType& posterior){ double sum=0; for(unsigned int i=0; i void StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage > ::SampleTractOrientation( vnl_random& randomgenerator, const ProbabilityDistributionImageType::PixelType& posterior, TractOrientationContainerType::ConstPointer orientations, TractOrientationContainerType::Element& choosendirection ){ double randomnum = randomgenerator.drand64(); int i=0; double cumsum=0; //will crash in the unlikely case that 0 was choosen as the randomnum while(cumsum < randomnum){ cumsum+=posterior[i]; i++; } choosendirection = orientations->GetElement(i-1); //std::cout<< "cumsum: " << cumsum< bool StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage > ::FiberExistenceTest( vnl_random& randomgenerator, typename InputWhiteMatterProbabilityImageType::ConstPointer wmpimage, typename InputWhiteMatterProbabilityImageType::IndexType index ){ double randomnum = randomgenerator.drand64(); if( randomnum < wmpimage->GetPixel( index ) ) return true; else return false; } //the seedindex is in continuous IJK coordinates template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage > void StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage > ::StochasticTractGeneration( typename InputDWIImageType::ConstPointer dwiimagePtr, typename InputWhiteMatterProbabilityImageType::ConstPointer wmpimagePtr, typename InputDWIImageType::IndexType seedindex, unsigned long randomseed, TractType::Pointer tract){ TractType::ContinuousIndexType cindex_curr = seedindex; typename InputDWIImageType::IndexType index_curr = {{0,0,0}}; ProbabilityDistributionImageType::PixelType prior_curr(this->m_SampleDirections->Size()); ProbabilityDistributionImageType::PixelType posterior_curr(this->m_SampleDirections->Size()); TractOrientationContainerType::Element v_curr(0,0,0); TractOrientationContainerType::Element v_prev(0,0,0); tract->Initialize(); vnl_random randomgenerator(randomseed); //std::cout<m_MaxTractLength; j++){ this->ProbabilisticallyInterpolate( randomgenerator, cindex_curr, index_curr ); if(!dwiimagePtr->GetLargestPossibleRegion().IsInside(index_curr)){ break; } if( FiberExistenceTest( randomgenerator, wmpimagePtr, index_curr ) ){ tract->AddVertex(cindex_curr); this->CalculatePrior( v_prev, this->m_SampleDirections, prior_curr); const ProbabilityDistributionImageType::PixelType& cachelikelihood_curr = this->AccessLikelihoodCache(index_curr); if( cachelikelihood_curr.GetSize() != 0){ //use the cached direction this->CalculatePosterior( cachelikelihood_curr, prior_curr, posterior_curr); } else{ //do the likelihood calculation and discard //std::cout<<"Cache Miss!\n"; ProbabilityDistributionImageType::PixelType likelihood_curr_temp(this->m_SampleDirections->Size()); this->CalculateLikelihood(static_cast< DWIVectorImageType::PixelType >( dwiimagePtr->GetPixel(index_curr)), this->m_SampleDirections, likelihood_curr_temp); this->CalculatePosterior( likelihood_curr_temp, prior_curr, posterior_curr); } this->SampleTractOrientation(randomgenerator, posterior_curr, this->m_SampleDirections, v_curr); //takes into account voxels of different sizes //converts from a step length of 1 mm to the corresponding length in IJK space const typename InputDWIImageType::SpacingType& spacing = dwiimagePtr->GetSpacing(); cindex_curr[0]+=v_curr[0]/spacing[0]; cindex_curr[1]+=v_curr[1]/spacing[1]; cindex_curr[2]+=v_curr[2]/spacing[2]; v_prev=v_curr; } else{ //fiber doesn't exist in this voxel //std::cout<<"Stopped Tracking: No Fiber in this Voxel\n"; break; } } } template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage > void StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage > ::GenerateTractContainerOutput(){ //allocate tractcontainer this->m_OutputTractContainer = TractContainerType::New(); this->UpdateGradientDirections(); this->UpdateTensorModelFittingMatrices(); this->m_TotalDelegatedTracts = 0; //calculate the number of voxels to cache from Megabyte memory size limit ProbabilityDistributionImageType::PixelType element(this->GetSampleDirections()->Size()); unsigned long elementsize = sizeof(ProbabilityDistributionImageType::PixelType) + sizeof(double)*element.Size(); this->m_MaxLikelihoodCacheElements = (this->m_MaxLikelihoodCacheSize*1048576)/elementsize; std::cout << "MaxLikelhoodCacheElements: " << this->m_MaxLikelihoodCacheElements << std::endl; //setup the multithreader StochasticTractGenerationCallbackStruct data; data.Filter = this; this->GetMultiThreader()->SetSingleMethod( StochasticTractGenerationCallback, &data ); this->GetMultiThreader()->SetNumberOfThreads(this->GetNumberOfThreads()); std::cout<<"Number of Threads: " << this->GetMultiThreader()->GetNumberOfThreads() << std::endl; //start the multithreaded execution this->GetMultiThreader()->SingleMethodExecute(); std::cout<< "CurrentLikelihoodCacheElements: " << this->m_CurrentLikelihoodCacheElements << std::endl; } template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage > void StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage > ::GenerateData(){ //Generate the tracts this->GenerateTractContainerOutput(); //allocate outputs this->AllocateOutputs(); //write tracts to output image this->TractContainerToConnectivityMap(this->m_OutputTractContainer); } template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage > ITK_THREAD_RETURN_TYPE StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage > ::StochasticTractGenerationCallback( void *arg ) { StochasticTractGenerationCallbackStruct* str= (StochasticTractGenerationCallbackStruct *) (((MultiThreader::ThreadInfoStruct *)(arg))->UserData); typename InputDWIImageType::ConstPointer inputDWIImagePtr = str->Filter->GetInput(); typename InputWhiteMatterProbabilityImageType::ConstPointer inputWMPImage = str->Filter->GetWhiteMatterProbabilityImage(); unsigned long randomseed=0; while(str->Filter->DelegateTract(randomseed)){ //std::cout<Filter->StochasticTractGeneration( inputDWIImagePtr, inputWMPImage, str->Filter->GetSeedIndex(), randomseed, tract); //only store tract if it is of nonzero length if( tract->GetVertexList()->Size() > 4 ){ //std::cout<<"Storing tract\n"; str->Filter->StoreTract(tract); } else{ //std::cout<<"Not Storing Tract\n"; } } return ITK_THREAD_RETURN_VALUE; } template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage > ProbabilityDistributionImageType::PixelType& StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage > ::AccessLikelihoodCache( typename InputDWIImageType::IndexType index ) { this->m_LikelihoodCacheMutexImagePtr->GetPixel(index).Lock(); ProbabilityDistributionImageType::PixelType& likelihood = m_LikelihoodCachePtr->GetPixel( index ); typename InputDWIImageType::ConstPointer inputDWIImagePtr = this->GetInput(); if( likelihood.GetSize() !=0){ //entry found in cache this->m_LikelihoodCacheMutexImagePtr->GetPixel(index).Unlock(); return likelihood; } //we need to lock m_CurrentLikelihoodCacheElements as well but not crucial right now else if( this->m_CurrentLikelihoodCacheElements < this->m_MaxLikelihoodCacheElements ){ //entry not found in cache but we have space to store it likelihood.SetSize(this->m_SampleDirections->Size()); this->CalculateLikelihood(static_cast< DWIVectorImageType::PixelType >( inputDWIImagePtr->GetPixel(index)), this->m_SampleDirections, likelihood); this->m_CurrentLikelihoodCacheElements++; this->m_LikelihoodCacheMutexImagePtr->GetPixel(index).Unlock(); return likelihood; } else{ //entry not found in cache and no space to store it this->m_LikelihoodCacheMutexImagePtr->GetPixel(index).Unlock(); return likelihood; } this->m_LikelihoodCacheMutexImagePtr->GetPixel(index).Unlock(); // dummy return likelihood; } template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage > bool StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage > ::DelegateTract(unsigned long& randomseed){ bool success = false; this->m_TotalDelegatedTractsMutex.Lock(); if(this->m_TotalDelegatedTracts < this->m_TotalTracts){ randomseed = this->m_RandomGenerator.lrand32(); this->m_TotalDelegatedTracts++; success = true; //a tract was successfully delegated } else success = false; //all tracts have been delegated this->m_TotalDelegatedTractsMutex.Unlock(); return success; } template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage > void StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage > ::TractContainerToConnectivityMap(TractContainerType::Pointer tractcontainer){ //zero the output image typename OutputConnectivityImageType::Pointer outputPtr = this->GetOutput(); outputPtr->FillBuffer(0); typedef PathIterator< OutputConnectivityImageType, TractType > OutputTractIteratorType; for(unsigned int i=0; iSize(); i++ ){ TractType::Pointer tract = tractcontainer->GetElement(i); //std::cout<< tract->EndOfInput() < void StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage > ::StoreTract(TractType::Pointer tract){ this->m_OutputTractContainerMutex.Lock(); this->m_OutputTractContainer->InsertElement( this->m_OutputTractContainer->Size(), tract); this->m_OutputTractContainerMutex.Unlock(); } template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage, class TOutputConnectivityImage > void StochasticTractographyFilter< TInputDWIImage, TInputWhiteMatterProbabilityImage, TOutputConnectivityImage > ::GenerateTensorImageOutput(void){ this->UpdateGradientDirections(); this->UpdateTensorModelFittingMatrices(); //allocate the tensor image this->m_OutputTensorImage = OutputTensorImageType::New(); m_OutputTensorImage->CopyInformation( this->GetInput() ); m_OutputTensorImage->SetBufferedRegion( this->GetInput()->GetBufferedRegion() ); m_OutputTensorImage->SetRequestedRegion( this->GetInput()->GetRequestedRegion() ); m_OutputTensorImage->Allocate(); //define an iterator for the input and output images typedef itk::ImageRegionConstIterator< InputDWIImageType > DWIImageIteratorType; typedef itk::ImageRegionIterator< OutputTensorImageType > TensorImageIteratorType; DWIImageIteratorType inputDWIit( this->GetInput(), m_OutputTensorImage->GetRequestedRegion() ); TensorImageIteratorType outputtensorit ( m_OutputTensorImage, m_OutputTensorImage->GetRequestedRegion() ); unsigned int N = this->m_TransformedGradients->Size(); TensorModelParamType tensormodelparams( 0.0 ); vnl_diag_matrix< double > W(N,0); for(inputDWIit.GoToBegin(), outputtensorit.GoToBegin(); !outputtensorit.IsAtEnd(); ++inputDWIit, ++outputtensorit){ CalculateTensorModelParameters( inputDWIit.Get(), W, tensormodelparams); OutputTensorImageType::PixelType& D = outputtensorit.Value(); //set the tensor model parameters into a Diffusion tensor D(0,0) = tensormodelparams[1]; D(0,1) = tensormodelparams[4]; D(0,2) = tensormodelparams[5]; D(1,0) = tensormodelparams[4]; D(1,1) = tensormodelparams[2]; D(1,2) = tensormodelparams[6]; D(2,0) = tensormodelparams[5]; D(2,1) = tensormodelparams[6]; D(2,2) = tensormodelparams[3]; //std::cout< #include #include // Qmitk #include "QmitkStochasticFiberTrackingView.h" #include "QmitkStdMultiWidget.h" // Qt #include // MITK #include #include // VTK #include #include #include #include #include #include const std::string QmitkStochasticFiberTrackingView::VIEW_ID = "org.mitk.views.stochasticfibertracking"; const std::string id_DataManager = "org.mitk.views.datamanager"; using namespace berry; QmitkStochasticFiberTrackingView::QmitkStochasticFiberTrackingView() : QmitkFunctionality() , m_Controls( 0 ) , m_MultiWidget( NULL ) - , m_DiffusionImage( NULL ) , m_SeedRoi( NULL ) + , m_DiffusionImage( NULL ) { } // Destructor QmitkStochasticFiberTrackingView::~QmitkStochasticFiberTrackingView() { } void QmitkStochasticFiberTrackingView::CreateQtPartControl( QWidget *parent ) { if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkStochasticFiberTrackingViewControls; m_Controls->setupUi( parent ); connect( m_Controls->commandLinkButton, SIGNAL(clicked()), this, SLOT(DoFiberTracking()) ); connect( m_Controls->m_SeedsPerVoxelSlider, SIGNAL(valueChanged(int)), this, SLOT(OnSeedsPerVoxelChanged(int)) ); connect( m_Controls->m_MaxCacheSizeSlider, SIGNAL(valueChanged(int)), this, SLOT(OnMaxCacheSizeChanged(int)) ); connect( m_Controls->m_MaxTractLengthSlider, SIGNAL(valueChanged(int)), this, SLOT(OnMaxTractLengthChanged(int)) ); } } void QmitkStochasticFiberTrackingView::OnSeedsPerVoxelChanged(int value) { m_Controls->m_SeedsPerVoxelLabel->setText(QString("Seeds per Voxel: ")+QString::number(value)); } void QmitkStochasticFiberTrackingView::OnMaxTractLengthChanged(int value) { m_Controls->m_MaxTractLengthLabel->setText(QString("Max. Tract Length: ")+QString::number(value)); } void QmitkStochasticFiberTrackingView::OnMaxCacheSizeChanged(int value) { m_Controls->m_MaxCacheSizeLabel->setText(QString("Max. Cache Size: ")+QString::number(value)+"GB"); } void QmitkStochasticFiberTrackingView::StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget) { m_MultiWidget = &stdMultiWidget; } void QmitkStochasticFiberTrackingView::StdMultiWidgetNotAvailable() { m_MultiWidget = NULL; } void QmitkStochasticFiberTrackingView::OnSelectionChanged( std::vector nodes ) { m_DiffusionImageNode = NULL; m_DiffusionImage = NULL; m_SeedRoi = NULL; m_Controls->m_DiffusionImageLabel->setText("mandatory"); m_Controls->m_RoiImageLabel->setText("mandatory"); for( std::vector::iterator it = nodes.begin(); it != nodes.end(); ++it ) { mitk::DataNode::Pointer node = *it; if( node.IsNotNull() && dynamic_cast(node->GetData()) ) { if( dynamic_cast*>(node->GetData()) ) { m_DiffusionImageNode = node; m_DiffusionImage = dynamic_cast*>(node->GetData()); m_Controls->m_DiffusionImageLabel->setText(node->GetName().c_str()); } else { bool isBinary = false; node->GetPropertyValue("binary", isBinary); if (isBinary) { m_SeedRoi = dynamic_cast(node->GetData()); m_Controls->m_RoiImageLabel->setText(node->GetName().c_str()); } } } } if(m_DiffusionImage.IsNotNull() && m_SeedRoi.IsNotNull()) { m_Controls->m_InputData->setTitle("Input Data"); m_Controls->commandLinkButton->setEnabled(true); } else { m_Controls->m_InputData->setTitle("Please Select Input Data"); m_Controls->commandLinkButton->setEnabled(false); } } void QmitkStochasticFiberTrackingView::DoFiberTracking() { typedef itk::VectorImage< short int, 3 > DWIVectorImageType; typedef itk::Image< float, 3 > FloatImageType; typedef itk::Image< unsigned int, 3 > CImageType; typedef itk::StochasticTractographyFilter< DWIVectorImageType, FloatImageType, CImageType > TrackingFilterType; typedef itk::DTITubeSpatialObject<3> DTITubeType; typedef itk::DTITubeSpatialObjectPoint<3> DTITubePointType; typedef itk::SceneSpatialObject<3> SceneSpatialObjectType; /* get Gradients/Direction of dwi */ itk::VectorContainer< unsigned int, vnl_vector_fixed >::Pointer Pdir = m_DiffusionImage->GetDirections(); /* bValueContainer, Container includes b-values according to corresponding gradient-direction*/ TrackingFilterType::bValueContainerType::Pointer vecCont = TrackingFilterType::bValueContainerType::New(); /* for each gradient set b-Value; for 0-gradient set b-value eq. 0 */ for ( int i=0; i<(int)Pdir->size(); ++i) { vnl_vector_fixed valsGrad = Pdir->at(i); if (valsGrad.get(0) == 0 && valsGrad.get(1) == 0 && valsGrad.get(2) == 0) { //set 0-Gradient to bValue 0 vecCont->InsertElement(i,0); }else{ vecCont->InsertElement(i,m_DiffusionImage->GetB_Value()); } } /* define measurement frame (identity-matrix 3x3) */ TrackingFilterType::MeasurementFrameType measurement_frame = m_DiffusionImage->GetMeasurementFrame(); /* generate white matterImage (dummy?)*/ FloatImageType::Pointer wmImage = FloatImageType::New(); wmImage->SetSpacing( m_DiffusionImage->GetVectorImage()->GetSpacing() ); wmImage->SetOrigin( m_DiffusionImage->GetVectorImage()->GetOrigin() ); wmImage->SetDirection( m_DiffusionImage->GetVectorImage()->GetDirection() ); wmImage->SetLargestPossibleRegion( m_DiffusionImage->GetVectorImage()->GetLargestPossibleRegion() ); wmImage->SetBufferedRegion( wmImage->GetLargestPossibleRegion() ); wmImage->SetRequestedRegion( wmImage->GetLargestPossibleRegion() ); wmImage->Allocate(); itk::ImageRegionIterator ot(wmImage, wmImage->GetLargestPossibleRegion() ); while (!ot.IsAtEnd()) { ot.Set(1); ++ot; } /* init TractographyFilter */ TrackingFilterType::Pointer trackingFilter = TrackingFilterType::New(); trackingFilter->SetInput(m_DiffusionImage->GetVectorImage().GetPointer()); trackingFilter->SetbValues(vecCont); trackingFilter->SetGradients(Pdir); trackingFilter->SetMeasurementFrame(measurement_frame); trackingFilter->SetWhiteMatterProbabilityImage(wmImage); trackingFilter->SetTotalTracts(m_Controls->m_SeedsPerVoxelSlider->value()); trackingFilter->SetMaxLikelihoodCacheSize(m_Controls->m_MaxCacheSizeSlider->value()*1000); trackingFilter->SetMaxTractLength(m_Controls->m_MaxTractLengthSlider->value()); //itk::Image< char, 3 > mitk::ImageToItk< itk::Image< unsigned char, 3 > >::Pointer binaryImageToItk1 = mitk::ImageToItk< itk::Image< unsigned char, 3 > >::New(); binaryImageToItk1->SetInput( m_SeedRoi ); binaryImageToItk1->Update(); vtkSmartPointer vPoints = vtkSmartPointer::New(); vtkSmartPointer vCellArray = vtkSmartPointer::New(); itk::ImageRegionConstIterator< BinaryImageType > it(binaryImageToItk1->GetOutput(), binaryImageToItk1->GetOutput()->GetRequestedRegion()); - it.Begin(); + it.GoToBegin(); mitk::Geometry3D* geom = m_DiffusionImage->GetGeometry(); while(!it.IsAtEnd()) { itk::ImageConstIterator::PixelType tmpPxValue = it.Get(); if(tmpPxValue != 0){ mitk::Point3D point; itk::ImageRegionConstIterator< BinaryImageType >::IndexType seedIdx = it.GetIndex(); trackingFilter->SetSeedIndex(seedIdx); trackingFilter->Update(); + /* get results from Filter */ /* write each single tract into member container */ TrackingFilterType::TractContainerType::Pointer container_tmp = trackingFilter->GetOutputTractContainer(); TrackingFilterType::TractContainerType::Iterator elIt = container_tmp->Begin(); TrackingFilterType::TractContainerType::Iterator end = container_tmp->End(); bool addTract = true; while( elIt != end ){ TrackingFilterType::TractContainerType::Element tract = elIt.Value(); TrackingFilterType::TractContainerType::Element::ObjectType::VertexListType::ConstPointer vertexlist = tract->GetVertexList(); vtkSmartPointer vPolyLine = vtkSmartPointer::New(); for( int j=0; j<(int)vertexlist->Size(); j++) { TrackingFilterType::TractContainerType::Element::ObjectType::VertexListType::Element vertex = vertexlist->GetElement(j); mitk::Point3D index; index[0] = (float)vertex[0]; index[1] = (float)vertex[1]; index[2] = (float)vertex[2]; if (geom->IsIndexInside(index)) { geom->IndexToWorld(index, point); vtkIdType id = vPoints->InsertNextPoint(point.GetDataPointer()); vPolyLine->GetPointIds()->InsertNextId(id); } else { addTract = false; break; } } if (addTract) vCellArray->InsertNextCell(vPolyLine); ++elIt; } } ++it; } vtkSmartPointer fiberPolyData = vtkSmartPointer::New(); fiberPolyData->SetPoints(vPoints); fiberPolyData->SetLines(vCellArray); mitk::FiberBundleX::Pointer fib = mitk::FiberBundleX::New(fiberPolyData); mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); fbNode->SetData(fib); QString name("FiberBundle_"); name += m_DiffusionImageNode->GetName().c_str(); name += "_Probabilistic"; fbNode->SetName(name.toStdString()); fbNode->SetVisibility(true); GetDataStorage()->Add(fbNode, m_DiffusionImageNode); }