diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingDataHandler.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingDataHandler.h index b70da654b0..48fd7cc884 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingDataHandler.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingDataHandler.h @@ -1,340 +1,340 @@ /*=================================================================== 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 _TrackingDataHandler #define _TrackingDataHandler #include #include #include #include #include #include #include #include namespace mitk { /** * \brief Abstract class for tracking handler. A tracking handler deals with determining the next progression direction of a streamline fiber. There are different handlers for tensor images, peak images, ... */ class MITKFIBERTRACKING_EXPORT TrackingDataHandler { public: enum MODE { DETERMINISTIC, PROBABILISTIC }; TrackingDataHandler(); virtual ~TrackingDataHandler(){} typedef boost::mt19937 BoostRngType; typedef itk::Image ItkUcharImgType; typedef itk::Image ItkShortImgType; typedef itk::Image ItkFloatImgType; typedef itk::Image ItkDoubleImgType; typedef vnl_vector_fixed< float, 3 > TrackingDirectionType; - virtual TrackingDirectionType ProposeDirection(itk::Point& pos, std::deque< TrackingDirectionType >& olddirs, itk::Index<3>& oldIndex) = 0; ///< predicts next progression direction at the given position + virtual TrackingDirectionType ProposeDirection(const itk::Point& pos, std::deque< TrackingDirectionType >& olddirs, itk::Index<3>& oldIndex) = 0; ///< predicts next progression direction at the given position virtual void InitForTracking() = 0; virtual itk::Vector GetSpacing() = 0; virtual itk::Point GetOrigin() = 0; virtual itk::Matrix GetDirection() = 0; virtual itk::ImageRegion<3> GetLargestPossibleRegion() = 0; virtual void SetMode(MODE m) = 0; MODE GetMode(){ return m_Mode; } void SetAngularThreshold( float a ){ m_AngularThreshold = a; } void SetInterpolate( bool interpolate ){ m_Interpolate = interpolate; } void SetFlipX( bool f ){ m_FlipX = f; } void SetFlipY( bool f ){ m_FlipY = f; } void SetFlipZ( bool f ){ m_FlipZ = f; } protected: float m_AngularThreshold; bool m_Interpolate; bool m_FlipX; bool m_FlipY; bool m_FlipZ; MODE m_Mode; BoostRngType m_Rng; template< class TPixelType > TPixelType GetImageValue(itk::Point itkP, itk::Image* image, vnl_vector_fixed& interpWeights){ // transform physical point to index coordinates itk::Index<3> idx; itk::ContinuousIndex< float, 3> cIdx; image->TransformPhysicalPointToIndex(itkP, idx); image->TransformPhysicalPointToContinuousIndex(itkP, cIdx); TPixelType pix = 0.0; if ( image->GetLargestPossibleRegion().IsInside(idx) ) { pix = image->GetPixel(idx); } else return pix; float frac_x = cIdx[0] - idx[0]; float frac_y = cIdx[1] - idx[1]; float frac_z = cIdx[2] - idx[2]; if (frac_x<0) { idx[0] -= 1; frac_x += 1; } if (frac_y<0) { idx[1] -= 1; frac_y += 1; } if (frac_z<0) { idx[2] -= 1; frac_z += 1; } frac_x = 1-frac_x; frac_y = 1-frac_y; frac_z = 1-frac_z; // int coordinates inside image? if (idx[0] >= 0 && idx[0] < image->GetLargestPossibleRegion().GetSize(0)-1 && idx[1] >= 0 && idx[1] < image->GetLargestPossibleRegion().GetSize(1)-1 && idx[2] >= 0 && idx[2] < image->GetLargestPossibleRegion().GetSize(2)-1) { // trilinear interpolation interpWeights[0] = ( frac_x)*( frac_y)*( frac_z); interpWeights[1] = (1-frac_x)*( frac_y)*( frac_z); interpWeights[2] = ( frac_x)*(1-frac_y)*( frac_z); interpWeights[3] = ( frac_x)*( frac_y)*(1-frac_z); interpWeights[4] = (1-frac_x)*(1-frac_y)*( frac_z); interpWeights[5] = ( frac_x)*(1-frac_y)*(1-frac_z); interpWeights[6] = (1-frac_x)*( frac_y)*(1-frac_z); interpWeights[7] = (1-frac_x)*(1-frac_y)*(1-frac_z); pix = image->GetPixel(idx) * interpWeights[0]; typename itk::Image::IndexType tmpIdx = idx; tmpIdx[0]++; pix += image->GetPixel(tmpIdx) * interpWeights[1]; tmpIdx = idx; tmpIdx[1]++; pix += image->GetPixel(tmpIdx) * interpWeights[2]; tmpIdx = idx; tmpIdx[2]++; pix += image->GetPixel(tmpIdx) * interpWeights[3]; tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++; pix += image->GetPixel(tmpIdx) * interpWeights[4]; tmpIdx = idx; tmpIdx[1]++; tmpIdx[2]++; pix += image->GetPixel(tmpIdx) * interpWeights[5]; tmpIdx = idx; tmpIdx[2]++; tmpIdx[0]++; pix += image->GetPixel(tmpIdx) * interpWeights[6]; tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++; tmpIdx[2]++; pix += image->GetPixel(tmpIdx) * interpWeights[7]; } if (pix!=pix) mitkThrow() << "nan values in image!"; return pix; } template< class TPixelType > TPixelType GetImageValue(itk::Point itkP, itk::Image* image, bool interpolate){ // transform physical point to index coordinates itk::Index<3> idx; itk::ContinuousIndex< float, 3> cIdx; image->TransformPhysicalPointToIndex(itkP, idx); image->TransformPhysicalPointToContinuousIndex(itkP, cIdx); TPixelType pix = 0.0; if ( image->GetLargestPossibleRegion().IsInside(idx) ) { pix = image->GetPixel(idx); if (!interpolate) return pix; } else return pix; float frac_x = cIdx[0] - idx[0]; float frac_y = cIdx[1] - idx[1]; float frac_z = cIdx[2] - idx[2]; if (frac_x<0) { idx[0] -= 1; frac_x += 1; } if (frac_y<0) { idx[1] -= 1; frac_y += 1; } if (frac_z<0) { idx[2] -= 1; frac_z += 1; } frac_x = 1-frac_x; frac_y = 1-frac_y; frac_z = 1-frac_z; // int coordinates inside image? if (idx[0] >= 0 && idx[0] < static_cast(image->GetLargestPossibleRegion().GetSize(0) - 1) && idx[1] >= 0 && idx[1] < static_cast(image->GetLargestPossibleRegion().GetSize(1) - 1) && idx[2] >= 0 && idx[2] < static_cast(image->GetLargestPossibleRegion().GetSize(2) - 1)) { // trilinear interpolation vnl_vector_fixed interpWeights; interpWeights[0] = ( frac_x)*( frac_y)*( frac_z); interpWeights[1] = (1-frac_x)*( frac_y)*( frac_z); interpWeights[2] = ( frac_x)*(1-frac_y)*( frac_z); interpWeights[3] = ( frac_x)*( frac_y)*(1-frac_z); interpWeights[4] = (1-frac_x)*(1-frac_y)*( frac_z); interpWeights[5] = ( frac_x)*(1-frac_y)*(1-frac_z); interpWeights[6] = (1-frac_x)*( frac_y)*(1-frac_z); interpWeights[7] = (1-frac_x)*(1-frac_y)*(1-frac_z); pix = image->GetPixel(idx) * interpWeights[0]; typename itk::Image::IndexType tmpIdx = idx; tmpIdx[0]++; pix += image->GetPixel(tmpIdx) * interpWeights[1]; tmpIdx = idx; tmpIdx[1]++; pix += image->GetPixel(tmpIdx) * interpWeights[2]; tmpIdx = idx; tmpIdx[2]++; pix += image->GetPixel(tmpIdx) * interpWeights[3]; tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++; pix += image->GetPixel(tmpIdx) * interpWeights[4]; tmpIdx = idx; tmpIdx[1]++; tmpIdx[2]++; pix += image->GetPixel(tmpIdx) * interpWeights[5]; tmpIdx = idx; tmpIdx[2]++; tmpIdx[0]++; pix += image->GetPixel(tmpIdx) * interpWeights[6]; tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++; tmpIdx[2]++; pix += image->GetPixel(tmpIdx) * interpWeights[7]; } if (pix!=pix) mitkThrow() << "nan values in image!"; return pix; } template< class TPixelType, int components > itk::Vector< TPixelType, components > GetImageValue(itk::Point itkP, itk::Image, 3>* image, bool interpolate){ // transform physical point to index coordinates itk::Index<3> idx; itk::ContinuousIndex< float, 3> cIdx; image->TransformPhysicalPointToIndex(itkP, idx); image->TransformPhysicalPointToContinuousIndex(itkP, cIdx); itk::Vector< TPixelType, components > pix = 0.0; if ( image->GetLargestPossibleRegion().IsInside(idx) ) { pix = image->GetPixel(idx); if (!interpolate) return pix; } else return pix; float frac_x = cIdx[0] - idx[0]; float frac_y = cIdx[1] - idx[1]; float frac_z = cIdx[2] - idx[2]; if (frac_x<0) { idx[0] -= 1; frac_x += 1; } if (frac_y<0) { idx[1] -= 1; frac_y += 1; } if (frac_z<0) { idx[2] -= 1; frac_z += 1; } frac_x = 1-frac_x; frac_y = 1-frac_y; frac_z = 1-frac_z; // int coordinates inside image? if (idx[0] >= 0 && idx[0] < static_cast(image->GetLargestPossibleRegion().GetSize(0) - 1) && idx[1] >= 0 && idx[1] < static_cast(image->GetLargestPossibleRegion().GetSize(1) - 1) && idx[2] >= 0 && idx[2] < static_cast(image->GetLargestPossibleRegion().GetSize(2) - 1)) { // trilinear interpolation vnl_vector_fixed interpWeights; interpWeights[0] = ( frac_x)*( frac_y)*( frac_z); interpWeights[1] = (1-frac_x)*( frac_y)*( frac_z); interpWeights[2] = ( frac_x)*(1-frac_y)*( frac_z); interpWeights[3] = ( frac_x)*( frac_y)*(1-frac_z); interpWeights[4] = (1-frac_x)*(1-frac_y)*( frac_z); interpWeights[5] = ( frac_x)*(1-frac_y)*(1-frac_z); interpWeights[6] = (1-frac_x)*( frac_y)*(1-frac_z); interpWeights[7] = (1-frac_x)*(1-frac_y)*(1-frac_z); pix = image->GetPixel(idx) * interpWeights[0]; typename itk::Image, 3>::IndexType tmpIdx = idx; tmpIdx[0]++; pix += image->GetPixel(tmpIdx) * interpWeights[1]; tmpIdx = idx; tmpIdx[1]++; pix += image->GetPixel(tmpIdx) * interpWeights[2]; tmpIdx = idx; tmpIdx[2]++; pix += image->GetPixel(tmpIdx) * interpWeights[3]; tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++; pix += image->GetPixel(tmpIdx) * interpWeights[4]; tmpIdx = idx; tmpIdx[1]++; tmpIdx[2]++; pix += image->GetPixel(tmpIdx) * interpWeights[5]; tmpIdx = idx; tmpIdx[2]++; tmpIdx[0]++; pix += image->GetPixel(tmpIdx) * interpWeights[6]; tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++; tmpIdx[2]++; pix += image->GetPixel(tmpIdx) * interpWeights[7]; } if (pix!=pix) mitkThrow() << "nan values in image!"; return pix; } }; } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.cpp index b406bec820..1480937425 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.cpp @@ -1,313 +1,313 @@ /*=================================================================== 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 "mitkTrackingHandlerOdf.h" #include #include #include namespace mitk { TrackingHandlerOdf::TrackingHandlerOdf() : m_GfaThreshold(0.1) , m_OdfPower(4) , m_SecondOrder(false) , m_MinMaxNormalize(true) { } TrackingHandlerOdf::~TrackingHandlerOdf() { } void TrackingHandlerOdf::InitForTracking() { MITK_INFO << "Initializing ODF tracker."; m_OdfHemisphereIndices.clear(); m_OdfReducedIndices.clear(); itk::OrientationDistributionFunction< float, QBALL_ODFSIZE > odf; vnl_vector_fixed ref; ref.fill(0); ref[0]=1; for (int i=0; i0) m_OdfHemisphereIndices.push_back(i); m_OdfFloatDirs.set_size(m_OdfHemisphereIndices.size(), 3); for (unsigned int i=0; i GfaFilterType; GfaFilterType::Pointer gfaFilter = GfaFilterType::New(); gfaFilter->SetInput(m_OdfImage); gfaFilter->SetComputationMethod(GfaFilterType::GFA_STANDARD); gfaFilter->Update(); m_GfaImage = gfaFilter->GetOutput(); } m_WorkingOdfImage = ItkOdfImageType::New(); m_WorkingOdfImage->SetSpacing( m_OdfImage->GetSpacing() ); m_WorkingOdfImage->SetOrigin( m_OdfImage->GetOrigin() ); m_WorkingOdfImage->SetDirection( m_OdfImage->GetDirection() ); m_WorkingOdfImage->SetRegions( m_OdfImage->GetLargestPossibleRegion() ); m_WorkingOdfImage->Allocate(); MITK_INFO << "Rescaling ODFs."; typedef itk::ImageRegionIterator< ItkOdfImageType > OdfIteratorType; OdfIteratorType it(m_OdfImage, m_OdfImage->GetLargestPossibleRegion() ); OdfIteratorType wit(m_WorkingOdfImage, m_WorkingOdfImage->GetLargestPossibleRegion() ); it.GoToBegin(); wit.GoToBegin(); while( !it.IsAtEnd() ) { ItkOdfImageType::PixelType odf_values = it.Get(); ItkOdfImageType::PixelType wodf_values = wit.Get(); for (int i=0; i=max) max = wodf_values[i]; else if (wodf_values[i]0.0) { wodf_values -= min; wodf_values /= max; } } for (int i=0; i TrackingHandlerOdf::GetSecondOrderProbabilities(itk::Point& itkP, vnl_vector< float >& angles, vnl_vector< float >& probs) +vnl_vector< float > TrackingHandlerOdf::GetSecondOrderProbabilities(const itk::Point& itkP, vnl_vector< float >& angles, vnl_vector< float >& probs) { vnl_vector< float > out_probs; out_probs.set_size(m_OdfHemisphereIndices.size()); out_probs.fill(0.0); float out_probs_sum = 0; int c = 0; for (unsigned int j=0; j=m_AngularThreshold) continue; vnl_vector_fixed d = m_OdfFloatDirs.get_row(j); itk::Point pos; pos[0] = itkP[0] + d[0]; pos[1] = itkP[1] + d[1]; pos[2] = itkP[2] + d[2]; ItkOdfImageType::PixelType odf_values = GetImageValue(pos, m_WorkingOdfImage, m_Interpolate); vnl_vector< float > new_angles = m_OdfFloatDirs*d; float probs_sum = 0; vnl_vector< float > new_probs; new_probs.set_size(m_OdfHemisphereIndices.size()); for (unsigned int i=0; i=m_AngularThreshold) { new_probs[i] = odf_values[m_OdfHemisphereIndices[i]]; probs_sum += new_probs[i]; } else new_probs[i] = 0; } if (probs_sum>0.0001) new_probs /= probs_sum; for (unsigned int i=0; i0.0001) out_probs /= out_probs_sum; return out_probs; } bool TrackingHandlerOdf::MinMaxNormalize() const { return m_MinMaxNormalize; } void TrackingHandlerOdf::SetMinMaxNormalize(bool MinMaxNormalize) { m_MinMaxNormalize = MinMaxNormalize; } void TrackingHandlerOdf::SetSecondOrder(bool SecondOrder) { m_SecondOrder = SecondOrder; } -vnl_vector_fixed TrackingHandlerOdf::ProposeDirection(itk::Point& pos, std::deque >& olddirs, itk::Index<3>& oldIndex) +vnl_vector_fixed TrackingHandlerOdf::ProposeDirection(const itk::Point& pos, std::deque >& olddirs, itk::Index<3>& oldIndex) { vnl_vector_fixed output_direction; output_direction.fill(0); itk::Index<3> idx; m_WorkingOdfImage->TransformPhysicalPointToIndex(pos, idx); if ( !m_WorkingOdfImage->GetLargestPossibleRegion().IsInside(idx) ) return output_direction; float gfa = GetImageValue(pos, m_GfaImage, m_Interpolate); if (gfa last_dir; if (!olddirs.empty()) last_dir = olddirs.back(); if (!m_Interpolate && oldIndex==idx) return last_dir; ItkOdfImageType::PixelType odf_values = GetImageValue(pos, m_WorkingOdfImage, m_Interpolate); float max = 0; int max_idx_v = -1; int max_idx_d = -1; int c = 0; for (int i : m_OdfHemisphereIndices) { if (odf_values[i]<0) odf_values[i] = 0; if (odf_values[i]>=max) { max = odf_values[i]; max_idx_v = i; max_idx_d = c; } c++; } if (max_idx_v>=0 && (olddirs.empty() || last_dir.magnitude()<=0.5) ) // no previous direction, so return principal diffusion direction { output_direction = m_OdfFloatDirs.get_row(max_idx_d); float odf_val = odf_values[max_idx_v]; return output_direction * odf_val; } else if (max_idx_v<0.0001) { return output_direction; } if (m_FlipX) last_dir[0] *= -1; if (m_FlipY) last_dir[1] *= -1; if (m_FlipZ) last_dir[2] *= -1; vnl_vector< float > angles = m_OdfFloatDirs*last_dir; vnl_vector< float > probs; probs.set_size(m_OdfHemisphereIndices.size()); float probs_sum = 0; for (unsigned int i=0; i d = m_OdfFloatDirs.get_row(i); if (angle<0) // make sure we don't walk backwards d *= -1; output_direction += odf_val*d; } else if (m_Mode==MODE::PROBABILISTIC) { probs[i] = odf_val; probs_sum += probs[i]; } } if (m_Mode==MODE::PROBABILISTIC && probs_sum>0.0001) { probs /= probs_sum; if (m_SecondOrder) probs = GetSecondOrderProbabilities(pos, angles, probs); boost::random::discrete_distribution dist(probs.begin(), probs.end()); int sampled_idx = 0; #pragma omp critical { boost::random::variate_generator> sampler(m_Rng, dist); sampled_idx = sampler(); } output_direction = m_OdfFloatDirs.get_row(sampled_idx); if (angles[sampled_idx]<0) // make sure we don't walk backwards output_direction *= -1; output_direction *= probs[sampled_idx]; } if (m_FlipX) output_direction[0] *= -1; if (m_FlipY) output_direction[1] *= -1; if (m_FlipZ) output_direction[2] *= -1; return output_direction; } int TrackingHandlerOdf::OdfPower() const { return m_OdfPower; } void TrackingHandlerOdf::SetOdfPower(int OdfPower) { m_OdfPower = OdfPower; } } diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.h index 78f943fd0d..d26baf8edb 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.h @@ -1,86 +1,86 @@ /*=================================================================== 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 _TrackingHandlerOdf #define _TrackingHandlerOdf #include "mitkTrackingDataHandler.h" #include #include #include namespace mitk { /** * \brief Enables streamline tracking on tensor images. Supports multi tensor tracking by adding multiple tensor images. */ class MITKFIBERTRACKING_EXPORT TrackingHandlerOdf : public TrackingDataHandler { public: TrackingHandlerOdf(); ~TrackingHandlerOdf(); typedef itk::DiffusionTensor3D TensorType; typedef itk::Image< itk::Vector< float, QBALL_ODFSIZE >, 3 > ItkOdfImageType; typedef itk::Image< vnl_vector_fixed, 3> ItkPDImgType; void InitForTracking(); ///< calls InputDataValidForTracking() and creates feature images - vnl_vector_fixed ProposeDirection(itk::Point& pos, std::deque< vnl_vector_fixed >& olddirs, itk::Index<3>& oldIndex); ///< predicts next progression direction at the given position + vnl_vector_fixed ProposeDirection(const itk::Point& pos, std::deque< vnl_vector_fixed >& olddirs, itk::Index<3>& oldIndex); ///< predicts next progression direction at the given position void SetGfaThreshold(float gfaThreshold){ m_GfaThreshold = gfaThreshold; } void SetOdfImage( ItkOdfImageType::Pointer img ){ m_OdfImage = img; } void SetGfaImage( ItkFloatImgType::Pointer img ){ m_GfaImage = img; } void SetMode( MODE m ){ m_Mode = m; } ItkUcharImgType::SpacingType GetSpacing(){ return m_OdfImage->GetSpacing(); } itk::Point GetOrigin(){ return m_OdfImage->GetOrigin(); } ItkUcharImgType::DirectionType GetDirection(){ return m_OdfImage->GetDirection(); } ItkUcharImgType::RegionType GetLargestPossibleRegion(){ return m_OdfImage->GetLargestPossibleRegion(); } int OdfPower() const; void SetOdfPower(int OdfPower); void SetSecondOrder(bool SecondOrder); bool MinMaxNormalize() const; void SetMinMaxNormalize(bool MinMaxNormalize); protected: - vnl_vector< float > GetSecondOrderProbabilities(itk::Point& itkP, vnl_vector< float >& angles, vnl_vector< float >& probs); + vnl_vector< float > GetSecondOrderProbabilities(const itk::Point& itkP, vnl_vector< float >& angles, vnl_vector< float >& probs); float m_GfaThreshold; ItkFloatImgType::Pointer m_GfaImage; ///< GFA image used to determine streamline termination. ItkOdfImageType::Pointer m_OdfImage; ///< Input odf image. ItkOdfImageType::Pointer m_WorkingOdfImage; ///< Modified odf image. std::vector< int > m_OdfHemisphereIndices; vnl_matrix< float > m_OdfFloatDirs; int m_OdfPower; bool m_SecondOrder; bool m_MinMaxNormalize; std::vector< int > m_OdfReducedIndices; }; } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.cpp index c71276bfc1..b8ac12b6c0 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.cpp @@ -1,271 +1,271 @@ /*=================================================================== 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 "mitkTrackingHandlerPeaks.h" namespace mitk { TrackingHandlerPeaks::TrackingHandlerPeaks() : m_PeakThreshold(0.1) , m_ApplyDirectionMatrix(false) { } TrackingHandlerPeaks::~TrackingHandlerPeaks() { } void TrackingHandlerPeaks::InitForTracking() { MITK_INFO << "Initializing peak tracker."; itk::Vector spacing4 = m_PeakImage->GetSpacing(); itk::Point origin4 = m_PeakImage->GetOrigin(); itk::Matrix direction4 = m_PeakImage->GetDirection(); itk::ImageRegion<4> imageRegion4 = m_PeakImage->GetLargestPossibleRegion(); spacing3[0] = spacing4[0]; spacing3[1] = spacing4[1]; spacing3[2] = spacing4[2]; origin3[0] = origin4[0]; origin3[1] = origin4[1]; origin3[2] = origin4[2]; for (int r=0; r<3; r++) for (int c=0; c<3; c++) { direction3[r][c] = direction4[r][c]; m_FloatImageRotation[r][c] = direction4[r][c]; } imageRegion3.SetSize(0, imageRegion4.GetSize()[0]); imageRegion3.SetSize(1, imageRegion4.GetSize()[1]); imageRegion3.SetSize(2, imageRegion4.GetSize()[2]); m_DummyImage = ItkUcharImgType::New(); m_DummyImage->SetSpacing( spacing3 ); m_DummyImage->SetOrigin( origin3 ); m_DummyImage->SetDirection( direction3 ); m_DummyImage->SetRegions( imageRegion3 ); m_DummyImage->Allocate(); m_DummyImage->FillBuffer(0.0); m_NumDirs = imageRegion4.GetSize(3)/3; } vnl_vector_fixed TrackingHandlerPeaks::GetMatchingDirection(itk::Index<3> idx3, vnl_vector_fixed& oldDir) { vnl_vector_fixed out_dir; out_dir.fill(0); float angle = 0; float mag = oldDir.magnitude(); if (magmitk::eps) { oldDir[0] = out_dir[0]; oldDir[1] = out_dir[1]; oldDir[2] = out_dir[2]; break; } } // if you didn't find a non-zero random direction, take first non-zero direction you find for (int i=0; imitk::eps) { oldDir[0] = out_dir[0]; oldDir[1] = out_dir[1]; oldDir[2] = out_dir[2]; break; } } } else { for (int i=0; i dir = GetDirection(idx3, i); mag = dir.magnitude(); if (mag>mitk::eps) dir.normalize(); float a = dot_product(dir, oldDir); if (fabs(a)>angle) { angle = fabs(a); if (a<0) out_dir = -dir; else out_dir = dir; out_dir *= mag; out_dir *= angle; // shrink contribution of direction if is less parallel to previous direction } } } return out_dir; } vnl_vector_fixed TrackingHandlerPeaks::GetDirection(itk::Index<3> idx3, int dirIdx) { vnl_vector_fixed dir; dir.fill(0.0); if ( !m_DummyImage->GetLargestPossibleRegion().IsInside(idx3) ) return dir; PeakImgType::IndexType idx4; idx4.SetElement(0,idx3[0]); idx4.SetElement(1,idx3[1]); idx4.SetElement(2,idx3[2]); for (int k=0; k<3; k++) { idx4.SetElement(3, dirIdx*3 + k); dir[k] = m_PeakImage->GetPixel(idx4); } if (m_FlipX) dir[0] *= -1; if (m_FlipY) dir[1] *= -1; if (m_FlipZ) dir[2] *= -1; if (m_ApplyDirectionMatrix) dir = m_FloatImageRotation*dir; return dir; } vnl_vector_fixed TrackingHandlerPeaks::GetDirection(itk::Point itkP, bool interpolate, vnl_vector_fixed oldDir){ // transform physical point to index coordinates itk::Index<3> idx3; itk::ContinuousIndex< float, 3> cIdx; m_DummyImage->TransformPhysicalPointToIndex(itkP, idx3); m_DummyImage->TransformPhysicalPointToContinuousIndex(itkP, cIdx); vnl_vector_fixed dir; dir.fill(0.0); if ( !m_DummyImage->GetLargestPossibleRegion().IsInside(idx3) ) return dir; if (interpolate) { float frac_x = cIdx[0] - idx3[0]; float frac_y = cIdx[1] - idx3[1]; float frac_z = cIdx[2] - idx3[2]; if (frac_x<0) { idx3[0] -= 1; frac_x += 1; } if (frac_y<0) { idx3[1] -= 1; frac_y += 1; } if (frac_z<0) { idx3[2] -= 1; frac_z += 1; } frac_x = 1-frac_x; frac_y = 1-frac_y; frac_z = 1-frac_z; // int coordinates inside image? if (idx3[0] >= 0 && idx3[0] < static_cast(m_DummyImage->GetLargestPossibleRegion().GetSize(0) - 1) && idx3[1] >= 0 && idx3[1] < static_cast(m_DummyImage->GetLargestPossibleRegion().GetSize(1) - 1) && idx3[2] >= 0 && idx3[2] < static_cast(m_DummyImage->GetLargestPossibleRegion().GetSize(2) - 1)) { // trilinear interpolation vnl_vector_fixed interpWeights; interpWeights[0] = ( frac_x)*( frac_y)*( frac_z); interpWeights[1] = (1-frac_x)*( frac_y)*( frac_z); interpWeights[2] = ( frac_x)*(1-frac_y)*( frac_z); interpWeights[3] = ( frac_x)*( frac_y)*(1-frac_z); interpWeights[4] = (1-frac_x)*(1-frac_y)*( frac_z); interpWeights[5] = ( frac_x)*(1-frac_y)*(1-frac_z); interpWeights[6] = (1-frac_x)*( frac_y)*(1-frac_z); interpWeights[7] = (1-frac_x)*(1-frac_y)*(1-frac_z); dir = GetMatchingDirection(idx3, oldDir) * interpWeights[0]; itk::Index<3> tmpIdx = idx3; tmpIdx[0]++; dir += GetMatchingDirection(tmpIdx, oldDir) * interpWeights[1]; tmpIdx = idx3; tmpIdx[1]++; dir += GetMatchingDirection(tmpIdx, oldDir) * interpWeights[2]; tmpIdx = idx3; tmpIdx[2]++; dir += GetMatchingDirection(tmpIdx, oldDir) * interpWeights[3]; tmpIdx = idx3; tmpIdx[0]++; tmpIdx[1]++; dir += GetMatchingDirection(tmpIdx, oldDir) * interpWeights[4]; tmpIdx = idx3; tmpIdx[1]++; tmpIdx[2]++; dir += GetMatchingDirection(tmpIdx, oldDir) * interpWeights[5]; tmpIdx = idx3; tmpIdx[2]++; tmpIdx[0]++; dir += GetMatchingDirection(tmpIdx, oldDir) * interpWeights[6]; tmpIdx = idx3; tmpIdx[0]++; tmpIdx[1]++; tmpIdx[2]++; dir += GetMatchingDirection(tmpIdx, oldDir) * interpWeights[7]; } } else dir = GetMatchingDirection(idx3, oldDir); return dir; } -vnl_vector_fixed TrackingHandlerPeaks::ProposeDirection(itk::Point& pos, std::deque >& olddirs, itk::Index<3>& oldIndex) +vnl_vector_fixed TrackingHandlerPeaks::ProposeDirection(const itk::Point& pos, std::deque >& olddirs, itk::Index<3>& oldIndex) { // CHECK: wann wird wo normalisiert vnl_vector_fixed output_direction; output_direction.fill(0); itk::Index<3> index; m_DummyImage->TransformPhysicalPointToIndex(pos, index); vnl_vector_fixed oldDir = olddirs.back(); float old_mag = oldDir.magnitude(); if (!m_Interpolate && oldIndex==index) return oldDir; output_direction = GetDirection(pos, m_Interpolate, oldDir); float mag = output_direction.magnitude(); if (mag>=m_PeakThreshold) { output_direction.normalize(); float a = 1; if (old_mag>0.5) a = dot_product(output_direction, oldDir); if (a>m_AngularThreshold) output_direction *= mag; else output_direction.fill(0); } else output_direction.fill(0); return output_direction; } } diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.h index 17207aa842..d5ee97f495 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.h @@ -1,83 +1,83 @@ /*=================================================================== 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 _TrackingHandlerPeaks #define _TrackingHandlerPeaks #include "mitkTrackingDataHandler.h" #include #include namespace mitk { /** * \brief Enables deterministic streamline tracking on MRtrix style peak images (4D float images) */ class MITKFIBERTRACKING_EXPORT TrackingHandlerPeaks : public TrackingDataHandler { public: TrackingHandlerPeaks(); ~TrackingHandlerPeaks(); typedef itk::Image< float, 4 > PeakImgType; ///< MRtrix peak image type void InitForTracking(); ///< calls InputDataValidForTracking() and creates feature images - vnl_vector_fixed ProposeDirection(itk::Point& pos, std::deque< vnl_vector_fixed >& olddirs, itk::Index<3>& oldIndex); ///< predicts next progression direction at the given position + vnl_vector_fixed ProposeDirection(const itk::Point& pos, std::deque< vnl_vector_fixed >& olddirs, itk::Index<3>& oldIndex); ///< predicts next progression direction at the given position void SetPeakThreshold(float thr){ m_PeakThreshold = thr; } void SetPeakImage( PeakImgType::Pointer image ){ m_PeakImage = image; } void SetApplyDirectionMatrix( bool applyDirectionMatrix ){ m_ApplyDirectionMatrix = applyDirectionMatrix; } itk::Vector GetSpacing(){ return spacing3; } itk::Point GetOrigin(){ return origin3; } itk::Matrix GetDirection(){ return direction3; } itk::ImageRegion<3> GetLargestPossibleRegion(){ return imageRegion3; } void SetMode( MODE m ) { if (m==MODE::DETERMINISTIC) m_Mode = m; else mitkThrow() << "Peak tracker is only implemented for deterministic mode."; } protected: vnl_vector_fixed GetDirection(itk::Point itkP, bool interpolate, vnl_vector_fixed oldDir); vnl_vector_fixed GetMatchingDirection(itk::Index<3> idx3, vnl_vector_fixed& oldDir); vnl_vector_fixed GetDirection(itk::Index<3> idx3, int dirIdx); PeakImgType::Pointer m_PeakImage; float m_PeakThreshold; int m_NumDirs; itk::Vector spacing3; itk::Point origin3; itk::Matrix direction3; itk::ImageRegion<3> imageRegion3; vnl_matrix_fixed m_FloatImageRotation; ItkUcharImgType::Pointer m_DummyImage; bool m_ApplyDirectionMatrix; }; } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.cpp index 8252cbcbfd..cc3b261fb1 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.cpp @@ -1,946 +1,946 @@ /*=================================================================== 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 _TrackingForestHandler_cpp #define _TrackingForestHandler_cpp #include "mitkTrackingHandlerRandomForest.h" #include #include namespace mitk { template< int ShOrder, int NumberOfSignalFeatures > TrackingHandlerRandomForest< ShOrder, NumberOfSignalFeatures >::TrackingHandlerRandomForest() : m_WmSampleDistance(-1) , m_NumTrees(30) , m_MaxTreeDepth(25) , m_SampleFraction(1.0) , m_NumberOfSamples(0) , m_GmSamplesPerVoxel(-1) , m_NumPreviousDirections(1) , m_ZeroDirWmFeatures(true) , m_BidirectionalFiberSampling(false) , m_MaxNumWmSamples(-1) { vnl_vector_fixed ref; ref.fill(0); ref[0]=1; itk::OrientationDistributionFunction< float, 200 > odf; m_DirectionContainer.clear(); for (unsigned int i = 0; i odf_dir; odf_dir[0] = odf.GetDirection(i)[0]; odf_dir[1] = odf.GetDirection(i)[1]; odf_dir[2] = odf.GetDirection(i)[2]; if (dot_product(ref, odf_dir)>0) // only used directions on one hemisphere m_DirectionContainer.push_back(odf_dir); // store indices for later mapping the classifier output to the actual direction } m_OdfFloatDirs.set_size(m_DirectionContainer.size(), 3); for (unsigned int i=0; i TrackingHandlerRandomForest< ShOrder, NumberOfSignalFeatures >::~TrackingHandlerRandomForest() { } template< int ShOrder, int NumberOfSignalFeatures > typename TrackingHandlerRandomForest< ShOrder, NumberOfSignalFeatures >::DwiFeatureImageType::PixelType TrackingHandlerRandomForest< ShOrder, NumberOfSignalFeatures >::GetDwiFeaturesAtPosition(itk::Point itkP, typename DwiFeatureImageType::Pointer image, bool interpolate) { // transform physical point to index coordinates itk::Index<3> idx; itk::ContinuousIndex< float, 3> cIdx; image->TransformPhysicalPointToIndex(itkP, idx); image->TransformPhysicalPointToContinuousIndex(itkP, cIdx); typename DwiFeatureImageType::PixelType pix; pix.Fill(0.0); if ( image->GetLargestPossibleRegion().IsInside(idx) ) { pix = image->GetPixel(idx); if (!interpolate) return pix; } else return pix; float frac_x = cIdx[0] - idx[0]; float frac_y = cIdx[1] - idx[1]; float frac_z = cIdx[2] - idx[2]; if (frac_x<0) { idx[0] -= 1; frac_x += 1; } if (frac_y<0) { idx[1] -= 1; frac_y += 1; } if (frac_z<0) { idx[2] -= 1; frac_z += 1; } frac_x = 1-frac_x; frac_y = 1-frac_y; frac_z = 1-frac_z; // int coordinates inside image? if (idx[0] >= 0 && idx[0] < static_cast(image->GetLargestPossibleRegion().GetSize(0) - 1) && idx[1] >= 0 && idx[1] < static_cast(image->GetLargestPossibleRegion().GetSize(1) - 1) && idx[2] >= 0 && idx[2] < static_cast(image->GetLargestPossibleRegion().GetSize(2) - 1)) { // trilinear interpolation vnl_vector_fixed interpWeights; interpWeights[0] = ( frac_x)*( frac_y)*( frac_z); interpWeights[1] = (1-frac_x)*( frac_y)*( frac_z); interpWeights[2] = ( frac_x)*(1-frac_y)*( frac_z); interpWeights[3] = ( frac_x)*( frac_y)*(1-frac_z); interpWeights[4] = (1-frac_x)*(1-frac_y)*( frac_z); interpWeights[5] = ( frac_x)*(1-frac_y)*(1-frac_z); interpWeights[6] = (1-frac_x)*( frac_y)*(1-frac_z); interpWeights[7] = (1-frac_x)*(1-frac_y)*(1-frac_z); pix = image->GetPixel(idx) * interpWeights[0]; typename DwiFeatureImageType::IndexType tmpIdx = idx; tmpIdx[0]++; pix += image->GetPixel(tmpIdx) * interpWeights[1]; tmpIdx = idx; tmpIdx[1]++; pix += image->GetPixel(tmpIdx) * interpWeights[2]; tmpIdx = idx; tmpIdx[2]++; pix += image->GetPixel(tmpIdx) * interpWeights[3]; tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++; pix += image->GetPixel(tmpIdx) * interpWeights[4]; tmpIdx = idx; tmpIdx[1]++; tmpIdx[2]++; pix += image->GetPixel(tmpIdx) * interpWeights[5]; tmpIdx = idx; tmpIdx[2]++; tmpIdx[0]++; pix += image->GetPixel(tmpIdx) * interpWeights[6]; tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++; tmpIdx[2]++; pix += image->GetPixel(tmpIdx) * interpWeights[7]; } return pix; } template< int ShOrder, int NumberOfSignalFeatures > void TrackingHandlerRandomForest< ShOrder, NumberOfSignalFeatures >::InputDataValidForTracking() { if (m_InputDwis.empty()) mitkThrow() << "No diffusion-weighted images set!"; if (!IsForestValid()) mitkThrow() << "No or invalid random forest detected!"; } template< int ShOrder, int NumberOfSignalFeatures> template typename std::enable_if< NumberOfSignalFeatures <=99, T >::type TrackingHandlerRandomForest< ShOrder, NumberOfSignalFeatures >::InitDwiImageFeatures(mitk::Image::Pointer mitk_dwi) { MITK_INFO << "Calculating spherical harmonics features"; typedef itk::AnalyticalDiffusionQballReconstructionImageFilter InterpolationFilterType; typename InterpolationFilterType::Pointer filter = InterpolationFilterType::New(); filter->SetGradientImage( mitk::DiffusionPropertyHelper::GetOriginalGradientContainer(mitk_dwi), mitk::DiffusionPropertyHelper::GetItkVectorImage(mitk_dwi) ); filter->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(mitk_dwi)); filter->SetLambda(0.006); filter->SetNormalizationMethod(InterpolationFilterType::QBAR_RAW_SIGNAL); filter->Update(); m_DwiFeatureImages.push_back(filter->GetCoefficientImage()); return true; } template< int ShOrder, int NumberOfSignalFeatures> template typename std::enable_if< NumberOfSignalFeatures >=100, T >::type TrackingHandlerRandomForest< ShOrder, NumberOfSignalFeatures >::InitDwiImageFeatures(mitk::Image::Pointer mitk_dwi) { MITK_INFO << "Interpolating raw dwi signal features"; typedef itk::AnalyticalDiffusionQballReconstructionImageFilter InterpolationFilterType; typename InterpolationFilterType::Pointer filter = InterpolationFilterType::New(); filter->SetGradientImage( mitk::DiffusionPropertyHelper::GetOriginalGradientContainer(mitk_dwi), mitk::DiffusionPropertyHelper::GetItkVectorImage(mitk_dwi) ); filter->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(mitk_dwi)); filter->SetLambda(0.006); filter->SetNormalizationMethod(InterpolationFilterType::QBAR_RAW_SIGNAL); filter->Update(); typename DwiFeatureImageType::Pointer dwiFeatureImage = DwiFeatureImageType::New(); dwiFeatureImage->SetSpacing(filter->GetOutput()->GetSpacing()); dwiFeatureImage->SetOrigin(filter->GetOutput()->GetOrigin()); dwiFeatureImage->SetDirection(filter->GetOutput()->GetDirection()); dwiFeatureImage->SetLargestPossibleRegion(filter->GetOutput()->GetLargestPossibleRegion()); dwiFeatureImage->SetBufferedRegion(filter->GetOutput()->GetLargestPossibleRegion()); dwiFeatureImage->SetRequestedRegion(filter->GetOutput()->GetLargestPossibleRegion()); dwiFeatureImage->Allocate(); // get signal values and store them in the feature image vnl_vector_fixed ref; ref.fill(0); ref[0]=1; itk::OrientationDistributionFunction< float, 2*NumberOfSignalFeatures > odf; itk::ImageRegionIterator< typename InterpolationFilterType::OutputImageType > it(filter->GetOutput(), filter->GetOutput()->GetLargestPossibleRegion()); while(!it.IsAtEnd()) { typename DwiFeatureImageType::PixelType pix; int f = 0; for (unsigned int i = 0; i0) // only used directions on one hemisphere { pix[f] = it.Get()[i]; f++; } } dwiFeatureImage->SetPixel(it.GetIndex(), pix); ++it; } m_DwiFeatureImages.push_back(dwiFeatureImage); return true; } template< int ShOrder, int NumberOfSignalFeatures > void TrackingHandlerRandomForest< ShOrder, NumberOfSignalFeatures >::InitForTracking() { MITK_INFO << "Initializing random forest tracker."; InputDataValidForTracking(); m_DwiFeatureImages.clear(); InitDwiImageFeatures<>(m_InputDwis.at(0)); } template< int ShOrder, int NumberOfSignalFeatures > -vnl_vector_fixed TrackingHandlerRandomForest< ShOrder, NumberOfSignalFeatures >::ProposeDirection(itk::Point& pos, std::deque >& olddirs, itk::Index<3>& oldIndex) +vnl_vector_fixed TrackingHandlerRandomForest< ShOrder, NumberOfSignalFeatures >::ProposeDirection(const itk::Point& pos, std::deque >& olddirs, itk::Index<3>& oldIndex) { vnl_vector_fixed output_direction; output_direction.fill(0); itk::Index<3> idx; m_DwiFeatureImages.at(0)->TransformPhysicalPointToIndex(pos, idx); bool check_last_dir = false; vnl_vector_fixed last_dir; if (!olddirs.empty()) { last_dir = olddirs.back(); if (last_dir.magnitude()>0.5) check_last_dir = true; } if (!m_Interpolate && oldIndex==idx) return last_dir; // store feature pixel values in a vigra data type vigra::MultiArray<2, float> featureData = vigra::MultiArray<2, float>( vigra::Shape2(1,m_Forest->feature_count()) ); typename DwiFeatureImageType::PixelType dwiFeaturePixel = GetDwiFeaturesAtPosition(pos, m_DwiFeatureImages.at(0), m_Interpolate); for (unsigned int f=0; f direction_matrix = m_DwiFeatureImages.at(0)->GetDirection().GetVnlMatrix(); vnl_matrix_fixed inverse_direction_matrix = m_DwiFeatureImages.at(0)->GetInverseDirection().GetVnlMatrix(); // append normalized previous direction(s) to feature vector int i = 0; vnl_vector_fixed ref; ref.fill(0); ref[0]=1; for (auto d : olddirs) { vnl_vector_fixed tempD; tempD[0] = d[0]; tempD[1] = d[1]; tempD[2] = d[2]; if (m_FlipX) tempD[0] *= -1; if (m_FlipY) tempD[1] *= -1; if (m_FlipZ) tempD[2] *= -1; tempD = inverse_direction_matrix * tempD; last_dir[0] = tempD[0]; last_dir[1] = tempD[1]; last_dir[2] = tempD[2]; int c = 0; for (int f=NumberOfSignalFeatures+3*i; f0) { int c = 0; for (auto img : m_AdditionalFeatureImages.at(0)) { float v = GetImageValue(pos, img, false); featureData(0,NumberOfSignalFeatures+m_NumPreviousDirections*3+c) = v; c++; } } // perform classification vigra::MultiArray<2, float> probs(vigra::Shape2(1, m_Forest->class_count())); m_Forest->predictProbabilities(featureData, probs); vnl_vector< float > angles = m_OdfFloatDirs*last_dir; vnl_vector< float > probs2; probs2.set_size(m_DirectionContainer.size()); probs2.fill(0.0); // used for probabilistic direction sampling float probs_sum = 0; float pNonFib = 0; // probability that we left the white matter float w = 0; // weight of the predicted direction for (int i=0; iclass_count(); i++) // for each class (number of possible directions + out-of-wm class) { if (probs(0,i)>0) // if probability of respective class is 0, do nothing { // get label of class (does not correspond to the loop variable i) unsigned int classLabel = 0; m_Forest->ext_param_.to_classlabel(i, classLabel); if (classLabel=m_AngularThreshold) probs2[classLabel] = probs(0,i); probs_sum += probs2[classLabel]; } else if (m_Mode==MODE::DETERMINISTIC) { vnl_vector_fixed d = m_DirectionContainer.at(classLabel); // get direction vector assiciated with the respective direction index if (check_last_dir) // do we have a previous streamline direction or did we just start? { // TODO: check if hard curvature threshold is necessary. // alternatively try square of dot product as weight. // TODO: check if additional weighting with dot product as directional prior is necessary. are there alternatives on the classification level? float dot = angles[classLabel]; // claculate angle between the candidate direction vector and the previous streamline direction if (fabs(dot)>=m_AngularThreshold) // is angle between the directions smaller than our hard threshold? { if (dot<0) // make sure we don't walk backwards d *= -1; float w_i = probs(0,i)*fabs(dot); output_direction += w_i*d; // weight contribution to output direction with its probability and the angular deviation from the previous direction w += w_i; // increase output weight of the final direction } } else { output_direction += probs(0,i)*d; w += probs(0,i); } } } else pNonFib += probs(0,i); // probability that we are not in the white matter anymore } } if (m_Mode==MODE::PROBABILISTIC && pNonFib<0.5) { probs2 /= probs_sum; boost::random::discrete_distribution dist(probs2.begin(), probs2.end()); int sampled_idx = 0; #pragma omp critical { boost::random::variate_generator> sampler(m_Rng, dist); sampled_idx = sampler(); } output_direction = m_DirectionContainer.at(sampled_idx); w = probs2[sampled_idx]; if (check_last_dir && angles[sampled_idx]<0) // make sure we don't walk backwards output_direction *= -1; } // if we did not find a suitable direction, make sure that we return (0,0,0) if (pNonFib>w && w>0) output_direction.fill(0.0); else { vnl_vector_fixed tempD; tempD[0] = output_direction[0]; tempD[1] = output_direction[1]; tempD[2] = output_direction[2]; tempD = direction_matrix * tempD; output_direction[0] = tempD[0]; output_direction[1] = tempD[1]; output_direction[2] = tempD[2]; if (m_FlipX) output_direction[0] *= -1; if (m_FlipY) output_direction[1] *= -1; if (m_FlipZ) output_direction[2] *= -1; } return output_direction * w; } template< int ShOrder, int NumberOfSignalFeatures > void TrackingHandlerRandomForest< ShOrder, NumberOfSignalFeatures >::StartTraining() { m_StartTime = std::chrono::system_clock::now(); InputDataValidForTraining(); InitForTraining(); CalculateTrainingSamples(); MITK_INFO << "Maximum tree depths: " << m_MaxTreeDepth; MITK_INFO << "Sample fraction per tree: " << m_SampleFraction; MITK_INFO << "Number of trees: " << m_NumTrees; DefaultSplitType splitter; splitter.UsePointBasedWeights(true); splitter.SetWeights(m_Weights); splitter.UseRandomSplit(false); splitter.SetPrecision(mitk::eps); splitter.SetMaximumTreeDepth(m_MaxTreeDepth); std::vector< std::shared_ptr< vigra::RandomForest > > trees; int count = 0; #pragma omp parallel for for (int i = 0; i < m_NumTrees; ++i) { std::shared_ptr< vigra::RandomForest > lrf = std::make_shared< vigra::RandomForest >(); lrf->set_options().use_stratification(vigra::RF_NONE); // How the data should be made equal lrf->set_options().sample_with_replacement(true); // if sampled with replacement or not lrf->set_options().samples_per_tree(m_SampleFraction); // Fraction of samples that are used to train a tree lrf->set_options().tree_count(1); // Number of trees that are calculated; lrf->set_options().min_split_node_size(5); // Minimum number of datapoints that must be in a node lrf->ext_param_.max_tree_depth = m_MaxTreeDepth; lrf->learn(m_FeatureData, m_LabelData,vigra::rf::visitors::VisitorBase(),splitter); #pragma omp critical { count++; MITK_INFO << "Tree " << count << " finished training."; trees.push_back(lrf); } } for (int i = 1; i < m_NumTrees; ++i) trees.at(0)->trees_.push_back(trees.at(i)->trees_[0]); m_Forest = trees.at(0); m_Forest->options_.tree_count_ = m_NumTrees; MITK_INFO << "Training finsihed"; m_EndTime = std::chrono::system_clock::now(); std::chrono::hours hh = std::chrono::duration_cast(m_EndTime - m_StartTime); std::chrono::minutes mm = std::chrono::duration_cast(m_EndTime - m_StartTime); mm %= 60; MITK_INFO << "Training took " << hh.count() << "h and " << mm.count() << "m"; } template< int ShOrder, int NumberOfSignalFeatures > void TrackingHandlerRandomForest< ShOrder, NumberOfSignalFeatures >::InputDataValidForTraining() { if (m_InputDwis.empty()) mitkThrow() << "No diffusion-weighted images set!"; if (m_Tractograms.empty()) mitkThrow() << "No tractograms set!"; if (m_InputDwis.size()!=m_Tractograms.size()) mitkThrow() << "Unequal number of diffusion-weighted images and tractograms detected!"; } template< int ShOrder, int NumberOfSignalFeatures > bool TrackingHandlerRandomForest< ShOrder, NumberOfSignalFeatures >::IsForestValid() { int additional_features = 0; if (m_AdditionalFeatureImages.size()>0) additional_features = m_AdditionalFeatureImages.at(0).size(); if (!m_Forest) MITK_ERROR << "Forest could not be read!"; else { if (m_Forest->tree_count()<=0) MITK_ERROR << "Forest contains no trees!"; if ( m_Forest->feature_count()!=(NumberOfSignalFeatures+3*m_NumPreviousDirections+additional_features) ) MITK_ERROR << "Wrong number of features in forest: got " << m_Forest->feature_count() << ", expected " << (NumberOfSignalFeatures+3*m_NumPreviousDirections+additional_features); } if(m_Forest && m_Forest->tree_count()>0 && m_Forest->feature_count()==(NumberOfSignalFeatures+3*m_NumPreviousDirections+additional_features)) return true; return false; } template< int ShOrder, int NumberOfSignalFeatures > void TrackingHandlerRandomForest< ShOrder, NumberOfSignalFeatures >::InitForTraining() { typedef itk::AnalyticalDiffusionQballReconstructionImageFilter InterpolationFilterType; MITK_INFO << "Spherical signal interpolation and sampling ..."; for (unsigned int i=0; i(m_InputDwis.at(i)); if (i>=m_AdditionalFeatureImages.size()) { m_AdditionalFeatureImages.push_back(std::vector< ItkFloatImgType::Pointer >()); } if (i>=m_FiberVolumeModImages.size()) { ItkFloatImgType::Pointer img = ItkFloatImgType::New(); img->SetSpacing( m_DwiFeatureImages.at(i)->GetSpacing() ); img->SetOrigin( m_DwiFeatureImages.at(i)->GetOrigin() ); img->SetDirection( m_DwiFeatureImages.at(i)->GetDirection() ); img->SetLargestPossibleRegion( m_DwiFeatureImages.at(i)->GetLargestPossibleRegion() ); img->SetBufferedRegion( m_DwiFeatureImages.at(i)->GetLargestPossibleRegion() ); img->SetRequestedRegion( m_DwiFeatureImages.at(i)->GetLargestPossibleRegion() ); img->Allocate(); img->FillBuffer(1); m_FiberVolumeModImages.push_back(img); } if (m_FiberVolumeModImages.at(i)==nullptr) { m_FiberVolumeModImages.at(i) = ItkFloatImgType::New(); m_FiberVolumeModImages.at(i)->SetSpacing( m_DwiFeatureImages.at(i)->GetSpacing() ); m_FiberVolumeModImages.at(i)->SetOrigin( m_DwiFeatureImages.at(i)->GetOrigin() ); m_FiberVolumeModImages.at(i)->SetDirection( m_DwiFeatureImages.at(i)->GetDirection() ); m_FiberVolumeModImages.at(i)->SetLargestPossibleRegion( m_DwiFeatureImages.at(i)->GetLargestPossibleRegion() ); m_FiberVolumeModImages.at(i)->SetBufferedRegion( m_DwiFeatureImages.at(i)->GetLargestPossibleRegion() ); m_FiberVolumeModImages.at(i)->SetRequestedRegion( m_DwiFeatureImages.at(i)->GetLargestPossibleRegion() ); m_FiberVolumeModImages.at(i)->Allocate(); m_FiberVolumeModImages.at(i)->FillBuffer(1); } if (i>=m_MaskImages.size()) { ItkUcharImgType::Pointer newMask = ItkUcharImgType::New(); newMask->SetSpacing( m_DwiFeatureImages.at(i)->GetSpacing() ); newMask->SetOrigin( m_DwiFeatureImages.at(i)->GetOrigin() ); newMask->SetDirection( m_DwiFeatureImages.at(i)->GetDirection() ); newMask->SetLargestPossibleRegion( m_DwiFeatureImages.at(i)->GetLargestPossibleRegion() ); newMask->SetBufferedRegion( m_DwiFeatureImages.at(i)->GetLargestPossibleRegion() ); newMask->SetRequestedRegion( m_DwiFeatureImages.at(i)->GetLargestPossibleRegion() ); newMask->Allocate(); newMask->FillBuffer(1); m_MaskImages.push_back(newMask); } if (m_MaskImages.at(i)==nullptr) { m_MaskImages.at(i) = ItkUcharImgType::New(); m_MaskImages.at(i)->SetSpacing( m_DwiFeatureImages.at(i)->GetSpacing() ); m_MaskImages.at(i)->SetOrigin( m_DwiFeatureImages.at(i)->GetOrigin() ); m_MaskImages.at(i)->SetDirection( m_DwiFeatureImages.at(i)->GetDirection() ); m_MaskImages.at(i)->SetLargestPossibleRegion( m_DwiFeatureImages.at(i)->GetLargestPossibleRegion() ); m_MaskImages.at(i)->SetBufferedRegion( m_DwiFeatureImages.at(i)->GetLargestPossibleRegion() ); m_MaskImages.at(i)->SetRequestedRegion( m_DwiFeatureImages.at(i)->GetLargestPossibleRegion() ); m_MaskImages.at(i)->Allocate(); m_MaskImages.at(i)->FillBuffer(1); } } MITK_INFO << "Resampling fibers and calculating number of samples ..."; m_NumberOfSamples = 0; m_SampleUsage.clear(); for (unsigned int t=0; t::Pointer env = itk::TractDensityImageFilter< ItkUcharImgType >::New(); env->SetFiberBundle(m_Tractograms.at(t)); env->SetInputImage(mask); env->SetBinaryOutput(true); env->SetUseImageGeometry(true); env->Update(); wmmask = env->GetOutput(); if (t>=m_WhiteMatterImages.size()) m_WhiteMatterImages.push_back(wmmask); else m_WhiteMatterImages.at(t) = wmmask; } // Calculate white-matter samples if (m_WmSampleDistance<0) { typename DwiFeatureImageType::Pointer image = m_DwiFeatureImages.at(t); float minSpacing = 1; if(image->GetSpacing()[0]GetSpacing()[1] && image->GetSpacing()[0]GetSpacing()[2]) minSpacing = image->GetSpacing()[0]; else if (image->GetSpacing()[1] < image->GetSpacing()[2]) minSpacing = image->GetSpacing()[1]; else minSpacing = image->GetSpacing()[2]; m_WmSampleDistance = minSpacing*0.5; } m_Tractograms.at(t)->ResampleLinear(m_WmSampleDistance); unsigned int wmSamples = m_Tractograms.at(t)->GetNumberOfPoints()-2*m_Tractograms.at(t)->GetNumFibers(); if (m_BidirectionalFiberSampling) wmSamples *= 2; if (m_ZeroDirWmFeatures) wmSamples *= (m_NumPreviousDirections+1); MITK_INFO << "White matter samples available: " << wmSamples; // upper limit for samples if (m_MaxNumWmSamples>0 && wmSamples>m_MaxNumWmSamples) { if ((float)m_MaxNumWmSamples/wmSamples > 0.8) { m_SampleUsage.push_back(std::vector(wmSamples, true)); m_NumberOfSamples += wmSamples; } else { m_SampleUsage.push_back(std::vector(wmSamples, false)); m_NumberOfSamples += m_MaxNumWmSamples; wmSamples = m_MaxNumWmSamples; MITK_INFO << "Limiting white matter samples to: " << m_MaxNumWmSamples; itk::Statistics::MersenneTwisterRandomVariateGenerator::Pointer randgen = itk::Statistics::MersenneTwisterRandomVariateGenerator::New(); randgen->SetSeed(); unsigned int c = 0; while (cGetIntegerVariate(m_MaxNumWmSamples-1); if (m_SampleUsage[t][idx]==false) { m_SampleUsage[t][idx]=true; c++; } } } } else { m_SampleUsage.push_back(std::vector(wmSamples, true)); m_NumberOfSamples += wmSamples; } // calculate gray-matter samples itk::ImageRegionConstIterator it(wmmask, wmmask->GetLargestPossibleRegion()); int OUTOFWM = 0; while(!it.IsAtEnd()) { if (it.Get()==0 && mask->GetPixel(it.GetIndex())>0) OUTOFWM++; ++it; } MITK_INFO << "Non-white matter voxels: " << OUTOFWM; if (m_GmSamplesPerVoxel>0) { m_GmSamples.push_back(m_GmSamplesPerVoxel); m_NumberOfSamples += m_GmSamplesPerVoxel*OUTOFWM; } else if (OUTOFWM>0) { int gm_per_voxel = 0.5+(float)wmSamples/(float)OUTOFWM; if (gm_per_voxel<=0) gm_per_voxel = 1; m_GmSamples.push_back(gm_per_voxel); m_NumberOfSamples += m_GmSamples.back()*OUTOFWM; MITK_INFO << "Non-white matter samples per voxel: " << m_GmSamples.back(); } else { m_GmSamples.push_back(0); } MITK_INFO << "Non-white matter samples: " << m_GmSamples.back()*OUTOFWM; } MITK_INFO << "Number of samples: " << m_NumberOfSamples; } template< int ShOrder, int NumberOfSignalFeatures > void TrackingHandlerRandomForest< ShOrder, NumberOfSignalFeatures >::CalculateTrainingSamples() { vnl_vector_fixed ref; ref.fill(0); ref[0]=1; m_FeatureData.reshape( vigra::Shape2(m_NumberOfSamples, NumberOfSignalFeatures+m_NumPreviousDirections*3+m_AdditionalFeatureImages.at(0).size()) ); m_LabelData.reshape( vigra::Shape2(m_NumberOfSamples,1) ); m_Weights.reshape( vigra::Shape2(m_NumberOfSamples,1) ); MITK_INFO << "Number of features: " << m_FeatureData.shape(1); itk::Statistics::MersenneTwisterRandomVariateGenerator::Pointer m_RandGen = itk::Statistics::MersenneTwisterRandomVariateGenerator::New(); m_RandGen->SetSeed(); MITK_INFO << "Creating training data ..."; unsigned int sampleCounter = 0; for (unsigned int t=0; t it(wmMask, wmMask->GetLargestPossibleRegion()); while(!it.IsAtEnd()) { if (it.Get()==0 && (mask.IsNull() || (mask.IsNotNull() && mask->GetPixel(it.GetIndex())>0))) { typename DwiFeatureImageType::PixelType pix = image->GetPixel(it.GetIndex()); // random directions for (unsigned int i=0; iGetIntegerVariate(m_NumPreviousDirections); // how many directions should be zero? for (unsigned int i=0; i probe; if (iGetVariate()*2-1; probe[1] = m_RandGen->GetVariate()*2-1; probe[2] = m_RandGen->GetVariate()*2-1; probe.normalize(); if (dot_product(ref, probe)<0) probe *= -1; } for (unsigned int f=NumberOfSignalFeatures+3*i; f itkP; image->TransformIndexToPhysicalPoint(it.GetIndex(), itkP); float v = GetImageValue(itkP, img, false); m_FeatureData(sampleCounter,NumberOfSignalFeatures+m_NumPreviousDirections*3+add_feat_c) = v; add_feat_c++; } // label and sample weight m_LabelData(sampleCounter,0) = m_DirectionContainer.size(); m_Weights(sampleCounter,0) = 1.0; sampleCounter++; } } ++it; } unsigned int num_gm_samples = sampleCounter; // white matter samples mitk::FiberBundle::Pointer fib = m_Tractograms.at(t); vtkSmartPointer< vtkPolyData > polyData = fib->GetFiberPolyData(); vnl_vector_fixed zero_dir; zero_dir.fill(0.0); for (int i=0; iGetNumFibers(); i++) { vtkCell* cell = polyData->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); float fiber_weight = fib->GetFiberWeight(i); for (int n = 0; n<=m_NumPreviousDirections; n++) { if (!m_ZeroDirWmFeatures) n = m_NumPreviousDirections; for (bool reverse : {false, true}) { for (int j=1; j itkP1, itkP2; int num_nonzero_dirs = m_NumPreviousDirections; if (!reverse) num_nonzero_dirs = std::min(n, j); else num_nonzero_dirs = std::min(n, numPoints-j-1); vnl_vector_fixed dir; // zero directions for (unsigned int k=0; kGetPoint(j-n_idx); itkP1[0] = p[0]; itkP1[1] = p[1]; itkP1[2] = p[2]; p = points->GetPoint(j-n_idx+1); itkP2[0] = p[0]; itkP2[1] = p[1]; itkP2[2] = p[2]; } else { p = points->GetPoint(j+n_idx); itkP1[0] = p[0]; itkP1[1] = p[1]; itkP1[2] = p[2]; p = points->GetPoint(j+n_idx-1); itkP2[0] = p[0]; itkP2[1] = p[1]; itkP2[2] = p[2]; } dir[0]=itkP1[0]-itkP2[0]; dir[1]=itkP1[1]-itkP2[1]; dir[2]=itkP1[2]-itkP2[2]; if (dir.magnitude()<0.0001) mitkThrow() << "streamline error!"; dir.normalize(); if (dir[0]!=dir[0] || dir[1]!=dir[1] || dir[2]!=dir[2]) mitkThrow() << "ERROR: NaN direction!"; if (dot_product(ref, dir)<0) dir *= -1; int c = 0; for (unsigned int f=NumberOfSignalFeatures+3*(k+m_NumPreviousDirections-num_nonzero_dirs); fGetPoint(j); itkP1[0] = p[0]; itkP1[1] = p[1]; itkP1[2] = p[2]; if (reverse) { p = points->GetPoint(j-1); itkP2[0] = p[0]; itkP2[1] = p[1]; itkP2[2] = p[2]; } else { p = points->GetPoint(j+1); itkP2[0] = p[0]; itkP2[1] = p[1]; itkP2[2] = p[2]; } dir[0]=itkP2[0]-itkP1[0]; dir[1]=itkP2[1]-itkP1[1]; dir[2]=itkP2[2]-itkP1[2]; if (dir.magnitude()<0.0001) mitkThrow() << "streamline error!"; dir.normalize(); if (dir[0]!=dir[0] || dir[1]!=dir[1] || dir[2]!=dir[2]) mitkThrow() << "ERROR: NaN direction!"; if (dot_product(ref, dir)<0) dir *= -1; // image features float volume_mod = GetImageValue(itkP1, fiber_folume, false); // diffusion signal features typename DwiFeatureImageType::PixelType pix = GetDwiFeaturesAtPosition(itkP1, image, m_Interpolate); for (unsigned int f=0; f(itkP1, img, false); add_feat_c++; m_FeatureData(sampleCounter,NumberOfSignalFeatures+2+add_feat_c) = v; } // set label values float angle = 0; float m = dir.magnitude(); if (m>0.0001) { int l = 0; for (auto d : m_DirectionContainer) { float a = fabs(dot_product(dir, d)); if (a>angle) { m_LabelData(sampleCounter,0) = l; m_Weights(sampleCounter,0) = fiber_weight*volume_mod; angle = a; } l++; } } sampleCounter++; } if (!m_BidirectionalFiberSampling) // don't sample fibers backward break; } } } } m_Tractograms.clear(); MITK_INFO << "done"; } template< int ShOrder, int NumberOfSignalFeatures > void TrackingHandlerRandomForest< ShOrder, NumberOfSignalFeatures >::SaveForest(std::string forestFile) { MITK_INFO << "Saving forest to " << forestFile; if (IsForestValid()) { vigra::rf_export_HDF5( *m_Forest, forestFile, "" ); MITK_INFO << "Forest saved successfully."; } else MITK_INFO << "Forest invalid! Could not be saved."; } template< int ShOrder, int NumberOfSignalFeatures > void TrackingHandlerRandomForest< ShOrder, NumberOfSignalFeatures >::LoadForest(std::string forestFile) { MITK_INFO << "Loading forest from " << forestFile; m_Forest = std::make_shared< vigra::RandomForest >(); vigra::rf_import_HDF5( *m_Forest, forestFile); } } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.h index 85e0dc5f7a..4a77723084 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.h @@ -1,160 +1,160 @@ /*=================================================================== 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 _TrackingForestHandler #define _TrackingForestHandler #include #include #include #include #include #include #include #undef DIFFERENCE #define VIGRA_STATIC_LIB #include #include #include #include #include #include #define _USE_MATH_DEFINES #include namespace mitk { /** * \brief Manages random forests for fiber tractography. The preparation of the features from the inputa data and the training process are handled here. The data preprocessing and actual prediction for the tracking process is also performed here. The tracking itself is performed in MLBSTrackingFilter. */ template< int ShOrder, int NumberOfSignalFeatures > class TrackingHandlerRandomForest : public TrackingDataHandler { public: TrackingHandlerRandomForest(); ~TrackingHandlerRandomForest(); typedef itk::Image< itk::Vector< float, NumberOfSignalFeatures > , 3 > DwiFeatureImageType; typedef mitk::ThresholdSplit >,int,vigra::ClassificationTag> DefaultSplitType; void SetDwis( std::vector< Image::Pointer > images ){ m_InputDwis = images; } void AddDwi( Image::Pointer img ){ m_InputDwis.push_back(img); } void SetTractograms( std::vector< FiberBundle::Pointer > tractograms ) { m_Tractograms.clear(); for (auto fib : tractograms) m_Tractograms.push_back(fib); } void SetMaskImages( std::vector< ItkUcharImgType::Pointer > images ){ m_MaskImages = images; } void SetWhiteMatterImages( std::vector< ItkUcharImgType::Pointer > images ){ m_WhiteMatterImages = images; } void SetFiberVolumeModImages( std::vector< ItkFloatImgType::Pointer > images ){ m_FiberVolumeModImages = images; } void SetAdditionalFeatureImages( std::vector< std::vector< ItkFloatImgType::Pointer > > images ){ m_AdditionalFeatureImages = images; } void StartTraining(); void SaveForest(std::string forestFile); void LoadForest(std::string forestFile); void SetMode( MODE m ){ m_Mode = m; } void SetMaxNumWmSamples(int num){ m_MaxNumWmSamples=num; } void SetNumPreviousDirections( int num ){ m_NumPreviousDirections=num; } void SetNumTrees(int num){ m_NumTrees = num; } void SetMaxTreeDepth(int depth){ m_MaxTreeDepth = depth; } void SetStepSize(float step){ m_WmSampleDistance = step; } void SetGrayMatterSamplesPerVoxel(int samples){ m_GmSamplesPerVoxel = samples; } void SetSampleFraction(float fraction){ m_SampleFraction = fraction; } void SetBidirectionalFiberSampling(bool val) { m_BidirectionalFiberSampling = val; } void SetZeroDirWmFeatures(bool val) { m_ZeroDirWmFeatures = val; } std::shared_ptr< vigra::RandomForest > GetForest(){ return m_Forest; } void InitForTracking(); ///< calls InputDataValidForTracking() and creates feature images - vnl_vector_fixed ProposeDirection(itk::Point& pos, std::deque< vnl_vector_fixed >& olddirs, itk::Index<3>& oldIndex); ///< predicts next progression direction at the given position + vnl_vector_fixed ProposeDirection(const itk::Point& pos, std::deque< vnl_vector_fixed >& olddirs, itk::Index<3>& oldIndex); ///< predicts next progression direction at the given position bool IsForestValid(); ///< true is forest is not null, has more than 0 trees and the correct number of features (NumberOfSignalFeatures + 3) ItkUcharImgType::SpacingType GetSpacing(){ return m_DwiFeatureImages.at(0)->GetSpacing(); } itk::Point GetOrigin(){ return m_DwiFeatureImages.at(0)->GetOrigin(); } ItkUcharImgType::DirectionType GetDirection(){ return m_DwiFeatureImages.at(0)->GetDirection(); } ItkUcharImgType::RegionType GetLargestPossibleRegion(){ return m_DwiFeatureImages.at(0)->GetLargestPossibleRegion(); } protected: void InputDataValidForTracking(); ///< check if raw data is set and tracking forest is valid template typename std::enable_if::type InitDwiImageFeatures(mitk::Image::Pointer mitk_dwi); template typename std::enable_if= 100, T>::type InitDwiImageFeatures(mitk::Image::Pointer mitk_dwi); void InputDataValidForTraining(); ///< Check if everything is tehere for training (raw datasets, fiber tracts) void InitForTraining(); ///< Generate masks if necessary, resample fibers, spherically interpolate raw DWIs void CalculateTrainingSamples(); ///< Calculate GM and WM features using the interpolated raw data, the WM masks and the fibers typename DwiFeatureImageType::PixelType GetDwiFeaturesAtPosition(itk::Point itkP, typename DwiFeatureImageType::Pointer image, bool interpolate); ///< get trilinearly interpolated raw image values at given world position std::vector< Image::Pointer > m_InputDwis; ///< original input DWI data std::shared_ptr< vigra::RandomForest > m_Forest; ///< random forest classifier std::chrono::time_point m_StartTime; std::chrono::time_point m_EndTime; std::vector< typename DwiFeatureImageType::Pointer > m_DwiFeatureImages; std::vector< std::vector< ItkFloatImgType::Pointer > > m_AdditionalFeatureImages; std::vector< ItkFloatImgType::Pointer > m_FiberVolumeModImages; ///< used to correct the fiber density std::vector< FiberBundle::Pointer > m_Tractograms; ///< training tractograms std::vector< ItkUcharImgType::Pointer > m_MaskImages; ///< binary mask images to constrain training to a certain area (e.g. brain mask) std::vector< ItkUcharImgType::Pointer > m_WhiteMatterImages; ///< defines white matter voxels. if not set, theses mask images are automatically generated from the input tractograms float m_WmSampleDistance; ///< deterines the number of white matter samples (distance of sampling points on each fiber). int m_NumTrees; ///< number of trees in random forest int m_MaxTreeDepth; ///< limits the tree depth float m_SampleFraction; ///< fraction of samples used to train each tree unsigned int m_NumberOfSamples; ///< stores overall number of samples used for training std::vector< unsigned int > m_GmSamples; ///< number of gray matter samples int m_GmSamplesPerVoxel; ///< number of gray matter samplees per voxel. if -1, then the number is automatically chosen to gain an overall number of GM samples close to the number of WM samples. vigra::MultiArray<2, float> m_FeatureData; ///< vigra container for training features unsigned int m_NumPreviousDirections; ///< How many "old" directions should be used as classification features? // only for tracking vigra::MultiArray<2, float> m_LabelData; ///< vigra container for training labels vigra::MultiArray<2, double> m_Weights; ///< vigra container for training sample weights std::vector< vnl_vector_fixed > m_DirectionContainer; vnl_matrix< float > m_OdfFloatDirs; bool m_BidirectionalFiberSampling; bool m_ZeroDirWmFeatures; int m_MaxNumWmSamples; std::vector< std::vector< bool > > m_SampleUsage; vnl_matrix_fixed m_ImageDirection; vnl_matrix_fixed m_ImageDirectionInverse; }; } #include "mitkTrackingHandlerRandomForest.cpp" #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.cpp index 255dba0ace..fef06bc2d2 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.cpp @@ -1,344 +1,344 @@ /*=================================================================== 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 "mitkTrackingHandlerTensor.h" namespace mitk { TrackingHandlerTensor::TrackingHandlerTensor() : m_InterpolateTensors(true) , m_FaThreshold(0.1) { } TrackingHandlerTensor::~TrackingHandlerTensor() { } void TrackingHandlerTensor::InitForTracking() { MITK_INFO << "Initializing tensor tracker."; m_NumberOfInputs = m_TensorImages.size(); for (int i=0; iSetSpacing( m_TensorImages.at(0)->GetSpacing() ); pdImage->SetOrigin( m_TensorImages.at(0)->GetOrigin() ); pdImage->SetDirection( m_TensorImages.at(0)->GetDirection() ); pdImage->SetRegions( m_TensorImages.at(0)->GetLargestPossibleRegion() ); pdImage->Allocate(); m_PdImage.push_back(pdImage); ItkDoubleImgType::Pointer emaxImage = ItkDoubleImgType::New(); emaxImage->SetSpacing( m_TensorImages.at(0)->GetSpacing() ); emaxImage->SetOrigin( m_TensorImages.at(0)->GetOrigin() ); emaxImage->SetDirection( m_TensorImages.at(0)->GetDirection() ); emaxImage->SetRegions( m_TensorImages.at(0)->GetLargestPossibleRegion() ); emaxImage->Allocate(); emaxImage->FillBuffer(1.0); m_EmaxImage.push_back(emaxImage); } bool useUserFaImage = true; if (m_FaImage.IsNull()) { m_FaImage = ItkFloatImgType::New(); m_FaImage->SetSpacing( m_TensorImages.at(0)->GetSpacing() ); m_FaImage->SetOrigin( m_TensorImages.at(0)->GetOrigin() ); m_FaImage->SetDirection( m_TensorImages.at(0)->GetDirection() ); m_FaImage->SetRegions( m_TensorImages.at(0)->GetLargestPossibleRegion() ); m_FaImage->Allocate(); m_FaImage->FillBuffer(0.0); useUserFaImage = false; } typedef itk::DiffusionTensor3D TensorType; TensorType::EigenValuesArrayType eigenvalues; TensorType::EigenVectorsMatrixType eigenvectors; for (unsigned int x=0; xGetLargestPossibleRegion().GetSize()[0]; x++) for (unsigned int y=0; yGetLargestPossibleRegion().GetSize()[1]; y++) for (unsigned int z=0; zGetLargestPossibleRegion().GetSize()[2]; z++) { ItkTensorImageType::IndexType index; index[0] = x; index[1] = y; index[2] = z; for (int i=0; iGetPixel(index); vnl_vector_fixed dir; tensor.ComputeEigenAnalysis(eigenvalues, eigenvectors); dir[0] = eigenvectors(2, 0); dir[1] = eigenvectors(2, 1); dir[2] = eigenvectors(2, 2); if (dir.magnitude()>mitk::eps) dir.normalize(); else dir.fill(0.0); m_PdImage.at(i)->SetPixel(index, dir); if (!useUserFaImage) m_FaImage->SetPixel(index, m_FaImage->GetPixel(index)+tensor.GetFractionalAnisotropy()); m_EmaxImage.at(i)->SetPixel(index, 2/eigenvalues[2]); } if (!useUserFaImage) m_FaImage->SetPixel(index, m_FaImage->GetPixel(index)/m_NumberOfInputs); } if (m_F+m_G>1.0) { float temp = m_F+m_G; m_F /= temp; m_G /= temp; } } vnl_vector_fixed TrackingHandlerTensor::GetMatchingDirection(itk::Index<3> idx, vnl_vector_fixed& oldDir, int& image_num) { vnl_vector_fixed out_dir; out_dir.fill(0); float angle = 0; float mag = oldDir.magnitude(); if (magGetPixel(idx); if (out_dir.magnitude()>0.5) { image_num = i; oldDir[0] = out_dir[0]; oldDir[1] = out_dir[1]; oldDir[2] = out_dir[2]; break; } } } else { for (unsigned int i=0; i dir = m_PdImage.at(i)->GetPixel(idx); float a = dot_product(dir, oldDir); if (fabs(a)>angle) { image_num = i; angle = fabs(a); if (a<0) out_dir = -dir; else out_dir = dir; out_dir *= angle; // shrink contribution of direction if is less parallel to previous direction } } } return out_dir; } vnl_vector_fixed TrackingHandlerTensor::GetDirection(itk::Point itkP, vnl_vector_fixed oldDir, TensorType& tensor) { // transform physical point to index coordinates itk::Index<3> idx; itk::ContinuousIndex< float, 3> cIdx; m_FaImage->TransformPhysicalPointToIndex(itkP, idx); m_FaImage->TransformPhysicalPointToContinuousIndex(itkP, cIdx); vnl_vector_fixed dir; dir.fill(0.0); if ( !m_FaImage->GetLargestPossibleRegion().IsInside(idx) ) return dir; int image_num = -1; if (!m_Interpolate) { dir = GetMatchingDirection(idx, oldDir, image_num); if (image_num>=0) tensor = m_TensorImages[image_num]->GetPixel(idx) * m_EmaxImage[image_num]->GetPixel(idx); } else { float frac_x = cIdx[0] - idx[0]; float frac_y = cIdx[1] - idx[1]; float frac_z = cIdx[2] - idx[2]; if (frac_x<0) { idx[0] -= 1; frac_x += 1; } if (frac_y<0) { idx[1] -= 1; frac_y += 1; } if (frac_z<0) { idx[2] -= 1; frac_z += 1; } frac_x = 1-frac_x; frac_y = 1-frac_y; frac_z = 1-frac_z; // int coordinates inside image? if (idx[0] >= 0 && idx[0] < static_cast(m_FaImage->GetLargestPossibleRegion().GetSize(0) - 1) && idx[1] >= 0 && idx[1] < static_cast(m_FaImage->GetLargestPossibleRegion().GetSize(1) - 1) && idx[2] >= 0 && idx[2] < static_cast(m_FaImage->GetLargestPossibleRegion().GetSize(2) - 1)) { // trilinear interpolation vnl_vector_fixed interpWeights; interpWeights[0] = ( frac_x)*( frac_y)*( frac_z); interpWeights[1] = (1-frac_x)*( frac_y)*( frac_z); interpWeights[2] = ( frac_x)*(1-frac_y)*( frac_z); interpWeights[3] = ( frac_x)*( frac_y)*(1-frac_z); interpWeights[4] = (1-frac_x)*(1-frac_y)*( frac_z); interpWeights[5] = ( frac_x)*(1-frac_y)*(1-frac_z); interpWeights[6] = (1-frac_x)*( frac_y)*(1-frac_z); interpWeights[7] = (1-frac_x)*(1-frac_y)*(1-frac_z); dir = GetMatchingDirection(idx, oldDir, image_num) * interpWeights[0]; if (image_num>=0) tensor += m_TensorImages[image_num]->GetPixel(idx) * m_EmaxImage[image_num]->GetPixel(idx) * interpWeights[0]; itk::Index<3> tmpIdx = idx; tmpIdx[0]++; dir += GetMatchingDirection(tmpIdx, oldDir, image_num) * interpWeights[1]; if (image_num>=0) tensor += m_TensorImages[image_num]->GetPixel(tmpIdx) * m_EmaxImage[image_num]->GetPixel(tmpIdx) * interpWeights[1]; tmpIdx = idx; tmpIdx[1]++; dir += GetMatchingDirection(tmpIdx, oldDir, image_num) * interpWeights[2]; if (image_num>=0) tensor += m_TensorImages[image_num]->GetPixel(tmpIdx) * m_EmaxImage[image_num]->GetPixel(tmpIdx) * interpWeights[2]; tmpIdx = idx; tmpIdx[2]++; dir += GetMatchingDirection(tmpIdx, oldDir, image_num) * interpWeights[3]; if (image_num>=0) tensor += m_TensorImages[image_num]->GetPixel(tmpIdx) * m_EmaxImage[image_num]->GetPixel(tmpIdx) * interpWeights[3]; tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++; dir += GetMatchingDirection(tmpIdx, oldDir, image_num) * interpWeights[4]; if (image_num>=0) tensor += m_TensorImages[image_num]->GetPixel(tmpIdx) * m_EmaxImage[image_num]->GetPixel(tmpIdx) * interpWeights[4]; tmpIdx = idx; tmpIdx[1]++; tmpIdx[2]++; dir += GetMatchingDirection(tmpIdx, oldDir, image_num) * interpWeights[5]; if (image_num>=0) tensor += m_TensorImages[image_num]->GetPixel(tmpIdx) * m_EmaxImage[image_num]->GetPixel(tmpIdx) * interpWeights[5]; tmpIdx = idx; tmpIdx[2]++; tmpIdx[0]++; dir += GetMatchingDirection(tmpIdx, oldDir, image_num) * interpWeights[6]; if (image_num>=0) tensor += m_TensorImages[image_num]->GetPixel(tmpIdx) * m_EmaxImage[image_num]->GetPixel(tmpIdx) * interpWeights[6]; tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++; tmpIdx[2]++; dir += GetMatchingDirection(tmpIdx, oldDir, image_num) * interpWeights[7]; if (image_num>=0) tensor += m_TensorImages[image_num]->GetPixel(tmpIdx) * m_EmaxImage[image_num]->GetPixel(tmpIdx) * interpWeights[7]; } } return dir; } vnl_vector_fixed TrackingHandlerTensor::GetLargestEigenvector(TensorType& tensor) { vnl_vector_fixed dir; TensorType::EigenValuesArrayType eigenvalues; TensorType::EigenVectorsMatrixType eigenvectors; tensor.ComputeEigenAnalysis(eigenvalues, eigenvectors); dir[0] = eigenvectors(2, 0); dir[1] = eigenvectors(2, 1); dir[2] = eigenvectors(2, 2); if (dir.magnitude() TrackingHandlerTensor::ProposeDirection(itk::Point& pos, std::deque >& olddirs, itk::Index<3>& oldIndex) +vnl_vector_fixed TrackingHandlerTensor::ProposeDirection(const itk::Point& pos, std::deque >& olddirs, itk::Index<3>& oldIndex) { vnl_vector_fixed output_direction; output_direction.fill(0); TensorType tensor; tensor.Fill(0); try { itk::Index<3> index; m_TensorImages.at(0)->TransformPhysicalPointToIndex(pos, index); float fa = GetImageValue(pos, m_FaImage, m_Interpolate); if (fa oldDir = olddirs.back(); if (m_FlipX) oldDir[0] *= -1; if (m_FlipY) oldDir[1] *= -1; if (m_FlipZ) oldDir[2] *= -1; float old_mag = oldDir.magnitude(); if (!m_Interpolate && oldIndex==index) return oldDir; output_direction = GetDirection(pos, oldDir, tensor); float mag = output_direction.magnitude(); if (mag>=mitk::eps) { output_direction.normalize(); if (old_mag>0.5 && m_G>mitk::eps) // TEND tracking { output_direction[0] = m_F*output_direction[0] + (1-m_F)*( (1-m_G)*oldDir[0] + m_G*(tensor[0]*oldDir[0] + tensor[1]*oldDir[1] + tensor[2]*oldDir[2])); output_direction[1] = m_F*output_direction[1] + (1-m_F)*( (1-m_G)*oldDir[1] + m_G*(tensor[1]*oldDir[0] + tensor[3]*oldDir[1] + tensor[4]*oldDir[2])); output_direction[2] = m_F*output_direction[2] + (1-m_F)*( (1-m_G)*oldDir[2] + m_G*(tensor[2]*oldDir[0] + tensor[4]*oldDir[1] + tensor[5]*oldDir[2])); output_direction.normalize(); } float a = 1; if (old_mag>0.5) a = dot_product(output_direction, oldDir); if (a>m_AngularThreshold) output_direction *= mag; else output_direction.fill(0); } else output_direction.fill(0); } catch(...) { } if (m_FlipX) output_direction[0] *= -1; if (m_FlipY) output_direction[1] *= -1; if (m_FlipZ) output_direction[2] *= -1; return output_direction; } } diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.h index 01ca4e5405..863e9b9452 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.h @@ -1,87 +1,87 @@ /*=================================================================== 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 _TrackingHandlerTensor #define _TrackingHandlerTensor #include "mitkTrackingDataHandler.h" #include #include namespace mitk { /** * \brief Enables streamline tracking on tensor images. Supports multi tensor tracking by adding multiple tensor images. */ class MITKFIBERTRACKING_EXPORT TrackingHandlerTensor : public TrackingDataHandler { public: TrackingHandlerTensor(); ~TrackingHandlerTensor(); typedef itk::DiffusionTensor3D TensorType; typedef itk::Image< TensorType, 3 > ItkTensorImageType; typedef itk::Image< vnl_vector_fixed, 3> ItkPDImgType; void InitForTracking(); ///< calls InputDataValidForTracking() and creates feature images - vnl_vector_fixed ProposeDirection(itk::Point& pos, std::deque< vnl_vector_fixed >& olddirs, itk::Index<3>& oldIndex); ///< predicts next progression direction at the given position + vnl_vector_fixed ProposeDirection(const itk::Point& pos, std::deque< vnl_vector_fixed >& olddirs, itk::Index<3>& oldIndex); ///< predicts next progression direction at the given position void SetF(float f){ m_F = f; } void SetG(float g){ m_G = g; } void SetFaThreshold(float FaThreshold){ m_FaThreshold = FaThreshold; } void AddTensorImage( ItkTensorImageType::Pointer img ){ m_TensorImages.push_back(img); } void ClearTensorImages(){ m_TensorImages.clear(); } void SetFaImage( ItkFloatImgType::Pointer img ){ m_FaImage = img; } void SetInterpolateTensors( bool interpolateTensors ){ m_InterpolateTensors = interpolateTensors; } void SetMode( MODE m ) { if (m==MODE::DETERMINISTIC) m_Mode = m; else mitkThrow() << "Tensor tracker is only implemented for deterministic mode."; } ItkUcharImgType::SpacingType GetSpacing(){ return m_FaImage->GetSpacing(); } itk::Point GetOrigin(){ return m_FaImage->GetOrigin(); } ItkUcharImgType::DirectionType GetDirection(){ return m_FaImage->GetDirection(); } ItkUcharImgType::RegionType GetLargestPossibleRegion(){ return m_FaImage->GetLargestPossibleRegion(); } protected: vnl_vector_fixed GetMatchingDirection(itk::Index<3> idx, vnl_vector_fixed& oldDir, int& image_num); vnl_vector_fixed GetDirection(itk::Point itkP, vnl_vector_fixed oldDir, TensorType& tensor); vnl_vector_fixed GetLargestEigenvector(TensorType& tensor); float m_FaThreshold; float m_F; float m_G; std::vector< ItkDoubleImgType::Pointer > m_EmaxImage; ///< Stores largest eigenvalues per voxel (one for each tensor) ItkFloatImgType::Pointer m_FaImage; ///< FA image used to determine streamline termination. std::vector< ItkPDImgType::Pointer > m_PdImage; ///< Stores principal direction of each tensor in each voxel. std::vector< ItkTensorImageType::Pointer > m_TensorImages; ///< Input tensor images. For multi tensor tracking provide multiple tensor images. bool m_InterpolateTensors; ///< If false, then the peaks are interpolated. Otherwiese, The tensors are interpolated. int m_NumberOfInputs; }; } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp index 7ef5f29891..3cc79f4a67 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp @@ -1,993 +1,993 @@ /*=================================================================== 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 #include #include #include #include "itkStreamlineTrackingFilter.h" #include #include #include #include #include "itkPointShell.h" #define _USE_MATH_DEFINES #include namespace itk { StreamlineTrackingFilter ::StreamlineTrackingFilter() - : m_PauseTracking(false) - , m_AbortTracking(false) - , m_FiberPolyData(nullptr) - , m_Points(nullptr) - , m_Cells(nullptr) - , m_AngularThreshold(-1) - , m_StepSize(0) - , m_MaxLength(10000) - , m_MinTractLength(20.0) - , m_MaxTractLength(400.0) - , m_SeedsPerVoxel(1) - , m_RandomSampling(false) - , m_SamplingDistance(-1) - , m_DeflectionMod(1.0) - , m_OnlyForwardSamples(true) - , m_UseStopVotes(true) - , m_NumberOfSamples(30) - , m_NumPreviousDirections(1) - , m_StoppingRegions(nullptr) - , m_SeedImage(nullptr) - , m_MaskImage(nullptr) - , m_AposterioriCurvCheck(false) - , m_AvoidStop(true) - , m_DemoMode(false) - , m_SeedOnlyGm(false) - , m_ControlGmEndings(false) - , m_WmLabel(3) // mrtrix 5ttseg labels - , m_GmLabel(1) // mrtrix 5ttseg labels - , m_StepSizeVox(-1) - , m_SamplingDistanceVox(-1) - , m_AngularThresholdDeg(-1) - , m_MaxNumTracts(-1) - , m_Random(true) - , m_Verbose(true) - , m_UseOutputProbabilityMap(false) + : m_PauseTracking(false) + , m_AbortTracking(false) + , m_FiberPolyData(nullptr) + , m_Points(nullptr) + , m_Cells(nullptr) + , m_AngularThreshold(-1) + , m_StepSize(0) + , m_MaxLength(10000) + , m_MinTractLength(20.0) + , m_MaxTractLength(400.0) + , m_SeedsPerVoxel(1) + , m_RandomSampling(false) + , m_SamplingDistance(-1) + , m_DeflectionMod(1.0) + , m_OnlyForwardSamples(true) + , m_UseStopVotes(true) + , m_NumberOfSamples(30) + , m_NumPreviousDirections(1) + , m_StoppingRegions(nullptr) + , m_SeedImage(nullptr) + , m_MaskImage(nullptr) + , m_AposterioriCurvCheck(false) + , m_AvoidStop(true) + , m_DemoMode(false) + , m_SeedOnlyGm(false) + , m_ControlGmEndings(false) + , m_WmLabel(3) // mrtrix 5ttseg labels + , m_GmLabel(1) // mrtrix 5ttseg labels + , m_StepSizeVox(-1) + , m_SamplingDistanceVox(-1) + , m_AngularThresholdDeg(-1) + , m_MaxNumTracts(-1) + , m_Random(true) + , m_Verbose(true) + , m_UseOutputProbabilityMap(false) { - this->SetNumberOfRequiredInputs(0); + this->SetNumberOfRequiredInputs(0); } void StreamlineTrackingFilter::BeforeTracking() { - std::cout.setf(std::ios::boolalpha); - - m_TrackingHandler->InitForTracking(); - - m_FiberPolyData = PolyDataType::New(); - m_Points = vtkSmartPointer< vtkPoints >::New(); - m_Cells = vtkSmartPointer< vtkCellArray >::New(); - - itk::Vector< double, 3 > imageSpacing = m_TrackingHandler->GetSpacing(); - - float minSpacing; - if(imageSpacing[0]InitForTracking(); + m_FiberPolyData = PolyDataType::New(); + m_Points = vtkSmartPointer< vtkPoints >::New(); + m_Cells = vtkSmartPointer< vtkCellArray >::New(); + + itk::Vector< double, 3 > imageSpacing = m_TrackingHandler->GetSpacing(); + + float minSpacing; + if(imageSpacing[0]SetAngularThreshold(m_AngularThreshold); + + if (m_SamplingDistanceVoxGetNumberOfThreads(); i++) + { + PolyDataType poly = PolyDataType::New(); + m_PolyDataContainer.push_back(poly); + } - if (m_StepSizeVoxSetSpacing(imageSpacing); + m_OutputProbabilityMap->SetOrigin(m_TrackingHandler->GetOrigin()); + m_OutputProbabilityMap->SetDirection(m_TrackingHandler->GetDirection()); + m_OutputProbabilityMap->SetRegions(m_TrackingHandler->GetLargestPossibleRegion()); + m_OutputProbabilityMap->Allocate(); + m_OutputProbabilityMap->FillBuffer(0); + } - if (m_AngularThresholdDeg<0) - { - if (m_StepSize/minSpacing<=1.0) - m_AngularThreshold = std::cos( 0.5 * M_PI * m_StepSize/minSpacing ); - else - m_AngularThreshold = std::cos( 0.5 * M_PI ); - } - else - m_AngularThreshold = std::cos( m_AngularThresholdDeg*M_PI/180.0 ); - m_TrackingHandler->SetAngularThreshold(m_AngularThreshold); + if (m_StoppingRegions.IsNull()) + { + m_StoppingRegions = ItkUcharImgType::New(); + m_StoppingRegions->SetSpacing( imageSpacing ); + m_StoppingRegions->SetOrigin( m_TrackingHandler->GetOrigin() ); + m_StoppingRegions->SetDirection( m_TrackingHandler->GetDirection() ); + m_StoppingRegions->SetRegions( m_TrackingHandler->GetLargestPossibleRegion() ); + m_StoppingRegions->Allocate(); + m_StoppingRegions->FillBuffer(0); + } + else + std::cout << "StreamlineTracking - Using stopping region image" << std::endl; + + if (m_SeedImage.IsNull()) + { + m_SeedImage = ItkUcharImgType::New(); + m_SeedImage->SetSpacing( imageSpacing ); + m_SeedImage->SetOrigin( m_TrackingHandler->GetOrigin() ); + m_SeedImage->SetDirection( m_TrackingHandler->GetDirection() ); + m_SeedImage->SetRegions( m_TrackingHandler->GetLargestPossibleRegion() ); + m_SeedImage->Allocate(); + m_SeedImage->FillBuffer(1); + } + else + std::cout << "StreamlineTracking - Using seed image" << std::endl; + + if (m_MaskImage.IsNull()) + { + // initialize mask image + m_MaskImage = ItkUcharImgType::New(); + m_MaskImage->SetSpacing( imageSpacing ); + m_MaskImage->SetOrigin( m_TrackingHandler->GetOrigin() ); + m_MaskImage->SetDirection( m_TrackingHandler->GetDirection() ); + m_MaskImage->SetRegions( m_TrackingHandler->GetLargestPossibleRegion() ); + m_MaskImage->Allocate(); + m_MaskImage->FillBuffer(1); + } + else + std::cout << "StreamlineTracking - Using mask image" << std::endl; - if (m_SamplingDistanceVoxGetNumberOfThreads(); i++) + if (m_TissueImage.IsNull()) + { + if (m_SeedOnlyGm) { - PolyDataType poly = PolyDataType::New(); - m_PolyDataContainer.push_back(poly); + MITK_WARN << "StreamlineTracking - Cannot seed in gray matter. No tissue type image set."; + m_SeedOnlyGm = false; } - if (m_UseOutputProbabilityMap) - { - m_OutputProbabilityMap = ItkDoubleImgType::New(); - m_OutputProbabilityMap->SetSpacing(imageSpacing); - m_OutputProbabilityMap->SetOrigin(m_TrackingHandler->GetOrigin()); - m_OutputProbabilityMap->SetDirection(m_TrackingHandler->GetDirection()); - m_OutputProbabilityMap->SetRegions(m_TrackingHandler->GetLargestPossibleRegion()); - m_OutputProbabilityMap->Allocate(); - m_OutputProbabilityMap->FillBuffer(0); - } - - if (m_StoppingRegions.IsNull()) + if (m_ControlGmEndings) { - m_StoppingRegions = ItkUcharImgType::New(); - m_StoppingRegions->SetSpacing( imageSpacing ); - m_StoppingRegions->SetOrigin( m_TrackingHandler->GetOrigin() ); - m_StoppingRegions->SetDirection( m_TrackingHandler->GetDirection() ); - m_StoppingRegions->SetRegions( m_TrackingHandler->GetLargestPossibleRegion() ); - m_StoppingRegions->Allocate(); - m_StoppingRegions->FillBuffer(0); + MITK_WARN << "StreamlineTracking - Cannot control gray matter endings. No tissue type image set."; + m_ControlGmEndings = false; } + } + else + { + if (m_ControlGmEndings) + m_SeedOnlyGm = true; + if (m_ControlGmEndings || m_SeedOnlyGm) + std::cout << "StreamlineTracking - Using tissue image" << std::endl; else - std::cout << "StreamlineTracking - Using stopping region image" << std::endl; + MITK_WARN << "StreamlineTracking - Tissue image set but no gray matter seeding or fiber endpoint-control enabled" << std::endl; + } - if (m_SeedImage.IsNull()) - { - m_SeedImage = ItkUcharImgType::New(); - m_SeedImage->SetSpacing( imageSpacing ); - m_SeedImage->SetOrigin( m_TrackingHandler->GetOrigin() ); - m_SeedImage->SetDirection( m_TrackingHandler->GetDirection() ); - m_SeedImage->SetRegions( m_TrackingHandler->GetLargestPossibleRegion() ); - m_SeedImage->Allocate(); - m_SeedImage->FillBuffer(1); - } - else - std::cout << "StreamlineTracking - Using seed image" << std::endl; + m_BuildFibersReady = 0; + m_BuildFibersFinished = false; + m_Tractogram.clear(); + m_SamplingPointset = mitk::PointSet::New(); + m_AlternativePointset = mitk::PointSet::New(); + m_StopVotePointset = mitk::PointSet::New(); + m_StartTime = std::chrono::system_clock::now(); + + if (m_SeedOnlyGm && m_ControlGmEndings) + InitGrayMatterEndings(); + + if (m_DemoMode) + omp_set_num_threads(1); + omp_set_num_threads(1); + + if (m_TrackingHandler->GetMode()==mitk::TrackingDataHandler::MODE::DETERMINISTIC) + std::cout << "StreamlineTracking - Mode: deterministic" << std::endl; + else if(m_TrackingHandler->GetMode()==mitk::TrackingDataHandler::MODE::PROBABILISTIC) + std::cout << "StreamlineTracking - Mode: probabilistic" << std::endl; + else + std::cout << "StreamlineTracking - Mode: ???" << std::endl; + + std::cout << "StreamlineTracking - Angular threshold: " << m_AngularThreshold << " (" << 180*std::acos( m_AngularThreshold )/M_PI << "°)" << std::endl; + std::cout << "StreamlineTracking - Stepsize: " << m_StepSize << "mm (" << m_StepSize/minSpacing << "*vox)" << std::endl; + std::cout << "StreamlineTracking - Seeds per voxel: " << m_SeedsPerVoxel << std::endl; + std::cout << "StreamlineTracking - Max. tract length: " << m_MaxTractLength << "mm" << std::endl; + std::cout << "StreamlineTracking - Min. tract length: " << m_MinTractLength << "mm" << std::endl; + std::cout << "StreamlineTracking - Max. num. tracts: " << m_MaxNumTracts << std::endl; + + std::cout << "StreamlineTracking - Num. neighborhood samples: " << m_NumberOfSamples << std::endl; + std::cout << "StreamlineTracking - Max. sampling distance: " << m_SamplingDistance << "mm (" << m_SamplingDistance/minSpacing << "*vox)" << std::endl; + std::cout << "StreamlineTracking - Deflection modifier: " << m_DeflectionMod << std::endl; + + std::cout << "StreamlineTracking - Use stop votes: " << m_UseStopVotes << std::endl; + std::cout << "StreamlineTracking - Only frontal samples: " << m_OnlyForwardSamples << std::endl; + + if (m_DemoMode) + { + std::cout << "StreamlineTracking - Running in demo mode"; + std::cout << "StreamlineTracking - Starting streamline tracking using 1 thread" << std::endl; + } + else + std::cout << "StreamlineTracking - Starting streamline tracking using " << omp_get_max_threads() << " threads" << std::endl; +} - if (m_MaskImage.IsNull()) - { - // initialize mask image - m_MaskImage = ItkUcharImgType::New(); - m_MaskImage->SetSpacing( imageSpacing ); - m_MaskImage->SetOrigin( m_TrackingHandler->GetOrigin() ); - m_MaskImage->SetDirection( m_TrackingHandler->GetDirection() ); - m_MaskImage->SetRegions( m_TrackingHandler->GetLargestPossibleRegion() ); - m_MaskImage->Allocate(); - m_MaskImage->FillBuffer(1); - } - else - std::cout << "StreamlineTracking - Using mask image" << std::endl; - if (m_TissueImage.IsNull()) - { - if (m_SeedOnlyGm) - { - MITK_WARN << "StreamlineTracking - Cannot seed in gray matter. No tissue type image set."; - m_SeedOnlyGm = false; - } - - if (m_ControlGmEndings) - { - MITK_WARN << "StreamlineTracking - Cannot control gray matter endings. No tissue type image set."; - m_ControlGmEndings = false; - } - } - else +void StreamlineTrackingFilter::InitGrayMatterEndings() +{ + m_TrackingHandler->SetAngularThreshold(0); + m_GmStubs.clear(); + if (m_TissueImage.IsNotNull()) + { + std::cout << "StreamlineTracking - initializing GM endings" << std::endl; + ImageRegionConstIterator< ItkUcharImgType > it(m_TissueImage, m_TissueImage->GetLargestPossibleRegion() ); + it.GoToBegin(); + + vnl_vector_fixed d1; d1.fill(0.0); + std::deque< vnl_vector_fixed > olddirs; + while (olddirs.size()GetMode()==mitk::TrackingDataHandler::MODE::DETERMINISTIC) - std::cout << "StreamlineTracking - Mode: deterministic" << std::endl; - else if(m_TrackingHandler->GetMode()==mitk::TrackingDataHandler::MODE::PROBABILISTIC) - std::cout << "StreamlineTracking - Mode: probabilistic" << std::endl; - else - std::cout << "StreamlineTracking - Mode: ???" << std::endl; - - std::cout << "StreamlineTracking - Angular threshold: " << m_AngularThreshold << " (" << 180*std::acos( m_AngularThreshold )/M_PI << "°)" << std::endl; - std::cout << "StreamlineTracking - Stepsize: " << m_StepSize << "mm (" << m_StepSize/minSpacing << "*vox)" << std::endl; - std::cout << "StreamlineTracking - Seeds per voxel: " << m_SeedsPerVoxel << std::endl; - std::cout << "StreamlineTracking - Max. tract length: " << m_MaxTractLength << "mm" << std::endl; - std::cout << "StreamlineTracking - Min. tract length: " << m_MinTractLength << "mm" << std::endl; - std::cout << "StreamlineTracking - Max. num. tracts: " << m_MaxNumTracts << std::endl; - - std::cout << "StreamlineTracking - Num. neighborhood samples: " << m_NumberOfSamples << std::endl; - std::cout << "StreamlineTracking - Max. sampling distance: " << m_SamplingDistance << "mm (" << m_SamplingDistance/minSpacing << "*vox)" << std::endl; - std::cout << "StreamlineTracking - Deflection modifier: " << m_DeflectionMod << std::endl; + if (it.Value()==m_GmLabel) + { + ItkUcharImgType::IndexType s_idx = it.GetIndex(); + itk::ContinuousIndex start; + m_TissueImage->TransformIndexToPhysicalPoint(s_idx, start); + itk::Point wm_p; + float max = -1; + FiberType fib; - std::cout << "StreamlineTracking - Use stop votes: " << m_UseStopVotes << std::endl; - std::cout << "StreamlineTracking - Only frontal samples: " << m_OnlyForwardSamples << std::endl; + for (int x : {-1,0,1}) + for (int y : {-1,0,1}) + for (int z : {-1,0,1}) + { + if (x==y && y==z) + continue; - if (m_DemoMode) - { - std::cout << "StreamlineTracking - Running in demo mode"; - std::cout << "StreamlineTracking - Starting streamline tracking using 1 thread" << std::endl; - } - else - std::cout << "StreamlineTracking - Starting streamline tracking using " << omp_get_max_threads() << " threads" << std::endl; -} + ItkUcharImgType::IndexType e_idx; + e_idx[0] = s_idx[0] + x; + e_idx[1] = s_idx[1] + y; + e_idx[2] = s_idx[2] + z; + if ( !m_TissueImage->GetLargestPossibleRegion().IsInside(e_idx) || m_TissueImage->GetPixel(e_idx)!=m_WmLabel ) + continue; -void StreamlineTrackingFilter::InitGrayMatterEndings() -{ - m_TrackingHandler->SetAngularThreshold(0); - m_GmStubs.clear(); - if (m_TissueImage.IsNotNull()) - { - std::cout << "StreamlineTracking - initializing GM endings" << std::endl; - ImageRegionConstIterator< ItkUcharImgType > it(m_TissueImage, m_TissueImage->GetLargestPossibleRegion() ); - it.GoToBegin(); + itk::ContinuousIndex end; + m_TissueImage->TransformIndexToPhysicalPoint(e_idx, end); - vnl_vector_fixed d1; d1.fill(0.0); - std::deque< vnl_vector_fixed > olddirs; - while (olddirs.size()ProposeDirection(end, olddirs, s_idx); + if (d1.magnitude()<0.0001) + continue; + d1.normalize(); + + vnl_vector_fixed< float, 3 > d2; + d2[0] = end[0] - start[0]; + d2[1] = end[1] - start[1]; + d2[2] = end[2] - start[2]; + d2.normalize(); + float a = fabs(dot_product(d1,d2)); + if (a>max) + { + max = a; + wm_p = end; + } + } - while( !it.IsAtEnd() ) + if (max>=0) { - if (it.Value()==m_GmLabel) - { - ItkUcharImgType::IndexType s_idx = it.GetIndex(); - itk::ContinuousIndex start; - m_TissueImage->TransformIndexToPhysicalPoint(s_idx, start); - itk::Point wm_p; - float max = -1; - FiberType fib; - - for (int x : {-1,0,1}) - for (int y : {-1,0,1}) - for (int z : {-1,0,1}) - { - if (x==y && y==z) - continue; - - ItkUcharImgType::IndexType e_idx; - e_idx[0] = s_idx[0] + x; - e_idx[1] = s_idx[1] + y; - e_idx[2] = s_idx[2] + z; - - if ( !m_TissueImage->GetLargestPossibleRegion().IsInside(e_idx) || m_TissueImage->GetPixel(e_idx)!=m_WmLabel ) - continue; - - itk::ContinuousIndex end; - m_TissueImage->TransformIndexToPhysicalPoint(e_idx, end); - - d1 = m_TrackingHandler->ProposeDirection(end, olddirs, s_idx); - if (d1.magnitude()<0.0001) - continue; - d1.normalize(); - - vnl_vector_fixed< float, 3 > d2; - d2[0] = end[0] - start[0]; - d2[1] = end[1] - start[1]; - d2[2] = end[2] - start[2]; - d2.normalize(); - float a = fabs(dot_product(d1,d2)); - if (a>max) - { - max = a; - wm_p = end; - } - } - - if (max>=0) - { - fib.push_back(start); - fib.push_back(wm_p); - m_GmStubs.push_back(fib); - } - } - ++it; + fib.push_back(start); + fib.push_back(wm_p); + m_GmStubs.push_back(fib); } + } + ++it; } - m_TrackingHandler->SetAngularThreshold(m_AngularThreshold); + } + m_TrackingHandler->SetAngularThreshold(m_AngularThreshold); } void StreamlineTrackingFilter::CalculateNewPosition(itk::Point& pos, vnl_vector_fixed& dir) { - pos[0] += dir[0]*m_StepSize; - pos[1] += dir[1]*m_StepSize; - pos[2] += dir[2]*m_StepSize; + pos[0] += dir[0]*m_StepSize; + pos[1] += dir[1]*m_StepSize; + pos[2] += dir[2]*m_StepSize; } bool StreamlineTrackingFilter -::IsValidPosition(itk::Point &pos) +::IsValidPosition(const itk::Point &pos) { - ItkUcharImgType::IndexType idx; - m_MaskImage->TransformPhysicalPointToIndex(pos, idx); - if (!m_MaskImage->GetLargestPossibleRegion().IsInside(idx) || m_MaskImage->GetPixel(idx)==0) - return false; + ItkUcharImgType::IndexType idx; + m_MaskImage->TransformPhysicalPointToIndex(pos, idx); + if (!m_MaskImage->GetLargestPossibleRegion().IsInside(idx) || m_MaskImage->GetPixel(idx)==0) + return false; - return true; + return true; } bool StreamlineTrackingFilter -::IsInGm(itk::Point &pos) +::IsInGm(const itk::Point &pos) { - if (m_TissueImage.IsNull()) - return true; + if (m_TissueImage.IsNull()) + return true; - ItkUcharImgType::IndexType idx; - m_TissueImage->TransformPhysicalPointToIndex(pos, idx); - if (m_TissueImage->GetLargestPossibleRegion().IsInside(idx) && m_TissueImage->GetPixel(idx) == m_GmLabel) - return true; + ItkUcharImgType::IndexType idx; + m_TissueImage->TransformPhysicalPointToIndex(pos, idx); + if (m_TissueImage->GetLargestPossibleRegion().IsInside(idx) && m_TissueImage->GetPixel(idx) == m_GmLabel) + return true; - return false; + return false; } float StreamlineTrackingFilter::GetRandDouble(float min, float max) { - return (float)(rand()%((int)(10000*(max-min))) + 10000*min)/10000; + return (float)(rand()%((int)(10000*(max-min))) + 10000*min)/10000; } std::vector< vnl_vector_fixed > StreamlineTrackingFilter::CreateDirections(int NPoints) { - std::vector< vnl_vector_fixed > pointshell; + std::vector< vnl_vector_fixed > pointshell; - if (NPoints<2) - return pointshell; + if (NPoints<2) + return pointshell; - std::vector< float > theta; theta.resize(NPoints); + std::vector< float > theta; theta.resize(NPoints); - std::vector< float > phi; phi.resize(NPoints); + std::vector< float > phi; phi.resize(NPoints); - float C = sqrt(4*M_PI); + float C = sqrt(4*M_PI); - phi[0] = 0.0; - phi[NPoints-1] = 0.0; + phi[0] = 0.0; + phi[NPoints-1] = 0.0; - for(int i=0; i0 && i0 && i d; - d[0] = cos(theta[i]) * cos(phi[i]); - d[1] = cos(theta[i]) * sin(phi[i]); - d[2] = sin(theta[i]); - pointshell.push_back(d); - } + for(int i=0; i d; + d[0] = cos(theta[i]) * cos(phi[i]); + d[1] = cos(theta[i]) * sin(phi[i]); + d[2] = sin(theta[i]); + pointshell.push_back(d); + } - return pointshell; + return pointshell; } vnl_vector_fixed StreamlineTrackingFilter::GetNewDirection(itk::Point &pos, std::deque >& olddirs, itk::Index<3> &oldIndex) { - if (m_DemoMode) - { - m_SamplingPointset->Clear(); - m_AlternativePointset->Clear(); - m_StopVotePointset->Clear(); - } - vnl_vector_fixed direction; direction.fill(0); + if (m_DemoMode) + { + m_SamplingPointset->Clear(); + m_AlternativePointset->Clear(); + m_StopVotePointset->Clear(); + } + vnl_vector_fixed direction; direction.fill(0); - ItkUcharImgType::IndexType idx; - m_MaskImage->TransformPhysicalPointToIndex(pos, idx); + ItkUcharImgType::IndexType idx; + m_MaskImage->TransformPhysicalPointToIndex(pos, idx); - if (IsValidPosition(pos)) + if (IsValidPosition(pos)) + { + if (m_StoppingRegions->GetPixel(idx)>0) + return direction; + direction = m_TrackingHandler->ProposeDirection(pos, olddirs, oldIndex); // get direction proposal at current streamline position + } + else + return direction; + + vnl_vector_fixed olddir = olddirs.back(); + std::vector< vnl_vector_fixed > probeVecs = CreateDirections(m_NumberOfSamples); + itk::Point sample_pos; + int alternatives = 1; + int stop_votes = 0; + int possible_stop_votes = 0; + for (int i=0; i d; + bool is_stop_voter = false; + if (m_RandomSampling) { - if (m_StoppingRegions->GetPixel(idx)>0) - return direction; - direction = m_TrackingHandler->ProposeDirection(pos, olddirs, oldIndex); // get direction proposal at current streamline position + d[0] = GetRandDouble(); + d[1] = GetRandDouble(); + d[2] = GetRandDouble(); + d.normalize(); + d *= GetRandDouble(0,m_SamplingDistance); } else - return direction; - - vnl_vector_fixed olddir = olddirs.back(); - std::vector< vnl_vector_fixed > probeVecs = CreateDirections(m_NumberOfSamples); - itk::Point sample_pos; - int alternatives = 1; - int stop_votes = 0; - int possible_stop_votes = 0; - for (int i=0; i d; - bool is_stop_voter = false; - if (m_RandomSampling) - { - d[0] = GetRandDouble(); - d[1] = GetRandDouble(); - d[2] = GetRandDouble(); - d.normalize(); - d *= GetRandDouble(0,m_SamplingDistance); - } - else - { - d = probeVecs.at(i); - float dot = dot_product(d, olddir); - if (m_UseStopVotes && dot>0.7) - { - is_stop_voter = true; - possible_stop_votes++; - } - else if (m_OnlyForwardSamples && dot<0) - continue; - d *= m_SamplingDistance; - } - - sample_pos[0] = pos[0] + d[0]; - sample_pos[1] = pos[1] + d[1]; - sample_pos[2] = pos[2] + d[2]; + d = probeVecs.at(i); + float dot = dot_product(d, olddir); + if (m_UseStopVotes && dot>0.7) + { + is_stop_voter = true; + possible_stop_votes++; + } + else if (m_OnlyForwardSamples && dot<0) + continue; + d *= m_SamplingDistance; + } - vnl_vector_fixed tempDir; tempDir.fill(0.0); - if (IsValidPosition(sample_pos)) - tempDir = m_TrackingHandler->ProposeDirection(sample_pos, olddirs, oldIndex); // sample neighborhood - if (tempDir.magnitude()>mitk::eps) - { - direction += tempDir; + sample_pos[0] = pos[0] + d[0]; + sample_pos[1] = pos[1] + d[1]; + sample_pos[2] = pos[2] + d[2]; - if(m_DemoMode) - m_SamplingPointset->InsertPoint(i, sample_pos); - } - else if (m_AvoidStop && olddir.magnitude()>0.5) // out of white matter - { - if (is_stop_voter) - stop_votes++; - if (m_DemoMode) - m_StopVotePointset->InsertPoint(i, sample_pos); - - float dot = dot_product(d, olddir); - if (dot >= 0.0) // in front of plane defined by pos and olddir - d = -d + 2*dot*olddir; // reflect - else - d = -d; // invert - - // look a bit further into the other direction - sample_pos[0] = pos[0] + d[0]; - sample_pos[1] = pos[1] + d[1]; - sample_pos[2] = pos[2] + d[2]; - alternatives++; - vnl_vector_fixed tempDir; tempDir.fill(0.0); - if (IsValidPosition(sample_pos)) - tempDir = m_TrackingHandler->ProposeDirection(sample_pos, olddirs, oldIndex); // sample neighborhood - - if (tempDir.magnitude()>mitk::eps) // are we back in the white matter? - { - direction += d * m_DeflectionMod; // go into the direction of the white matter - direction += tempDir; // go into the direction of the white matter direction at this location + vnl_vector_fixed tempDir; tempDir.fill(0.0); + if (IsValidPosition(sample_pos)) + tempDir = m_TrackingHandler->ProposeDirection(sample_pos, olddirs, oldIndex); // sample neighborhood + if (tempDir.magnitude()>mitk::eps) + { + direction += tempDir; - if(m_DemoMode) - m_AlternativePointset->InsertPoint(alternatives, sample_pos); - } - else - { - if (m_DemoMode) - m_StopVotePointset->InsertPoint(i, sample_pos); - } - } - else - { - if (m_DemoMode) - m_StopVotePointset->InsertPoint(i, sample_pos); + if(m_DemoMode) + m_SamplingPointset->InsertPoint(i, sample_pos); + } + else if (m_AvoidStop && olddir.magnitude()>0.5) // out of white matter + { + if (is_stop_voter) + stop_votes++; + if (m_DemoMode) + m_StopVotePointset->InsertPoint(i, sample_pos); + + float dot = dot_product(d, olddir); + if (dot >= 0.0) // in front of plane defined by pos and olddir + d = -d + 2*dot*olddir; // reflect + else + d = -d; // invert + + // look a bit further into the other direction + sample_pos[0] = pos[0] + d[0]; + sample_pos[1] = pos[1] + d[1]; + sample_pos[2] = pos[2] + d[2]; + alternatives++; + vnl_vector_fixed tempDir; tempDir.fill(0.0); + if (IsValidPosition(sample_pos)) + tempDir = m_TrackingHandler->ProposeDirection(sample_pos, olddirs, oldIndex); // sample neighborhood + + if (tempDir.magnitude()>mitk::eps) // are we back in the white matter? + { + direction += d * m_DeflectionMod; // go into the direction of the white matter + direction += tempDir; // go into the direction of the white matter direction at this location + + if(m_DemoMode) + m_AlternativePointset->InsertPoint(alternatives, sample_pos); + } + else + { + if (m_DemoMode) + m_StopVotePointset->InsertPoint(i, sample_pos); + } + } + else + { + if (m_DemoMode) + m_StopVotePointset->InsertPoint(i, sample_pos); - if (is_stop_voter) - stop_votes++; - } + if (is_stop_voter) + stop_votes++; } + } - if (direction.magnitude()>0.001 && (possible_stop_votes==0 || (float)stop_votes/possible_stop_votes<0.5) ) - direction.normalize(); - else - direction.fill(0); + if (direction.magnitude()>0.001 && (possible_stop_votes==0 || (float)stop_votes/possible_stop_votes<0.5) ) + direction.normalize(); + else + direction.fill(0); - return direction; + return direction; } float StreamlineTrackingFilter::FollowStreamline(itk::Point pos, vnl_vector_fixed dir, FiberType* fib, float tractLength, bool front) { - vnl_vector_fixed zero_dir; zero_dir.fill(0.0); - std::deque< vnl_vector_fixed > last_dirs; - for (unsigned int i=0; i zero_dir; zero_dir.fill(0.0); + std::deque< vnl_vector_fixed > last_dirs; + for (unsigned int i=0; iTransformPhysicalPointToIndex(pos, oldIndex); + for (int step=0; step< m_MaxLength/2; step++) + { + ItkUcharImgType::IndexType oldIndex; + m_StoppingRegions->TransformPhysicalPointToIndex(pos, oldIndex); - // get new position - CalculateNewPosition(pos, dir); + // get new position + CalculateNewPosition(pos, dir); - // is new position inside of image and mask - if (m_AbortTracking) // if not end streamline - { - return tractLength; - } - else // if yes, add new point to streamline + // is new position inside of image and mask + if (m_AbortTracking) // if not end streamline + { + return tractLength; + } + else // if yes, add new point to streamline + { + tractLength += m_StepSize; + if (front) + fib->push_front(pos); + else + fib->push_back(pos); + + if (m_AposterioriCurvCheck) + { + int curv = CheckCurvature(fib, front); + if (curv>0) { - tractLength += m_StepSize; + tractLength -= m_StepSize*curv; + while (curv>0) + { if (front) - fib->push_front(pos); + fib->pop_front(); else - fib->push_back(pos); - - if (m_AposterioriCurvCheck) - { - int curv = CheckCurvature(fib, front); - if (curv>0) - { - tractLength -= m_StepSize*curv; - while (curv>0) - { - if (front) - fib->pop_front(); - else - fib->pop_back(); - curv--; - } - return tractLength; - } - } - - if (tractLength>m_MaxTractLength) - return tractLength; + fib->pop_back(); + curv--; + } + return tractLength; } + } + + if (tractLength>m_MaxTractLength) + return tractLength; + } #pragma omp critical - if (m_DemoMode && !m_UseOutputProbabilityMap) // CHECK: warum sind die samplingpunkte der streamline in der visualisierung immer einen schritt voras? - { - m_BuildFibersReady++; - m_Tractogram.push_back(*fib); - BuildFibers(true); - m_Stop = true; + if (m_DemoMode && !m_UseOutputProbabilityMap) // CHECK: warum sind die samplingpunkte der streamline in der visualisierung immer einen schritt voras? + { + m_BuildFibersReady++; + m_Tractogram.push_back(*fib); + BuildFibers(true); + m_Stop = true; - while (m_Stop){ - } - } + while (m_Stop){ + } + } - dir.normalize(); - last_dirs.push_back(dir); - if (last_dirs.size()>m_NumPreviousDirections) - last_dirs.pop_front(); - dir = GetNewDirection(pos, last_dirs, oldIndex); + dir.normalize(); + last_dirs.push_back(dir); + if (last_dirs.size()>m_NumPreviousDirections) + last_dirs.pop_front(); + dir = GetNewDirection(pos, last_dirs, oldIndex); - while (m_PauseTracking){} + while (m_PauseTracking){} - if (dir.magnitude()<0.0001) - return tractLength; - } - return tractLength; + if (dir.magnitude()<0.0001) + return tractLength; + } + return tractLength; } int StreamlineTrackingFilter::CheckCurvature(FiberType* fib, bool front) { - float m_Distance = 5; - if (fib->size()<3) - return 0; - - float dist = 0; - std::vector< vnl_vector_fixed< float, 3 > > vectors; - vnl_vector_fixed< float, 3 > meanV; meanV.fill(0); - float dev = 0; - - if (front) + float m_Distance = 5; + if (fib->size()<3) + return 0; + + float dist = 0; + std::vector< vnl_vector_fixed< float, 3 > > vectors; + vnl_vector_fixed< float, 3 > meanV; meanV.fill(0); + float dev = 0; + + if (front) + { + int c=0; + while(distsize()-1) { - int c=0; - while(distsize()-1) - { - itk::Point p1 = fib->at(c); - itk::Point p2 = fib->at(c+1); - - vnl_vector_fixed< float, 3 > v; - v[0] = p2[0]-p1[0]; - v[1] = p2[1]-p1[1]; - v[2] = p2[2]-p1[2]; - dist += v.magnitude(); - v.normalize(); - vectors.push_back(v); - if (c==0) - meanV += v; - c++; - } - } - else - { - int c=fib->size()-1; - while(dist0) - { - itk::Point p1 = fib->at(c); - itk::Point p2 = fib->at(c-1); - - vnl_vector_fixed< float, 3 > v; - v[0] = p2[0]-p1[0]; - v[1] = p2[1]-p1[1]; - v[2] = p2[2]-p1[2]; - dist += v.magnitude(); - v.normalize(); - vectors.push_back(v); - if (c==fib->size()-1) - meanV += v; - c--; - } + itk::Point p1 = fib->at(c); + itk::Point p2 = fib->at(c+1); + + vnl_vector_fixed< float, 3 > v; + v[0] = p2[0]-p1[0]; + v[1] = p2[1]-p1[1]; + v[2] = p2[2]-p1[2]; + dist += v.magnitude(); + v.normalize(); + vectors.push_back(v); + if (c==0) + meanV += v; + c++; } - meanV.normalize(); - - for (int c=0; csize()-1; + while(dist0) { - float angle = dot_product(meanV, vectors.at(c)); - if (angle>1.0) - angle = 1.0; - if (angle<-1.0) - angle = -1.0; - dev += acos(angle)*180/M_PI; + itk::Point p1 = fib->at(c); + itk::Point p2 = fib->at(c-1); + + vnl_vector_fixed< float, 3 > v; + v[0] = p2[0]-p1[0]; + v[1] = p2[1]-p1[1]; + v[2] = p2[2]-p1[2]; + dist += v.magnitude(); + v.normalize(); + vectors.push_back(v); + if (c==fib->size()-1) + meanV += v; + c--; } - if (vectors.size()>0) - dev /= vectors.size(); + } + meanV.normalize(); + + for (int c=0; c1.0) + angle = 1.0; + if (angle<-1.0) + angle = -1.0; + dev += acos(angle)*180/M_PI; + } + if (vectors.size()>0) + dev /= vectors.size(); - if (dev<30) - return 0; - else - return vectors.size(); + if (dev<30) + return 0; + else + return vectors.size(); } void StreamlineTrackingFilter::GetSeedPointsFromSeedImage() { - MITK_INFO << "Calculating seed points."; - m_SeedPoints.clear(); - - if (!m_ControlGmEndings) - { - typedef ImageRegionConstIterator< ItkUcharImgType > MaskIteratorType; - MaskIteratorType sit(m_SeedImage, m_SeedImage->GetLargestPossibleRegion()); - sit.GoToBegin(); - - while (!sit.IsAtEnd()) - { - if (sit.Value()>0) - { - ItkUcharImgType::IndexType index = sit.GetIndex(); - itk::ContinuousIndex start; - start[0] = index[0]; - start[1] = index[1]; - start[2] = index[2]; - itk::Point worldPos; - m_SeedImage->TransformContinuousIndexToPhysicalPoint(start, worldPos); - - if (IsValidPosition(worldPos) && (!m_SeedOnlyGm || IsInGm(worldPos))) - { - m_SeedPoints.push_back(worldPos); - for (int s = 1; s < m_SeedsPerVoxel; s++) - { - start[0] = index[0] + GetRandDouble(-0.5, 0.5); - start[1] = index[1] + GetRandDouble(-0.5, 0.5); - start[2] = index[2] + GetRandDouble(-0.5, 0.5); - - itk::Point worldPos; - m_SeedImage->TransformContinuousIndexToPhysicalPoint(start, worldPos); - m_SeedPoints.push_back(worldPos); - } - } - } - ++sit; - } - } - else - { - for (auto s : m_GmStubs) - m_SeedPoints.push_back(s[1]); - } + MITK_INFO << "Calculating seed points."; + m_SeedPoints.clear(); + + if (!m_ControlGmEndings) + { + typedef ImageRegionConstIterator< ItkUcharImgType > MaskIteratorType; + MaskIteratorType sit(m_SeedImage, m_SeedImage->GetLargestPossibleRegion()); + sit.GoToBegin(); + + while (!sit.IsAtEnd()) + { + if (sit.Value()>0) + { + ItkUcharImgType::IndexType index = sit.GetIndex(); + itk::ContinuousIndex start; + start[0] = index[0]; + start[1] = index[1]; + start[2] = index[2]; + itk::Point worldPos; + m_SeedImage->TransformContinuousIndexToPhysicalPoint(start, worldPos); + + if (IsValidPosition(worldPos) && (!m_SeedOnlyGm || IsInGm(worldPos))) + { + m_SeedPoints.push_back(worldPos); + for (int s = 1; s < m_SeedsPerVoxel; s++) + { + start[0] = index[0] + GetRandDouble(-0.5, 0.5); + start[1] = index[1] + GetRandDouble(-0.5, 0.5); + start[2] = index[2] + GetRandDouble(-0.5, 0.5); + + itk::Point worldPos; + m_SeedImage->TransformContinuousIndexToPhysicalPoint(start, worldPos); + m_SeedPoints.push_back(worldPos); + } + } + } + ++sit; + } + } + else + { + for (auto s : m_GmStubs) + m_SeedPoints.push_back(s[1]); + } } void StreamlineTrackingFilter::GenerateData() { - this->BeforeTracking(); - - if (m_Random) - { - std::srand(std::time(0)); - std::random_shuffle(m_SeedPoints.begin(), m_SeedPoints.end()); - } + this->BeforeTracking(); + + if (m_Random) + { + std::srand(std::time(0)); + std::random_shuffle(m_SeedPoints.begin(), m_SeedPoints.end()); + } - bool stop = false; - unsigned int current_tracts = 0; - int num_seeds = m_SeedPoints.size(); - itk::Index<3> zeroIndex; zeroIndex.Fill(0); - int progress = 0; - int i = 0; - int print_interval = num_seeds/100; - if (print_interval<100) - m_Verbose=false; + bool stop = false; + unsigned int current_tracts = 0; + int num_seeds = m_SeedPoints.size(); + itk::Index<3> zeroIndex; zeroIndex.Fill(0); + int progress = 0; + int i = 0; + int print_interval = num_seeds/100; + if (print_interval<100) + m_Verbose=false; #pragma omp parallel - while (i=num_seeds || stop) - continue; - else if (m_Verbose && i%print_interval==0) + if (temp_i>=num_seeds || stop) + continue; + else if (m_Verbose && i%print_interval==0) #pragma omp critical - { - progress += print_interval; - std::cout << " \r"; - if (m_MaxNumTracts>0) - std::cout << "Tried: " << progress << "/" << num_seeds << " | Accepted: " << current_tracts << "/" << m_MaxNumTracts << '\r'; - else - std::cout << "Tried: " << progress << "/" << num_seeds << " | Accepted: " << current_tracts << '\r'; - cout.flush(); - } + { + progress += print_interval; + std::cout << " \r"; + if (m_MaxNumTracts>0) + std::cout << "Tried: " << progress << "/" << num_seeds << " | Accepted: " << current_tracts << "/" << m_MaxNumTracts << '\r'; + else + std::cout << "Tried: " << progress << "/" << num_seeds << " | Accepted: " << current_tracts << '\r'; + cout.flush(); + } - itk::Point worldPos = m_SeedPoints.at(temp_i); - FiberType fib; - float tractLength = 0; - unsigned int counter = 0; + const itk::Point worldPos = m_SeedPoints.at(temp_i); + FiberType fib; + float tractLength = 0; + unsigned int counter = 0; - // get starting direction - vnl_vector_fixed dir; dir.fill(0.0); - std::deque< vnl_vector_fixed > olddirs; - while (olddirs.size() dir; dir.fill(0.0); + std::deque< vnl_vector_fixed > olddirs; + while (olddirs.size() gm_start_dir; - if (m_ControlGmEndings) - { - gm_start_dir[0] = m_GmStubs[temp_i][1][0] - m_GmStubs[temp_i][0][0]; - gm_start_dir[1] = m_GmStubs[temp_i][1][1] - m_GmStubs[temp_i][0][1]; - gm_start_dir[2] = m_GmStubs[temp_i][1][2] - m_GmStubs[temp_i][0][2]; - gm_start_dir.normalize(); - olddirs.pop_back(); - olddirs.push_back(gm_start_dir); - } + vnl_vector_fixed< float, 3 > gm_start_dir; + if (m_ControlGmEndings) + { + gm_start_dir[0] = m_GmStubs[temp_i][1][0] - m_GmStubs[temp_i][0][0]; + gm_start_dir[1] = m_GmStubs[temp_i][1][1] - m_GmStubs[temp_i][0][1]; + gm_start_dir[2] = m_GmStubs[temp_i][1][2] - m_GmStubs[temp_i][0][2]; + gm_start_dir.normalize(); + olddirs.pop_back(); + olddirs.push_back(gm_start_dir); + } - if (IsValidPosition(worldPos)) - dir = m_TrackingHandler->ProposeDirection(worldPos, olddirs, zeroIndex); + if (IsValidPosition(worldPos)) + dir = m_TrackingHandler->ProposeDirection(worldPos, olddirs, zeroIndex); - if (dir.magnitude()>0.0001) + if (dir.magnitude()>0.0001) + { + if (m_ControlGmEndings) + { + float a = dot_product(gm_start_dir, dir); + if (a<0) + dir = -dir; + } + + // forward tracking + tractLength = FollowStreamline(worldPos, dir, &fib, 0, false); + fib.push_front(worldPos); + + if (m_ControlGmEndings) + { + fib.push_front(m_GmStubs[temp_i][0]); + CheckFiberForGmEnding(&fib); + } + else + { + // backward tracking (only if we don't explicitely start in the GM) + tractLength = FollowStreamline(worldPos, -dir, &fib, tractLength, true); + if (m_ControlGmEndings) { - if (m_ControlGmEndings) - { - float a = dot_product(gm_start_dir, dir); - if (a<0) - dir = -dir; - } - - // forward tracking - tractLength = FollowStreamline(worldPos, dir, &fib, 0, false); - fib.push_front(worldPos); - - if (m_ControlGmEndings) - { - fib.push_front(m_GmStubs[temp_i][0]); - CheckFiberForGmEnding(&fib); - } - else - { - // backward tracking (only if we don't explicitely start in the GM) - tractLength = FollowStreamline(worldPos, -dir, &fib, tractLength, true); - if (m_ControlGmEndings) - { - CheckFiberForGmEnding(&fib); - std::reverse(fib.begin(),fib.end()); - CheckFiberForGmEnding(&fib); - } - } - counter = fib.size(); + CheckFiberForGmEnding(&fib); + std::reverse(fib.begin(),fib.end()); + CheckFiberForGmEnding(&fib); + } + } + counter = fib.size(); #pragma omp critical - if (tractLength>=m_MinTractLength && counter>=2) - { - if (!stop) - { - if (!m_UseOutputProbabilityMap) - m_Tractogram.push_back(fib); - else - FiberToProbmap(&fib); - current_tracts++; - } - if (m_MaxNumTracts > 0 && current_tracts>=static_cast(m_MaxNumTracts)) - { - if (!stop) - { - std::cout << " \r"; - MITK_INFO << "Reconstructed maximum number of tracts (" << current_tracts << "). Stopping tractography."; - } - stop = true; - } - } + if (tractLength>=m_MinTractLength && counter>=2) + { + if (!stop) + { + if (!m_UseOutputProbabilityMap) + m_Tractogram.push_back(fib); + else + FiberToProbmap(&fib); + current_tracts++; } + if (m_MaxNumTracts > 0 && current_tracts>=static_cast(m_MaxNumTracts)) + { + if (!stop) + { + std::cout << " \r"; + MITK_INFO << "Reconstructed maximum number of tracts (" << current_tracts << "). Stopping tractography."; + } + stop = true; + } + } } + } - this->AfterTracking(); + this->AfterTracking(); } void StreamlineTrackingFilter::CheckFiberForGmEnding(FiberType* fib) { - if (m_TissueImage.IsNull()) - return; - - // first check if the current fibe rendpoint is located inside of the white matter - // if not, remove last fiber point and repeat - bool in_wm = false; - while (!in_wm && fib->size()>2) - { - ItkUcharImgType::IndexType idx; - m_TissueImage->TransformPhysicalPointToIndex(fib->back(), idx); - if (m_TissueImage->GetPixel(idx)==m_WmLabel) - in_wm = true; - else - fib->pop_back(); - } - if (fib->size()<3 || !in_wm) - { - fib->clear(); - return; - } + if (m_TissueImage.IsNull()) + return; + + // first check if the current fibe rendpoint is located inside of the white matter + // if not, remove last fiber point and repeat + bool in_wm = false; + while (!in_wm && fib->size()>2) + { + ItkUcharImgType::IndexType idx; + m_TissueImage->TransformPhysicalPointToIndex(fib->back(), idx); + if (m_TissueImage->GetPixel(idx)==m_WmLabel) + in_wm = true; + else + fib->pop_back(); + } + if (fib->size()<3 || !in_wm) + { + fib->clear(); + return; + } - // get fiber direction at end point - vnl_vector_fixed< float, 3 > d1; - d1[0] = fib->back()[0] - fib->at(fib->size()-2)[0]; - d1[1] = fib->back()[1] - fib->at(fib->size()-2)[1]; - d1[2] = fib->back()[2] - fib->at(fib->size()-2)[2]; - d1.normalize(); - - // find closest gray matter voxel - ItkUcharImgType::IndexType s_idx; - m_TissueImage->TransformPhysicalPointToIndex(fib->back(), s_idx); - itk::Point gm_endp; - float max = -1; - - for (int x : {-1,0,1}) - for (int y : {-1,0,1}) - for (int z : {-1,0,1}) - { - if (x==y && y==z) - continue; - - ItkUcharImgType::IndexType e_idx; - e_idx[0] = s_idx[0] + x; - e_idx[1] = s_idx[1] + y; - e_idx[2] = s_idx[2] + z; - - if ( !m_TissueImage->GetLargestPossibleRegion().IsInside(e_idx) || m_TissueImage->GetPixel(e_idx)!=m_GmLabel ) - continue; - - itk::ContinuousIndex end; - m_TissueImage->TransformIndexToPhysicalPoint(e_idx, end); - vnl_vector_fixed< float, 3 > d2; - d2[0] = end[0] - fib->back()[0]; - d2[1] = end[1] - fib->back()[1]; - d2[2] = end[2] - fib->back()[2]; - d2.normalize(); - float a = dot_product(d1,d2); - if (a>max) - { - max = a; - gm_endp = end; - } - } + // get fiber direction at end point + vnl_vector_fixed< float, 3 > d1; + d1[0] = fib->back()[0] - fib->at(fib->size()-2)[0]; + d1[1] = fib->back()[1] - fib->at(fib->size()-2)[1]; + d1[2] = fib->back()[2] - fib->at(fib->size()-2)[2]; + d1.normalize(); + + // find closest gray matter voxel + ItkUcharImgType::IndexType s_idx; + m_TissueImage->TransformPhysicalPointToIndex(fib->back(), s_idx); + itk::Point gm_endp; + float max = -1; + + for (int x : {-1,0,1}) + for (int y : {-1,0,1}) + for (int z : {-1,0,1}) + { + if (x==y && y==z) + continue; + + ItkUcharImgType::IndexType e_idx; + e_idx[0] = s_idx[0] + x; + e_idx[1] = s_idx[1] + y; + e_idx[2] = s_idx[2] + z; + + if ( !m_TissueImage->GetLargestPossibleRegion().IsInside(e_idx) || m_TissueImage->GetPixel(e_idx)!=m_GmLabel ) + continue; + + itk::ContinuousIndex end; + m_TissueImage->TransformIndexToPhysicalPoint(e_idx, end); + vnl_vector_fixed< float, 3 > d2; + d2[0] = end[0] - fib->back()[0]; + d2[1] = end[1] - fib->back()[1]; + d2[2] = end[2] - fib->back()[2]; + d2.normalize(); + float a = dot_product(d1,d2); + if (a>max) + { + max = a; + gm_endp = end; + } + } - if (max>=0) - fib->push_back(gm_endp); - else // no gray matter enpoint found -> delete fiber - fib->clear(); + if (max>=0) + fib->push_back(gm_endp); + else // no gray matter enpoint found -> delete fiber + fib->clear(); } void StreamlineTrackingFilter::FiberToProbmap(FiberType* fib) { ItkDoubleImgType::IndexType last_idx; last_idx.Fill(0); for (auto p : *fib) - { + { ItkDoubleImgType::IndexType idx; - m_OutputProbabilityMap->TransformPhysicalPointToIndex(p, idx); + m_OutputProbabilityMap->TransformPhysicalPointToIndex(p, idx); if (idx != last_idx) { if (m_OutputProbabilityMap->GetLargestPossibleRegion().IsInside(idx)) - m_OutputProbabilityMap->SetPixel(idx, m_OutputProbabilityMap->GetPixel(idx)+1); - last_idx = idx; - } + m_OutputProbabilityMap->SetPixel(idx, m_OutputProbabilityMap->GetPixel(idx)+1); + last_idx = idx; + } } } void StreamlineTrackingFilter::BuildFibers(bool check) { - if (m_BuildFibersReady::New(); - vtkSmartPointer vNewLines = vtkSmartPointer::New(); - vtkSmartPointer vNewPoints = vtkSmartPointer::New(); - - for (int i=0; i::New(); + vtkSmartPointer vNewLines = vtkSmartPointer::New(); + vtkSmartPointer vNewPoints = vtkSmartPointer::New(); + + for (int i=0; i container = vtkSmartPointer::New(); + FiberType fib = m_Tractogram.at(i); + for (FiberType::iterator it = fib.begin(); it!=fib.end(); ++it) { - vtkSmartPointer container = vtkSmartPointer::New(); - FiberType fib = m_Tractogram.at(i); - for (FiberType::iterator it = fib.begin(); it!=fib.end(); ++it) - { - vtkIdType id = vNewPoints->InsertNextPoint((*it).GetDataPointer()); - container->GetPointIds()->InsertNextId(id); - } - vNewLines->InsertNextCell(container); + vtkIdType id = vNewPoints->InsertNextPoint((*it).GetDataPointer()); + container->GetPointIds()->InsertNextId(id); } + vNewLines->InsertNextCell(container); + } - if (check) - for (int i=0; iSetPoints(vNewPoints); - m_FiberPolyData->SetLines(vNewLines); - m_BuildFibersFinished = true; + m_FiberPolyData->SetPoints(vNewPoints); + m_FiberPolyData->SetLines(vNewLines); + m_BuildFibersFinished = true; } void StreamlineTrackingFilter::AfterTracking() { if (m_Verbose) - std::cout << " \r"; - if (!m_UseOutputProbabilityMap) - { - MITK_INFO << "Reconstructed " << m_Tractogram.size() << " fibers."; - MITK_INFO << "Generating polydata "; - BuildFibers(false); - } + std::cout << " \r"; + if (!m_UseOutputProbabilityMap) + { + MITK_INFO << "Reconstructed " << m_Tractogram.size() << " fibers."; + MITK_INFO << "Generating polydata "; + BuildFibers(false); + } MITK_INFO << "done"; m_EndTime = std::chrono::system_clock::now(); std::chrono::hours hh = std::chrono::duration_cast(m_EndTime - m_StartTime); std::chrono::minutes mm = std::chrono::duration_cast(m_EndTime - m_StartTime); std::chrono::seconds ss = std::chrono::duration_cast(m_EndTime - m_StartTime); mm %= 60; ss %= 60; MITK_INFO << "Tracking took " << hh.count() << "h, " << mm.count() << "m and " << ss.count() << "s"; - m_SeedPoints.clear(); + m_SeedPoints.clear(); } } diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.h index 43461d9825..e3b7c75927 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.h @@ -1,211 +1,211 @@ /*=================================================================== 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 __itkMLBSTrackingFilter_h_ #define __itkMLBSTrackingFilter_h_ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include namespace itk{ /** * \brief Performs streamline tracking on the input image. Depending on the tracking handler this can be a tensor, peak or machine learning based tracking. */ class MITKFIBERTRACKING_EXPORT StreamlineTrackingFilter : public ProcessObject { public: - typedef StreamlineTrackingFilter Self; - typedef SmartPointer Pointer; - typedef SmartPointer ConstPointer; - typedef ProcessObject Superclass; - - /** Method for creation through the object factory. */ - itkFactorylessNewMacro(Self) - itkCloneMacro(Self) - - /** Runtime information support. */ - itkTypeMacro(MLBSTrackingFilter, ImageToImageFilter) - - typedef itk::Image ItkUcharImgType; - typedef itk::Image ItkDoubleImgType; - typedef itk::Image ItkFloatImgType; - typedef vtkSmartPointer< vtkPolyData > PolyDataType; - - typedef std::deque< itk::Point > FiberType; - typedef std::vector< FiberType > BundleType; - - volatile bool m_PauseTracking; - bool m_AbortTracking; - bool m_BuildFibersFinished; - int m_BuildFibersReady; - volatile bool m_Stop; - mitk::PointSet::Pointer m_SamplingPointset; - mitk::PointSet::Pointer m_StopVotePointset; - mitk::PointSet::Pointer m_AlternativePointset; - - void SetStepSize(float v) ///< Integration step size in voxels, default is 0.5 * voxel - { m_StepSizeVox = v; } - void SetAngularThreshold(float v) ///< Angular threshold per step (in degree), default is 90deg x stepsize - { m_AngularThresholdDeg = v; } - void SetSamplingDistance(float v) ///< Maximum distance of sampling points in voxels, default is 0.25 * voxel - { m_SamplingDistanceVox = v; } - - itkGetMacro( OutputProbabilityMap, ItkDoubleImgType::Pointer) ///< Output probability map - itkGetMacro( FiberPolyData, PolyDataType ) ///< Output fibers - itkGetMacro( UseOutputProbabilityMap, bool) - - itkSetMacro( SeedImage, ItkUcharImgType::Pointer) ///< Seeds are only placed inside of this mask. - itkSetMacro( MaskImage, ItkUcharImgType::Pointer) ///< Tracking is only performed inside of this mask image. - itkSetMacro( TissueImage, ItkUcharImgType::Pointer) ///< - itkSetMacro( SeedsPerVoxel, int) ///< One seed placed in the center of each voxel or multiple seeds randomly placed inside each voxel. - itkSetMacro( MinTractLength, float ) ///< Shorter tracts are discarded. - itkSetMacro( MaxTractLength, float ) ///< Streamline progression stops if tract is longer than specified. - - itkSetMacro( UseStopVotes, bool ) ///< Frontal sampling points can vote for stopping the streamline even if the remaining sampling points keep pushing - itkSetMacro( OnlyForwardSamples, bool ) ///< Don't use sampling points behind the current position in progression direction - itkSetMacro( DeflectionMod, float ) ///< Deflection distance modifier - itkSetMacro( StoppingRegions, ItkUcharImgType::Pointer) ///< Streamlines entering a stopping region will stop immediately - itkSetMacro( DemoMode, bool ) - itkSetMacro( SeedOnlyGm, bool ) ///< place seed points only in the gray matter - itkSetMacro( ControlGmEndings, bool ) ///< - itkSetMacro( NumberOfSamples, unsigned int ) ///< Number of neighborhood sampling points - itkSetMacro( AposterioriCurvCheck, bool ) ///< Checks fiber curvature (angular deviation across 5mm) is larger than 30°. If yes, the streamline progression is stopped. - itkSetMacro( AvoidStop, bool ) ///< Use additional sampling points to avoid premature streamline termination - itkSetMacro( RandomSampling, bool ) ///< If true, the sampling points are distributed randomly around the current position, not sphericall in the specified sampling distance. - itkSetMacro( NumPreviousDirections, unsigned int ) ///< How many "old" steps do we want to consider in our decision where to go next? - itkSetMacro( MaxNumTracts, unsigned int ) ///< Tracking is stopped if the maximum number of tracts is exceeded - itkSetMacro( Random, bool ) ///< If true, seedpoints are shuffled randomly before tracking - itkSetMacro( Verbose, bool ) ///< If true, output tracking progress (might be slower) - itkSetMacro( UseOutputProbabilityMap, bool) ///< If true, no tractogram but a probability map is created as output. - itkSetMacro( SeedPoints, std::vector< itk::Point >) ///< Use manually defined points in physical space as seed points instead of seed image - - void SetTrackingHandler( mitk::TrackingDataHandler* h ) ///< - { - m_TrackingHandler = h; - } - - virtual void Update() override{ - this->GenerateData(); - } + typedef StreamlineTrackingFilter Self; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + typedef ProcessObject Superclass; + + /** Method for creation through the object factory. */ + itkFactorylessNewMacro(Self) + itkCloneMacro(Self) + + /** Runtime information support. */ + itkTypeMacro(MLBSTrackingFilter, ImageToImageFilter) + + typedef itk::Image ItkUcharImgType; + typedef itk::Image ItkDoubleImgType; + typedef itk::Image ItkFloatImgType; + typedef vtkSmartPointer< vtkPolyData > PolyDataType; + + typedef std::deque< itk::Point > FiberType; + typedef std::vector< FiberType > BundleType; + + volatile bool m_PauseTracking; + bool m_AbortTracking; + bool m_BuildFibersFinished; + int m_BuildFibersReady; + volatile bool m_Stop; + mitk::PointSet::Pointer m_SamplingPointset; + mitk::PointSet::Pointer m_StopVotePointset; + mitk::PointSet::Pointer m_AlternativePointset; + + void SetStepSize(float v) ///< Integration step size in voxels, default is 0.5 * voxel + { m_StepSizeVox = v; } + void SetAngularThreshold(float v) ///< Angular threshold per step (in degree), default is 90deg x stepsize + { m_AngularThresholdDeg = v; } + void SetSamplingDistance(float v) ///< Maximum distance of sampling points in voxels, default is 0.25 * voxel + { m_SamplingDistanceVox = v; } + + itkGetMacro( OutputProbabilityMap, ItkDoubleImgType::Pointer) ///< Output probability map + itkGetMacro( FiberPolyData, PolyDataType ) ///< Output fibers + itkGetMacro( UseOutputProbabilityMap, bool) + + itkSetMacro( SeedImage, ItkUcharImgType::Pointer) ///< Seeds are only placed inside of this mask. + itkSetMacro( MaskImage, ItkUcharImgType::Pointer) ///< Tracking is only performed inside of this mask image. + itkSetMacro( TissueImage, ItkUcharImgType::Pointer) ///< + itkSetMacro( SeedsPerVoxel, int) ///< One seed placed in the center of each voxel or multiple seeds randomly placed inside each voxel. + itkSetMacro( MinTractLength, float ) ///< Shorter tracts are discarded. + itkSetMacro( MaxTractLength, float ) ///< Streamline progression stops if tract is longer than specified. + + itkSetMacro( UseStopVotes, bool ) ///< Frontal sampling points can vote for stopping the streamline even if the remaining sampling points keep pushing + itkSetMacro( OnlyForwardSamples, bool ) ///< Don't use sampling points behind the current position in progression direction + itkSetMacro( DeflectionMod, float ) ///< Deflection distance modifier + itkSetMacro( StoppingRegions, ItkUcharImgType::Pointer) ///< Streamlines entering a stopping region will stop immediately + itkSetMacro( DemoMode, bool ) + itkSetMacro( SeedOnlyGm, bool ) ///< place seed points only in the gray matter + itkSetMacro( ControlGmEndings, bool ) ///< + itkSetMacro( NumberOfSamples, unsigned int ) ///< Number of neighborhood sampling points + itkSetMacro( AposterioriCurvCheck, bool ) ///< Checks fiber curvature (angular deviation across 5mm) is larger than 30°. If yes, the streamline progression is stopped. + itkSetMacro( AvoidStop, bool ) ///< Use additional sampling points to avoid premature streamline termination + itkSetMacro( RandomSampling, bool ) ///< If true, the sampling points are distributed randomly around the current position, not sphericall in the specified sampling distance. + itkSetMacro( NumPreviousDirections, unsigned int ) ///< How many "old" steps do we want to consider in our decision where to go next? + itkSetMacro( MaxNumTracts, unsigned int ) ///< Tracking is stopped if the maximum number of tracts is exceeded + itkSetMacro( Random, bool ) ///< If true, seedpoints are shuffled randomly before tracking + itkSetMacro( Verbose, bool ) ///< If true, output tracking progress (might be slower) + itkSetMacro( UseOutputProbabilityMap, bool) ///< If true, no tractogram but a probability map is created as output. + itkSetMacro( SeedPoints, std::vector< itk::Point >) ///< Use manually defined points in physical space as seed points instead of seed image + + void SetTrackingHandler( mitk::TrackingDataHandler* h ) ///< + { + m_TrackingHandler = h; + } + + virtual void Update() override{ + this->GenerateData(); + } protected: - void GenerateData() override; - - StreamlineTrackingFilter(); - ~StreamlineTrackingFilter() {} - - void InitGrayMatterEndings(); - void CheckFiberForGmEnding(FiberType* fib); - void FiberToProbmap(FiberType* fib); - - void GetSeedPointsFromSeedImage(); - void CalculateNewPosition(itk::Point& pos, vnl_vector_fixed& dir); ///< Calculate next integration step. - float FollowStreamline(itk::Point pos, vnl_vector_fixed dir, FiberType* fib, float tractLength, bool front); ///< Start streamline in one direction. - bool IsValidPosition(itk::Point& pos); ///< Are we outside of the mask image? - bool IsInGm(itk::Point &pos); - vnl_vector_fixed GetNewDirection(itk::Point& pos, std::deque< vnl_vector_fixed >& olddirs, itk::Index<3>& oldIndex); ///< Determine new direction by sample voting at the current position taking the last progression direction into account. - - float GetRandDouble(float min=-1, float max=1); - std::vector< vnl_vector_fixed > CreateDirections(int NPoints); - - void BeforeTracking(); - void AfterTracking(); - - PolyDataType m_FiberPolyData; - vtkSmartPointer m_Points; - vtkSmartPointer m_Cells; - BundleType m_Tractogram; - BundleType m_GmStubs; - - float m_AngularThresholdDeg; - float m_StepSizeVox; - float m_SamplingDistanceVox; - - float m_AngularThreshold; - float m_StepSize; - int m_MaxLength; - float m_MinTractLength; - float m_MaxTractLength; - int m_SeedsPerVoxel; - bool m_RandomSampling; - float m_SamplingDistance; - float m_DeflectionMod; - bool m_OnlyForwardSamples; - bool m_UseStopVotes; - unsigned int m_NumberOfSamples; - unsigned int m_NumPreviousDirections; - int m_WmLabel; - int m_GmLabel; - bool m_SeedOnlyGm; - bool m_ControlGmEndings; - int m_MaxNumTracts; - - ItkUcharImgType::Pointer m_StoppingRegions; - ItkUcharImgType::Pointer m_SeedImage; - ItkUcharImgType::Pointer m_MaskImage; - ItkUcharImgType::Pointer m_TissueImage; - ItkDoubleImgType::Pointer m_OutputProbabilityMap; - - bool m_Verbose; - bool m_AposterioriCurvCheck; - bool m_AvoidStop; - bool m_DemoMode; - bool m_Random; - bool m_UseOutputProbabilityMap; - std::vector< itk::Point > m_SeedPoints; - - void BuildFibers(bool check); - int CheckCurvature(FiberType* fib, bool front); - - // decision forest - mitk::TrackingDataHandler* m_TrackingHandler; - std::vector< PolyDataType > m_PolyDataContainer; - - std::chrono::time_point m_StartTime; - std::chrono::time_point m_EndTime; + void GenerateData() override; + + StreamlineTrackingFilter(); + ~StreamlineTrackingFilter() {} + + void InitGrayMatterEndings(); + void CheckFiberForGmEnding(FiberType* fib); + void FiberToProbmap(FiberType* fib); + + void GetSeedPointsFromSeedImage(); + void CalculateNewPosition(itk::Point& pos, vnl_vector_fixed& dir); ///< Calculate next integration step. + float FollowStreamline(itk::Point start_pos, vnl_vector_fixed dir, FiberType* fib, float tractLength, bool front); ///< Start streamline in one direction. + bool IsValidPosition(const itk::Point& pos); ///< Are we outside of the mask image? + bool IsInGm(const itk::Point &pos); + vnl_vector_fixed GetNewDirection(itk::Point& pos, std::deque< vnl_vector_fixed >& olddirs, itk::Index<3>& oldIndex); ///< Determine new direction by sample voting at the current position taking the last progression direction into account. + + float GetRandDouble(float min=-1, float max=1); + std::vector< vnl_vector_fixed > CreateDirections(int NPoints); + + void BeforeTracking(); + void AfterTracking(); + + PolyDataType m_FiberPolyData; + vtkSmartPointer m_Points; + vtkSmartPointer m_Cells; + BundleType m_Tractogram; + BundleType m_GmStubs; + + float m_AngularThresholdDeg; + float m_StepSizeVox; + float m_SamplingDistanceVox; + + float m_AngularThreshold; + float m_StepSize; + int m_MaxLength; + float m_MinTractLength; + float m_MaxTractLength; + int m_SeedsPerVoxel; + bool m_RandomSampling; + float m_SamplingDistance; + float m_DeflectionMod; + bool m_OnlyForwardSamples; + bool m_UseStopVotes; + unsigned int m_NumberOfSamples; + unsigned int m_NumPreviousDirections; + int m_WmLabel; + int m_GmLabel; + bool m_SeedOnlyGm; + bool m_ControlGmEndings; + int m_MaxNumTracts; + + ItkUcharImgType::Pointer m_StoppingRegions; + ItkUcharImgType::Pointer m_SeedImage; + ItkUcharImgType::Pointer m_MaskImage; + ItkUcharImgType::Pointer m_TissueImage; + ItkDoubleImgType::Pointer m_OutputProbabilityMap; + + bool m_Verbose; + bool m_AposterioriCurvCheck; + bool m_AvoidStop; + bool m_DemoMode; + bool m_Random; + bool m_UseOutputProbabilityMap; + std::vector< itk::Point > m_SeedPoints; + + void BuildFibers(bool check); + int CheckCurvature(FiberType* fib, bool front); + + // decision forest + mitk::TrackingDataHandler* m_TrackingHandler; + std::vector< PolyDataType > m_PolyDataContainer; + + std::chrono::time_point m_StartTime; + std::chrono::time_point m_EndTime; private: }; } //#ifndef ITK_MANUAL_INSTANTIATION //#include "itkMLBSTrackingFilter.cpp" //#endif #endif //__itkMLBSTrackingFilter_h_ diff --git a/Modules/DiffusionImaging/cmdapps/StreamlineTracking.cpp b/Modules/DiffusionImaging/cmdapps/StreamlineTracking.cpp index 5da594fb9b..b7c6b752a7 100755 --- a/Modules/DiffusionImaging/cmdapps/StreamlineTracking.cpp +++ b/Modules/DiffusionImaging/cmdapps/StreamlineTracking.cpp @@ -1,431 +1,431 @@ /*=================================================================== 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 #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include using namespace std; const int numOdfSamples = 200; typedef itk::Image< itk::Vector< float, numOdfSamples > , 3 > SampledShImageType; /*! \brief */ int main(int argc, char* argv[]) { - mitkCommandLineParser parser; - - parser.setTitle("Streamline Tractography"); - parser.setCategory("Fiber Tracking and Processing Methods"); - parser.setDescription("Perform streamline tractography"); - parser.setContributor("MIC"); - - // parameters fo all methods - parser.setArgumentPrefix("--", "-"); - parser.addArgument("input", "i", mitkCommandLineParser::StringList, "Input:", "input image (multiple possible for 'Tensor' algorithm)", us::Any(), false); - parser.addArgument("algorithm", "a", mitkCommandLineParser::String, "Algorithm:", "which algorithm to use (Peaks, Tensor, DetODF, ProbODF, DetRF, ProbRF)", us::Any(), false); - parser.addArgument("out", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output fiberbundle/probability map", us::Any(), false); - - parser.addArgument("stop_mask", "", mitkCommandLineParser::String, "Stop image:", "streamlines entering the binary mask will stop immediately", us::Any()); - parser.addArgument("tracking_mask", "", mitkCommandLineParser::String, "Mask image:", "restrict tractography with a binary mask image", us::Any()); - parser.addArgument("seed_mask", "", mitkCommandLineParser::String, "Seed image:", "binary mask image defining seed voxels", us::Any()); - parser.addArgument("tissue_image", "", mitkCommandLineParser::String, "Tissue type image:", "image with tissue type labels (WM=3, GM=1)", us::Any()); - - parser.addArgument("cutoff", "", mitkCommandLineParser::Float, "Cutoff:", "set the FA, GFA or Peak amplitude cutoff for terminating tracks", 0.1); - parser.addArgument("step_size", "", mitkCommandLineParser::Float, "Step size:", "step size (in voxels)", 0.5); - parser.addArgument("angular_threshold", "", mitkCommandLineParser::Float, "Angular threshold:", "angular threshold between two successive steps, (default: 90° * step_size)"); - parser.addArgument("seeds", "", mitkCommandLineParser::Int, "Seeds per voxel:", "number of seed points per voxel", 1); - parser.addArgument("seed_gm", "", mitkCommandLineParser::Bool, "Seed only GM:", "Seed only in gray matter (requires tissue type image --tissue_image)"); - parser.addArgument("control_gm_endings", "", mitkCommandLineParser::Bool, "Control GM endings:", "Seed perpendicular to gray matter and enforce endings inside of the gray matter (requires tissue type image --tissue_image)"); - parser.addArgument("max_tracts", "", mitkCommandLineParser::Int, "Max. number of tracts:", "tractography is stopped if the reconstructed number of tracts is exceeded.", -1); - - parser.addArgument("num_samples", "", mitkCommandLineParser::Int, "Num. neighborhood samples:", "number of neighborhood samples that are use to determine the next progression direction", 0); - parser.addArgument("sampling_distance", "", mitkCommandLineParser::Float, "Sampling distance:", "distance of neighborhood sampling points (in voxels)", 0.25); - parser.addArgument("use_stop_votes", "", mitkCommandLineParser::Bool, "Use stop votes:", "use stop votes"); - parser.addArgument("use_only_forward_samples", "", mitkCommandLineParser::Bool, "Use only forward samples:", "use only forward samples"); - parser.addArgument("output_prob_map", "", mitkCommandLineParser::Bool, "Output probability map:", "output probability map instead of tractogram"); - - parser.addArgument("no_interpolation", "", mitkCommandLineParser::Bool, "Don't interpolate:", "don't interpolate image values"); - parser.addArgument("flip_x", "", mitkCommandLineParser::Bool, "Flip X:", "multiply x-coordinate of direction proposal by -1"); - parser.addArgument("flip_y", "", mitkCommandLineParser::Bool, "Flip Y:", "multiply y-coordinate of direction proposal by -1"); - parser.addArgument("flip_z", "", mitkCommandLineParser::Bool, "Flip Z:", "multiply z-coordinate of direction proposal by -1"); - //parser.addArgument("apply_image_rotation", "", mitkCommandLineParser::Bool, "Apply image rotation:", "applies image rotation to image peaks (only for 'Peaks' algorithm). This is necessary in some cases, e.g. if the peaks were obtained with MRtrix."); - - parser.addArgument("compress", "", mitkCommandLineParser::Float, "Compress:", "Compress output fibers using the given error threshold (in mm)"); - parser.addArgument("additional_images", "", mitkCommandLineParser::StringList, "Additional images:", "specify a list of float images that hold additional information (FA, GFA, additional Features)", us::Any()); - - // parameters for ODF based tractography - parser.addArgument("no_odf_normalization", "", mitkCommandLineParser::Bool, "Don't min-max normalize ODFs:", "No min-max normalization of ODFs"); - - // parameters for random forest based tractography - parser.addArgument("forest", "", mitkCommandLineParser::String, "Forest:", "input random forest (HDF5 file)", us::Any()); - parser.addArgument("use_sh_features", "", mitkCommandLineParser::Bool, "Use SH features:", "use SH features"); - - // parameters for tensor tractography - parser.addArgument("tend_f", "", mitkCommandLineParser::Float, "Weight f", "Weighting factor between first eigenvector (f=1 equals FACT tracking) and input vector dependent direction (f=0).", 1.0); - parser.addArgument("tend_g", "", mitkCommandLineParser::Float, "Weight g", "Weighting factor between input vector (g=0) and tensor deflection (g=1 equals TEND tracking)", 0.0); - - - map parsedArgs = parser.parseArguments(argc, argv); - if (parsedArgs.size()==0) - return EXIT_FAILURE; - - mitkCommandLineParser::StringContainerType input_files = us::any_cast(parsedArgs["input"]); - string outFile = us::any_cast(parsedArgs["out"]); - string algorithm = us::any_cast(parsedArgs["algorithm"]); - -// mitk::LoggingBackend::SetLogFile( (outFile + ".log").c_str() ); -// MITK_INFO << "LOG"; - - bool interpolate = true; - if (parsedArgs.count("no_interpolation")) - interpolate = !us::any_cast(parsedArgs["no_interpolation"]); - - bool use_sh_features = false; - if (parsedArgs.count("use_sh_features")) - use_sh_features = us::any_cast(parsedArgs["use_sh_features"]); - - bool seed_gm = false; - if (parsedArgs.count("seed_gm")) - seed_gm = us::any_cast(parsedArgs["seed_gm"]); - - bool control_gm_endings = false; - if (parsedArgs.count("control_gm_endings")) - control_gm_endings = us::any_cast(parsedArgs["control_gm_endings"]); - - bool use_stop_votes = false; - if (parsedArgs.count("use_stop_votes")) - use_stop_votes = us::any_cast(parsedArgs["use_stop_votes"]); - - bool use_only_forward_samples = false; - if (parsedArgs.count("use_only_forward_samples")) - use_only_forward_samples = us::any_cast(parsedArgs["use_only_forward_samples"]); - - bool output_prob_map = false; - if (parsedArgs.count("output_prob_map")) - output_prob_map = us::any_cast(parsedArgs["output_prob_map"]); - - bool flip_x = false; - if (parsedArgs.count("flip_x")) - flip_x = us::any_cast(parsedArgs["flip_x"]); - bool flip_y = false; - if (parsedArgs.count("flip_y")) - flip_y = us::any_cast(parsedArgs["flip_y"]); - bool flip_z = false; - if (parsedArgs.count("flip_z")) - flip_z = us::any_cast(parsedArgs["flip_z"]); - - bool apply_image_rotation = false; - if (parsedArgs.count("apply_image_rotation")) - apply_image_rotation = us::any_cast(parsedArgs["apply_image_rotation"]); - - float compress = -1; - if (parsedArgs.count("compress")) - compress = us::any_cast(parsedArgs["compress"]); - - string forestFile; - if (parsedArgs.count("forest")) - forestFile = us::any_cast(parsedArgs["forest"]); - - string maskFile = ""; - if (parsedArgs.count("tracking_mask")) - maskFile = us::any_cast(parsedArgs["tracking_mask"]); - - string seedFile = ""; - if (parsedArgs.count("seed_mask")) - seedFile = us::any_cast(parsedArgs["seed_mask"]); - - string stopFile = ""; - if (parsedArgs.count("stop_mask")) - stopFile = us::any_cast(parsedArgs["stop_mask"]); - - string tissueFile = ""; - if (parsedArgs.count("tissue_image")) - tissueFile = us::any_cast(parsedArgs["tissue_image"]); - - float cutoff = 0.1; - if (parsedArgs.count("cutoff")) - cutoff = us::any_cast(parsedArgs["cutoff"]); - - float stepsize = -1; - if (parsedArgs.count("step_size")) - stepsize = us::any_cast(parsedArgs["step_size"]); - - float sampling_distance = -1; - if (parsedArgs.count("sampling_distance")) - sampling_distance = us::any_cast(parsedArgs["sampling_distance"]); - - int num_samples = 0; - if (parsedArgs.count("num_samples")) - num_samples = us::any_cast(parsedArgs["num_samples"]); - - int seeds = 1; - if (parsedArgs.count("seeds")) - seeds = us::any_cast(parsedArgs["seeds"]); - - float tend_f = 1; - if (parsedArgs.count("tend_f")) - tend_f = us::any_cast(parsedArgs["tend_f"]); - - float tend_g = 0; - if (parsedArgs.count("tend_g")) - tend_g = us::any_cast(parsedArgs["tend_g"]); - - float angular_threshold = -1; - if (parsedArgs.count("angular_threshold")) - angular_threshold = us::any_cast(parsedArgs["angular_threshold"]); - - bool normalize_odfs = true; - if (parsedArgs.count("no_odf_normalization")) - normalize_odfs = !us::any_cast(parsedArgs["no_odf_normalization"]); - - unsigned int max_tracts = -1; - if (parsedArgs.count("max_tracts")) - max_tracts = us::any_cast(parsedArgs["max_tracts"]); - - - // LOAD DATASETS - - mitkCommandLineParser::StringContainerType addFiles; - if (parsedArgs.count("additional_images")) - addFiles = us::any_cast(parsedArgs["additional_images"]); - - typedef itk::Image ItkUcharImgType; - - MITK_INFO << "loading input"; - std::vector< mitk::Image::Pointer > input_images; - for (unsigned int i=0; i parsedArgs = parser.parseArguments(argc, argv); + if (parsedArgs.size()==0) + return EXIT_FAILURE; + + mitkCommandLineParser::StringContainerType input_files = us::any_cast(parsedArgs["input"]); + string outFile = us::any_cast(parsedArgs["out"]); + string algorithm = us::any_cast(parsedArgs["algorithm"]); + + // mitk::LoggingBackend::SetLogFile( (outFile + ".log").c_str() ); + // MITK_INFO << "LOG"; + + bool interpolate = true; + if (parsedArgs.count("no_interpolation")) + interpolate = !us::any_cast(parsedArgs["no_interpolation"]); + + bool use_sh_features = false; + if (parsedArgs.count("use_sh_features")) + use_sh_features = us::any_cast(parsedArgs["use_sh_features"]); + + bool seed_gm = false; + if (parsedArgs.count("seed_gm")) + seed_gm = us::any_cast(parsedArgs["seed_gm"]); + + bool control_gm_endings = false; + if (parsedArgs.count("control_gm_endings")) + control_gm_endings = us::any_cast(parsedArgs["control_gm_endings"]); + + bool use_stop_votes = false; + if (parsedArgs.count("use_stop_votes")) + use_stop_votes = us::any_cast(parsedArgs["use_stop_votes"]); + + bool use_only_forward_samples = false; + if (parsedArgs.count("use_only_forward_samples")) + use_only_forward_samples = us::any_cast(parsedArgs["use_only_forward_samples"]); + + bool output_prob_map = false; + if (parsedArgs.count("output_prob_map")) + output_prob_map = us::any_cast(parsedArgs["output_prob_map"]); + + bool flip_x = false; + if (parsedArgs.count("flip_x")) + flip_x = us::any_cast(parsedArgs["flip_x"]); + bool flip_y = false; + if (parsedArgs.count("flip_y")) + flip_y = us::any_cast(parsedArgs["flip_y"]); + bool flip_z = false; + if (parsedArgs.count("flip_z")) + flip_z = us::any_cast(parsedArgs["flip_z"]); + + bool apply_image_rotation = false; + if (parsedArgs.count("apply_image_rotation")) + apply_image_rotation = us::any_cast(parsedArgs["apply_image_rotation"]); + + float compress = -1; + if (parsedArgs.count("compress")) + compress = us::any_cast(parsedArgs["compress"]); + + string forestFile; + if (parsedArgs.count("forest")) + forestFile = us::any_cast(parsedArgs["forest"]); + + string maskFile = ""; + if (parsedArgs.count("tracking_mask")) + maskFile = us::any_cast(parsedArgs["tracking_mask"]); + + string seedFile = ""; + if (parsedArgs.count("seed_mask")) + seedFile = us::any_cast(parsedArgs["seed_mask"]); + + string stopFile = ""; + if (parsedArgs.count("stop_mask")) + stopFile = us::any_cast(parsedArgs["stop_mask"]); + + string tissueFile = ""; + if (parsedArgs.count("tissue_image")) + tissueFile = us::any_cast(parsedArgs["tissue_image"]); + + float cutoff = 0.1; + if (parsedArgs.count("cutoff")) + cutoff = us::any_cast(parsedArgs["cutoff"]); + + float stepsize = -1; + if (parsedArgs.count("step_size")) + stepsize = us::any_cast(parsedArgs["step_size"]); + + float sampling_distance = -1; + if (parsedArgs.count("sampling_distance")) + sampling_distance = us::any_cast(parsedArgs["sampling_distance"]); + + int num_samples = 0; + if (parsedArgs.count("num_samples")) + num_samples = us::any_cast(parsedArgs["num_samples"]); + + int seeds = 1; + if (parsedArgs.count("seeds")) + seeds = us::any_cast(parsedArgs["seeds"]); + + float tend_f = 1; + if (parsedArgs.count("tend_f")) + tend_f = us::any_cast(parsedArgs["tend_f"]); + + float tend_g = 0; + if (parsedArgs.count("tend_g")) + tend_g = us::any_cast(parsedArgs["tend_g"]); + + float angular_threshold = -1; + if (parsedArgs.count("angular_threshold")) + angular_threshold = us::any_cast(parsedArgs["angular_threshold"]); + + bool normalize_odfs = true; + if (parsedArgs.count("no_odf_normalization")) + normalize_odfs = !us::any_cast(parsedArgs["no_odf_normalization"]); + + unsigned int max_tracts = -1; + if (parsedArgs.count("max_tracts")) + max_tracts = us::any_cast(parsedArgs["max_tracts"]); + + + // LOAD DATASETS + + mitkCommandLineParser::StringContainerType addFiles; + if (parsedArgs.count("additional_images")) + addFiles = us::any_cast(parsedArgs["additional_images"]); + + typedef itk::Image ItkUcharImgType; + + MITK_INFO << "loading input"; + std::vector< mitk::Image::Pointer > input_images; + for (unsigned int i=0; i(mitk::IOUtil::LoadDataNode(input_files.at(i))->GetData()); + input_images.push_back(mitkImage); + } + + ItkUcharImgType::Pointer mask; + if (!maskFile.empty()) + { + MITK_INFO << "loading mask image"; + mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::LoadImage(maskFile).GetPointer()); + mask = ItkUcharImgType::New(); + mitk::CastToItkImage(img, mask); + } + + ItkUcharImgType::Pointer seed; + if (!seedFile.empty()) + { + MITK_INFO << "loading seed image"; + mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::LoadImage(seedFile).GetPointer()); + seed = ItkUcharImgType::New(); + mitk::CastToItkImage(img, seed); + } + + ItkUcharImgType::Pointer stop; + if (!stopFile.empty()) + { + MITK_INFO << "loading stop image"; + mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::LoadImage(stopFile).GetPointer()); + stop = ItkUcharImgType::New(); + mitk::CastToItkImage(img, stop); + } + + + ItkUcharImgType::Pointer tissue; + if (!tissueFile.empty()) + { + MITK_INFO << "loading tissue image"; + mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::LoadImage(tissueFile).GetPointer()); + tissue = ItkUcharImgType::New(); + mitk::CastToItkImage(img, tissue); + } + + MITK_INFO << "loading additional images"; + typedef itk::Image ItkFloatImgType; + std::vector< std::vector< ItkFloatImgType::Pointer > > addImages; + addImages.push_back(std::vector< ItkFloatImgType::Pointer >()); + for (auto file : addFiles) + { + mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::LoadImage(file).GetPointer()); + ItkFloatImgType::Pointer itkimg = ItkFloatImgType::New(); + mitk::CastToItkImage(img, itkimg); + addImages.at(0).push_back(itkimg); + } + + // ////////////////////////////////////////////////////////////////// + // omp_set_num_threads(1); + + typedef itk::StreamlineTrackingFilter TrackerType; + TrackerType::Pointer tracker = TrackerType::New(); + + mitk::TrackingDataHandler* handler; + if (algorithm == "DetRF" || algorithm == "ProbRF") + { + if (use_sh_features) { - mitk::Image::Pointer mitkImage = dynamic_cast(mitk::IOUtil::LoadDataNode(input_files.at(i))->GetData()); - input_images.push_back(mitkImage); + handler = new mitk::TrackingHandlerRandomForest<6,28>(); + dynamic_cast*>(handler)->LoadForest(forestFile); + dynamic_cast*>(handler)->AddDwi(input_images.at(0)); + dynamic_cast*>(handler)->SetAdditionalFeatureImages(addImages); } - - ItkUcharImgType::Pointer mask; - if (!maskFile.empty()) - { - MITK_INFO << "loading mask image"; - mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::LoadImage(maskFile).GetPointer()); - mask = ItkUcharImgType::New(); - mitk::CastToItkImage(img, mask); - } - - ItkUcharImgType::Pointer seed; - if (!seedFile.empty()) - { - MITK_INFO << "loading seed image"; - mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::LoadImage(seedFile).GetPointer()); - seed = ItkUcharImgType::New(); - mitk::CastToItkImage(img, seed); - } - - ItkUcharImgType::Pointer stop; - if (!stopFile.empty()) - { - MITK_INFO << "loading stop image"; - mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::LoadImage(stopFile).GetPointer()); - stop = ItkUcharImgType::New(); - mitk::CastToItkImage(img, stop); - } - - - ItkUcharImgType::Pointer tissue; - if (!tissueFile.empty()) + else { - MITK_INFO << "loading tissue image"; - mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::LoadImage(tissueFile).GetPointer()); - tissue = ItkUcharImgType::New(); - mitk::CastToItkImage(img, tissue); + handler = new mitk::TrackingHandlerRandomForest<6,100>(); + dynamic_cast*>(handler)->LoadForest(forestFile); + dynamic_cast*>(handler)->AddDwi(input_images.at(0)); + dynamic_cast*>(handler)->SetAdditionalFeatureImages(addImages); } - MITK_INFO << "loading additional images"; - typedef itk::Image ItkFloatImgType; - std::vector< std::vector< ItkFloatImgType::Pointer > > addImages; - addImages.push_back(std::vector< ItkFloatImgType::Pointer >()); - for (auto file : addFiles) + if (algorithm == "ProbRF") + handler->SetMode(mitk::TrackingDataHandler::MODE::PROBABILISTIC); + } + else if (algorithm == "Peaks") + { + handler = new mitk::TrackingHandlerPeaks(); + + typedef mitk::ImageToItk< mitk::TrackingHandlerPeaks::PeakImgType > CasterType; + CasterType::Pointer caster = CasterType::New(); + caster->SetInput(input_images.at(0)); + caster->Update(); + mitk::TrackingHandlerPeaks::PeakImgType::Pointer itkImg = caster->GetOutput(); + + dynamic_cast(handler)->SetPeakImage(itkImg); + dynamic_cast(handler)->SetApplyDirectionMatrix(apply_image_rotation); + dynamic_cast(handler)->SetPeakThreshold(cutoff); + } + else if (algorithm == "Tensor") + { + handler = new mitk::TrackingHandlerTensor(); + + for (auto input_image : input_images) { - mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::LoadImage(file).GetPointer()); - ItkFloatImgType::Pointer itkimg = ItkFloatImgType::New(); - mitk::CastToItkImage(img, itkimg); - addImages.at(0).push_back(itkimg); + typedef mitk::ImageToItk< mitk::TrackingHandlerTensor::ItkTensorImageType > CasterType; + CasterType::Pointer caster = CasterType::New(); + caster->SetInput(input_image); + caster->Update(); + mitk::TrackingHandlerTensor::ItkTensorImageType::Pointer itkImg = caster->GetOutput(); + dynamic_cast(handler)->AddTensorImage(itkImg); } -// ////////////////////////////////////////////////////////////////// -// omp_set_num_threads(1); + dynamic_cast(handler)->SetFaThreshold(cutoff); + dynamic_cast(handler)->SetF(tend_f); + dynamic_cast(handler)->SetG(tend_g); - typedef itk::StreamlineTrackingFilter TrackerType; - TrackerType::Pointer tracker = TrackerType::New(); + if (addImages.size()>0) + dynamic_cast(handler)->SetFaImage(addImages.at(0).at(0)); + } + else if (algorithm == "DetODF" || algorithm == "ProbODF") + { + handler = new mitk::TrackingHandlerOdf(); - mitk::TrackingDataHandler* handler; - if (algorithm == "DetRF" || algorithm == "ProbRF") + for (auto input_image : input_images) { - if (use_sh_features) - { - handler = new mitk::TrackingHandlerRandomForest<6,28>(); - dynamic_cast*>(handler)->LoadForest(forestFile); - dynamic_cast*>(handler)->AddDwi(input_images.at(0)); - dynamic_cast*>(handler)->SetAdditionalFeatureImages(addImages); - } - else - { - handler = new mitk::TrackingHandlerRandomForest<6,100>(); - dynamic_cast*>(handler)->LoadForest(forestFile); - dynamic_cast*>(handler)->AddDwi(input_images.at(0)); - dynamic_cast*>(handler)->SetAdditionalFeatureImages(addImages); - } - - if (algorithm == "ProbRF") - handler->SetMode(mitk::TrackingDataHandler::MODE::PROBABILISTIC); + typedef mitk::ImageToItk< mitk::TrackingHandlerOdf::ItkOdfImageType > CasterType; + CasterType::Pointer caster = CasterType::New(); + caster->SetInput(input_image); + caster->Update(); + mitk::TrackingHandlerOdf::ItkOdfImageType::Pointer itkImg = caster->GetOutput(); + dynamic_cast(handler)->SetOdfImage(itkImg); } - else if (algorithm == "Peaks") - { - handler = new mitk::TrackingHandlerPeaks(); - typedef mitk::ImageToItk< mitk::TrackingHandlerPeaks::PeakImgType > CasterType; - CasterType::Pointer caster = CasterType::New(); - caster->SetInput(input_images.at(0)); - caster->Update(); - mitk::TrackingHandlerPeaks::PeakImgType::Pointer itkImg = caster->GetOutput(); - - dynamic_cast(handler)->SetPeakImage(itkImg); - dynamic_cast(handler)->SetApplyDirectionMatrix(apply_image_rotation); - dynamic_cast(handler)->SetPeakThreshold(cutoff); - } - else if (algorithm == "Tensor") - { - handler = new mitk::TrackingHandlerTensor(); - - for (auto input_image : input_images) - { - typedef mitk::ImageToItk< mitk::TrackingHandlerTensor::ItkTensorImageType > CasterType; - CasterType::Pointer caster = CasterType::New(); - caster->SetInput(input_image); - caster->Update(); - mitk::TrackingHandlerTensor::ItkTensorImageType::Pointer itkImg = caster->GetOutput(); - dynamic_cast(handler)->AddTensorImage(itkImg); - } - - dynamic_cast(handler)->SetFaThreshold(cutoff); - dynamic_cast(handler)->SetF(tend_f); - dynamic_cast(handler)->SetG(tend_g); - - if (addImages.size()>0) - dynamic_cast(handler)->SetFaImage(addImages.at(0).at(0)); - } - else if (algorithm == "DetODF" || algorithm == "ProbODF") - { - handler = new mitk::TrackingHandlerOdf(); - - for (auto input_image : input_images) - { - typedef mitk::ImageToItk< mitk::TrackingHandlerOdf::ItkOdfImageType > CasterType; - CasterType::Pointer caster = CasterType::New(); - caster->SetInput(input_image); - caster->Update(); - mitk::TrackingHandlerOdf::ItkOdfImageType::Pointer itkImg = caster->GetOutput(); - dynamic_cast(handler)->SetOdfImage(itkImg); - } - - dynamic_cast(handler)->SetGfaThreshold(cutoff); - if (algorithm == "ProbODF") - dynamic_cast(handler)->SetMode(mitk::TrackingHandlerOdf::MODE::PROBABILISTIC); - - if (addImages.size()>0) - dynamic_cast(handler)->SetGfaImage(addImages.at(0).at(0)); - - dynamic_cast(handler)->SetMinMaxNormalize(normalize_odfs); - } - handler->SetInterpolate(interpolate); - handler->SetFlipX(flip_x); - handler->SetFlipY(flip_y); - handler->SetFlipZ(flip_z); - - MITK_INFO << "Tractography algorithm: " << algorithm; - - tracker->SetNumberOfSamples(num_samples); - tracker->SetAngularThreshold(angular_threshold); - tracker->SetMaskImage(mask); - tracker->SetSeedImage(seed); - tracker->SetStoppingRegions(stop); - tracker->SetSeedsPerVoxel(seeds); - tracker->SetStepSize(stepsize); - tracker->SetSamplingDistance(sampling_distance); - tracker->SetUseStopVotes(use_stop_votes); - tracker->SetOnlyForwardSamples(use_only_forward_samples); - tracker->SetAposterioriCurvCheck(false); - tracker->SetTissueImage(tissue); - tracker->SetSeedOnlyGm(seed_gm); - tracker->SetControlGmEndings(control_gm_endings); - tracker->SetMaxNumTracts(max_tracts); - tracker->SetTrackingHandler(handler); - tracker->SetUseOutputProbabilityMap(output_prob_map); - tracker->Update(); - - std::string ext = itksys::SystemTools::GetFilenameExtension(outFile); - - if (!output_prob_map) - { - vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); - mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); - if (compress > 0) - outFib->Compress(compress); - - if (ext != ".fib" && ext != ".trk" && ext != ".tck") - outFile += ".fib"; - - mitk::IOUtil::Save(outFib, outFile); - } - else - { - TrackerType::ItkDoubleImgType::Pointer outImg = tracker->GetOutputProbabilityMap(); - mitk::Image::Pointer img = mitk::Image::New(); - img->InitializeByItk(outImg.GetPointer()); - img->SetVolume(outImg->GetBufferPointer()); - - if (ext != ".nii" && ext != ".nii.gz" && ext != ".nrrd") - outFile += ".nii.gz"; - - mitk::IOUtil::SaveBaseData(img, outFile); - } - - delete handler; - - return EXIT_SUCCESS; + dynamic_cast(handler)->SetGfaThreshold(cutoff); + if (algorithm == "ProbODF") + dynamic_cast(handler)->SetMode(mitk::TrackingHandlerOdf::MODE::PROBABILISTIC); + + if (addImages.size()>0) + dynamic_cast(handler)->SetGfaImage(addImages.at(0).at(0)); + + dynamic_cast(handler)->SetMinMaxNormalize(normalize_odfs); + } + handler->SetInterpolate(interpolate); + handler->SetFlipX(flip_x); + handler->SetFlipY(flip_y); + handler->SetFlipZ(flip_z); + + MITK_INFO << "Tractography algorithm: " << algorithm; + + tracker->SetNumberOfSamples(num_samples); + tracker->SetAngularThreshold(angular_threshold); + tracker->SetMaskImage(mask); + tracker->SetSeedImage(seed); + tracker->SetStoppingRegions(stop); + tracker->SetSeedsPerVoxel(seeds); + tracker->SetStepSize(stepsize); + tracker->SetSamplingDistance(sampling_distance); + tracker->SetUseStopVotes(use_stop_votes); + tracker->SetOnlyForwardSamples(use_only_forward_samples); + tracker->SetAposterioriCurvCheck(false); + tracker->SetTissueImage(tissue); + tracker->SetSeedOnlyGm(seed_gm); + tracker->SetControlGmEndings(control_gm_endings); + tracker->SetMaxNumTracts(max_tracts); + tracker->SetTrackingHandler(handler); + tracker->SetUseOutputProbabilityMap(output_prob_map); + tracker->Update(); + + std::string ext = itksys::SystemTools::GetFilenameExtension(outFile); + + if (!output_prob_map) + { + vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); + mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); + if (compress > 0) + outFib->Compress(compress); + + if (ext != ".fib" && ext != ".trk" && ext != ".tck") + outFile += ".fib"; + + mitk::IOUtil::Save(outFib, outFile); + } + else + { + TrackerType::ItkDoubleImgType::Pointer outImg = tracker->GetOutputProbabilityMap(); + mitk::Image::Pointer img = mitk::Image::New(); + img->InitializeByItk(outImg.GetPointer()); + img->SetVolume(outImg->GetBufferPointer()); + + if (ext != ".nii" && ext != ".nii.gz" && ext != ".nrrd") + outFile += ".nii.gz"; + + mitk::IOUtil::SaveBaseData(img, outFile); + } + + delete handler; + + return EXIT_SUCCESS; } 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 455762f504..1ba7a4f1f3 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,441 +1,557 @@ /*=================================================================== 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 // VTK +#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; QmitkStreamlineTrackingView::QmitkStreamlineTrackingView() : m_Controls(nullptr) { } // Destructor QmitkStreamlineTrackingView::~QmitkStreamlineTrackingView() { } void QmitkStreamlineTrackingView::CreateQtPartControl( QWidget *parent ) { - if ( !m_Controls ) + 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_StopImageBox->SetDataStorage(this->GetDataStorage()); + m_Controls->m_TissueImageBox->SetDataStorage(this->GetDataStorage()); + + mitk::TNodePredicateDataType::Pointer isImagePredicate = 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_FaImageBox->SetPredicate( mitk::NodePredicateAnd::New(isNotABinaryImagePredicate, dimensionPredicate) ); + m_Controls->m_FaImageBox->SetZeroEntryText("--"); + m_Controls->m_SeedImageBox->SetPredicate( mitk::NodePredicateAnd::New(isBinaryPredicate, dimensionPredicate) ); + m_Controls->m_SeedImageBox->SetZeroEntryText("--"); + m_Controls->m_MaskImageBox->SetPredicate( mitk::NodePredicateAnd::New(isBinaryPredicate, dimensionPredicate) ); + m_Controls->m_MaskImageBox->SetZeroEntryText("--"); + m_Controls->m_StopImageBox->SetPredicate( mitk::NodePredicateAnd::New(isBinaryPredicate, dimensionPredicate) ); + m_Controls->m_StopImageBox->SetZeroEntryText("--"); + m_Controls->m_TissueImageBox->SetPredicate( mitk::NodePredicateAnd::New(isNotABinaryImagePredicate, dimensionPredicate) ); + m_Controls->m_TissueImageBox->SetZeroEntryText("--"); + + connect( m_Controls->commandLinkButton, SIGNAL(clicked()), this, SLOT(DoFiberTracking()) ); + connect( m_Controls->m_InteractiveBox, SIGNAL(stateChanged(int)), this, SLOT(ToggleInteractive()) ); + connect( m_Controls->m_TissueImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui()) ); + connect( m_Controls->m_ModeBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui()) ); + + m_FirstTensorProbRun = true; + } + + UpdateGui(); +} + +void QmitkStreamlineTrackingView::ToggleInteractive() +{ + UpdateGui(); + + mitk::IRenderWindowPart* renderWindow = this->GetRenderWindowPart(); + + if ( m_Controls->m_InteractiveBox->isChecked() ) + { + m_InteractiveNode = mitk::DataNode::New(); + + QString name("InteractiveFib"); + m_InteractiveNode->SetName(name.toStdString()); + GetDataStorage()->Add(m_InteractiveNode); + + 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); + + if (renderWindow) { - // 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_StopImageBox->SetDataStorage(this->GetDataStorage()); - m_Controls->m_TissueImageBox->SetDataStorage(this->GetDataStorage()); - - mitk::TNodePredicateDataType::Pointer isImagePredicate = 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_FaImageBox->SetPredicate( mitk::NodePredicateAnd::New(isNotABinaryImagePredicate, dimensionPredicate) ); - m_Controls->m_FaImageBox->SetZeroEntryText("--"); - m_Controls->m_SeedImageBox->SetPredicate( mitk::NodePredicateAnd::New(isBinaryPredicate, dimensionPredicate) ); - m_Controls->m_SeedImageBox->SetZeroEntryText("--"); - m_Controls->m_MaskImageBox->SetPredicate( mitk::NodePredicateAnd::New(isBinaryPredicate, dimensionPredicate) ); - m_Controls->m_MaskImageBox->SetZeroEntryText("--"); - m_Controls->m_StopImageBox->SetPredicate( mitk::NodePredicateAnd::New(isBinaryPredicate, dimensionPredicate) ); - m_Controls->m_StopImageBox->SetZeroEntryText("--"); - m_Controls->m_TissueImageBox->SetPredicate( mitk::NodePredicateAnd::New(isNotABinaryImagePredicate, dimensionPredicate) ); - m_Controls->m_TissueImageBox->SetZeroEntryText("--"); - - connect( m_Controls->commandLinkButton, SIGNAL(clicked()), this, SLOT(DoFiberTracking()) ); - connect( m_Controls->m_TissueImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui()) ); - connect( m_Controls->m_ModeBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui()) ); - - m_FirstTensorProbRun = true; + { + mitk::SliceNavigationController* slicer = renderWindow->GetQmitkRenderWindow(QString("axial"))->GetSliceNavigationController(); + itk::ReceptorMemberCommand::Pointer command = itk::ReceptorMemberCommand::New(); + command->SetCallbackFunction( this, &QmitkStreamlineTrackingView::OnSliceChanged ); + m_SliceObserverTag1 = slicer->AddObserver( mitk::SliceNavigationController::GeometrySliceEvent(nullptr, 0), command ); + } + + { + mitk::SliceNavigationController* slicer = renderWindow->GetQmitkRenderWindow(QString("sagittal"))->GetSliceNavigationController(); + itk::ReceptorMemberCommand::Pointer command = itk::ReceptorMemberCommand::New(); + command->SetCallbackFunction( this, &QmitkStreamlineTrackingView::OnSliceChanged ); + m_SliceObserverTag2 = slicer->AddObserver( mitk::SliceNavigationController::GeometrySliceEvent(nullptr, 0), command ); + } + + { + mitk::SliceNavigationController* slicer = renderWindow->GetQmitkRenderWindow(QString("coronal"))->GetSliceNavigationController(); + itk::ReceptorMemberCommand::Pointer command = itk::ReceptorMemberCommand::New(); + command->SetCallbackFunction( this, &QmitkStreamlineTrackingView::OnSliceChanged ); + m_SliceObserverTag3 = slicer->AddObserver( mitk::SliceNavigationController::GeometrySliceEvent(nullptr, 0), command ); + } } + } + else + { + m_InteractiveNode = nullptr; + m_InteractivePointSetNode = nullptr; + + mitk::SliceNavigationController* slicer = renderWindow->GetQmitkRenderWindow(QString("axial"))->GetSliceNavigationController(); + slicer->RemoveObserver(m_SliceObserverTag1); + slicer = renderWindow->GetQmitkRenderWindow(QString("sagittal"))->GetSliceNavigationController(); + slicer->RemoveObserver(m_SliceObserverTag2); + slicer = renderWindow->GetQmitkRenderWindow(QString("coronal"))->GetSliceNavigationController(); + slicer->RemoveObserver(m_SliceObserverTag3); + } +} - UpdateGui(); +void QmitkStreamlineTrackingView::OnSliceChanged(const itk::EventObject& /*e*/) +{ + 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); + } + + DoFiberTracking(); } void QmitkStreamlineTrackingView::SetFocus() { } void QmitkStreamlineTrackingView::OnSelectionChanged( berry::IWorkbenchPart::Pointer part, const QList& nodes ) { - m_InputImageNodes.clear(); - m_InputImages.clear(); + m_InputImageNodes.clear(); + m_InputImages.clear(); - for( auto node : nodes ) - { + for( auto node : nodes ) + { - if( node.IsNotNull() && dynamic_cast(node->GetData()) ) + 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())); + } + else if ( dynamic_cast(node->GetData()) ) + { + m_InputImageNodes.push_back(node); + m_InputImages.push_back(dynamic_cast(node->GetData())); + } + else + { + mitk::Image* img = dynamic_cast(node->GetData()); + if (img!=nullptr) { - if( dynamic_cast(node->GetData()) ) - { - m_InputImageNodes.push_back(node); - m_InputImages.push_back(dynamic_cast(node->GetData())); - } - else if ( dynamic_cast(node->GetData()) ) - { - m_InputImageNodes.push_back(node); - m_InputImages.push_back(dynamic_cast(node->GetData())); - } - 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())); - } - } - } + 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())); + } } + } } + } - UpdateGui(); + UpdateGui(); } void QmitkStreamlineTrackingView::UpdateGui() { - m_Controls->m_TensorImageLabel->setText("mandatory"); - - m_Controls->m_fBox->setVisible(false); - m_Controls->m_fLabel->setVisible(false); - m_Controls->m_gBox->setVisible(false); - m_Controls->m_gLabel->setVisible(false); - m_Controls->m_FaImageBox->setVisible(false); - m_Controls->mFaImageLabel->setVisible(false); - m_Controls->m_NormalizeODFsBox->setVisible(false); - - if (m_Controls->m_TissueImageBox->GetSelectedNode().IsNotNull()) - m_Controls->m_SeedGmBox->setVisible(true); + m_Controls->m_TensorImageLabel->setText("mandatory"); + + m_Controls->m_fBox->setVisible(false); + m_Controls->m_fLabel->setVisible(false); + m_Controls->m_gBox->setVisible(false); + m_Controls->m_gLabel->setVisible(false); + m_Controls->m_FaImageBox->setVisible(false); + m_Controls->mFaImageLabel->setVisible(false); + m_Controls->m_NormalizeODFsBox->setVisible(false); + + if (m_Controls->m_TissueImageBox->GetSelectedNode().IsNotNull()) + m_Controls->m_SeedGmBox->setVisible(true); + else + m_Controls->m_SeedGmBox->setVisible(false); + + if(!m_InputImageNodes.empty()) + { + if (m_InputImageNodes.size()>1) + m_Controls->m_TensorImageLabel->setText(m_InputImageNodes.size()+" images selected"); else - m_Controls->m_SeedGmBox->setVisible(false); + m_Controls->m_TensorImageLabel->setText(m_InputImageNodes.at(0)->GetName().c_str()); + m_Controls->m_InputData->setTitle("Input Data"); + m_Controls->commandLinkButton->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); + m_Controls->m_InteractiveBox->setEnabled(true); - if(!m_InputImageNodes.empty()) + if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { - if (m_InputImageNodes.size()>1) - m_Controls->m_TensorImageLabel->setText(m_InputImageNodes.size()+" images selected"); - else - m_Controls->m_TensorImageLabel->setText(m_InputImageNodes.at(0)->GetName().c_str()); - m_Controls->m_InputData->setTitle("Input Data"); - m_Controls->commandLinkButton->setEnabled(true); - - if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) - { - m_Controls->m_fBox->setVisible(true); - m_Controls->m_fLabel->setVisible(true); - m_Controls->m_gBox->setVisible(true); - m_Controls->m_gLabel->setVisible(true); - m_Controls->mFaImageLabel->setVisible(true); - m_Controls->m_FaImageBox->setVisible(true); - - if (m_Controls->m_ModeBox->currentIndex()==1) - m_Controls->m_NormalizeODFsBox->setVisible(true); - } - else if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) - { - m_Controls->mFaImageLabel->setVisible(true); - m_Controls->m_FaImageBox->setVisible(true); - m_Controls->m_NormalizeODFsBox->setVisible(true); - } - else - { - - } + m_Controls->m_fBox->setVisible(true); + m_Controls->m_fLabel->setVisible(true); + m_Controls->m_gBox->setVisible(true); + m_Controls->m_gLabel->setVisible(true); + m_Controls->mFaImageLabel->setVisible(true); + m_Controls->m_FaImageBox->setVisible(true); + + if (m_Controls->m_ModeBox->currentIndex()==1) + m_Controls->m_NormalizeODFsBox->setVisible(true); + } + else if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) + { + m_Controls->mFaImageLabel->setVisible(true); + m_Controls->m_FaImageBox->setVisible(true); + m_Controls->m_NormalizeODFsBox->setVisible(true); } else { - m_Controls->m_InputData->setTitle("Please Select Input Data"); - m_Controls->commandLinkButton->setEnabled(false); + } + } + else + { + m_Controls->m_InputData->setTitle("Please Select Input Data"); + m_Controls->commandLinkButton->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); + m_Controls->m_InteractiveBox->setEnabled(false); + } } void QmitkStreamlineTrackingView::DoFiberTracking() { - if (m_InputImages.empty()) - return; + if (m_InputImages.empty()) + return; - mitk::TrackingDataHandler* trackingHandler; + mitk::TrackingDataHandler* trackingHandler; - if( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) - { - typedef itk::Image< itk::DiffusionTensor3D, 3> TensorImageType; - typedef mitk::ImageToItk CasterType; + if( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) + { + typedef itk::Image< itk::DiffusionTensor3D, 3> TensorImageType; + typedef mitk::ImageToItk CasterType; - 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!"); - return; - } - - if (m_FirstTensorProbRun) - { - QMessageBox::information(nullptr, "Information", "Internally calculating ODF from tensor image and performing probabilistic ODF tractography. Please keep the state of the ODF normalization box (see advanced parameters) in mind. TEND parameters are ignored."); - m_FirstTensorProbRun = false; - } - - TensorImageType::Pointer itkImg = TensorImageType::New(); - mitk::CastToItkImage(m_InputImages.at(0), itkImg); - - typedef itk::TensorImageToQBallImageFilter< float, float > FilterType; - FilterType::Pointer filter = FilterType::New(); - filter->SetInput( itkImg ); - filter->Update(); - - typedef mitk::ImageToItk< mitk::TrackingHandlerOdf::ItkOdfImageType > CasterType; - trackingHandler = new mitk::TrackingHandlerOdf(); - dynamic_cast(trackingHandler)->SetOdfImage(filter->GetOutput()); - dynamic_cast(trackingHandler)->SetGfaThreshold(m_Controls->m_ScalarThresholdBox->value()); - dynamic_cast(trackingHandler)->SetMinMaxNormalize(m_Controls->m_NormalizeODFsBox->isChecked()); - dynamic_cast(trackingHandler)->SetOdfPower(1); - - if (m_Controls->m_FaImageBox->GetSelectedNode().IsNotNull()) - { - ItkFloatImageType::Pointer itkImg = ItkFloatImageType::New(); - mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageBox->GetSelectedNode()->GetData()), itkImg); - - dynamic_cast(trackingHandler)->SetGfaImage(itkImg); - } - } - else - { - trackingHandler = new mitk::TrackingHandlerTensor(); - for (int i=0; i<(int)m_InputImages.size(); i++) - { - TensorImageType::Pointer itkImg = TensorImageType::New(); - mitk::CastToItkImage(m_InputImages.at(i), itkImg); - - dynamic_cast(trackingHandler)->AddTensorImage(itkImg); - } - - if (m_Controls->m_FaImageBox->GetSelectedNode().IsNotNull()) - { - ItkFloatImageType::Pointer itkImg = ItkFloatImageType::New(); - mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageBox->GetSelectedNode()->GetData()), itkImg); - - dynamic_cast(trackingHandler)->SetFaImage(itkImg); - } - - dynamic_cast(trackingHandler)->SetFaThreshold(m_Controls->m_ScalarThresholdBox->value()); - dynamic_cast(trackingHandler)->SetF((float)m_Controls->m_fBox->value()); - dynamic_cast(trackingHandler)->SetG((float)m_Controls->m_gBox->value()); - } - } - else if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) + if (m_Controls->m_ModeBox->currentIndex()==1) { - typedef mitk::ImageToItk< mitk::TrackingHandlerOdf::ItkOdfImageType > CasterType; - trackingHandler = new mitk::TrackingHandlerOdf(); - mitk::TrackingHandlerOdf::ItkOdfImageType::Pointer itkImg = mitk::TrackingHandlerOdf::ItkOdfImageType::New(); - mitk::CastToItkImage(m_InputImages.at(0), itkImg); - dynamic_cast(trackingHandler)->SetOdfImage(itkImg); - dynamic_cast(trackingHandler)->SetGfaThreshold(m_Controls->m_ScalarThresholdBox->value()); - dynamic_cast(trackingHandler)->SetMinMaxNormalize(m_Controls->m_NormalizeODFsBox->isChecked()); - - if (m_Controls->m_FaImageBox->GetSelectedNode().IsNotNull()) - { - ItkFloatImageType::Pointer itkImg = ItkFloatImageType::New(); - mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageBox->GetSelectedNode()->GetData()), itkImg); - - dynamic_cast(trackingHandler)->SetGfaImage(itkImg); - } + if (m_InputImages.size()>1) + { + QMessageBox::information(nullptr, "Information", "Probabilistic tensor tractography is only implemented for single-tensor mode!"); + return; + } + + if (m_FirstTensorProbRun) + { + QMessageBox::information(nullptr, "Information", "Internally calculating ODF from tensor image and performing probabilistic ODF tractography. Please keep the state of the ODF normalization box (see advanced parameters) in mind. TEND parameters are ignored."); + m_FirstTensorProbRun = false; + } + + TensorImageType::Pointer itkImg = TensorImageType::New(); + mitk::CastToItkImage(m_InputImages.at(0), itkImg); + + typedef itk::TensorImageToQBallImageFilter< float, float > FilterType; + FilterType::Pointer filter = FilterType::New(); + filter->SetInput( itkImg ); + filter->Update(); + + typedef mitk::ImageToItk< mitk::TrackingHandlerOdf::ItkOdfImageType > CasterType; + trackingHandler = new mitk::TrackingHandlerOdf(); + dynamic_cast(trackingHandler)->SetOdfImage(filter->GetOutput()); + dynamic_cast(trackingHandler)->SetGfaThreshold(m_Controls->m_ScalarThresholdBox->value()); + dynamic_cast(trackingHandler)->SetMinMaxNormalize(m_Controls->m_NormalizeODFsBox->isChecked()); + dynamic_cast(trackingHandler)->SetOdfPower(1); + + if (m_Controls->m_FaImageBox->GetSelectedNode().IsNotNull()) + { + ItkFloatImageType::Pointer itkImg = ItkFloatImageType::New(); + mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageBox->GetSelectedNode()->GetData()), itkImg); + + dynamic_cast(trackingHandler)->SetGfaImage(itkImg); + } } else { - if (m_Controls->m_ModeBox->currentIndex()==1) - { - QMessageBox::information(nullptr, "Information", "Probabilstic tractography is only implementedfor ODF images."); - return; - } - try { - typedef mitk::ImageToItk< mitk::TrackingHandlerPeaks::PeakImgType > CasterType; - CasterType::Pointer caster = CasterType::New(); - caster->SetInput(m_InputImages.at(0)); - caster->Update(); - mitk::TrackingHandlerPeaks::PeakImgType::Pointer itkImg = caster->GetOutput(); - trackingHandler = new mitk::TrackingHandlerPeaks(); - dynamic_cast(trackingHandler)->SetPeakImage(itkImg); - dynamic_cast(trackingHandler)->SetPeakThreshold(m_Controls->m_ScalarThresholdBox->value()); - } - catch(...) - { - return; - } + trackingHandler = new mitk::TrackingHandlerTensor(); + for (int i=0; i<(int)m_InputImages.size(); i++) + { + TensorImageType::Pointer itkImg = TensorImageType::New(); + mitk::CastToItkImage(m_InputImages.at(i), itkImg); + + dynamic_cast(trackingHandler)->AddTensorImage(itkImg); + } + + if (m_Controls->m_FaImageBox->GetSelectedNode().IsNotNull()) + { + ItkFloatImageType::Pointer itkImg = ItkFloatImageType::New(); + mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageBox->GetSelectedNode()->GetData()), itkImg); + + dynamic_cast(trackingHandler)->SetFaImage(itkImg); + } + + dynamic_cast(trackingHandler)->SetFaThreshold(m_Controls->m_ScalarThresholdBox->value()); + dynamic_cast(trackingHandler)->SetF((float)m_Controls->m_fBox->value()); + dynamic_cast(trackingHandler)->SetG((float)m_Controls->m_gBox->value()); } - - trackingHandler->SetFlipX(m_Controls->m_FlipXBox->isChecked()); - trackingHandler->SetFlipY(m_Controls->m_FlipYBox->isChecked()); - trackingHandler->SetFlipZ(m_Controls->m_FlipZBox->isChecked()); - trackingHandler->SetInterpolate(m_Controls->m_InterpolationBox->isChecked()); - switch (m_Controls->m_ModeBox->currentIndex()) + } + else if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) + { + typedef mitk::ImageToItk< mitk::TrackingHandlerOdf::ItkOdfImageType > CasterType; + trackingHandler = new mitk::TrackingHandlerOdf(); + mitk::TrackingHandlerOdf::ItkOdfImageType::Pointer itkImg = mitk::TrackingHandlerOdf::ItkOdfImageType::New(); + mitk::CastToItkImage(m_InputImages.at(0), itkImg); + dynamic_cast(trackingHandler)->SetOdfImage(itkImg); + dynamic_cast(trackingHandler)->SetGfaThreshold(m_Controls->m_ScalarThresholdBox->value()); + dynamic_cast(trackingHandler)->SetMinMaxNormalize(m_Controls->m_NormalizeODFsBox->isChecked()); + + if (m_Controls->m_FaImageBox->GetSelectedNode().IsNotNull()) { - case 0: - trackingHandler->SetMode(mitk::TrackingDataHandler::MODE::DETERMINISTIC); - break; - case 1: - trackingHandler->SetMode(mitk::TrackingDataHandler::MODE::PROBABILISTIC); - break; - default: - trackingHandler->SetMode(mitk::TrackingDataHandler::MODE::DETERMINISTIC); - } + ItkFloatImageType::Pointer itkImg = ItkFloatImageType::New(); + mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageBox->GetSelectedNode()->GetData()), itkImg); - typedef itk::StreamlineTrackingFilter TrackerType; - TrackerType::Pointer tracker = TrackerType::New(); - - if (m_Controls->m_SeedImageBox->GetSelectedNode().IsNotNull()) + dynamic_cast(trackingHandler)->SetGfaImage(itkImg); + } + } + else + { + if (m_Controls->m_ModeBox->currentIndex()==1) { - ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); - mitk::CastToItkImage(dynamic_cast(m_Controls->m_SeedImageBox->GetSelectedNode()->GetData()), mask); - tracker->SetSeedImage(mask); + QMessageBox::information(nullptr, "Information", "Probabilstic tractography is only implementedfor ODF images."); + return; } - - if (m_Controls->m_MaskImageBox->GetSelectedNode().IsNotNull()) + try { + typedef mitk::ImageToItk< mitk::TrackingHandlerPeaks::PeakImgType > CasterType; + CasterType::Pointer caster = CasterType::New(); + caster->SetInput(m_InputImages.at(0)); + caster->Update(); + mitk::TrackingHandlerPeaks::PeakImgType::Pointer itkImg = caster->GetOutput(); + trackingHandler = new mitk::TrackingHandlerPeaks(); + dynamic_cast(trackingHandler)->SetPeakImage(itkImg); + dynamic_cast(trackingHandler)->SetPeakThreshold(m_Controls->m_ScalarThresholdBox->value()); + } + catch(...) { - ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); - mitk::CastToItkImage(dynamic_cast(m_Controls->m_MaskImageBox->GetSelectedNode()->GetData()), mask); - tracker->SetMaskImage(mask); + return; } - - if (m_Controls->m_StopImageBox->GetSelectedNode().IsNotNull()) + } + + trackingHandler->SetFlipX(m_Controls->m_FlipXBox->isChecked()); + trackingHandler->SetFlipY(m_Controls->m_FlipYBox->isChecked()); + trackingHandler->SetFlipZ(m_Controls->m_FlipZBox->isChecked()); + trackingHandler->SetInterpolate(m_Controls->m_InterpolationBox->isChecked()); + switch (m_Controls->m_ModeBox->currentIndex()) + { + case 0: + trackingHandler->SetMode(mitk::TrackingDataHandler::MODE::DETERMINISTIC); + break; + case 1: + trackingHandler->SetMode(mitk::TrackingDataHandler::MODE::PROBABILISTIC); + break; + default: + trackingHandler->SetMode(mitk::TrackingDataHandler::MODE::DETERMINISTIC); + } + + typedef itk::StreamlineTrackingFilter TrackerType; + TrackerType::Pointer tracker = TrackerType::New(); + + if (m_Controls->m_InteractiveBox->isChecked()) + { + tracker->SetSeedPoints(m_SeedPoints); + } + else if (m_Controls->m_SeedImageBox->GetSelectedNode().IsNotNull()) + { + ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); + mitk::CastToItkImage(dynamic_cast(m_Controls->m_SeedImageBox->GetSelectedNode()->GetData()), mask); + tracker->SetSeedImage(mask); + } + + if (m_Controls->m_MaskImageBox->GetSelectedNode().IsNotNull()) + { + ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); + mitk::CastToItkImage(dynamic_cast(m_Controls->m_MaskImageBox->GetSelectedNode()->GetData()), mask); + tracker->SetMaskImage(mask); + } + + if (m_Controls->m_StopImageBox->GetSelectedNode().IsNotNull()) + { + ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); + mitk::CastToItkImage(dynamic_cast(m_Controls->m_StopImageBox->GetSelectedNode()->GetData()), mask); + tracker->SetStoppingRegions(mask); + } + + if (m_Controls->m_TissueImageBox->GetSelectedNode().IsNotNull()) + { + ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); + mitk::CastToItkImage(dynamic_cast(m_Controls->m_TissueImageBox->GetSelectedNode()->GetData()), mask); + tracker->SetTissueImage(mask); + tracker->SetSeedOnlyGm(m_Controls->m_SeedGmBox->isChecked()); + } + + tracker->SetSeedsPerVoxel(m_Controls->m_SeedsPerVoxelBox->value()); + tracker->SetStepSize(m_Controls->m_StepSizeBox->value()); + //tracker->SetSamplingDistance(0.7); + tracker->SetUseStopVotes(false); + tracker->SetOnlyForwardSamples(false); + tracker->SetAposterioriCurvCheck(false); + tracker->SetMaxNumTracts(m_Controls->m_NumFibersBox->value()); + tracker->SetNumberOfSamples(m_Controls->m_NumSamplesBox->value()); + tracker->SetTrackingHandler(trackingHandler); + tracker->SetAngularThreshold(m_Controls->m_AngularThresholdBox->value()); + tracker->SetMinTractLength(m_Controls->m_MinTractLengthBox->value()); + tracker->SetUseOutputProbabilityMap(m_Controls->m_OutputProbMap->isChecked() && !m_Controls->m_InteractiveBox->isChecked()); + tracker->Update(); + + delete trackingHandler; + + if (m_Controls->m_InteractiveBox->isChecked()) + { + vtkSmartPointer fiberBundle = tracker->GetFiberPolyData(); + if (m_InteractiveNode.IsNull()) { - ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); - mitk::CastToItkImage(dynamic_cast(m_Controls->m_StopImageBox->GetSelectedNode()->GetData()), mask); - tracker->SetStoppingRegions(mask); + m_InteractiveNode = mitk::DataNode::New(); + QString name("InteractiveFib"); + m_InteractiveNode->SetName(name.toStdString()); + GetDataStorage()->Add(m_InteractiveNode); } - if (m_Controls->m_TissueImageBox->GetSelectedNode().IsNotNull()) + mitk::FiberBundle::Pointer fib = mitk::FiberBundle::New(fiberBundle); + fib->SetReferenceGeometry(dynamic_cast(m_InputImageNodes.at(0)->GetData())->GetGeometry()); + if (m_Controls->m_ResampleFibersBox->isChecked()) + fib->Compress(m_Controls->m_FiberErrorBox->value()); + + m_InteractiveNode->SetData(fib); + } + else if (!tracker->GetUseOutputProbabilityMap()) + { + vtkSmartPointer fiberBundle = tracker->GetFiberPolyData(); + if (fiberBundle->GetNumberOfLines() == 0) { - ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); - mitk::CastToItkImage(dynamic_cast(m_Controls->m_TissueImageBox->GetSelectedNode()->GetData()), mask); - tracker->SetTissueImage(mask); - tracker->SetSeedOnlyGm(m_Controls->m_SeedGmBox->isChecked()); + QMessageBox warnBox; + warnBox.setWindowTitle("Warning"); + warnBox.setText("No fiberbundle was generated!"); + warnBox.setDetailedText("No fibers were generated using the parameters: \n\n" + m_Controls->m_FaThresholdLabel->text() + "\n" + m_Controls->m_AngularThresholdLabel->text() + "\n" + m_Controls->m_fLabel->text() + "\n" + m_Controls->m_gLabel->text() + "\n" + m_Controls->m_StepsizeLabel->text() + "\n" + m_Controls->m_MinTractLengthLabel->text() + "\n" + m_Controls->m_SeedsPerVoxelLabel->text() + "\n\nPlease check your parametersettings."); + warnBox.setIcon(QMessageBox::Warning); + warnBox.exec(); + return; } - - tracker->SetSeedsPerVoxel(m_Controls->m_SeedsPerVoxelBox->value()); - tracker->SetStepSize(m_Controls->m_StepSizeBox->value()); - //tracker->SetSamplingDistance(0.7); - tracker->SetUseStopVotes(false); - tracker->SetOnlyForwardSamples(false); - tracker->SetAposterioriCurvCheck(false); - tracker->SetMaxNumTracts(m_Controls->m_NumFibersBox->value()); - tracker->SetNumberOfSamples(m_Controls->m_NumSamplesBox->value()); - tracker->SetTrackingHandler(trackingHandler); - tracker->SetAngularThreshold(m_Controls->m_AngularThresholdBox->value()); - tracker->SetMinTractLength(m_Controls->m_MinTractLengthBox->value()); - tracker->SetUseOutputProbabilityMap(m_Controls->m_OutputProbMap->isChecked()); - tracker->Update(); - - delete trackingHandler; - - if (!tracker->GetUseOutputProbabilityMap()) - { - vtkSmartPointer fiberBundle = tracker->GetFiberPolyData(); - if (fiberBundle->GetNumberOfLines() == 0) - { - QMessageBox warnBox; - warnBox.setWindowTitle("Warning"); - warnBox.setText("No fiberbundle was generated!"); - warnBox.setDetailedText("No fibers were generated using the parameters: \n\n" + m_Controls->m_FaThresholdLabel->text() + "\n" + m_Controls->m_AngularThresholdLabel->text() + "\n" + m_Controls->m_fLabel->text() + "\n" + m_Controls->m_gLabel->text() + "\n" + m_Controls->m_StepsizeLabel->text() + "\n" + m_Controls->m_MinTractLengthLabel->text() + "\n" + m_Controls->m_SeedsPerVoxelLabel->text() + "\n\nPlease check your parametersettings."); - warnBox.setIcon(QMessageBox::Warning); - warnBox.exec(); - return; - } - mitk::FiberBundle::Pointer fib = mitk::FiberBundle::New(fiberBundle); - fib->SetReferenceGeometry(dynamic_cast(m_InputImageNodes.at(0)->GetData())->GetGeometry()); - if (m_Controls->m_ResampleFibersBox->isChecked()) - fib->Compress(m_Controls->m_FiberErrorBox->value()); - fib->ColorFibersByOrientation(); - - mitk::DataNode::Pointer node = mitk::DataNode::New(); - node->SetData(fib); - QString name("FiberBundle_"); - name += m_InputImageNodes.at(0)->GetName().c_str(); - name += "_Streamline"; - node->SetName(name.toStdString()); - node->SetVisibility(true); - GetDataStorage()->Add(node, m_InputImageNodes.at(0)); - } - else - { - TrackerType::ItkDoubleImgType::Pointer outImg = tracker->GetOutputProbabilityMap(); - mitk::Image::Pointer img = mitk::Image::New(); - img->InitializeByItk(outImg.GetPointer()); - img->SetVolume(outImg->GetBufferPointer()); - - mitk::DataNode::Pointer node = mitk::DataNode::New(); - node->SetData(img); - QString name("ProbabilityMap_"); - name += m_InputImageNodes.at(0)->GetName().c_str(); - node->SetName(name.toStdString()); - node->SetVisibility(true); - - 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_InputImageNodes.at(0)); - } + mitk::FiberBundle::Pointer fib = mitk::FiberBundle::New(fiberBundle); + fib->SetReferenceGeometry(dynamic_cast(m_InputImageNodes.at(0)->GetData())->GetGeometry()); + if (m_Controls->m_ResampleFibersBox->isChecked()) + fib->Compress(m_Controls->m_FiberErrorBox->value()); + fib->ColorFibersByOrientation(); + + mitk::DataNode::Pointer node = mitk::DataNode::New(); + node->SetData(fib); + QString name("FiberBundle_"); + name += m_InputImageNodes.at(0)->GetName().c_str(); + name += "_Streamline"; + node->SetName(name.toStdString()); + GetDataStorage()->Add(node, m_InputImageNodes.at(0)); + } + else + { + TrackerType::ItkDoubleImgType::Pointer outImg = tracker->GetOutputProbabilityMap(); + mitk::Image::Pointer img = mitk::Image::New(); + img->InitializeByItk(outImg.GetPointer()); + img->SetVolume(outImg->GetBufferPointer()); + + mitk::DataNode::Pointer node = mitk::DataNode::New(); + node->SetData(img); + QString name("ProbabilityMap_"); + name += m_InputImageNodes.at(0)->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_InputImageNodes.at(0)); + } } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.h b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.h index ac09642ea2..1da8700f21 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.h @@ -1,85 +1,99 @@ /*=================================================================== 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 QmitkStreamlineTrackingView_h #define QmitkStreamlineTrackingView_h #include #include "ui_QmitkStreamlineTrackingViewControls.h" #include #include +#include #include #include #include #include #include #include +#include +#include +#include /*! \brief View for tensor based deterministic streamline fiber tracking. */ class QmitkStreamlineTrackingView : public QmitkAbstractView { - // this is needed for all Qt objects that should have a Qt meta-object - // (everything that derives from QObject and wants to have signal/slots) - Q_OBJECT + // this is needed for all Qt objects that should have a Qt meta-object + // (everything that derives from QObject and wants to have signal/slots) + Q_OBJECT public: - static const std::string VIEW_ID; + static const std::string VIEW_ID; - typedef itk::Image< unsigned char, 3 > ItkUCharImageType; - typedef itk::Image< float, 3 > ItkFloatImageType; + typedef itk::Image< unsigned char, 3 > ItkUCharImageType; + typedef itk::Image< float, 3 > ItkFloatImageType; - QmitkStreamlineTrackingView(); - virtual ~QmitkStreamlineTrackingView(); + QmitkStreamlineTrackingView(); + virtual ~QmitkStreamlineTrackingView(); - virtual void CreateQtPartControl(QWidget *parent) override; + virtual void CreateQtPartControl(QWidget *parent) override; - /// - /// Sets the focus to an internal widget. - /// - virtual void SetFocus() override; + /// + /// Sets the focus to an internal widget. + /// + virtual void SetFocus() override; protected slots: - void DoFiberTracking(); ///< start fiber tracking - void UpdateGui(); + void DoFiberTracking(); ///< start fiber tracking + void UpdateGui(); + void ToggleInteractive(); protected: - /// \brief called by QmitkAbstractView when DataManager's selection has changed - virtual void OnSelectionChanged(berry::IWorkbenchPart::Pointer part, const QList& nodes) override; + /// \brief called by QmitkAbstractView when DataManager's selection has changed + virtual void OnSelectionChanged(berry::IWorkbenchPart::Pointer part, const QList& nodes) override; - Ui::QmitkStreamlineTrackingViewControls* m_Controls; + Ui::QmitkStreamlineTrackingViewControls* m_Controls; protected slots: private: - std::vector< mitk::DataNode::Pointer > m_InputImageNodes; ///< input images - std::vector< mitk::Image::Pointer > m_InputImages; ///< input image datanode - bool m_FirstTensorProbRun; + int m_SliceObserverTag1; + int m_SliceObserverTag2; + int m_SliceObserverTag3; + std::vector< itk::Point > m_SeedPoints; + mitk::DataNode::Pointer m_InteractiveNode; + mitk::DataNode::Pointer m_InteractivePointSetNode; + + std::vector< mitk::DataNode::Pointer > m_InputImageNodes; ///< input images + std::vector< mitk::Image::Pointer > m_InputImages; ///< input image datanode + bool m_FirstTensorProbRun; + + void OnSliceChanged(const itk::EventObject& e); }; #endif // _QMITKFIBERTRACKINGVIEW_H_INCLUDED 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 d0d2631389..96b98a3cdf 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,685 +1,765 @@ QmitkStreamlineTrackingViewControls 0 0 358 - 733 + 779 0 0 QmitkTemplate 0 0 3 3 + + + + Please Select Input Data + + + + + + Only track insida mask area. + + + Mask Image: + + + + + + + Fibers that enter a region defined in this image will stop immediately. + + + QComboBox::AdjustToMinimumContentsLength + + + + - + + + + + + + + If an image is selected, the stopping criterion is not calculated from the input image but instead the selected image is used. + + + QComboBox::AdjustToMinimumContentsLength + + + + - + + + + + + + + Binary seed ROI. If not specified, the whole image area is seeded. + + + Seed Image: + + + + + + + <html><head/><body><p><span style=" color:#ff0000;">mandatory</span></p></body></html> + + + true + + + + + + + Seed points are only placed inside the mask image. + + + QComboBox::AdjustToMinimumContentsLength + + + + - + + + + + + + + Only track insida mask area. + + + Stop Image: + + + + + + + Input DTI + + + Tensor/Peak/ODF Image: + + + + + + + Only track insida mask area. + + + FA/GFA Image: + + + + + + + Tractography is only performed inside the mask image. Fibers that leave the mask image are stopped. + + + QComboBox::AdjustToMinimumContentsLength + + + + - + + + + + + + + Only track insida mask area. + + + Tissue Image: + + + + + + + Tissue label image needed for gray matter seeding (WM=3, GM=1). Use e.g. MRtrix 5ttgen to generate such a label image. + + + QComboBox::AdjustToMinimumContentsLength + + + + - + + + + + + + false Start Tractography - - - - Qt::Vertical - - - QSizePolicy::Expanding + + + + + 0 + 0 + - - - 20 - 220 - + + Parameters - + + + + + + + + Mode: + + + + + + + Toggle between deterministic and probabilistic tractography. Some modes might not be available for all types of tractography. + + + + Deterministic + + + + + Probabilistic + + + + + + + + + + + Seeds per Voxel: + + + + + + + Number of seed points placed in each voxel. + + + 1 + + + 9999999 + + + + + + + + + + Max. num. fibers: + + + + + + + Tractography is stopped after the desired number of fibers is reached, even before all seed points are processed. + + + 999999999 + + + + + + + + + + Cutoff: + + + + + + + Threshold on peak magnitude, FA, GFA, ... + + + 5 + + + 1.000000000000000 + + + 0.100000000000000 + + + 0.100000000000000 + + + + + + + Qt::Horizontal + + + QSizePolicy::Fixed + + + + 200 + 0 + + + + + + - + Advanced Paramaters Step Size: Step size (in voxels) 2 0.010000000000000 10.000000000000000 0.100000000000000 0.500000000000000 Angular Threshold: Default: 90° * step_size -1 90 1 -1 Min. Tract Length: Shorter fibers are discarded. 1 999.000000000000000 1.000000000000000 20.000000000000000 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: 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 g: 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 Flip directions: 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 Neighborhood Samples: Number of neighborhood samples that are used to determine the next fiber progression direction. 50 Resample 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 If false, nearest neighbor interpolation is used. Enable Trilinear Interpolation true Requires tissue image. Enable gray matter seeding false Min-max normalize ODFs true Output probability map instead of tractogram. Output probability map false - - - - Please Select Input Data + + + + Qt::Vertical - - - - - Only track insida mask area. - - - Mask Image: - - - - - - - Fibers that enter a region defined in this image will stop immediately. - - - QComboBox::AdjustToMinimumContentsLength - - - - - - - - - - - - - If an image is selected, the stopping criterion is not calculated from the input image but instead the selected image is used. - - - QComboBox::AdjustToMinimumContentsLength - - - - - - - - - - - - - Binary seed ROI. If not specified, the whole image area is seeded. - - - Seed Image: - - - - - - - <html><head/><body><p><span style=" color:#ff0000;">mandatory</span></p></body></html> - - - true - - - - - - - Seed points are only placed inside the mask image. - - - QComboBox::AdjustToMinimumContentsLength - - - - - - - - - - - - - Only track insida mask area. - - - Stop Image: - - - - - - - Input DTI - - - Tensor/Peak/ODF Image: - - - - - - - Only track insida mask area. - - - FA/GFA Image: - - - - - - - Tractography is only performed inside the mask image. Fibers that leave the mask image are stopped. - - - QComboBox::AdjustToMinimumContentsLength - - - - - - - - - - - - - Only track insida mask area. - - - Tissue Image: - - - - - - - Tissue label image needed for gray matter seeding (WM=3, GM=1). Use e.g. MRtrix 5ttgen to generate such a label image. - - - QComboBox::AdjustToMinimumContentsLength - - - - - - - - - - - + + QSizePolicy::Expanding + + + + 20 + 220 + + + - + 0 0 - Parameters + Interactive Tracking - - - - - - - - Mode: - - - - - - - Toggle between deterministic and probabilistic tractography. Some modes might not be available for all types of tractography. - - - - Deterministic - - - - - Probabilistic - - - - + - + - Seeds per Voxel: + Radius: - - - - Number of seed points placed in each voxel. - - - 1 + + + + false - - 9999999 + + Enable interactive tracking - + - Max. num. fibers: + Num. Seeds: - + - Tractography is stopped after the desired number of fibers is reached, even before all seed points are processed. + Number of seed points normally distributed around selected position. - - 999999999 + + 1 - - - - - - + + 9999999 - - Cutoff: + + 50 - - + + - Threshold on peak magnitude, FA, GFA, ... + Seedpoints are normally distributed within a sphere centered at the selected position with the specified radius. - 5 + 2 - 1.000000000000000 + 50.000000000000000 0.100000000000000 - 0.100000000000000 + 2.500000000000000 - - - - Qt::Horizontal - - - QSizePolicy::Fixed - - - - 200 - 0 - - - - QmitkDataStorageComboBoxWithSelectNone QComboBox
QmitkDataStorageComboBoxWithSelectNone.h