diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStochasticTractographyFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStochasticTractographyFilter.h index c8770d5a56..d9f5378022 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStochasticTractographyFilter.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStochasticTractographyFilter.h @@ -1,305 +1,315 @@ /*=================================================================== 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) itkTypeMacro( StochasticTractographyFilter, 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) /** the maximum length of Tract **/ itkSetMacro( MaxTractLength, unsigned int ) itkGetMacro( MaxTractLength, unsigned int ) /** Set/Get bvalues **/ itkSetConstObjectMacro( bValues, bValueContainerType ) itkGetConstObjectMacro( bValues, bValueContainerType ) /** Set/Get of gradient directions **/ 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) //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 ); + void SetPrimaryInput(DataObject *input) + { + typename InputDWIImageType::Pointer dwiimagePtr; + try + { + dwiimagePtr = dynamic_cast( input ); + } + catch(itk::ExceptionObject &e) + { + MITK_INFO << "wrong image type: " << e.what(); + throw; + } + Superclass::SetPrimaryInput( input ); //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 ) /** Set/Get the list of directions to sample **/ itkSetConstObjectMacro( SampleDirections, TractOrientationContainerType ) itkGetConstObjectMacro( SampleDirections, TractOrientationContainerType ) /** Set/Get the Measurement Frame **/ 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 ) /** Get the Tracts that are generated **/ itkGetObjectMacro( OutputTractContainer, TractContainerType ) /** Get TensorImage **/ 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/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkStochasticFiberTrackingView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkStochasticFiberTrackingView.cpp index 1ec19523bc..65e522a80e 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkStochasticFiberTrackingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkStochasticFiberTrackingView.cpp @@ -1,289 +1,289 @@ /*=================================================================== 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 "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_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->SetPrimaryInput(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.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); }