diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingDataHandler.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingDataHandler.h index a4dbc28c1f..ea4f1b0a11 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingDataHandler.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingDataHandler.h @@ -1,379 +1,380 @@ /*=================================================================== 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 itk::Statistics::MersenneTwisterRandomVariateGenerator ItkRngType; 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(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 bool WorldToIndex(itk::Point& pos, itk::Index<3>& index) = 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; } void SetRandom( bool random ) { m_Random = random; if (!random) { m_Rng.seed(0); std::srand(0); m_RngItk->SetSeed(0); } else { m_Rng.seed(); m_RngItk->SetSeed(); std::srand(std::time(0)); } } double GetRandDouble(const double & a, const double & b) { return m_RngItk->GetUniformVariate(a, b); } protected: float m_AngularThreshold; bool m_Interpolate; bool m_FlipX; bool m_FlipY; bool m_FlipZ; MODE m_Mode; BoostRngType m_Rng; ItkRngType::Pointer m_RngItk; bool m_NeedsDataInit; bool m_Random; void DataModified() { m_NeedsDataInit = true; } 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) { pix = 0; MITK_WARN << "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) { pix = 0; MITK_WARN << "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; 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 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) { pix.Fill(0.0); 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 c368ef2ff5..f36ba50b80 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.cpp @@ -1,273 +1,279 @@ /*=================================================================== 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 #include #include namespace mitk { TrackingHandlerOdf::TrackingHandlerOdf() : m_GfaThreshold(0.2) , m_OdfThreshold(0.1) , m_SharpenOdfs(false) , m_NumProbSamples(1) { } TrackingHandlerOdf::~TrackingHandlerOdf() { } +bool TrackingHandlerOdf::WorldToIndex(itk::Point& pos, itk::Index<3>& index) +{ + m_OdfImage->TransformPhysicalPointToIndex(pos, index); + return m_OdfImage->GetLargestPossibleRegion().IsInside(index); +} + void TrackingHandlerOdf::InitForTracking() { MITK_INFO << "Initializing ODF tracker."; if (m_NeedsDataInit) { m_OdfHemisphereIndices.clear(); itk::OrientationDistributionFunction< float, ODF_SAMPLING_SIZE > 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_NeedsDataInit = false; } std::cout << "TrackingHandlerOdf - GFA threshold: " << m_GfaThreshold << std::endl; std::cout << "TrackingHandlerOdf - ODF threshold: " << m_OdfThreshold << std::endl; if (m_SharpenOdfs) std::cout << "TrackingHandlerOdf - Sharpening ODfs" << std::endl; } int TrackingHandlerOdf::SampleOdf(vnl_vector< float >& probs, vnl_vector< float >& angles) { boost::random::discrete_distribution dist(probs.begin(), probs.end()); int sampled_idx = 0; int max_sample_idx = -1; float max_prob = 0; int trials = 0; for (int i=0; i> sampler(m_Rng, dist); sampled_idx = sampler(); } if (probs[sampled_idx]>max_prob && probs[sampled_idx]>m_OdfThreshold && fabs(angles[sampled_idx])>=m_AngularThreshold) { max_prob = probs[sampled_idx]; max_sample_idx = sampled_idx; } else if ( (probs[sampled_idx]<=m_OdfThreshold || fabs(angles[sampled_idx]) 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_OdfImage->TransformPhysicalPointToIndex(pos, idx); if ( !m_OdfImage->GetLargestPossibleRegion().IsInside(idx) ) return output_direction; // check GFA threshold for termination 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_OdfImage, m_Interpolate); vnl_vector< float > probs; probs.set_size(m_OdfHemisphereIndices.size()); vnl_vector< float > angles; angles.set_size(m_OdfHemisphereIndices.size()); angles.fill(1.0); // Find ODF maximum and remove <0 values float max_odf_val = 0; float min_odf_val = 999; 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_odf_val) { max_odf_val = odf_values[i]; max_idx_d = c; } if (odf_values[i]0) { probs /= odf_sum; max_odf_val /= odf_sum; } } // no previous direction if (max_odf_val>m_OdfThreshold && (olddirs.empty() || last_dir.magnitude()<=0.5)) { if (m_Mode==MODE::DETERMINISTIC) // return maximum peak { output_direction = m_OdfFloatDirs.get_row(max_idx_d); return output_direction * max_odf_val; } else if (m_Mode==MODE::PROBABILISTIC) // sample from complete ODF { int max_sample_idx = SampleOdf(probs, angles); if (max_sample_idx>=0) output_direction = m_OdfFloatDirs.get_row(max_sample_idx) * probs[max_sample_idx]; return output_direction; } } else if (max_odf_val<=m_OdfThreshold) // return (0,0,0) { return output_direction; } // correct previous direction if (m_FlipX) last_dir[0] *= -1; if (m_FlipY) last_dir[1] *= -1; if (m_FlipZ) last_dir[2] *= -1; // calculate angles between previous direction and ODF directions angles = m_OdfFloatDirs*last_dir; float probs_sum = 0; float max_prob = 0; for (unsigned int i=0; imax_prob && odf_val>m_OdfThreshold) { // use maximum peak of the ODF weighted with the directional prior max_prob = odf_val; vnl_vector_fixed d = m_OdfFloatDirs.get_row(i); if (angle<0) d *= -1; output_direction = odf_val*d; } else if (m_Mode==MODE::PROBABILISTIC) { // update ODF probabilties with the ODF values pow(abs_angle, m_DirPriorPower) probs[i] = odf_val; probs_sum += probs[i]; } } // do probabilistic sampling if (m_Mode==MODE::PROBABILISTIC && probs_sum>0.0001) { int max_sample_idx = SampleOdf(probs, angles); if (max_sample_idx>=0) { output_direction = m_OdfFloatDirs.get_row(max_sample_idx); if (angles[max_sample_idx]<0) output_direction *= -1; output_direction *= probs[max_sample_idx]; } } // check hard angular threshold float mag = output_direction.magnitude(); if (mag>=0.0001) { output_direction.normalize(); float a = dot_product(output_direction, last_dir); if (a #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(); + TrackingHandlerOdf(); + ~TrackingHandlerOdf(); - typedef TensorImage::PixelType TensorType; - typedef OdfImage::ItkOdfImageType ItkOdfImageType; - typedef itk::Image< vnl_vector_fixed, 3> ItkPDImgType; + typedef TensorImage::PixelType TensorType; + typedef OdfImage::ItkOdfImageType ItkOdfImageType; + typedef itk::Image< vnl_vector_fixed, 3> ItkPDImgType; - void InitForTracking(); ///< calls InputDataValidForTracking() and creates feature images - 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 InitForTracking(); ///< calls InputDataValidForTracking() and creates feature images + 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 WorldToIndex(itk::Point& pos, itk::Index<3>& index); - void SetSharpenOdfs(bool doSharpen) { m_SharpenOdfs=doSharpen; } - void SetOdfThreshold(float odfThreshold){ m_OdfThreshold = odfThreshold; } - void SetGfaThreshold(float gfaThreshold){ m_GfaThreshold = gfaThreshold; } - void SetOdfImage( ItkOdfImageType::Pointer img ){ m_OdfImage = img; DataModified(); } - void SetGfaImage( ItkFloatImgType::Pointer img ){ m_GfaImage = img; DataModified(); } - void SetMode( MODE m ){ m_Mode = m; } + void SetSharpenOdfs(bool doSharpen) { m_SharpenOdfs=doSharpen; } + void SetOdfThreshold(float odfThreshold){ m_OdfThreshold = odfThreshold; } + void SetGfaThreshold(float gfaThreshold){ m_GfaThreshold = gfaThreshold; } + void SetOdfImage( ItkOdfImageType::Pointer img ){ m_OdfImage = img; DataModified(); } + void SetGfaImage( ItkFloatImgType::Pointer img ){ m_GfaImage = img; DataModified(); } + 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(); } + 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 SetNumProbSamples(int NumProbSamples); + int OdfPower() const; + void SetNumProbSamples(int NumProbSamples); protected: - int SampleOdf(vnl_vector< float >& probs, vnl_vector< float >& angles); - - float m_GfaThreshold; - float m_OdfThreshold; - bool m_SharpenOdfs; - 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_NumProbSamples; + int SampleOdf(vnl_vector< float >& probs, vnl_vector< float >& angles); + + float m_GfaThreshold; + float m_OdfThreshold; + bool m_SharpenOdfs; + 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_NumProbSamples; }; } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.cpp index 111bf7baf6..34e609b0ca 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.cpp @@ -1,287 +1,293 @@ /*=================================================================== 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() { } +bool TrackingHandlerPeaks::WorldToIndex(itk::Point& pos, itk::Index<3>& index) +{ + m_DummyImage->TransformPhysicalPointToIndex(pos, index); + return m_DummyImage->GetLargestPossibleRegion().IsInside(index); +} + void TrackingHandlerPeaks::InitForTracking() { MITK_INFO << "Initializing peak tracker."; if (m_NeedsDataInit) { 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; m_NeedsDataInit = false; } std::cout << "TrackingHandlerPeaks - Peak threshold: " << m_PeakThreshold << std::endl; } 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 (magGetIntegerVariate(m_NumDirs-1); out_dir = GetDirection(idx3, i); if (out_dir.magnitude()>mitk::eps) { oldDir[0] = out_dir[0]; oldDir[1] = out_dir[1]; oldDir[2] = out_dir[2]; found = true; break; } } } if (!found) { // 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(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 3f00a26c13..c1d5be6bf5 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.h @@ -1,83 +1,84 @@ /*=================================================================== 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(); + TrackingHandlerPeaks(); + ~TrackingHandlerPeaks(); - typedef itk::Image< float, 4 > PeakImgType; ///< MRtrix peak image type + typedef itk::Image< float, 4 > PeakImgType; ///< MRtrix peak image type - void InitForTracking(); ///< calls InputDataValidForTracking() and creates feature images - 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 InitForTracking(); ///< calls InputDataValidForTracking() and creates feature images + 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 WorldToIndex(itk::Point& pos, itk::Index<3>& index); - void SetPeakThreshold(float thr){ m_PeakThreshold = thr; } - void SetPeakImage( PeakImgType::Pointer image ){ m_PeakImage = image; DataModified(); } - void SetApplyDirectionMatrix( bool applyDirectionMatrix ){ m_ApplyDirectionMatrix = applyDirectionMatrix; } + void SetPeakThreshold(float thr){ m_PeakThreshold = thr; } + void SetPeakImage( PeakImgType::Pointer image ){ m_PeakImage = image; DataModified(); } + 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."; - } + 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); + 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::ConstPointer m_PeakImage; - float m_PeakThreshold; - int m_NumDirs; + PeakImgType::ConstPointer 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; + itk::Vector spacing3; + itk::Point origin3; + itk::Matrix direction3; + itk::ImageRegion<3> imageRegion3; + vnl_matrix_fixed m_FloatImageRotation; - ItkUcharImgType::Pointer m_DummyImage; + ItkUcharImgType::Pointer m_DummyImage; - bool m_ApplyDirectionMatrix; + bool m_ApplyDirectionMatrix; }; } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.cpp index ad9edb8b68..0206c82f0a 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.cpp @@ -1,932 +1,939 @@ /*=================================================================== 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_BidirectionalFiberSampling(false) , m_ZeroDirWmFeatures(true) , 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 > +bool TrackingHandlerRandomForest< ShOrder, NumberOfSignalFeatures >::WorldToIndex(itk::Point& pos, itk::Index<3>& index) +{ + m_DwiFeatureImages.at(0)->TransformPhysicalPointToIndex(pos, index); + return m_DwiFeatureImages.at(0)->GetLargestPossibleRegion().IsInside(index); +} + 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->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(mitk_dwi)); filter->SetGradientImage( mitk::DiffusionPropertyHelper::GetOriginalGradientContainer(mitk_dwi), mitk::DiffusionPropertyHelper::GetItkVectorImage(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->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(mitk_dwi)); filter->SetGradientImage( mitk::DiffusionPropertyHelper::GetOriginalGradientContainer(mitk_dwi), mitk::DiffusionPropertyHelper::GetItkVectorImage(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."; if (m_NeedsDataInit) { InputDataValidForTracking(); m_DwiFeatureImages.clear(); InitDwiImageFeatures<>(m_InputDwis.at(0)); m_NeedsDataInit = false; } } template< int ShOrder, int NumberOfSignalFeatures > 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->GetNumFeatures()) ); 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->GetNumClasses())); 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; iGetNumClasses(); 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 = m_Forest->IndexToClassLabel(i); if (classLabel 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? { if (abs_angle>=m_AngularThreshold) // is angle between the directions smaller than our hard threshold? { if (angle<0) // make sure we don't walk backwards d *= -1; float w_i = probs(0,i)*abs_angle; 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) { boost::random::discrete_distribution dist(probs2.begin(), probs2.end()); int sampled_idx = 0; for (int i=0; i<50; i++) // we allow 50 trials to exceed m_AngularThreshold { #pragma omp critical { boost::random::variate_generator> sampler(m_Rng, dist); sampled_idx = sampler(); } if ( probs2[sampled_idx]>0.1 && (!check_last_dir || (check_last_dir && fabs(angles[sampled_idx])>=m_AngularThreshold)) ) break; } 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]); std::shared_ptr< vigra::RandomForest > forest = trees.at(0); forest->options_.tree_count_ = m_NumTrees; m_Forest = mitk::TractographyForest::New(forest); 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_INFO << "No forest available!"; else { if (m_Forest->GetNumTrees() <= 0) MITK_ERROR << "Forest contains no trees!"; if ( m_Forest->GetNumFeatures() != static_cast(NumberOfSignalFeatures+3*m_NumPreviousDirections+additional_features) ) MITK_ERROR << "Wrong number of features in forest: got " << m_Forest->GetNumFeatures() << ", expected " << (NumberOfSignalFeatures+3*m_NumPreviousDirections+additional_features); } if(m_Forest && m_Forest->GetNumTrees()>0 && m_Forest->GetNumFeatures() == static_cast(NumberOfSignalFeatures+3*m_NumPreviousDirections+additional_features)) return true; return false; } template< int ShOrder, int NumberOfSignalFeatures > void TrackingHandlerRandomForest< ShOrder, NumberOfSignalFeatures >::InitForTraining() { 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); 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(); 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 (static_cast(i)GetVariate()*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 <= static_cast(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"; } } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.h index 0c3a37df5e..8595555f5c 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.h @@ -1,161 +1,162 @@ /*=================================================================== 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 #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; DataModified(); } void AddDwi( Image::Pointer img ){ m_InputDwis.push_back(img); DataModified(); } void SetTractograms( std::vector< FiberBundle::Pointer > tractograms ) { m_Tractograms.clear(); for (auto fib : tractograms) m_Tractograms.push_back(fib); DataModified(); } void SetMaskImages( std::vector< ItkUcharImgType::Pointer > images ){ m_MaskImages = images; DataModified(); } void SetWhiteMatterImages( std::vector< ItkUcharImgType::Pointer > images ){ m_WhiteMatterImages = images; DataModified(); } void SetFiberVolumeModImages( std::vector< ItkFloatImgType::Pointer > images ){ m_FiberVolumeModImages = images; DataModified(); } void SetAdditionalFeatureImages( std::vector< std::vector< ItkFloatImgType::Pointer > > images ){ m_AdditionalFeatureImages = images; DataModified(); } void StartTraining(); void SetMode( MODE m ){ m_Mode = m; } void SetForest(mitk::TractographyForest::Pointer forest){ m_Forest = forest; } 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 SetFiberSamplingStep(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; } void InitForTracking(); ///< calls InputDataValidForTracking() and creates feature images 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 WorldToIndex(itk::Point& pos, itk::Index<3>& index); bool IsForestValid(); ///< true is forest is not null, has more than 0 trees and the correct number of features (NumberOfSignalFeatures + 3) mitk::TractographyForest::Pointer GetForest(){ return m_Forest; } 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 mitk::TractographyForest::Pointer 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 36c0894ac8..148ca0dabf 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.cpp @@ -1,361 +1,367 @@ /*=================================================================== 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_FaThreshold(0.1) , m_F(1.0) , m_G(0.0) , m_InterpolateTensors(true) , m_NumberOfInputs(0) { } TrackingHandlerTensor::~TrackingHandlerTensor() { } void TrackingHandlerTensor::InitForTracking() { MITK_INFO << "Initializing tensor tracker."; if (m_NeedsDataInit) { 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; for (int x=0; x<(int)m_TensorImages.at(0)->GetLargestPossibleRegion().GetSize()[0]; x++) for (int y=0; y<(int)m_TensorImages.at(0)->GetLargestPossibleRegion().GetSize()[1]; y++) for (int z=0; z<(int)m_TensorImages.at(0)->GetLargestPossibleRegion().GetSize()[2]; z++) { ItkTensorImageType::IndexType index; index[0] = x; index[1] = y; index[2] = z; for (int i=0; i dir; tensor = m_TensorImages.at(i)->GetPixel(index); 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); } m_NeedsDataInit = false; } if (m_F+m_G>1.0) { float temp = m_F+m_G; m_F /= temp; m_G /= temp; } std::cout << "TrackingHandlerTensor - FA threshold: " << m_FaThreshold << std::endl; std::cout << "TrackingHandlerTensor - f: " << m_F << std::endl; std::cout << "TrackingHandlerTensor - g: " << m_G << std::endl; } 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; } +bool TrackingHandlerTensor::WorldToIndex(itk::Point& pos, itk::Index<3>& index) +{ + m_TensorImages.at(0)->TransformPhysicalPointToIndex(pos, index); + return m_TensorImages.at(0)->GetLargestPossibleRegion().IsInside(index); +} + 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(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 700481142f..6e4311d56b 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.h @@ -1,89 +1,90 @@ /*=================================================================== 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 #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 TensorImage::PixelType TensorType; - typedef TensorImage::ItkTensorImageType ItkTensorImageType; - typedef itk::Image< vnl_vector_fixed, 3> ItkPDImgType; - - - void InitForTracking(); ///< calls InputDataValidForTracking() and creates feature images - 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::ConstPointer img ){ m_TensorImages.push_back(img); DataModified(); } - void SetTensorImage( ItkTensorImageType::ConstPointer img ){ m_TensorImages.clear(); m_TensorImages.push_back(img); DataModified(); } - void ClearTensorImages(){ m_TensorImages.clear(); DataModified(); } - void SetFaImage( ItkFloatImgType::Pointer img ){ m_FaImage = img; DataModified(); } - 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(); } + TrackingHandlerTensor(); + ~TrackingHandlerTensor(); + + typedef TensorImage::PixelType TensorType; + typedef TensorImage::ItkTensorImageType ItkTensorImageType; + typedef itk::Image< vnl_vector_fixed, 3> ItkPDImgType; + + + void InitForTracking(); ///< calls InputDataValidForTracking() and creates feature images + 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 WorldToIndex(itk::Point& pos, itk::Index<3>& index); + + void SetF(float f){ m_F = f; } + void SetG(float g){ m_G = g; } + void SetFaThreshold(float FaThreshold){ m_FaThreshold = FaThreshold; } + void AddTensorImage( ItkTensorImageType::ConstPointer img ){ m_TensorImages.push_back(img); DataModified(); } + void SetTensorImage( ItkTensorImageType::ConstPointer img ){ m_TensorImages.clear(); m_TensorImages.push_back(img); DataModified(); } + void ClearTensorImages(){ m_TensorImages.clear(); DataModified(); } + void SetFaImage( ItkFloatImgType::Pointer img ){ m_FaImage = img; DataModified(); } + 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); + 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; + 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::ConstPointer > 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; + 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::ConstPointer > 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 597e842961..b4fc5f6c9d 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp @@ -1,1042 +1,1036 @@ /*=================================================================== 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" #include #include #define _USE_MATH_DEFINES #include namespace itk { StreamlineTrackingFilter ::StreamlineTrackingFilter() : m_PauseTracking(false) , m_AbortTracking(false) , m_BuildFibersFinished(false) , m_BuildFibersReady(0) , m_FiberPolyData(nullptr) , m_Points(nullptr) , m_Cells(nullptr) , m_StoppingRegions(nullptr) , m_TargetRegions(nullptr) , m_SeedImage(nullptr) , m_MaskImage(nullptr) , m_TissueImage(nullptr) , m_OutputProbabilityMap(nullptr) , m_AngularThresholdDeg(-1) , m_StepSizeVox(-1) , m_SamplingDistanceVox(-1) , m_AngularThreshold(-1) , m_StepSize(0) , m_MaxLength(10000) , m_MinTractLength(20.0) , m_MaxTractLength(400.0) , m_SeedsPerVoxel(1) , m_AvoidStop(true) , m_RandomSampling(false) , m_SamplingDistance(-1) , m_DeflectionMod(1.0) , m_OnlyForwardSamples(true) , m_UseStopVotes(true) , m_NumberOfSamples(30) , m_NumPreviousDirections(1) , m_WmLabel(3) // mrtrix 5ttseg labels , m_GmLabel(1) // mrtrix 5ttseg labels , m_SeedOnlyGm(false) , m_ControlGmEndings(false) , m_MaxNumTracts(-1) , m_Verbose(true) , m_AposterioriCurvCheck(false) , m_DemoMode(false) , m_Random(true) , m_UseOutputProbabilityMap(false) , m_CurrentTracts(0) , m_Progress(0) , m_StopTracking(false) { this->SetNumberOfRequiredInputs(0); } std::string StreamlineTrackingFilter::GetStatusText() { std::string status = "Seedpoints processed: " + boost::lexical_cast(m_Progress) + "/" + boost::lexical_cast(m_SeedPoints.size()); if (m_SeedPoints.size()>0) status += " (" + boost::lexical_cast(100*m_Progress/m_SeedPoints.size()) + "%)"; if (m_MaxNumTracts>0) status += "\nFibers accepted: " + boost::lexical_cast(m_CurrentTracts) + "/" + boost::lexical_cast(m_MaxNumTracts); else status += "\nFibers accepted: " + boost::lexical_cast(m_CurrentTracts); return status; } void StreamlineTrackingFilter::BeforeTracking() { m_StopTracking = false; m_TrackingHandler->SetRandom(m_Random); 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]SetAngularThreshold(m_AngularThreshold); if (m_SamplingDistanceVoxGetNumberOfThreads(); i++) { PolyDataType poly = PolyDataType::New(); m_PolyDataContainer.push_back(poly); } 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()) { 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_TargetRegions.IsNull()) { m_TargetRegions = ItkUintImgType::New(); m_TargetRegions->SetSpacing( imageSpacing ); m_TargetRegions->SetOrigin( m_TrackingHandler->GetOrigin() ); m_TargetRegions->SetDirection( m_TrackingHandler->GetDirection() ); m_TargetRegions->SetRegions( m_TrackingHandler->GetLargestPossibleRegion() ); m_TargetRegions->Allocate(); m_TargetRegions->FillBuffer(1); } else std::cout << "StreamlineTracking - Using target 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_SeedPoints.empty()) GetSeedPointsFromSeedImage(); else m_SeedOnlyGm = false; 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 { if (m_ControlGmEndings) m_SeedOnlyGm = true; if (m_ControlGmEndings || m_SeedOnlyGm) std::cout << "StreamlineTracking - Using tissue image" << std::endl; else MITK_WARN << "StreamlineTracking - Tissue image set but no gray matter seeding or fiber endpoint-control enabled" << 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); 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; } 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() 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; } } 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; } bool StreamlineTrackingFilter -::IsValidPosition(const itk::Point &pos) +::IsInsideMask(const itk::Point &pos, ItkUcharImgType::Pointer mask) { + if (mask.IsNull()) + return true; + ItkUcharImgType::IndexType idx; - m_MaskImage->TransformPhysicalPointToIndex(pos, idx); - if (!m_MaskImage->GetLargestPossibleRegion().IsInside(idx) || m_MaskImage->GetPixel(idx)==0) + mask->TransformPhysicalPointToIndex(pos, idx); + if (!mask->GetLargestPossibleRegion().IsInside(idx) || mask->GetPixel(idx)==0) return false; return true; } bool StreamlineTrackingFilter ::IsInGm(const itk::Point &pos) { 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; return false; } std::vector< vnl_vector_fixed > StreamlineTrackingFilter::CreateDirections(int NPoints) { std::vector< vnl_vector_fixed > pointshell; if (NPoints<2) return pointshell; std::vector< float > theta; theta.resize(NPoints); std::vector< float > phi; phi.resize(NPoints); float C = sqrt(4*M_PI); phi[0] = 0.0; phi[NPoints-1] = 0.0; for(int i=0; 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); } 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); - ItkUcharImgType::IndexType idx; - m_MaskImage->TransformPhysicalPointToIndex(pos, idx); - - if (IsValidPosition(pos)) - { - if (m_StoppingRegions->GetPixel(idx)>0) - return direction; + if (IsInsideMask(pos, m_MaskImage) && !IsInsideMask(pos, m_StoppingRegions)) 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 (unsigned int i=0; i d; bool is_stop_voter = false; if (m_Random && m_RandomSampling) { d[0] = m_TrackingHandler->GetRandDouble(-0.5, 0.5); d[1] = m_TrackingHandler->GetRandDouble(-0.5, 0.5); d[2] = m_TrackingHandler->GetRandDouble(-0.5, 0.5); d.normalize(); d *= m_TrackingHandler->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]; vnl_vector_fixed tempDir; tempDir.fill(0.0); - if (IsValidPosition(sample_pos)) + if (IsInsideMask(sample_pos, m_MaskImage)) tempDir = m_TrackingHandler->ProposeDirection(sample_pos, olddirs, oldIndex); // sample neighborhood if (tempDir.magnitude()>mitk::eps) { direction += tempDir; 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)) + if (IsInsideMask(sample_pos, m_MaskImage)) 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 (direction.magnitude()>0.001 && (possible_stop_votes==0 || (float)stop_votes/possible_stop_votes<0.5) ) direction.normalize(); else direction.fill(0); 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; iTransformPhysicalPointToIndex(pos, oldIndex); + itk::Index<3> oldIndex; + m_TrackingHandler->WorldToIndex(pos, oldIndex); // 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 { 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*curv; while (curv>0) { if (front) fib->pop_front(); else fib->pop_back(); curv--; } return tractLength; } } if (tractLength>m_MaxTractLength) return tractLength; } if (m_DemoMode && !m_UseOutputProbabilityMap) // CHECK: warum sind die samplingpunkte der streamline in der visualisierung immer einen schritt voras? { #pragma omp critical { m_BuildFibersReady++; m_Tractogram.push_back(*fib); BuildFibers(true); m_Stop = true; 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); while (m_PauseTracking){} 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) { 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==(int)fib->size()-1) meanV += v; c--; } } meanV.normalize(); for (unsigned 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(); } void StreamlineTrackingFilter::GetSeedPointsFromSeedImage() { MITK_INFO << "StreamlineTracking - 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))) + if (IsInsideMask(worldPos, m_MaskImage) && (!m_SeedOnlyGm || IsInGm(worldPos))) { m_SeedPoints.push_back(worldPos); for (int s = 1; s < m_SeedsPerVoxel; s++) { start[0] = index[0] + m_TrackingHandler->GetRandDouble(-0.5, 0.5); start[1] = index[1] + m_TrackingHandler->GetRandDouble(-0.5, 0.5); start[2] = index[2] + m_TrackingHandler->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::random_shuffle(m_SeedPoints.begin(), m_SeedPoints.end()); m_CurrentTracts = 0; int num_seeds = m_SeedPoints.size(); itk::Index<3> zeroIndex; zeroIndex.Fill(0); m_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 || m_StopTracking) continue; else if (m_Verbose && i%print_interval==0) #pragma omp critical { m_Progress += print_interval; std::cout << " \r"; if (m_MaxNumTracts>0) std::cout << "Tried: " << m_Progress << "/" << num_seeds << " | Accepted: " << m_CurrentTracts << "/" << m_MaxNumTracts << '\r'; else std::cout << "Tried: " << m_Progress << "/" << num_seeds << " | Accepted: " << m_CurrentTracts << '\r'; cout.flush(); } 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() 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)) + if (IsInsideMask(worldPos, m_MaskImage)) dir = m_TrackingHandler->ProposeDirection(worldPos, olddirs, zeroIndex); 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) { CheckFiberForGmEnding(&fib); std::reverse(fib.begin(),fib.end()); CheckFiberForGmEnding(&fib); } } counter = fib.size(); if (tractLength>=m_MinTractLength && counter>=2) { ItkUintImgType::IndexType idx_begin, idx_end; m_TargetRegions->TransformPhysicalPointToIndex(fib.front(), idx_begin); m_TargetRegions->TransformPhysicalPointToIndex(fib.back(), idx_end); #pragma omp critical if ( !m_TargetRegions->GetLargestPossibleRegion().IsInside(idx_end) || !m_TargetRegions->GetLargestPossibleRegion().IsInside(idx_begin) || (m_TargetRegions->GetPixel(idx_begin)>0 && m_TargetRegions->GetPixel(idx_end)==m_TargetRegions->GetPixel(idx_begin)) ) { if (!m_StopTracking) { if (!m_UseOutputProbabilityMap) m_Tractogram.push_back(fib); else FiberToProbmap(&fib); m_CurrentTracts++; } if (m_MaxNumTracts > 0 && m_CurrentTracts>=static_cast(m_MaxNumTracts)) { if (!m_StopTracking) { std::cout << " \r"; MITK_INFO << "Reconstructed maximum number of tracts (" << m_CurrentTracts << "). Stopping tractography."; } m_StopTracking = true; } } } } } 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; } // 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(); } 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); if (idx != last_idx) { if (m_OutputProbabilityMap->GetLargestPossibleRegion().IsInside(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 (unsigned int i=0; i 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); } if (check) for (int i=0; iSetPoints(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); } else { itk::RescaleIntensityImageFilter< ItkDoubleImgType, ItkDoubleImgType >::Pointer filter = itk::RescaleIntensityImageFilter< ItkDoubleImgType, ItkDoubleImgType >::New(); filter->SetInput(m_OutputProbabilityMap); filter->SetOutputMaximum(1.0); filter->SetOutputMinimum(0.0); filter->Update(); m_OutputProbabilityMap = filter->GetOutput(); } 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(); } } diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.h index 3f2ac4f526..249a989a14 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.h @@ -1,224 +1,224 @@ /*=================================================================== 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 ItkUintImgType; 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( TargetRegions, ItkUintImgType::Pointer) ///< Only streamline starting and ending in this mask are retained 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, 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( StopTracking, bool ) ///< Use manually defined points in physical space as seed points instead of seed image void SetSeedPoints( const std::vector< itk::Point >& sP) { m_SeedPoints = sP; } void SetTrackingHandler( mitk::TrackingDataHandler* h ) ///< { m_TrackingHandler = h; } virtual void Update() override{ this->GenerateData(); } std::string GetStatusText(); 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 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 IsInsideMask(const itk::Point& pos, ItkUcharImgType::Pointer mask); ///< 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. 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; ItkUcharImgType::Pointer m_StoppingRegions; ItkUintImgType::Pointer m_TargetRegions; ItkUcharImgType::Pointer m_SeedImage; ItkUcharImgType::Pointer m_MaskImage; ItkUcharImgType::Pointer m_TissueImage; ItkDoubleImgType::Pointer m_OutputProbabilityMap; 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_AvoidStop; 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; bool m_Verbose; bool m_AposterioriCurvCheck; bool m_DemoMode; bool m_Random; bool m_UseOutputProbabilityMap; std::vector< itk::Point > m_SeedPoints; unsigned int m_CurrentTracts; unsigned int m_Progress; bool m_StopTracking; 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/FiberTracking/Algorithms/itkTractClusteringFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractClusteringFilter.cpp index 8c94149c0e..b7d9c8e655 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractClusteringFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractClusteringFilter.cpp @@ -1,538 +1,542 @@ /*=================================================================== 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 "itkTractClusteringFilter.h" #define _USE_MATH_DEFINES #include #include namespace itk{ TractClusteringFilter::TractClusteringFilter() : m_NumPoints(12) , m_InCentroids(nullptr) , m_MinClusterSize(1) , m_MaxClusters(0) , m_Metric(Metric::MDF) , m_ScalarMap(nullptr) , m_Scale(0) , m_MergeDuplicateThreshold(-1) { } TractClusteringFilter::~TractClusteringFilter() { } std::vector TractClusteringFilter::GetOutClusters() const { return m_OutClusters; } std::vector TractClusteringFilter::GetOutCentroids() const { return m_OutCentroids; } void TractClusteringFilter::SetMetric(const Metric &Metric) { m_Metric = Metric; } std::vector TractClusteringFilter::GetOutTractograms() const { return m_OutTractograms; } void TractClusteringFilter::SetDistances(const std::vector &Distances) { m_Distances = Distances; } float TractClusteringFilter::CalcMDF(vnl_matrix& s, vnl_matrix& t, bool& flipped) { float d_direct = 0; float d_flipped = 0; for (unsigned int i=0; id_flipped) { flipped = true; return d_flipped/m_NumPoints; } flipped = false; return d_direct/m_NumPoints; } float TractClusteringFilter::CalcMDF_STD(vnl_matrix& s, vnl_matrix& t, bool& flipped) { float d_direct = 0; float d_flipped = 0; vnl_vector dists_d; dists_d.set_size(m_NumPoints); vnl_vector dists_f; dists_f.set_size(m_NumPoints); for (unsigned int i=0; id_flipped) { float d = d_flipped/m_NumPoints; dists_f -= d; d += dists_f.magnitude(); flipped = true; return d/2; } float d = d_direct/m_NumPoints; dists_d -= d; d += dists_d.magnitude(); flipped = false; return d/2; } float TractClusteringFilter::CalcMAX_MDF(vnl_matrix& s, vnl_matrix& t, bool& flipped) { float d_direct = 0; float d_flipped = 0; vnl_vector dists_d; dists_d.set_size(m_NumPoints); vnl_vector dists_f; dists_f.set_size(m_NumPoints); for (unsigned int i=0; id_flipped) { float d = dists_f.max_value(); flipped = true; return d; } float d = dists_d.max_value(); flipped = false; return d; } std::vector > TractClusteringFilter::ResampleFibers(mitk::FiberBundle::Pointer tractogram) { mitk::FiberBundle::Pointer temp_fib = tractogram->GetDeepCopy(); temp_fib->ResampleToNumPoints(m_NumPoints); std::vector< vnl_matrix > out_fib; for (int i=0; iGetFiberPolyData()->GetNumberOfCells(); i++) { vtkCell* cell = temp_fib->GetFiberPolyData()->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vnl_matrix streamline; if (m_ScalarMap.IsNull()) streamline.set_size(3, m_NumPoints); else streamline.set_size(4, m_NumPoints); streamline.fill(0.0); for (int j=0; jGetPoint(j, cand); if (m_ScalarMap.IsNull()) { vnl_vector_fixed< float, 3 > candV; candV[0]=cand[0]; candV[1]=cand[1]; candV[2]=cand[2]; streamline.set_column(j, candV); } else { vnl_vector_fixed< float, 4 > candV; candV[0]=cand[0]; candV[1]=cand[1]; candV[2]=cand[2]; candV[3]=0; itk::Point wp; wp[0]=cand[0]; wp[1]=cand[1]; wp[2]=cand[2]; itk::Index<3> idx; m_ScalarMap->TransformPhysicalPointToIndex(wp, idx); if (m_ScalarMap->GetLargestPossibleRegion().IsInside(idx)) candV[3]=m_ScalarMap->GetPixel(idx)*m_Scale; streamline.set_column(j, candV); } } out_fib.push_back(streamline); } return out_fib; } std::vector< TractClusteringFilter::Cluster > TractClusteringFilter::ClusterStep(std::vector< long > f_indices, std::vector distances) { float dist_thres = distances.back(); distances.pop_back(); std::vector< Cluster > C; int N = f_indices.size(); Cluster c1; c1.I.push_back(f_indices.at(0)); c1.h = T[f_indices.at(0)]; c1.n = 1; C.push_back(c1); if (f_indices.size()==1) return C; for (int i=1; i t = T.at(f_indices.at(i)); int min_cluster_index = -1; float min_cluster_distance = 99999; bool flip = false; for (unsigned int k=0; k v = C.at(k).h / C.at(k).n; bool f = false; float d = 0; if (m_Metric==Metric::MDF) d = CalcMDF(t, v, f); else if (m_Metric==Metric::MDF_STD) d = CalcMDF_STD(t, v, f); else if (m_Metric==Metric::MAX_MDF) d = CalcMAX_MDF(t, v, f); if (d=0 && min_cluster_distance outC; #pragma omp parallel for for (int c=0; c<(int)C.size(); c++) { std::vector< Cluster > tempC = ClusterStep(C.at(c).I, distances); AppendCluster(outC, tempC); } return outC; } else return C; } void TractClusteringFilter::AppendCluster(std::vector< Cluster >& a, std::vector< Cluster >&b) { for (auto c : b) a.push_back(c); } void TractClusteringFilter::MergeDuplicateClusters(std::vector< TractClusteringFilter::Cluster >& clusters) { if (m_MergeDuplicateThreshold<0) m_MergeDuplicateThreshold = m_Distances.at(0); bool found = true; MITK_INFO << "Merging duplicate clusters with distance threshold " << m_MergeDuplicateThreshold; while (found && m_MergeDuplicateThreshold>mitk::eps) { std::cout << " \r"; std::cout << "Number of clusters: " << clusters.size() << '\r'; cout.flush(); found = false; for (int k1=0; k1<(int)clusters.size(); ++k1) { Cluster c1 = clusters.at(k1); vnl_matrix t = c1.h / c1.n; - int min_cluster_index = -1; - float min_cluster_distance = 99999; - bool flip = false; + std::vector< int > merge_indices; + std::vector< bool > flip_indices; #pragma omp parallel for for (int k2=0; k2<(int)clusters.size(); ++k2) { if (k1!=k2) { Cluster c2 = clusters.at(k2); vnl_matrix v = c2.h / c2.n; bool f = false; float d = 0; if (m_Metric==Metric::MDF) d = CalcMDF(t, v, f); else if (m_Metric==Metric::MDF_STD) d = CalcMDF_STD(t, v, f); else if (m_Metric==Metric::MAX_MDF) d = CalcMAX_MDF(t, v, f); #pragma omp critical - if (d=0 && min_cluster_distance TractClusteringFilter::AddToKnownClusters(std::vector< long > f_indices, std::vector >& centroids) { float dist_thres = m_Distances.at(0); int N = f_indices.size(); std::vector< Cluster > C; vnl_matrix zero_h; zero_h.set_size(T.at(0).rows(), T.at(0).cols()); zero_h.fill(0.0); Cluster no_fit; no_fit.h = zero_h; for (unsigned int i=0; i t = T.at(f_indices.at(i)); int min_cluster_index = -1; float min_cluster_distance = 99999; bool flip = false; int c_idx = 0; for (vnl_matrix centroid : centroids) { bool f = false; float d = 0; if (m_Metric==Metric::MDF) d = CalcMDF(t, centroid, f); else if (m_Metric==Metric::MDF_STD) d = CalcMDF_STD(t, centroid, f); else if (m_Metric==Metric::MAX_MDF) d = CalcMAX_MDF(t, centroid, f); if (d=0 && min_cluster_distance f_indices; for (unsigned int i=0; i clusters; if (m_InCentroids.IsNull()) { MITK_INFO << "Clustering fibers"; clusters = ClusterStep(f_indices, m_Distances); MITK_INFO << "Number of clusters: " << clusters.size(); MergeDuplicateClusters(clusters); std::sort(clusters.begin(),clusters.end()); } else { MITK_INFO << "Clustering with input centroids"; std::vector > centroids = ResampleFibers(m_InCentroids); clusters = AddToKnownClusters(f_indices, centroids); no_match = clusters.back(); clusters.pop_back(); MITK_INFO << "Number of clusters: " << clusters.size(); MergeDuplicateClusters(clusters); } MITK_INFO << "Clustering finished"; int max = clusters.size()-1; if (m_MaxClusters>0 && clusters.size()-1>m_MaxClusters) max = m_MaxClusters; for (int i=clusters.size()-1; i>=0; --i) { Cluster c = clusters.at(i); if ( c.n>=(int)m_MinClusterSize && !(m_MaxClusters>0 && clusters.size()-i>=m_MaxClusters) ) { m_OutClusters.push_back(c); vtkSmartPointer weights = vtkSmartPointer::New(); vtkSmartPointer pTmp = m_Tractogram->GeneratePolyDataByIds(c.I, weights); mitk::FiberBundle::Pointer fib = mitk::FiberBundle::New(pTmp); if (max>0) fib->SetFiberWeights((float)i/max); m_OutTractograms.push_back(fib); // create centroid vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); vtkSmartPointer polyData = vtkSmartPointer::New(); vtkSmartPointer container = vtkSmartPointer::New(); vnl_matrix centroid_points = c.h / c.n; for (unsigned int j=0; jInsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); polyData->SetPoints(vtkNewPoints); polyData->SetLines(vtkNewCells); mitk::FiberBundle::Pointer centroid = mitk::FiberBundle::New(polyData); centroid->SetFiberColors(255, 255, 255); m_OutCentroids.push_back(centroid); } } MITK_INFO << "Final number of clusters: " << m_OutTractograms.size(); int w = 0; for (auto fib : m_OutTractograms) { if (m_OutTractograms.size()>1) { fib->SetFiberWeights((float)w/(m_OutTractograms.size()-1)); m_OutCentroids.at(w)->SetFiberWeights((float)w/(m_OutTractograms.size()-1)); } else { fib->SetFiberWeights(1); m_OutCentroids.at(w)->SetFiberWeights(1); } fib->ColorFibersByFiberWeights(false, false); ++w; } if (no_match.n>0) { vtkSmartPointer weights = vtkSmartPointer::New(); vtkSmartPointer pTmp = m_Tractogram->GeneratePolyDataByIds(no_match.I, weights); mitk::FiberBundle::Pointer fib = mitk::FiberBundle::New(pTmp); fib->SetFiberColors(0, 0, 0); m_OutTractograms.push_back(fib); } } }