diff --git a/Modules/DiffusionImaging/Tractography/itkStreamlineTrackingFilter.cpp b/Modules/DiffusionImaging/Tractography/itkStreamlineTrackingFilter.cpp index 5b13461545..18a89b9b8c 100644 --- a/Modules/DiffusionImaging/Tractography/itkStreamlineTrackingFilter.cpp +++ b/Modules/DiffusionImaging/Tractography/itkStreamlineTrackingFilter.cpp @@ -1,369 +1,422 @@ /*=================================================================== 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 __itkStreamlineTrackingFilter_txx #define __itkStreamlineTrackingFilter_txx #include #include #include #include "itkStreamlineTrackingFilter.h" #include #include #include #define _USE_MATH_DEFINES #include namespace itk { //#define QBALL_RECON_PI M_PI template< class TTensorPixelType, class TPDPixelType> StreamlineTrackingFilter< TTensorPixelType, TPDPixelType> ::StreamlineTrackingFilter(): m_FaThreshold(0.2), m_StepSize(1), m_MaxLength(10000), m_SeedsPerVoxel(1), m_AngularThreshold(0.7), m_F(1.0), m_G(0.0) { // 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 TTensorPixelType, class TPDPixelType> double StreamlineTrackingFilter< TTensorPixelType, TPDPixelType> ::RoundToNearest(double num) { return (num > 0.0) ? floor(num + 0.5) : ceil(num - 0.5); } template< class TTensorPixelType, class TPDPixelType> void StreamlineTrackingFilter< TTensorPixelType, TPDPixelType> ::BeforeThreadedGenerateData() { m_FiberPolyData = FiberPolyDataType::New(); m_Points = vtkPoints::New(); m_Cells = vtkCellArray::New(); typename InputImageType::Pointer inputImage = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) ); m_ImageSize.resize(3); m_ImageSize[0] = inputImage->GetLargestPossibleRegion().GetSize()[0]; m_ImageSize[1] = inputImage->GetLargestPossibleRegion().GetSize()[1]; m_ImageSize[2] = inputImage->GetLargestPossibleRegion().GetSize()[2]; m_ImageSpacing.resize(3); m_ImageSpacing[0] = inputImage->GetSpacing()[0]; m_ImageSpacing[1] = inputImage->GetSpacing()[1]; m_ImageSpacing[2] = inputImage->GetSpacing()[2]; if (m_StepSize<0.005) { float minSpacing; if(m_ImageSpacing[0]::New(); for (int i=0; iGetNumberOfThreads(); i++) { FiberPolyDataType poly = FiberPolyDataType::New(); m_PolyDataContainer->InsertElement(i, poly); } if (m_MaskImage.IsNull()) { - itk::Vector spacing = inputImage->GetSpacing(); - itk::Point origin = inputImage->GetOrigin(); - itk::Matrix direction = inputImage->GetDirection(); - ImageRegion<3> imageRegion = inputImage->GetLargestPossibleRegion(); - - // initialize crossings image + // initialize mask image m_MaskImage = ItkUcharImgType::New(); - m_MaskImage->SetSpacing( spacing ); - m_MaskImage->SetOrigin( origin ); - m_MaskImage->SetDirection( direction ); - m_MaskImage->SetRegions( imageRegion ); + m_MaskImage->SetSpacing( inputImage->GetSpacing() ); + m_MaskImage->SetOrigin( inputImage->GetOrigin() ); + m_MaskImage->SetDirection( inputImage->GetDirection() ); + m_MaskImage->SetRegions( inputImage->GetLargestPossibleRegion() ); m_MaskImage->Allocate(); m_MaskImage->FillBuffer(1); } + + m_FaImage = ItkFloatImgType::New(); + m_FaImage->SetSpacing( inputImage->GetSpacing() ); + m_FaImage->SetOrigin( inputImage->GetOrigin() ); + m_FaImage->SetDirection( inputImage->GetDirection() ); + m_FaImage->SetRegions( inputImage->GetLargestPossibleRegion() ); + m_FaImage->Allocate(); + m_FaImage->FillBuffer(0.0); + + m_PdImage = ItkPDImgType::New(); + m_PdImage->SetSpacing( inputImage->GetSpacing() ); + m_PdImage->SetOrigin( inputImage->GetOrigin() ); + m_PdImage->SetDirection( inputImage->GetDirection() ); + m_PdImage->SetRegions( inputImage->GetLargestPossibleRegion() ); + m_PdImage->Allocate(); + + m_EmaxImage = ItkFloatImgType::New(); + m_EmaxImage->SetSpacing( inputImage->GetSpacing() ); + m_EmaxImage->SetOrigin( inputImage->GetOrigin() ); + m_EmaxImage->SetDirection( inputImage->GetDirection() ); + m_EmaxImage->SetRegions( inputImage->GetLargestPossibleRegion() ); + m_EmaxImage->Allocate(); + m_EmaxImage->FillBuffer(1.0); + + typedef itk::DiffusionTensor3D TensorType; + typename TensorType::EigenValuesArrayType eigenvalues; + typename TensorType::EigenVectorsMatrixType eigenvectors; + for (int x=0; xGetPixel(index); + if(tensor.GetTrace()!=0 && tensor.GetFractionalAnisotropy()>m_FaThreshold) + { + vnl_vector_fixed dir; + tensor.ComputeEigenAnalysis(eigenvalues, eigenvectors); + dir[0] = eigenvectors(2, 0); + dir[1] = eigenvectors(2, 1); + dir[2] = eigenvectors(2, 2); + dir.normalize(); + m_PdImage->SetPixel(index, dir); + m_FaImage->SetPixel(index, tensor.GetFractionalAnisotropy()); + m_EmaxImage->SetPixel(index, 2/eigenvalues[2]); + } + } + std::cout << "StreamlineTrackingFilter: Angular threshold: " << m_AngularThreshold << std::endl; std::cout << "StreamlineTrackingFilter: FA threshold: " << m_FaThreshold << std::endl; std::cout << "StreamlineTrackingFilter: stepsize: " << m_StepSize << " mm" << std::endl; + std::cout << "StreamlineTrackingFilter: f: " << m_F << std::endl; + std::cout << "StreamlineTrackingFilter: g: " << m_G << std::endl; std::cout << "StreamlineTrackingFilter: starting streamline tracking" << std::endl; } template< class TTensorPixelType, class TPDPixelType> void StreamlineTrackingFilter< TTensorPixelType, TPDPixelType> ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId) { FiberPolyDataType poly = m_PolyDataContainer->GetElement(threadId); vtkSmartPointer Points = vtkPoints::New(); vtkSmartPointer Cells = vtkCellArray::New(); typedef itk::DiffusionTensor3D TensorType; typedef ImageRegionConstIterator< InputImageType > InputIteratorType; typedef ImageRegionConstIterator< ItkUcharImgType > MaskIteratorType; typedef typename InputImageType::PixelType InputTensorType; typename InputImageType::Pointer inputImage = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) ); InputIteratorType it(inputImage, outputRegionForThread ); MaskIteratorType mit(m_MaskImage, outputRegionForThread ); it.GoToBegin(); mit.GoToBegin(); while( !it.IsAtEnd() ) { if (mit.Value()==0) { ++mit; ++it; continue; } - typename TensorType::EigenValuesArrayType eigenvalues; - typename TensorType::EigenVectorsMatrixType eigenvectors; - for (int s=0; s container = vtkSmartPointer::New(); std::vector< vtkIdType > pointISs; typename InputImageType::IndexType index = it.GetIndex(); + typename InputImageType::IndexType indexOld; indexOld[0] = -1; indexOld[1] = -1; indexOld[2] = -1; itk::ContinuousIndex pos; itk::ContinuousIndex start; if (m_SeedsPerVoxel>1) { pos[0] = index[0]+(double)(rand()%99-49)/100; pos[1] = index[1]+(double)(rand()%99-49)/100; pos[2] = index[2]+(double)(rand()%99-49)/100; } else { pos[0] = index[0]; pos[1] = index[1]; pos[2] = index[2]; } start = pos; int step = 0; vnl_vector_fixed dirOld; dirOld.fill(0.0); // do forward tracking while (step < m_MaxLength) { ++step; index[0] = RoundToNearest(pos[0]); index[1] = RoundToNearest(pos[1]); index[2] = RoundToNearest(pos[2]); if (!inputImage->GetLargestPossibleRegion().IsInside(index)) break; typename InputImageType::PixelType tensor = inputImage->GetPixel(index); - if(tensor.GetTrace()!=0 && tensor.GetFractionalAnisotropy()>m_FaThreshold) + if(m_FaImage->GetPixel(index)>m_FaThreshold) { - tensor.ComputeEigenAnalysis(eigenvalues, eigenvectors); vnl_vector_fixed dir; - dir[0] = eigenvectors(2, 0); - dir[1] = eigenvectors(2, 1); - dir[2] = eigenvectors(2, 2); - dir.normalize(); - - if (!dirOld.is_zero()) + if (indexOld!=index) { - dir[0] = m_F*dir[0] + (1-m_F)*( (1-m_G)*dirOld[0] + m_G*(tensor[0]*dirOld[0] + tensor[1]*dirOld[1] + tensor[2]*dirOld[2])); - dir[1] = m_F*dir[1] + (1-m_F)*( (1-m_G)*dirOld[1] + m_G*(tensor[1]*dirOld[0] + tensor[3]*dirOld[1] + tensor[4]*dirOld[2])); - dir[2] = m_F*dir[2] + (1-m_F)*( (1-m_G)*dirOld[2] + m_G*(tensor[2]*dirOld[0] + tensor[4]*dirOld[1] + tensor[5]*dirOld[2])); - dir.normalize(); - - float angle = dot_product(dirOld, dir); - if (angle<0) - dir *= -1; - angle = fabs(dot_product(dirOld, dir)); - if (angleGetPixel(index); + + if (!dirOld.is_zero()) + { + float scale = m_EmaxImage->GetPixel(index); + dir[0] = m_F*dir[0] + (1-m_F)*( (1-m_G)*dirOld[0] + scale*m_G*(tensor[0]*dirOld[0] + tensor[1]*dirOld[1] + tensor[2]*dirOld[2])); + dir[1] = m_F*dir[1] + (1-m_F)*( (1-m_G)*dirOld[1] + scale*m_G*(tensor[1]*dirOld[0] + tensor[3]*dirOld[1] + tensor[4]*dirOld[2])); + dir[2] = m_F*dir[2] + (1-m_F)*( (1-m_G)*dirOld[2] + scale*m_G*(tensor[2]*dirOld[0] + tensor[4]*dirOld[1] + tensor[5]*dirOld[2])); + dir.normalize(); + + float angle = dot_product(dirOld, dir); + if (angle<0) + dir *= -1; + angle = dot_product(dirOld, dir); + if (angle worldPos; inputImage->TransformContinuousIndexToPhysicalPoint( pos, worldPos ); vtkIdType id = Points->InsertNextPoint(worldPos.GetDataPointer()); pointISs.push_back(id); counter++; pos[0] += dir[0]/m_ImageSpacing[0]; pos[1] += dir[1]/m_ImageSpacing[1]; pos[2] += dir[2]/m_ImageSpacing[2]; } } // insert reverse IDs while (!pointISs.empty()) { container->GetPointIds()->InsertNextId(pointISs.back()); pointISs.pop_back(); } // do backward tracking index = it.GetIndex(); + indexOld[0] = -1; indexOld[1] = -1; indexOld[2] = -1; pos = start; dirOld.fill(0.0); while (step < m_MaxLength) { ++step; index[0] = RoundToNearest(pos[0]); index[1] = RoundToNearest(pos[1]); index[2] = RoundToNearest(pos[2]); if (index[0] < 0 || index[0]>=m_ImageSize[0]) break; if (index[1] < 0 || index[1]>=m_ImageSize[1]) break; if (index[2] < 0 || index[2]>=m_ImageSize[2]) break; typename InputImageType::PixelType tensor = inputImage->GetPixel(index); - if(tensor.GetTrace()!=0 && tensor.GetFractionalAnisotropy()>m_FaThreshold) + if(m_FaImage->GetPixel(index)>m_FaThreshold) { - tensor.ComputeEigenAnalysis(eigenvalues, eigenvectors); - vnl_vector_fixed dir; - dir[0] = eigenvectors(2, 0); - dir[1] = eigenvectors(2, 1); - dir[2] = eigenvectors(2, 2); - dir.normalize(); - dir *= -1; // reverse direction - - if (!dirOld.is_zero()) + if (indexOld!=index) { - float angle = dot_product(dirOld, dir); - if (angle<0) - dir *= -1; - angle = fabs(dot_product(dirOld, dir)); - if (angle<0.7) - break; + dir = m_PdImage->GetPixel(index); + dir *= -1; // reverse direction + + if (!dirOld.is_zero()) + { + float scale = m_EmaxImage->GetPixel(index); + dir[0] = m_F*dir[0] + (1-m_F)*( (1-m_G)*dirOld[0] + scale*m_G*(tensor[0]*dirOld[0] + tensor[1]*dirOld[1] + tensor[2]*dirOld[2])); + dir[1] = m_F*dir[1] + (1-m_F)*( (1-m_G)*dirOld[1] + scale*m_G*(tensor[1]*dirOld[0] + tensor[3]*dirOld[1] + tensor[4]*dirOld[2])); + dir[2] = m_F*dir[2] + (1-m_F)*( (1-m_G)*dirOld[2] + scale*m_G*(tensor[2]*dirOld[0] + tensor[4]*dirOld[1] + tensor[5]*dirOld[2])); + dir.normalize(); + + float angle = dot_product(dirOld, dir); + if (angle<0) + dir *= -1; + angle = dot_product(dirOld, dir); + if (angle worldPos; inputImage->TransformContinuousIndexToPhysicalPoint( pos, worldPos ); vtkIdType id = Points->InsertNextPoint(worldPos.GetDataPointer()); container->GetPointIds()->InsertNextId(id); counter++; pos[0] += dir[0]/m_ImageSpacing[0]; pos[1] += dir[1]/m_ImageSpacing[1]; pos[2] += dir[2]/m_ImageSpacing[2]; } } if (counter>0) Cells->InsertNextCell(container); } ++mit; ++it; } poly->SetPoints(Points); poly->SetLines(Cells); std::cout << "Thread " << threadId << " finished tracking" << std::endl; } template< class TTensorPixelType, class TPDPixelType> vtkSmartPointer< vtkPolyData > StreamlineTrackingFilter< TTensorPixelType, TPDPixelType> ::AddPolyData(FiberPolyDataType poly1, FiberPolyDataType poly2) { vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); vtkSmartPointer vNewLines = poly1->GetLines(); vtkSmartPointer vNewPoints = poly1->GetPoints(); vtkSmartPointer vLines = poly2->GetLines(); vLines->InitTraversal(); for( int i=0; iGetNumberOfCells(); i++ ) { vtkIdType numPoints(0); vtkIdType* points(NULL); vLines->GetNextCell ( numPoints, points ); vtkSmartPointer container = vtkSmartPointer::New(); for( int j=0; jInsertNextPoint(poly2->GetPoint(points[j])); container->GetPointIds()->InsertNextId(id); } vNewLines->InsertNextCell(container); } // initialize polydata vNewPolyData->SetPoints(vNewPoints); vNewPolyData->SetLines(vNewLines); return vNewPolyData; } template< class TTensorPixelType, class TPDPixelType> void StreamlineTrackingFilter< TTensorPixelType, TPDPixelType> ::AfterThreadedGenerateData() { MITK_INFO << "Generating polydata "; m_FiberPolyData = m_PolyDataContainer->GetElement(0); for (int i=1; iGetNumberOfThreads(); i++) { m_FiberPolyData = AddPolyData(m_FiberPolyData, m_PolyDataContainer->GetElement(i)); } MITK_INFO << "done"; } template< class TTensorPixelType, class TPDPixelType> void StreamlineTrackingFilter< TTensorPixelType, TPDPixelType> ::PrintSelf(std::ostream& os, Indent indent) const { } } #endif // __itkDiffusionQballPrincipleDirectionsImageFilter_txx diff --git a/Modules/DiffusionImaging/Tractography/itkStreamlineTrackingFilter.h b/Modules/DiffusionImaging/Tractography/itkStreamlineTrackingFilter.h index dcdd0643a9..617fc1c3fc 100644 --- a/Modules/DiffusionImaging/Tractography/itkStreamlineTrackingFilter.h +++ b/Modules/DiffusionImaging/Tractography/itkStreamlineTrackingFilter.h @@ -1,116 +1,123 @@ /*=================================================================== 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. ===================================================================*/ /*=================================================================== This file is based heavily on a corresponding ITK filter. ===================================================================*/ #ifndef __itkStreamlineTrackingFilter_h_ #define __itkStreamlineTrackingFilter_h_ #include "MitkDiffusionImagingExports.h" #include #include #include #include #include #include #include #include #include namespace itk{ /** \class StreamlineTrackingFilter */ template< class TTensorPixelType, class TPDPixelType=double> class StreamlineTrackingFilter : public ImageToImageFilter< Image< DiffusionTensor3D, 3 >, Image< Vector< TPDPixelType, 3 >, 3 > > { public: typedef StreamlineTrackingFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ImageToImageFilter< Image< DiffusionTensor3D, 3 >, Image< Vector< TPDPixelType, 3 >, 3 > > Superclass; /** Method for creation through the object factory. */ itkNewMacro(Self) /** Runtime information support. */ itkTypeMacro(StreamlineTrackingFilter, ImageToImageFilter) - typedef TTensorPixelType TensorComponentType; - typedef TPDPixelType DirectionPixelType; - typedef typename Superclass::InputImageType InputImageType; - typedef typename Superclass::OutputImageType OutputImageType; - typedef typename Superclass::OutputImageRegionType OutputImageRegionType; - typedef itk::Image ItkUcharImgType; + typedef TTensorPixelType TensorComponentType; + typedef TPDPixelType DirectionPixelType; + typedef typename Superclass::InputImageType InputImageType; + typedef typename Superclass::OutputImageType OutputImageType; + typedef typename Superclass::OutputImageRegionType OutputImageRegionType; + typedef itk::Image ItkUcharImgType; + typedef itk::Image ItkFloatImgType; + typedef itk::Image< vnl_vector_fixed, 3> ItkPDImgType; typedef vtkSmartPointer< vtkPolyData > FiberPolyDataType; itkGetMacro( FiberPolyData, FiberPolyDataType ) itkSetMacro( MaskImage, ItkUcharImgType::Pointer) itkSetMacro( SeedsPerVoxel, int) itkSetMacro( FaThreshold, float) itkSetMacro( AngularThreshold, float) itkSetMacro( StepSize, float) itkSetMacro( F, float ) itkSetMacro( G, float ) protected: StreamlineTrackingFilter(); ~StreamlineTrackingFilter() {} void PrintSelf(std::ostream& os, Indent indent) const; double RoundToNearest(double num); void BeforeThreadedGenerateData(); void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, int threadId); void AfterThreadedGenerateData(); FiberPolyDataType AddPolyData(FiberPolyDataType poly1, FiberPolyDataType poly2); FiberPolyDataType m_FiberPolyData; vtkSmartPointer m_Points; vtkSmartPointer m_Cells; + + ItkFloatImgType::Pointer m_EmaxImage; + ItkFloatImgType::Pointer m_FaImage; + ItkPDImgType::Pointer m_PdImage; + float m_FaThreshold; float m_AngularThreshold; float m_StepSize; int m_MaxLength; int m_SeedsPerVoxel; float m_F; float m_G; std::vector< int > m_ImageSize; std::vector< float > m_ImageSpacing; ItkUcharImgType::Pointer m_MaskImage; itk::VectorContainer< int, FiberPolyDataType >::Pointer m_PolyDataContainer; private: }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkStreamlineTrackingFilter.cpp" #endif #endif //__itkStreamlineTrackingFilter_h_ diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkStreamlineTrackingViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkStreamlineTrackingViewControls.ui index 6e9462c529..fbdc9112e5 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkStreamlineTrackingViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkStreamlineTrackingViewControls.ui @@ -1,358 +1,358 @@ QmitkStreamlineTrackingViewControls 0 0 480 553 0 0 QmitkTemplate 3 3 0 Number of tracts started in each voxel of the seed ROI. Mori et al. Annals Neurology 1999 false Start Tracking Qt::Vertical QSizePolicy::Expanding 20 220 Parameters Minimum tract length in mm. Angular Threshold: 45° Minimum tract length in mm. Step Size: auto Minimum tract length in mm. f: 1 Number of tracts started in each voxel of the seed ROI. 1 10 Qt::Horizontal Fractional Anisotropy Threshold 0 - 180 + 90 45 Qt::Horizontal Minimum tract length in mm. 0 500 40 Qt::Horizontal Number of tracts started in each voxel of the seed ROI. Seeds per Voxel: 1 Qt::Horizontal QSizePolicy::Fixed 200 0 Fractional Anisotropy Threshold 0 100 20 Qt::Horizontal Weighting factor between first eigenvector (f=1 equals FACT tracking) and input vector dependent direction (f=0). 0 100 100 Qt::Horizontal Minimum tract length in mm. FA Threshold: 0.2 Minimum tract length in mm. Min. Tract Length: 40mm - Stepsize in mm (auto = 0.5*minimal spacing) + Stepsize in mm (auto = 0.1*minimal spacing) 0 100 0 Qt::Horizontal Minimum tract length in mm. - g: 1 + g: 0 Weighting factor between input vector (g=0) and tensor deflection (g=1 equals TEND tracking) 0 100 0 Qt::Horizontal Data Tensor Image: - Seed ROI Image: - Number of tracts started in each voxel of the seed ROI. Lazar et al. Human Brain Mapping 2003 Number of tracts started in each voxel of the seed ROI. Weinstein et al. Proceedings of IEEE Visualization 1999 commandLinkButton