diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingDataHandler.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingDataHandler.cpp index 94f25ae860..bae0ea1b8d 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingDataHandler.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingDataHandler.cpp @@ -1,34 +1,27 @@ /*=================================================================== 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 "mitkTrackingDataHandler.h" namespace mitk { TrackingDataHandler::TrackingDataHandler() - : m_AngularThreshold(0.7) - , m_Interpolate(true) - , m_FlipX(false) - , m_FlipY(false) - , m_FlipZ(false) - , m_Mode(MODE::DETERMINISTIC) - , m_RngItk(ItkRngType::New()) + : m_RngItk(ItkRngType::New()) , m_NeedsDataInit(true) - , m_Random(true) { } } diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingDataHandler.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingDataHandler.h index 886a8792b3..c64b09d514 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingDataHandler.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingDataHandler.h @@ -1,123 +1,103 @@ /*=================================================================== 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 #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; + typedef mitk::StreamlineTractographyParameters::MODE MODE; 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() const { return m_Mode; } - - void SetInterpolate( bool interpolate ){ m_Interpolate = interpolate; } - bool GetInterpolate() const { return m_Interpolate; } - void SetAngularThreshold( float a ){ m_AngularThreshold = a; } - 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 ) + void SetParameters(mitk::StreamlineTractographyParameters* parameters) { - m_Random = random; - if (!random) + m_Parameters = parameters; + + if (m_Parameters->m_FixRandomSeed) { m_Rng.seed(0); std::srand(0); m_RngItk->SetSeed(0); } else { m_Rng.seed(); m_RngItk->SetSeed(); std::srand(std::time(nullptr)); } } double GetRandDouble(const double & a, const double & b) { return m_RngItk->GetUniformVariate(a, b); } - bool GetFlipX() const { return m_FlipX; } - bool GetFlipY() const { return m_FlipY; } - bool GetFlipZ() const { return m_FlipZ; } 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; + mitk::StreamlineTractographyParameters* m_Parameters; void DataModified() { m_NeedsDataInit = true; } }; } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.cpp index 2fed13346f..4cd88df053 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.cpp @@ -1,294 +1,291 @@ /*=================================================================== 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) + : m_NumProbSamples(1) , m_OdfFromTensor(false) { m_GfaInterpolator = itk::LinearInterpolateImageFunction< itk::Image< float, 3 >, float >::New(); m_OdfInterpolator = itk::LinearInterpolateImageFunction< itk::Image< ItkOdfImageType::PixelType, 3 >, float >::New(); } 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; } m_GfaInterpolator->SetInputImage(m_GfaImage); m_OdfInterpolator->SetInputImage(m_OdfImage); - std::cout << "TrackingHandlerOdf - GFA threshold: " << m_GfaThreshold << std::endl; - std::cout << "TrackingHandlerOdf - ODF threshold: " << m_OdfThreshold << std::endl; - if (m_SharpenOdfs) + std::cout << "TrackingHandlerOdf - GFA threshold: " << m_Parameters->m_Cutoff << std::endl; + std::cout << "TrackingHandlerOdf - ODF threshold: " << m_Parameters->m_OdfCutoff << std::endl; + if (m_Parameters->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) + if (probs[sampled_idx]>max_prob && probs[sampled_idx]>m_Parameters->m_OdfCutoff && fabs(angles[sampled_idx])>=m_Parameters->m_AngularThreshold) { max_prob = probs[sampled_idx]; max_sample_idx = sampled_idx; } - else if ( (probs[sampled_idx]<=m_OdfThreshold || fabs(angles[sampled_idx])m_OdfCutoff || 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 = mitk::imv::GetImageValue(pos, m_Interpolate, m_GfaInterpolator); - if (gfa(pos, m_Parameters->m_InterpolateTractographyData, m_GfaInterpolator); + if (gfam_Cutoff) return output_direction; vnl_vector_fixed last_dir; if (!olddirs.empty()) last_dir = olddirs.back(); - if (!m_Interpolate && oldIndex==idx) + if (!m_Parameters->m_InterpolateTractographyData && oldIndex==idx) return last_dir; - ItkOdfImageType::PixelType odf_values = mitk::imv::GetImageValue(pos, m_Interpolate, m_OdfInterpolator); + ItkOdfImageType::PixelType odf_values = mitk::imv::GetImageValue(pos, m_Parameters->m_InterpolateTractographyData, m_OdfInterpolator); 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]m_SharpenOdfs) { // sharpen ODF probs -= min_odf_val; probs /= (max_odf_val-min_odf_val); for (unsigned int i=0; i0) { 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 (max_odf_val>m_Parameters->m_OdfCutoff && (olddirs.empty() || last_dir.magnitude()<=0.5)) { - if (m_Mode==MODE::DETERMINISTIC) // return maximum peak + if (m_Parameters->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 + else if (m_Parameters->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) + else if (max_odf_val<=m_Parameters->m_OdfCutoff) // return (0,0,0) { return output_direction; } // correct previous direction - if (m_FlipX) + if (m_Parameters->m_FlipX) last_dir[0] *= -1; - if (m_FlipY) + if (m_Parameters->m_FlipY) last_dir[1] *= -1; - if (m_FlipZ) + if (m_Parameters->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) + if (m_Parameters->m_Mode==MODE::DETERMINISTIC && odf_val>max_prob && odf_val>m_Parameters->m_OdfCutoff) { // 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) + else if (m_Parameters->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) + if (m_Parameters->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 (aget) output_direction.fill(0); } else output_direction.fill(0); - if (m_FlipX) + if (m_Parameters->m_FlipX) output_direction[0] *= -1; - if (m_FlipY) + if (m_Parameters->m_FlipY) output_direction[1] *= -1; - if (m_FlipZ) + if (m_Parameters->m_FlipZ) output_direction[2] *= -1; return output_direction; } void TrackingHandlerOdf::SetNumProbSamples(int NumProbSamples) { m_NumProbSamples = NumProbSamples; } } diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.h index 593c02520e..4a70fd4f0c 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.h @@ -1,87 +1,80 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef _TrackingHandlerOdf #define _TrackingHandlerOdf #include "mitkTrackingDataHandler.h" #include #include #include #include namespace mitk { /** * \brief Enables streamline tracking on sampled ODF images. */ class MITKFIBERTRACKING_EXPORT TrackingHandlerOdf : public TrackingDataHandler { public: TrackingHandlerOdf(); ~TrackingHandlerOdf() override; typedef OdfImage::ItkOdfImageType ItkOdfImageType; typedef itk::Image< vnl_vector_fixed, 3> ItkPDImgType; void InitForTracking() override; ///< calls InputDataValidForTracking() and creates feature images vnl_vector_fixed ProposeDirection(const itk::Point& pos, std::deque< vnl_vector_fixed >& olddirs, itk::Index<3>& oldIndex) override; ///< predicts next progression direction at the given position bool WorldToIndex(itk::Point& pos, itk::Index<3>& index) override; - 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 ) override{ m_Mode = m; } ItkUcharImgType::SpacingType GetSpacing() override{ return m_OdfImage->GetSpacing(); } itk::Point GetOrigin() override{ return m_OdfImage->GetOrigin(); } ItkUcharImgType::DirectionType GetDirection() override{ return m_OdfImage->GetDirection(); } ItkUcharImgType::RegionType GetLargestPossibleRegion() override{ return m_OdfImage->GetLargestPossibleRegion(); } int OdfPower() const; void SetNumProbSamples(int NumProbSamples); bool GetIsOdfFromTensor() const; void SetIsOdfFromTensor(bool OdfFromTensor); 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; bool m_OdfFromTensor; itk::LinearInterpolateImageFunction< itk::Image< float, 3 >, float >::Pointer m_GfaInterpolator; itk::LinearInterpolateImageFunction< itk::Image< ItkOdfImageType::PixelType, 3 >, float >::Pointer m_OdfInterpolator; }; } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.cpp index 5243db3e71..5264b4b03c 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.cpp @@ -1,295 +1,294 @@ /*=================================================================== 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) + : 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; + std::cout << "TrackingHandlerPeaks - Peak threshold: " << m_Parameters->m_Cutoff << 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 (magm_FixRandomSeed) { // try m_NumDirs times to get a non-zero random direction for (int j=0; jGetIntegerVariate(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) + if (m_Parameters->m_FlipX) dir[0] *= -1; - if (m_FlipY) + if (m_Parameters->m_FlipY) dir[1] *= -1; - if (m_FlipZ) + if (m_Parameters->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; oldDir.fill(0.0); if (!olddirs.empty()) oldDir = olddirs.back(); float old_mag = oldDir.magnitude(); - if (!m_Interpolate && oldIndex==index) + if (!m_Parameters->m_InterpolateTractographyData && oldIndex==index) return oldDir; - output_direction = GetDirection(pos, m_Interpolate, oldDir); + output_direction = GetDirection(pos, m_Parameters->m_InterpolateTractographyData, oldDir); float mag = output_direction.magnitude(); - if (mag>=m_PeakThreshold) + if (mag>=m_Parameters->m_Cutoff) { 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 42c7288777..98da81d71c 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.h @@ -1,84 +1,75 @@ /*=================================================================== 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() override; typedef itk::Image< float, 4 > PeakImgType; ///< MRtrix peak image type void InitForTracking() override; ///< calls InputDataValidForTracking() and creates feature images vnl_vector_fixed ProposeDirection(const itk::Point& pos, std::deque< vnl_vector_fixed >& olddirs, itk::Index<3>& oldIndex) override; ///< predicts next progression direction at the given position bool WorldToIndex(itk::Point& pos, itk::Index<3>& index) override; - 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() override{ return spacing3; } itk::Point GetOrigin() override{ return origin3; } itk::Matrix GetDirection() override{ return direction3; } itk::ImageRegion<3> GetLargestPossibleRegion() override{ return imageRegion3; } - void SetMode( MODE m ) override - { - if (m==MODE::DETERMINISTIC) - m_Mode = m; - else - mitkThrow() << "Peak tracker is only implemented for deterministic mode."; - } protected: vnl_vector_fixed GetDirection(itk::Point itkP, bool interpolate, vnl_vector_fixed oldDir); vnl_vector_fixed GetMatchingDirection(itk::Index<3> idx3, vnl_vector_fixed& oldDir); vnl_vector_fixed GetDirection(itk::Index<3> idx3, int dirIdx); PeakImgType::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; ItkUcharImgType::Pointer m_DummyImage; bool m_ApplyDirectionMatrix; }; } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.cpp index dfeee974d3..3b383ecd00 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.cpp @@ -1,900 +1,900 @@ /*=================================================================== 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 > 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)); // initialize interpolators m_DwiFeatureImageInterpolator = DwiFeatureImageInterpolatorType::New(); m_DwiFeatureImageInterpolator->SetInputImage(m_DwiFeatureImages.at(0)); m_AdditionalFeatureImageInterpolators.clear(); for (auto afi_vec : m_AdditionalFeatureImages) { std::vector< FloatImageInterpolatorType::Pointer > v; for (auto img : afi_vec) { FloatImageInterpolatorType::Pointer interp = FloatImageInterpolatorType::New(); interp->SetInputImage(img); v.push_back(interp); } m_AdditionalFeatureImageInterpolators.push_back(v); } 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) + if (!m_Parameters->m_InterpolateTractographyData && 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()) ); featureData.init(0.0); - typename DwiFeatureImageType::PixelType dwiFeaturePixel = mitk::imv::GetImageValue< typename DwiFeatureImageType::PixelType >(pos, m_Interpolate, m_DwiFeatureImageInterpolator); + typename DwiFeatureImageType::PixelType dwiFeaturePixel = mitk::imv::GetImageValue< typename DwiFeatureImageType::PixelType >(pos, m_Parameters->m_InterpolateTractographyData, m_DwiFeatureImageInterpolator); 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) + if (m_Parameters->m_FlipX) tempD[0] *= -1; - if (m_FlipY) + if (m_Parameters->m_FlipY) tempD[1] *= -1; - if (m_FlipZ) + if (m_Parameters->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 interpolator : m_AdditionalFeatureImageInterpolators.at(0)) { float v = mitk::imv::GetImageValue(pos, false, interpolator); 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 (classLabelm_Mode==MODE::PROBABILISTIC) { probs2[classLabel] = probs(0,i); if (check_last_dir) probs2[classLabel] *= abs_angle; probs_sum += probs2[classLabel]; } - else if (m_Mode==MODE::DETERMINISTIC) + else if (m_Parameters->m_Mode==MODE::DETERMINISTIC) { vnl_vector_fixed d = m_DirectionContainer.at(classLabel); // get direction vector assiciated with the respective direction index if (check_last_dir) // do we have a previous streamline direction or did we just start? { 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) + if (m_Parameters->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) + if (m_Parameters->m_FlipX) output_direction[0] *= -1; - if (m_FlipY) + if (m_Parameters->m_FlipY) output_direction[1] *= -1; - if (m_FlipZ) + if (m_Parameters->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); } } // initialize interpolators m_AdditionalFeatureImageInterpolators.clear(); for (auto afi_vec : m_AdditionalFeatureImages) { std::vector< FloatImageInterpolatorType::Pointer > v; for (auto img : afi_vec) { FloatImageInterpolatorType::Pointer interp = FloatImageInterpolatorType::New(); interp->SetInputImage(img); v.push_back(interp); } m_AdditionalFeatureImageInterpolators.push_back(v); } 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; tSetInputImage(fiber_folume); typename DwiFeatureImageType::Pointer image = m_DwiFeatureImages.at(t); typename DwiFeatureImageInterpolatorType::Pointer dwi_interp = DwiFeatureImageInterpolatorType::New(); dwi_interp->SetInputImage(image); ItkUcharImgType::Pointer wmMask = m_WhiteMatterImages.at(t); ItkUcharImgType::Pointer mask; if (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 = mitk::imv::GetImageValue(itkP, false, interpolator); 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 (unsigned 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 = mitk::imv::GetImageValue(itkP1, false, volume_interpolator); // diffusion signal features - typename DwiFeatureImageType::PixelType pix = mitk::imv::GetImageValue< typename DwiFeatureImageType::PixelType >(itkP1, m_Interpolate, dwi_interp); + typename DwiFeatureImageType::PixelType pix = mitk::imv::GetImageValue< typename DwiFeatureImageType::PixelType >(itkP1, m_Parameters->m_InterpolateTractographyData, dwi_interp); for (unsigned int f=0; f(itkP1, false, interpolator); 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 5b911770d4..4e28ed5e2c 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.h @@ -1,166 +1,165 @@ /*=================================================================== 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() override; typedef itk::Image< itk::Vector< float, NumberOfSignalFeatures > , 3 > DwiFeatureImageType; typedef itk::LinearInterpolateImageFunction< DwiFeatureImageType, float > DwiFeatureImageInterpolatorType; typedef itk::LinearInterpolateImageFunction< ItkFloatImgType, float > FloatImageInterpolatorType; 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 ) override{ m_Mode = m; } void SetForest(mitk::TractographyForest::Pointer forest){ m_Forest = forest; } void SetMaxNumWmSamples(int num){ m_MaxNumWmSamples=num; } void SetNumPreviousDirections( unsigned 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() override; ///< calls InputDataValidForTracking() and creates feature images vnl_vector_fixed ProposeDirection(const itk::Point& pos, std::deque< vnl_vector_fixed >& olddirs, itk::Index<3>& oldIndex) override; ///< predicts next progression direction at the given position bool WorldToIndex(itk::Point& pos, itk::Index<3>& index) override; 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() override{ return m_DwiFeatureImages.at(0)->GetSpacing(); } itk::Point GetOrigin() override{ return m_DwiFeatureImages.at(0)->GetOrigin(); } ItkUcharImgType::DirectionType GetDirection() override{ return m_DwiFeatureImages.at(0)->GetDirection(); } ItkUcharImgType::RegionType GetLargestPossibleRegion() override{ 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 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; typename DwiFeatureImageInterpolatorType::Pointer m_DwiFeatureImageInterpolator; std::vector< std::vector< FloatImageInterpolatorType::Pointer > > m_AdditionalFeatureImageInterpolators; }; } #include "mitkTrackingHandlerRandomForest.cpp" #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.cpp index 0581e41ec7..97a2d13e2b 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.cpp @@ -1,371 +1,368 @@ /*=================================================================== 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_InterpolateTensors(true) , m_NumberOfInputs(0) { m_FaInterpolator = itk::LinearInterpolateImageFunction< itk::Image< float, 3 >, float >::New(); } 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) + if (m_Parameters->m_F+m_Parameters->m_G>1.0) { - float temp = m_F+m_G; - m_F /= temp; - m_G /= temp; + float temp = m_Parameters->m_F+m_Parameters->m_G; + m_Parameters->m_F /= temp; + m_Parameters->m_G /= temp; } m_FaInterpolator->SetInputImage(m_FaImage); - 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; + std::cout << "TrackingHandlerTensor - FA threshold: " << m_Parameters->m_Cutoff << std::endl; + std::cout << "TrackingHandlerTensor - f: " << m_Parameters->m_F << std::endl; + std::cout << "TrackingHandlerTensor - g: " << m_Parameters->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) + if (!m_Parameters->m_InterpolateTractographyData) { 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 = mitk::imv::GetImageValue(pos, m_Interpolate, m_FaInterpolator); - if (fa(pos, m_Parameters->m_InterpolateTractographyData, m_FaInterpolator); + if (fam_Cutoff) return output_direction; vnl_vector_fixed oldDir; oldDir.fill(0.0); if (!olddirs.empty()) oldDir = olddirs.back(); - if (m_FlipX) + if (m_Parameters->m_FlipX) oldDir[0] *= -1; - if (m_FlipY) + if (m_Parameters->m_FlipY) oldDir[1] *= -1; - if (m_FlipZ) + if (m_Parameters->m_FlipZ) oldDir[2] *= -1; float old_mag = oldDir.magnitude(); - if (!m_Interpolate && oldIndex==index) + if (!m_Parameters->m_InterpolateTractographyData && 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 + if (old_mag>0.5 && m_Parameters->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[0] = m_Parameters->m_F*output_direction[0] + (1-m_Parameters->m_F)*( (1-m_Parameters->m_G)*oldDir[0] + m_Parameters->m_G*(tensor[0]*oldDir[0] + tensor[1]*oldDir[1] + tensor[2]*oldDir[2])); + output_direction[1] = m_Parameters->m_F*output_direction[1] + (1-m_Parameters->m_F)*( (1-m_Parameters->m_G)*oldDir[1] + m_Parameters->m_G*(tensor[1]*oldDir[0] + tensor[3]*oldDir[1] + tensor[4]*oldDir[2])); + output_direction[2] = m_Parameters->m_F*output_direction[2] + (1-m_Parameters->m_F)*( (1-m_Parameters->m_G)*oldDir[2] + m_Parameters->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) + if (a>=m_Parameters->GetAngularThreshold()) output_direction *= mag; else output_direction.fill(0); } else output_direction.fill(0); } catch(...) { } - if (m_FlipX) + if (m_Parameters->m_FlipX) output_direction[0] *= -1; - if (m_FlipY) + if (m_Parameters->m_FlipY) output_direction[1] *= -1; - if (m_FlipZ) + if (m_Parameters->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 5884ee8be1..6175bec856 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.h @@ -1,94 +1,80 @@ /*=================================================================== 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 #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() override; typedef TensorImage::PixelType TensorType; typedef TensorImage::ItkTensorImageType ItkTensorImageType; typedef itk::Image< vnl_vector_fixed, 3> ItkPDImgType; void InitForTracking() override; ///< calls InputDataValidForTracking() and creates feature images vnl_vector_fixed ProposeDirection(const itk::Point& pos, std::deque< vnl_vector_fixed >& olddirs, itk::Index<3>& oldIndex) override; ///< predicts next progression direction at the given position bool WorldToIndex(itk::Point& pos, itk::Index<3>& index) override; - 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 ) override - { - if (m==MODE::DETERMINISTIC) - m_Mode = m; - else - mitkThrow() << "Tensor tracker is only implemented for deterministic mode."; - } int GetNumTensorImages() const { return m_TensorImages.size(); } ItkUcharImgType::SpacingType GetSpacing() override{ return m_FaImage->GetSpacing(); } itk::Point GetOrigin() override{ return m_FaImage->GetOrigin(); } ItkUcharImgType::DirectionType GetDirection() override{ return m_FaImage->GetDirection(); } ItkUcharImgType::RegionType GetLargestPossibleRegion() override{ return m_FaImage->GetLargestPossibleRegion(); } protected: vnl_vector_fixed GetMatchingDirection(itk::Index<3> idx, vnl_vector_fixed& oldDir, int& image_num); vnl_vector_fixed GetDirection(itk::Point itkP, vnl_vector_fixed oldDir, TensorType& tensor); vnl_vector_fixed GetLargestEigenvector(TensorType& tensor); - float m_FaThreshold; - float m_F; - float m_G; - std::vector< ItkDoubleImgType::Pointer > m_EmaxImage; ///< Stores largest eigenvalues per voxel (one for each tensor) ItkFloatImgType::Pointer m_FaImage; ///< FA image used to determine streamline termination. std::vector< ItkPDImgType::Pointer > m_PdImage; ///< Stores principal direction of each tensor in each voxel. std::vector< ItkTensorImageType::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; itk::LinearInterpolateImageFunction< itk::Image< float, 3 >, float >::Pointer m_FaInterpolator; }; } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp index 7d76d2624c..b488f7f1e6 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp @@ -1,1025 +1,983 @@ /*=================================================================== 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 "itkPointShell.h" #include #include #include #include #include #include #include #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_ExclusionRegions(nullptr) , m_OutputProbabilityMap(nullptr) , m_MinVoxelSize(-1) - , 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_MaxNumTracts(-1) , m_Verbose(true) - , m_LoopCheck(-1) , m_DemoMode(false) - , m_Random(true) - , m_UseOutputProbabilityMap(false) , m_CurrentTracts(0) , m_Progress(0) , m_StopTracking(false) - , m_InterpolateMasks(true) - , m_TrialsPerSeed(10) - , m_EndpointConstraint(EndpointConstraints::NONE) - , m_IntroduceDirectionsFromPrior(true) - , m_TrackingPriorAsMask(true) - , m_TrackingPriorWeight(1.0) , m_TrackingPriorHandler(nullptr) { 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) + if (m_Parameters->m_MaxNumFibers>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->SetParameters(m_Parameters); 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(); if(imageSpacing[0](imageSpacing[0]); else if (imageSpacing[1] < imageSpacing[2]) m_MinVoxelSize = static_cast(imageSpacing[1]); else m_MinVoxelSize = static_cast(imageSpacing[2]); - if (m_StepSizeVox(mitk::eps)) - m_StepSize = 0.5f*m_MinVoxelSize; - else - m_StepSize = m_StepSizeVox*m_MinVoxelSize; - - if (m_AngularThresholdDeg<0) - { - if (m_StepSize/m_MinVoxelSize<=0.966f) // minimum 15° for automatic estimation - m_AngularThreshold = static_cast(std::cos( 0.5 * itk::Math::pi * static_cast(m_StepSize/m_MinVoxelSize) )); - else - m_AngularThreshold = static_cast(std::cos( 0.5 * itk::Math::pi * 0.966 )); - } - else - m_AngularThreshold = static_cast(std::cos( static_cast(m_AngularThresholdDeg)*itk::Math::pi/180.0 )); - m_TrackingHandler->SetAngularThreshold(m_AngularThreshold); - if (m_TrackingPriorHandler!=nullptr) { m_TrackingPriorHandler->SetRandom(m_Random); m_TrackingPriorHandler->InitForTracking(); m_TrackingPriorHandler->SetAngularThreshold(m_AngularThreshold); } if (m_SamplingDistanceVox(mitk::eps)) m_SamplingDistance = m_MinVoxelSize*0.25f; else m_SamplingDistance = m_SamplingDistanceVox * m_MinVoxelSize; m_PolyDataContainer.clear(); for (unsigned int i=0; iGetNumberOfThreads(); 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); } m_MaskInterpolator = itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::New(); m_StopInterpolator = itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::New(); m_SeedInterpolator = itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::New(); m_TargetInterpolator = itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::New(); m_ExclusionInterpolator = itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::New(); if (m_StoppingRegions.IsNull()) { m_StoppingRegions = ItkFloatImgType::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; m_StopInterpolator->SetInputImage(m_StoppingRegions); if (m_ExclusionRegions.IsNotNull()) { std::cout << "StreamlineTracking - Using exclusion region image" << std::endl; m_ExclusionInterpolator->SetInputImage(m_ExclusionRegions); } if (m_TargetRegions.IsNull()) { m_TargetImageSet = false; m_TargetRegions = ItkFloatImgType::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 { m_TargetImageSet = true; m_TargetInterpolator->SetInputImage(m_TargetRegions); std::cout << "StreamlineTracking - Using target region image" << std::endl; } if (m_SeedImage.IsNull()) { m_SeedImageSet = false; m_SeedImage = ItkFloatImgType::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 { m_SeedImageSet = true; std::cout << "StreamlineTracking - Using seed image" << std::endl; } m_SeedInterpolator->SetInputImage(m_SeedImage); if (m_MaskImage.IsNull()) { // initialize mask image m_MaskImage = ItkFloatImgType::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; m_MaskInterpolator->SetInputImage(m_MaskImage); // Autosettings for endpoint constraints if (m_EndpointConstraint==EndpointConstraints::NONE && m_TargetImageSet && m_SeedImageSet) { MITK_INFO << "No endpoint constraint chosen but seed and target image set --> setting constraint to EPS_IN_SEED_AND_TARGET"; m_EndpointConstraint = EndpointConstraints::EPS_IN_SEED_AND_TARGET; } else if (m_EndpointConstraint==EndpointConstraints::NONE && m_TargetImageSet) { MITK_INFO << "No endpoint constraint chosen but target image set --> setting constraint to EPS_IN_TARGET"; m_EndpointConstraint = EndpointConstraints::EPS_IN_TARGET; } // Check if endpoint constraints are valid FiberType test_fib; itk::Point p; p.Fill(0); test_fib.push_back(p); test_fib.push_back(p); IsValidFiber(&test_fib); if (m_SeedPoints.empty()) GetSeedPointsFromSeedImage(); 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_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; std::cout << "StreamlineTracking - Trials per seed: " << m_TrialsPerSeed << std::endl; } else std::cout << "StreamlineTracking - Mode: ???" << std::endl; if (m_EndpointConstraint==EndpointConstraints::NONE) std::cout << "StreamlineTracking - Endpoint constraint: NONE" << std::endl; else if (m_EndpointConstraint==EndpointConstraints::EPS_IN_TARGET) std::cout << "StreamlineTracking - Endpoint constraint: EPS_IN_TARGET" << std::endl; else if (m_EndpointConstraint==EndpointConstraints::EPS_IN_TARGET_LABELDIFF) std::cout << "StreamlineTracking - Endpoint constraint: EPS_IN_TARGET_LABELDIFF" << std::endl; else if (m_EndpointConstraint==EndpointConstraints::EPS_IN_SEED_AND_TARGET) std::cout << "StreamlineTracking - Endpoint constraint: EPS_IN_SEED_AND_TARGET" << std::endl; else if (m_EndpointConstraint==EndpointConstraints::MIN_ONE_EP_IN_TARGET) std::cout << "StreamlineTracking - Endpoint constraint: MIN_ONE_EP_IN_TARGET" << std::endl; else if (m_EndpointConstraint==EndpointConstraints::ONE_EP_IN_TARGET) std::cout << "StreamlineTracking - Endpoint constraint: ONE_EP_IN_TARGET" << std::endl; else if (m_EndpointConstraint==EndpointConstraints::NO_EP_IN_TARGET) std::cout << "StreamlineTracking - Endpoint constraint: NO_EP_IN_TARGET" << std::endl; std::cout << "StreamlineTracking - Angular threshold: " << m_AngularThreshold << " (" << 180*std::acos( static_cast(m_AngularThreshold) )/itk::Math::pi << "°)" << std::endl; std::cout << "StreamlineTracking - Stepsize: " << m_StepSize << "mm (" << m_StepSize/m_MinVoxelSize << "*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 - Loop check: " << m_LoopCheck << "°" << std::endl; std::cout << "StreamlineTracking - Num. neighborhood samples: " << m_NumberOfSamples << std::endl; std::cout << "StreamlineTracking - Max. sampling distance: " << m_SamplingDistance << "mm (" << m_SamplingDistance/m_MinVoxelSize << "*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_TrackingPriorHandler!=nullptr) std::cout << "StreamlineTracking - Using directional prior for tractography (w=" << m_TrackingPriorWeight << ")" << 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::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; } std::vector< vnl_vector_fixed > StreamlineTrackingFilter::CreateDirections(unsigned int NPoints) { std::vector< vnl_vector_fixed > pointshell; if (NPoints<2) return pointshell; std::vector< double > theta; theta.resize(NPoints); std::vector< double > phi; phi.resize(NPoints); auto C = sqrt(4*itk::Math::pi); phi[0] = 0.0; phi[NPoints-1] = 0.0; for(unsigned int i=0; i0 && i d; d[0] = static_cast(cos(theta[i]) * cos(phi[i])); d[1] = static_cast(cos(theta[i]) * sin(phi[i])); d[2] = static_cast(sin(theta[i])); pointshell.push_back(d); } return pointshell; } vnl_vector_fixed StreamlineTrackingFilter::GetNewDirection(const itk::Point &pos, std::deque >& olddirs, itk::Index<3> &oldIndex) { if (m_DemoMode) { m_SamplingPointset->Clear(); m_AlternativePointset->Clear(); m_StopVotePointset->Clear(); } vnl_vector_fixed direction; direction.fill(0); if (mitk::imv::IsInsideMask(pos, m_InterpolateMasks, m_MaskInterpolator) && !mitk::imv::IsInsideMask(pos, m_InterpolateMasks, m_StopInterpolator)) direction = m_TrackingHandler->ProposeDirection(pos, olddirs, oldIndex); // get direction proposal at current streamline position else return direction; int stop_votes = 0; int possible_stop_votes = 0; if (!olddirs.empty()) { vnl_vector_fixed olddir = olddirs.back(); std::vector< vnl_vector_fixed > probeVecs = CreateDirections(m_NumberOfSamples); itk::Point sample_pos; unsigned int alternatives = 1; for (unsigned int i=0; i d; bool is_stop_voter = false; if (m_Random && m_RandomSampling) { d[0] = static_cast(m_TrackingHandler->GetRandDouble(-0.5, 0.5)); d[1] = static_cast(m_TrackingHandler->GetRandDouble(-0.5, 0.5)); d[2] = static_cast(m_TrackingHandler->GetRandDouble(-0.5, 0.5)); d.normalize(); d *= static_cast(m_TrackingHandler->GetRandDouble(0, static_cast(m_SamplingDistance))); } else { d = probeVecs.at(i); float dot = dot_product(d, olddir); if (m_UseStopVotes && dot>0.7f) { 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 (mitk::imv::IsInsideMask(sample_pos, m_InterpolateMasks, m_MaskInterpolator)) tempDir = m_TrackingHandler->ProposeDirection(sample_pos, olddirs, oldIndex); // sample neighborhood if (tempDir.magnitude()>static_cast(mitk::eps)) { direction += tempDir; if(m_DemoMode) m_SamplingPointset->InsertPoint(i, sample_pos); } else if (m_AvoidStop && olddir.magnitude()>0.5f) // 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.0f) // 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 (mitk::imv::IsInsideMask(sample_pos, m_InterpolateMasks, m_MaskInterpolator)) tempDir = m_TrackingHandler->ProposeDirection(sample_pos, olddirs, oldIndex); // sample neighborhood if (tempDir.magnitude()>static_cast(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++; } } } bool valid = false; if (direction.magnitude()>0.001f && (possible_stop_votes==0 || static_cast(stop_votes)/possible_stop_votes<0.5f) ) { direction.normalize(); valid = true; } else direction.fill(0); if (m_TrackingPriorHandler!=nullptr && (m_IntroduceDirectionsFromPrior || valid)) { vnl_vector_fixed prior = m_TrackingPriorHandler->ProposeDirection(pos, olddirs, oldIndex); if (prior.magnitude()>0.001f) { prior.normalize(); if (dot_product(prior,direction)<0) prior *= -1; direction = (1.0f-m_TrackingPriorWeight) * direction + m_TrackingPriorWeight * prior; direction.normalize(); } else if (m_TrackingPriorAsMask) direction.fill(0.0); } return direction; } float StreamlineTrackingFilter::FollowStreamline(itk::Point pos, vnl_vector_fixed dir, FiberType* fib, DirectionContainer* container, float tractLength, bool front, bool &exclude) { vnl_vector_fixed zero_dir; zero_dir.fill(0.0); std::deque< vnl_vector_fixed > last_dirs; for (unsigned int i=0; i oldIndex; m_TrackingHandler->WorldToIndex(pos, oldIndex); // get new position CalculateNewPosition(pos, dir); if (m_ExclusionRegions.IsNotNull() && mitk::imv::IsInsideMask(pos, m_InterpolateMasks, m_ExclusionInterpolator)) { exclude = true; return tractLength; } if (m_AbortTracking) return tractLength; // if yes, add new point to streamline dir.normalize(); if (front) { fib->push_front(pos); container->push_front(dir); } else { fib->push_back(pos); container->push_back(dir); } tractLength += m_StepSize; if (m_LoopCheck>=0 && CheckCurvature(container, front)>m_LoopCheck) 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){ } } } 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.0001f) return tractLength; } return tractLength; } float StreamlineTrackingFilter::CheckCurvature(DirectionContainer* fib, bool front) { if (fib->size()<8) return 0; float m_Distance = std::max(m_MinVoxelSize*4, m_StepSize*8); 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(dist(fib->size())-1) { dist += m_StepSize; vnl_vector_fixed< float, 3 > v = fib->at(static_cast(c)); if (dot_product(v,meanV)<0) v = -v; vectors.push_back(v); meanV += v; c++; } } else { int c = static_cast(fib->size())-1; while(dist=0) { dist += m_StepSize; vnl_vector_fixed< float, 3 > v = fib->at(static_cast(c)); if (dot_product(v,meanV)<0) v = -v; vectors.push_back(v); meanV += v; c--; } } meanV.normalize(); for (unsigned int c=0; c1.0f) angle = 1.0; dev += acos(angle)*180.0f/static_cast(itk::Math::pi); } if (vectors.size()>0) dev /= vectors.size(); return dev; } void StreamlineTrackingFilter::SetTrackingPriorHandler(mitk::TrackingDataHandler *TrackingPriorHandler) { m_TrackingPriorHandler = TrackingPriorHandler; } void StreamlineTrackingFilter::GetSeedPointsFromSeedImage() { MITK_INFO << "StreamlineTracking - Calculating seed points."; m_SeedPoints.clear(); typedef ImageRegionConstIterator< ItkFloatImgType > MaskIteratorType; MaskIteratorType sit(m_SeedImage, m_SeedImage->GetLargestPossibleRegion()); sit.GoToBegin(); while (!sit.IsAtEnd()) { if (sit.Value()>0) { ItkFloatImgType::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 ( mitk::imv::IsInsideMask(worldPos, m_InterpolateMasks, m_MaskInterpolator) ) { m_SeedPoints.push_back(worldPos); for (int s = 1; s < m_SeedsPerVoxel; s++) { start[0] = index[0] + static_cast(m_TrackingHandler->GetRandDouble(-0.5, 0.5)); start[1] = index[1] + static_cast(m_TrackingHandler->GetRandDouble(-0.5, 0.5)); start[2] = index[2] + static_cast(m_TrackingHandler->GetRandDouble(-0.5, 0.5)); itk::Point worldPos; m_SeedImage->TransformContinuousIndexToPhysicalPoint(start, worldPos); m_SeedPoints.push_back(worldPos); } } } ++sit; } if (m_SeedPoints.empty()) mitkThrow() << "No valid seed point in seed image! Is your seed image registered with the image you are tracking on?"; } void StreamlineTrackingFilter::GenerateData() { this->BeforeTracking(); if (m_Random) std::random_shuffle(m_SeedPoints.begin(), m_SeedPoints.end()); m_CurrentTracts = 0; int num_seeds = static_cast(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 += static_cast(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(static_cast(temp_i)); for (unsigned int trials=0; trials dir; dir.fill(0.0); std::deque< vnl_vector_fixed > olddirs; dir = GetNewDirection(worldPos, olddirs, zeroIndex) * 0.5f; bool exclude = false; if (m_ExclusionRegions.IsNotNull() && mitk::imv::IsInsideMask(worldPos, m_InterpolateMasks, m_ExclusionInterpolator)) exclude = true; bool success = false; if (dir.magnitude()>0.0001f && !exclude) { // forward tracking tractLength = FollowStreamline(worldPos, dir, &fib, &direction_container, 0, false, exclude); fib.push_front(worldPos); // backward tracking if (!exclude) tractLength = FollowStreamline(worldPos, -dir, &fib, &direction_container, tractLength, true, exclude); counter = fib.size(); if (tractLength>=m_MinTractLength && counter>=2 && !exclude) { #pragma omp critical if ( IsValidFiber(&fib) ) { if (!m_StopTracking) { if (!m_UseOutputProbabilityMap) m_Tractogram.push_back(fib); else FiberToProbmap(&fib); m_CurrentTracts++; success = true; } 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; } } } } - if (success || m_TrackingHandler->GetMode()!=mitk::TrackingDataHandler::PROBABILISTIC) + if (success || m_TrackingHandler->GetMode()!=MODE::PROBABILISTIC) break; // we only try one seed point multiple times if we use a probabilistic tracker and have not found a valid streamline yet }// trials per seed }// seed points this->AfterTracking(); } bool StreamlineTrackingFilter::IsValidFiber(FiberType* fib) { if (m_EndpointConstraint==EndpointConstraints::NONE) { return true; } else if (m_EndpointConstraint==EndpointConstraints::EPS_IN_TARGET) { if (m_TargetImageSet) { if ( mitk::imv::IsInsideMask(fib->front(), m_InterpolateMasks, m_TargetInterpolator) && mitk::imv::IsInsideMask(fib->back(), m_InterpolateMasks, m_TargetInterpolator) ) return true; return false; } else mitkThrow() << "No target image set but endpoint constraint EPS_IN_TARGET chosen!"; } else if (m_EndpointConstraint==EndpointConstraints::EPS_IN_TARGET_LABELDIFF) { if (m_TargetImageSet) { float v1 = mitk::imv::GetImageValue(fib->front(), false, m_TargetInterpolator); float v2 = mitk::imv::GetImageValue(fib->back(), false, m_TargetInterpolator); if ( v1>0.0f && v2>0.0f && v1!=v2 ) return true; return false; } else mitkThrow() << "No target image set but endpoint constraint EPS_IN_TARGET_LABELDIFF chosen!"; } else if (m_EndpointConstraint==EndpointConstraints::EPS_IN_SEED_AND_TARGET) { if (m_TargetImageSet && m_SeedImageSet) { if ( mitk::imv::IsInsideMask(fib->front(), m_InterpolateMasks, m_SeedInterpolator) && mitk::imv::IsInsideMask(fib->back(), m_InterpolateMasks, m_TargetInterpolator) ) return true; if ( mitk::imv::IsInsideMask(fib->back(), m_InterpolateMasks, m_SeedInterpolator) && mitk::imv::IsInsideMask(fib->front(), m_InterpolateMasks, m_TargetInterpolator) ) return true; return false; } else mitkThrow() << "No target or seed image set but endpoint constraint EPS_IN_SEED_AND_TARGET chosen!"; } else if (m_EndpointConstraint==EndpointConstraints::MIN_ONE_EP_IN_TARGET) { if (m_TargetImageSet) { if ( mitk::imv::IsInsideMask(fib->front(), m_InterpolateMasks, m_TargetInterpolator) || mitk::imv::IsInsideMask(fib->back(), m_InterpolateMasks, m_TargetInterpolator) ) return true; return false; } else mitkThrow() << "No target image set but endpoint constraint MIN_ONE_EP_IN_TARGET chosen!"; } else if (m_EndpointConstraint==EndpointConstraints::ONE_EP_IN_TARGET) { if (m_TargetImageSet) { if ( mitk::imv::IsInsideMask(fib->front(), m_InterpolateMasks, m_TargetInterpolator) && !mitk::imv::IsInsideMask(fib->back(), m_InterpolateMasks, m_TargetInterpolator) ) return true; if ( !mitk::imv::IsInsideMask(fib->back(), m_InterpolateMasks, m_TargetInterpolator) && mitk::imv::IsInsideMask(fib->front(), m_InterpolateMasks, m_TargetInterpolator) ) return true; return false; } else mitkThrow() << "No target image set but endpoint constraint ONE_EP_IN_TARGET chosen!"; } else if (m_EndpointConstraint==EndpointConstraints::NO_EP_IN_TARGET) { if (m_TargetImageSet) { if ( mitk::imv::IsInsideMask(fib->front(), m_InterpolateMasks, m_TargetInterpolator) || mitk::imv::IsInsideMask(fib->back(), m_InterpolateMasks, m_TargetInterpolator) ) return false; return true; } else mitkThrow() << "No target image set but endpoint constraint NO_EP_IN_TARGET chosen!"; } return true; } 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(); } void StreamlineTrackingFilter::SetDicomProperties(mitk::FiberBundle::Pointer fib) { std::string model_code_value = "-"; std::string model_code_meaning = "-"; std::string algo_code_value = "-"; std::string algo_code_meaning = "-"; - if (m_TrackingHandler->GetMode()==mitk::TrackingDataHandler::DETERMINISTIC && dynamic_cast(m_TrackingHandler) && !m_TrackingHandler->GetInterpolate()) + if (m_TrackingHandler->GetMode()==MODE::DETERMINISTIC && dynamic_cast(m_TrackingHandler) && !m_TrackingHandler->GetInterpolate()) { algo_code_value = "sup181_ee04"; algo_code_meaning = "FACT"; } - else if (m_TrackingHandler->GetMode()==mitk::TrackingDataHandler::DETERMINISTIC) + else if (m_TrackingHandler->GetMode()==MODE::DETERMINISTIC) { algo_code_value = "sup181_ee01"; algo_code_meaning = "Deterministic"; } - else if (m_TrackingHandler->GetMode()==mitk::TrackingDataHandler::PROBABILISTIC) + else if (m_TrackingHandler->GetMode()==MODE::PROBABILISTIC) { algo_code_value = "sup181_ee02"; algo_code_meaning = "Probabilistic"; } if (dynamic_cast(m_TrackingHandler) || (dynamic_cast(m_TrackingHandler) && dynamic_cast(m_TrackingHandler)->GetIsOdfFromTensor() ) ) { if ( dynamic_cast(m_TrackingHandler) && dynamic_cast(m_TrackingHandler)->GetNumTensorImages()>1 ) { model_code_value = "sup181_bb02"; model_code_meaning = "Multi Tensor"; } else { model_code_value = "sup181_bb01"; model_code_meaning = "Single Tensor"; } } else if (dynamic_cast*>(m_TrackingHandler) || dynamic_cast*>(m_TrackingHandler)) { model_code_value = "sup181_bb03"; model_code_meaning = "Model Free"; } else if (dynamic_cast(m_TrackingHandler)) { model_code_value = "-"; model_code_meaning = "ODF"; } else if (dynamic_cast(m_TrackingHandler)) { model_code_value = "-"; model_code_meaning = "Peaks"; } fib->SetProperty("DICOM.anatomy.value", mitk::StringProperty::New("T-A0095")); fib->SetProperty("DICOM.anatomy.meaning", mitk::StringProperty::New("White matter of brain and spinal cord")); fib->SetProperty("DICOM.algo_code.value", mitk::StringProperty::New(algo_code_value)); fib->SetProperty("DICOM.algo_code.meaning", mitk::StringProperty::New(algo_code_meaning)); fib->SetProperty("DICOM.model_code.value", mitk::StringProperty::New(model_code_value)); fib->SetProperty("DICOM.model_code.meaning", mitk::StringProperty::New(model_code_meaning)); } } diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.h index 3817102215..004c7c7524 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.h @@ -1,258 +1,194 @@ /*=================================================================== 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 #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: - enum EndpointConstraints { - NONE, ///< No constraints on endpoint locations - EPS_IN_TARGET, ///< Both EPs are required to be located in the target image - EPS_IN_TARGET_LABELDIFF, ///< Both EPs are required to be located in the target image and the image values at the respective position needs to be distinct - EPS_IN_SEED_AND_TARGET, ///< One EP is required to be located in the seed image and one in the target image - MIN_ONE_EP_IN_TARGET, ///< At least one EP is required to be located in the target image - ONE_EP_IN_TARGET, ///< Exactly one EP is required to be located in the target image - NO_EP_IN_TARGET ///< No EP is allowed to be located in the target image - }; - 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< vnl_vector_fixed > DirectionContainer; 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; } - void SetDicomProperties(mitk::FiberBundle::Pointer fib); itkGetMacro( OutputProbabilityMap, ItkDoubleImgType::Pointer) ///< Output probability map itkGetMacro( FiberPolyData, PolyDataType ) ///< Output fibers - itkGetMacro( UseOutputProbabilityMap, bool) itkGetMacro( MinVoxelSize, float) - itkGetMacro( EndpointConstraint, EndpointConstraints) itkSetMacro( SeedImage, ItkFloatImgType::Pointer) ///< Seeds are only placed inside of this mask. itkSetMacro( MaskImage, ItkFloatImgType::Pointer) ///< Tracking is only performed inside of this mask image. itkSetMacro( ExclusionRegions, ItkFloatImgType::Pointer)///< Fibers passing any of the ROIs in this image are discarded. - 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( EndpointConstraint, EndpointConstraints) ///< Determines what fibers are accepted based on their endpoint location - - 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, ItkFloatImgType::Pointer) ///< Streamlines entering a stopping region will stop immediately itkSetMacro( TargetRegions, ItkFloatImgType::Pointer) ///< Only streamline starting and ending in this mask are retained itkSetMacro( DemoMode, bool ) - itkSetMacro( NumberOfSamples, unsigned int ) ///< Number of neighborhood sampling points - itkSetMacro( LoopCheck, float ) ///< 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 ) - itkSetMacro( InterpolateMasks, bool ) - itkSetMacro( TrialsPerSeed, unsigned int ) ///< When using probabilistic tractography, each seed point is used N times until a valid streamline that is compliant with all thresholds etc. is found - itkSetMacro( TrackingPriorWeight, float) ///< Weight between prior and data [0-1]. One mean tracking only on the prior peaks, zero only on the data. - itkSetMacro( TrackingPriorAsMask, bool) ///< If true, data directions in voxels where prior directions are invalid are set to zero - itkSetMacro( IntroduceDirectionsFromPrior, bool) ///< If false, prior voxels with invalid data voxel are ignored ///< 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; } void Update() override{ this->GenerateData(); } std::string GetStatusText(); void SetTrackingPriorHandler(mitk::TrackingDataHandler *TrackingPriorHandler); protected: void GenerateData() override; StreamlineTrackingFilter(); ~StreamlineTrackingFilter() override {} bool IsValidFiber(FiberType* fib); ///< Check endpoints 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, DirectionContainer* container, float tractLength, bool front, bool& exclude); ///< Start streamline in one direction. vnl_vector_fixed GetNewDirection(const 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(unsigned int NPoints); void BeforeTracking(); void AfterTracking(); PolyDataType m_FiberPolyData; vtkSmartPointer m_Points; vtkSmartPointer m_Cells; BundleType m_Tractogram; BundleType m_GmStubs; ItkFloatImgType::Pointer m_StoppingRegions; ItkFloatImgType::Pointer m_TargetRegions; ItkFloatImgType::Pointer m_SeedImage; ItkFloatImgType::Pointer m_MaskImage; ItkFloatImgType::Pointer m_ExclusionRegions; ItkDoubleImgType::Pointer m_OutputProbabilityMap; float m_MinVoxelSize; - 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_MaxNumTracts; bool m_Verbose; - float m_LoopCheck; 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; - bool m_InterpolateMasks; - unsigned int m_TrialsPerSeed; - EndpointConstraints m_EndpointConstraint; void BuildFibers(bool check); float CheckCurvature(DirectionContainer *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; itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::Pointer m_MaskInterpolator; itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::Pointer m_StopInterpolator; itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::Pointer m_TargetInterpolator; itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::Pointer m_SeedInterpolator; itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::Pointer m_ExclusionInterpolator; bool m_SeedImageSet; bool m_TargetImageSet; + mitk::StreamlineTractographyParameters* m_Parameters; + // Directional priors - bool m_IntroduceDirectionsFromPrior; - bool m_TrackingPriorAsMask; - float m_TrackingPriorWeight; mitk::TrackingDataHandler* m_TrackingPriorHandler; private: }; } //#ifndef ITK_MANUAL_INSTANTIATION //#include "itkMLBSTrackingFilter.cpp" //#endif #endif //__itkMLBSTrackingFilter_h_ diff --git a/Modules/DiffusionImaging/FiberTracking/files.cmake b/Modules/DiffusionImaging/FiberTracking/files.cmake index 19ac812b48..39210be9cf 100644 --- a/Modules/DiffusionImaging/FiberTracking/files.cmake +++ b/Modules/DiffusionImaging/FiberTracking/files.cmake @@ -1,109 +1,111 @@ set(CPP_FILES mitkFiberTrackingModuleActivator.cpp ## IO datastructures IODataStructures/FiberBundle/mitkFiberBundle.cpp IODataStructures/FiberBundle/mitkTrackvis.cpp IODataStructures/PlanarFigureComposite/mitkPlanarFigureComposite.cpp IODataStructures/mitkTractographyForest.cpp IODataStructures/mitkFiberfoxParameters.cpp + IODataStructures/mitkStreamlineTractographyParameters.cpp # Interactions # Tractography Algorithms/GibbsTracking/mitkParticleGrid.cpp Algorithms/GibbsTracking/mitkMetropolisHastingsSampler.cpp Algorithms/GibbsTracking/mitkEnergyComputer.cpp Algorithms/GibbsTracking/mitkGibbsEnergyComputer.cpp Algorithms/GibbsTracking/mitkFiberBuilder.cpp Algorithms/GibbsTracking/mitkSphereInterpolator.cpp Algorithms/itkStreamlineTrackingFilter.cpp Algorithms/TrackingHandlers/mitkTrackingDataHandler.cpp Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.cpp Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.cpp Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.cpp ) set(H_FILES # DataStructures -> FiberBundle IODataStructures/FiberBundle/mitkFiberBundle.h IODataStructures/FiberBundle/mitkTrackvis.h IODataStructures/mitkFiberfoxParameters.h IODataStructures/mitkTractographyForest.h + IODataStructures/mitkStreamlineTractographyParameters.h # Algorithms Algorithms/itkTractDensityImageFilter.h Algorithms/itkTractsToFiberEndingsImageFilter.h Algorithms/itkTractsToRgbaImageFilter.h Algorithms/itkTractsToVectorImageFilter.h Algorithms/itkEvaluateDirectionImagesFilter.h Algorithms/itkEvaluateTractogramDirectionsFilter.h Algorithms/itkFiberCurvatureFilter.h Algorithms/itkFitFibersToImageFilter.h Algorithms/itkTractClusteringFilter.h Algorithms/itkTractDistanceFilter.h Algorithms/itkFiberExtractionFilter.h Algorithms/itkTdiToVolumeFractionFilter.h # Tractography Algorithms/TrackingHandlers/mitkTrackingDataHandler.h Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.h Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.h Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.h Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.h Algorithms/itkGibbsTrackingFilter.h Algorithms/itkStochasticTractographyFilter.h Algorithms/GibbsTracking/mitkParticle.h Algorithms/GibbsTracking/mitkParticleGrid.h Algorithms/GibbsTracking/mitkMetropolisHastingsSampler.h Algorithms/GibbsTracking/mitkSimpSamp.h Algorithms/GibbsTracking/mitkEnergyComputer.h Algorithms/GibbsTracking/mitkGibbsEnergyComputer.h Algorithms/GibbsTracking/mitkSphereInterpolator.h Algorithms/GibbsTracking/mitkFiberBuilder.h Algorithms/itkStreamlineTrackingFilter.h # Clustering Algorithms/ClusteringMetrics/mitkClusteringMetric.h Algorithms/ClusteringMetrics/mitkClusteringMetricEuclideanMean.h Algorithms/ClusteringMetrics/mitkClusteringMetricEuclideanMax.h Algorithms/ClusteringMetrics/mitkClusteringMetricEuclideanStd.h Algorithms/ClusteringMetrics/mitkClusteringMetricAnatomic.h Algorithms/ClusteringMetrics/mitkClusteringMetricScalarMap.h Algorithms/ClusteringMetrics/mitkClusteringMetricInnerAngles.h Algorithms/ClusteringMetrics/mitkClusteringMetricLength.h # Fiberfox Fiberfox/itkFibersFromPlanarFiguresFilter.h Fiberfox/itkTractsToDWIImageFilter.h Fiberfox/itkKspaceImageFilter.h Fiberfox/itkDftImageFilter.h Fiberfox/itkFieldmapGeneratorFilter.h Fiberfox/itkRandomPhantomFilter.h Fiberfox/SignalModels/mitkDiffusionSignalModel.h Fiberfox/SignalModels/mitkTensorModel.h Fiberfox/SignalModels/mitkBallModel.h Fiberfox/SignalModels/mitkDotModel.h Fiberfox/SignalModels/mitkAstroStickModel.h Fiberfox/SignalModels/mitkStickModel.h Fiberfox/SignalModels/mitkRawShModel.h Fiberfox/SignalModels/mitkDiffusionNoiseModel.h Fiberfox/SignalModels/mitkRicianNoiseModel.h Fiberfox/SignalModels/mitkChiSquareNoiseModel.h Fiberfox/Sequences/mitkAcquisitionType.h Fiberfox/Sequences/mitkSingleShotEpi.h Fiberfox/Sequences/mitkConventionalSpinEcho.h Fiberfox/Sequences/mitkFastSpinEcho.h ) set(RESOURCE_FILES # Binary directory resources FiberTrackingLUTBaryCoords.bin FiberTrackingLUTIndices.bin ) diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.cpp index ebb483faa3..98cdc3c58a 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.cpp @@ -1,1001 +1,1026 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ // Blueberry #include #include #include // Qmitk #include "QmitkStreamlineTrackingView.h" #include "QmitkStdMultiWidget.h" // Qt #include +#include // MITK #include #include #include #include #include #include #include #include #include #include #include #include #include #include // VTK #include #include #include #include #include #include #include #include #include #include const std::string QmitkStreamlineTrackingView::VIEW_ID = "org.mitk.views.streamlinetracking"; const std::string id_DataManager = "org.mitk.views.datamanager"; using namespace berry; QmitkStreamlineTrackingWorker::QmitkStreamlineTrackingWorker(QmitkStreamlineTrackingView* view) : m_View(view) { } void QmitkStreamlineTrackingWorker::run() { m_View->m_Tracker->Update(); m_View->m_TrackingThread.quit(); } QmitkStreamlineTrackingView::QmitkStreamlineTrackingView() : m_TrackingWorker(this) , m_Controls(nullptr) , m_FirstTensorProbRun(true) , m_FirstInteractiveRun(true) , m_TrackingHandler(nullptr) , m_ThreadIsRunning(false) , m_DeleteTrackingHandler(false) , m_Visible(false) , m_LastPrior(nullptr) , m_TrackingPriorHandler(nullptr) { m_TrackingWorker.moveToThread(&m_TrackingThread); connect(&m_TrackingThread, SIGNAL(started()), this, SLOT(BeforeThread())); connect(&m_TrackingThread, SIGNAL(started()), &m_TrackingWorker, SLOT(run())); connect(&m_TrackingThread, SIGNAL(finished()), this, SLOT(AfterThread())); m_TrackingTimer = new QTimer(this); } // Destructor QmitkStreamlineTrackingView::~QmitkStreamlineTrackingView() { if (m_Tracker.IsNull()) return; m_Tracker->SetStopTracking(true); m_TrackingThread.wait(); } void QmitkStreamlineTrackingView::CreateQtPartControl( QWidget *parent ) { if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkStreamlineTrackingViewControls; m_Controls->setupUi( parent ); m_Controls->m_FaImageSelectionWidget->SetDataStorage(this->GetDataStorage()); m_Controls->m_SeedImageSelectionWidget->SetDataStorage(this->GetDataStorage()); m_Controls->m_MaskImageSelectionWidget->SetDataStorage(this->GetDataStorage()); m_Controls->m_TargetImageSelectionWidget->SetDataStorage(this->GetDataStorage()); m_Controls->m_PriorImageSelectionWidget->SetDataStorage(this->GetDataStorage()); m_Controls->m_StopImageSelectionWidget->SetDataStorage(this->GetDataStorage()); m_Controls->m_ForestSelectionWidget->SetDataStorage(this->GetDataStorage()); m_Controls->m_ExclusionImageSelectionWidget->SetDataStorage(this->GetDataStorage()); mitk::TNodePredicateDataType::Pointer isPeakImagePredicate = mitk::TNodePredicateDataType::New(); mitk::TNodePredicateDataType::Pointer isImagePredicate = mitk::TNodePredicateDataType::New(); mitk::TNodePredicateDataType::Pointer isTractographyForest = mitk::TNodePredicateDataType::New(); mitk::NodePredicateProperty::Pointer isBinaryPredicate = mitk::NodePredicateProperty::New("binary", mitk::BoolProperty::New(true)); mitk::NodePredicateNot::Pointer isNotBinaryPredicate = mitk::NodePredicateNot::New( isBinaryPredicate ); mitk::NodePredicateAnd::Pointer isNotABinaryImagePredicate = mitk::NodePredicateAnd::New( isImagePredicate, isNotBinaryPredicate ); mitk::NodePredicateDimension::Pointer dimensionPredicate = mitk::NodePredicateDimension::New(3); m_Controls->m_ForestSelectionWidget->SetNodePredicate(isTractographyForest); m_Controls->m_FaImageSelectionWidget->SetNodePredicate( mitk::NodePredicateAnd::New(isNotABinaryImagePredicate, dimensionPredicate) ); m_Controls->m_FaImageSelectionWidget->SetEmptyInfo("--"); m_Controls->m_FaImageSelectionWidget->SetSelectionIsOptional(true); m_Controls->m_SeedImageSelectionWidget->SetNodePredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); m_Controls->m_SeedImageSelectionWidget->SetEmptyInfo("--"); m_Controls->m_SeedImageSelectionWidget->SetSelectionIsOptional(true); m_Controls->m_MaskImageSelectionWidget->SetNodePredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); m_Controls->m_MaskImageSelectionWidget->SetEmptyInfo("--"); m_Controls->m_MaskImageSelectionWidget->SetSelectionIsOptional(true); m_Controls->m_StopImageSelectionWidget->SetNodePredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); m_Controls->m_StopImageSelectionWidget->SetEmptyInfo("--"); m_Controls->m_StopImageSelectionWidget->SetSelectionIsOptional(true); m_Controls->m_TargetImageSelectionWidget->SetNodePredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); m_Controls->m_TargetImageSelectionWidget->SetEmptyInfo("--"); m_Controls->m_TargetImageSelectionWidget->SetSelectionIsOptional(true); m_Controls->m_PriorImageSelectionWidget->SetNodePredicate( isPeakImagePredicate ); m_Controls->m_PriorImageSelectionWidget->SetEmptyInfo("--"); m_Controls->m_PriorImageSelectionWidget->SetSelectionIsOptional(true); m_Controls->m_ExclusionImageSelectionWidget->SetNodePredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); m_Controls->m_ExclusionImageSelectionWidget->SetEmptyInfo("--"); m_Controls->m_ExclusionImageSelectionWidget->SetSelectionIsOptional(true); connect( m_TrackingTimer, SIGNAL(timeout()), this, SLOT(TimerUpdate()) ); + connect( m_Controls->m_SaveParametersButton, SIGNAL(clicked()), this, SLOT(SaveParameters()) ); connect( m_Controls->commandLinkButton_2, SIGNAL(clicked()), this, SLOT(StopTractography()) ); connect( m_Controls->commandLinkButton, SIGNAL(clicked()), this, SLOT(DoFiberTracking()) ); connect( m_Controls->m_InteractiveBox, SIGNAL(stateChanged(int)), this, SLOT(ToggleInteractive()) ); connect( m_Controls->m_ModeBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui()) ); connect( m_Controls->m_FaImageSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(DeleteTrackingHandler()) ); connect( m_Controls->m_ModeBox, SIGNAL(currentIndexChanged(int)), this, SLOT(DeleteTrackingHandler()) ); connect( m_Controls->m_OutputProbMap, SIGNAL(stateChanged(int)), this, SLOT(OutputStyleSwitched()) ); connect( m_Controls->m_SeedImageSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_ModeBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_StopImageSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_TargetImageSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_PriorImageSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_ExclusionImageSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_MaskImageSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FaImageSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_ForestSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(ForestSwitched()) ); connect( m_Controls->m_ForestSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SeedsPerVoxelBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_NumFibersBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_ScalarThresholdBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_OdfCutoffBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_StepSizeBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SamplingDistanceBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_AngularThresholdBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_MinTractLengthBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_fBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_gBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_NumSamplesBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SeedRadiusBox, SIGNAL(editingFinished()), this, SLOT(InteractiveSeedChanged()) ); connect( m_Controls->m_NumSeedsBox, SIGNAL(editingFinished()), this, SLOT(InteractiveSeedChanged()) ); connect( m_Controls->m_OutputProbMap, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SharpenOdfsBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_InterpolationBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_MaskInterpolationBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FlipXBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FlipYBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FlipZBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FrontalSamplesBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_StopVotesBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_LoopCheckBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_TrialsPerSeedBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_EpConstraintsBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); m_Controls->m_SeedsPerVoxelBox->editingFinished(); m_Controls->m_NumFibersBox->editingFinished(); m_Controls->m_ScalarThresholdBox->editingFinished(); m_Controls->m_OdfCutoffBox->editingFinished(); m_Controls->m_StepSizeBox->editingFinished(); m_Controls->m_SamplingDistanceBox->editingFinished(); m_Controls->m_AngularThresholdBox->editingFinished(); m_Controls->m_MinTractLengthBox->editingFinished(); m_Controls->m_fBox->editingFinished(); m_Controls->m_gBox->editingFinished(); m_Controls->m_NumSamplesBox->editingFinished(); m_Controls->m_SeedRadiusBox->editingFinished(); m_Controls->m_NumSeedsBox->editingFinished(); m_Controls->m_LoopCheckBox->editingFinished(); m_Controls->m_TrialsPerSeedBox->editingFinished(); StartStopTrackingGui(false); } + m_ParameterFile = QDir::currentPath()+"/param.stp"; UpdateGui(); } +mitk::StreamlineTractographyParameters QmitkStreamlineTrackingView::GetParametersFromGui() +{ + mitk::StreamlineTractographyParameters params; + params.m_InteractiveRadius = m_Controls->m_SeedRadiusBox->value(); + params.m_MaxNumFibers = m_Controls->m_NumFibersBox->value(); + return params; +} + +void QmitkStreamlineTrackingView::SaveParameters() +{ + QString filename = QFileDialog::getSaveFileName( + 0, + tr("Save Parameters"), + m_ParameterFile, + tr("Streamline Tractography Parameters (*.stp)") ); + + m_ParameterFile = filename; + + auto params = GetParametersFromGui(); + params.SaveParameters(m_ParameterFile.toStdString()); +} + void QmitkStreamlineTrackingView::StopTractography() { if (m_Tracker.IsNull()) return; m_Tracker->SetStopTracking(true); } void QmitkStreamlineTrackingView::TimerUpdate() { if (m_Tracker.IsNull()) return; QString status_text(m_Tracker->GetStatusText().c_str()); m_Controls->m_StatusTextBox->setText(status_text); } void QmitkStreamlineTrackingView::BeforeThread() { m_TrackingTimer->start(1000); } void QmitkStreamlineTrackingView::AfterThread() { m_TrackingTimer->stop(); if (!m_Tracker->GetUseOutputProbabilityMap()) { vtkSmartPointer fiberBundle = m_Tracker->GetFiberPolyData(); if (!m_Controls->m_InteractiveBox->isChecked() && fiberBundle->GetNumberOfLines() == 0) { QMessageBox warnBox; warnBox.setWindowTitle("Warning"); warnBox.setText("No fiberbundle was generated!"); warnBox.setDetailedText("No fibers were generated using the chosen parameters. Typical reasons are:\n\n- Cutoff too high. Some images feature very low FA/GFA/peak size. Try to lower this parameter.\n- Angular threshold too strict. Try to increase this parameter.\n- A small step sizes also means many steps to go wrong. Especially in the case of probabilistic tractography. Try to adjust the angular threshold."); warnBox.setIcon(QMessageBox::Warning); warnBox.exec(); if (m_InteractivePointSetNode.IsNotNull()) m_InteractivePointSetNode->SetProperty("color", mitk::ColorProperty::New(1,1,1)); StartStopTrackingGui(false); if (m_DeleteTrackingHandler) DeleteTrackingHandler(); UpdateGui(); return; } mitk::FiberBundle::Pointer fib = mitk::FiberBundle::New(fiberBundle); fib->SetReferenceGeometry(dynamic_cast(m_ParentNode->GetData())->GetGeometry()); if (m_Controls->m_ResampleFibersBox->isChecked() && fiberBundle->GetNumberOfLines()>0) fib->Compress(m_Controls->m_FiberErrorBox->value()); fib->ColorFibersByOrientation(); m_Tracker->SetDicomProperties(fib); if (m_Controls->m_InteractiveBox->isChecked()) { if (m_InteractiveNode.IsNull()) { m_InteractiveNode = mitk::DataNode::New(); QString name("Interactive"); m_InteractiveNode->SetName(name.toStdString()); GetDataStorage()->Add(m_InteractiveNode); } m_InteractiveNode->SetData(fib); m_InteractiveNode->SetFloatProperty("Fiber2DSliceThickness", m_Tracker->GetMinVoxelSize()/2); if (auto renderWindowPart = this->GetRenderWindowPart()) renderWindowPart->RequestUpdate(); } else { mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(fib); QString name("FiberBundle_"); name += m_ParentNode->GetName().c_str(); name += "_Streamline"; node->SetName(name.toStdString()); node->SetFloatProperty("Fiber2DSliceThickness", m_Tracker->GetMinVoxelSize()/2); GetDataStorage()->Add(node, m_ParentNode); } } else { TrackerType::ItkDoubleImgType::Pointer outImg = m_Tracker->GetOutputProbabilityMap(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); if (m_Controls->m_InteractiveBox->isChecked()) { if (m_InteractiveNode.IsNull()) { m_InteractiveNode = mitk::DataNode::New(); QString name("Interactive"); m_InteractiveNode->SetName(name.toStdString()); GetDataStorage()->Add(m_InteractiveNode); } m_InteractiveNode->SetData(img); mitk::LookupTable::Pointer lut = mitk::LookupTable::New(); lut->SetType(mitk::LookupTable::JET_TRANSPARENT); mitk::LookupTableProperty::Pointer lut_prop = mitk::LookupTableProperty::New(); lut_prop->SetLookupTable(lut); m_InteractiveNode->SetProperty("LookupTable", lut_prop); m_InteractiveNode->SetProperty("opacity", mitk::FloatProperty::New(0.5)); m_InteractiveNode->SetFloatProperty("Fiber2DSliceThickness", m_Tracker->GetMinVoxelSize()/2); if (auto renderWindowPart = this->GetRenderWindowPart()) renderWindowPart->RequestUpdate(); } else { mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(img); QString name("ProbabilityMap_"); name += m_ParentNode->GetName().c_str(); node->SetName(name.toStdString()); mitk::LookupTable::Pointer lut = mitk::LookupTable::New(); lut->SetType(mitk::LookupTable::JET_TRANSPARENT); mitk::LookupTableProperty::Pointer lut_prop = mitk::LookupTableProperty::New(); lut_prop->SetLookupTable(lut); node->SetProperty("LookupTable", lut_prop); node->SetProperty("opacity", mitk::FloatProperty::New(0.5)); GetDataStorage()->Add(node, m_ParentNode); } } if (m_InteractivePointSetNode.IsNotNull()) m_InteractivePointSetNode->SetProperty("color", mitk::ColorProperty::New(1,1,1)); StartStopTrackingGui(false); if (m_DeleteTrackingHandler) DeleteTrackingHandler(); UpdateGui(); } void QmitkStreamlineTrackingView::InteractiveSeedChanged(bool posChanged) { if(!CheckAndStoreLastParams(sender()) && !posChanged) return; if (m_ThreadIsRunning || !m_Visible) return; if (!posChanged && (!m_Controls->m_InteractiveBox->isChecked() || !m_Controls->m_ParamUpdateBox->isChecked()) ) return; std::srand(std::time(0)); m_SeedPoints.clear(); itk::Point world_pos = this->GetRenderWindowPart()->GetSelectedPosition(); m_SeedPoints.push_back(world_pos); float radius = m_Controls->m_SeedRadiusBox->value(); int num = m_Controls->m_NumSeedsBox->value(); mitk::PointSet::Pointer pointset = mitk::PointSet::New(); pointset->InsertPoint(0, world_pos); m_InteractivePointSetNode->SetProperty("pointsize", mitk::FloatProperty::New(radius*2)); m_InteractivePointSetNode->SetProperty("point 2D size", mitk::FloatProperty::New(radius*2)); m_InteractivePointSetNode->SetData(pointset); for (int i=1; i p; p[0] = rand()%1000-500; p[1] = rand()%1000-500; p[2] = rand()%1000-500; p.Normalize(); p *= radius; m_SeedPoints.push_back(world_pos+p); } m_InteractivePointSetNode->SetProperty("color", mitk::ColorProperty::New(1,0,0)); DoFiberTracking(); } bool QmitkStreamlineTrackingView::CheckAndStoreLastParams(QObject* obj) { if (obj!=nullptr) { std::string new_val = ""; if(qobject_cast(obj)!=nullptr) new_val = boost::lexical_cast(qobject_cast(obj)->value()); else if (qobject_cast(obj)!=nullptr) new_val = boost::lexical_cast(qobject_cast(obj)->value()); else return true; if (m_LastTractoParams.find(obj->objectName())==m_LastTractoParams.end()) { m_LastTractoParams[obj->objectName()] = new_val; return false; } else if (m_LastTractoParams.at(obj->objectName()) != new_val) { m_LastTractoParams[obj->objectName()] = new_val; return true; } else if (m_LastTractoParams.at(obj->objectName()) == new_val) return false; } return true; } void QmitkStreamlineTrackingView::OnParameterChanged() { UpdateGui(); if(!CheckAndStoreLastParams(sender())) return; if (m_Controls->m_InteractiveBox->isChecked() && m_Controls->m_ParamUpdateBox->isChecked()) DoFiberTracking(); } void QmitkStreamlineTrackingView::ToggleInteractive() { UpdateGui(); m_Controls->m_SeedsPerVoxelBox->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->m_SeedsPerVoxelLabel->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->m_SeedImageSelectionWidget->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->label_6->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); if ( m_Controls->m_InteractiveBox->isChecked() ) { if (m_FirstInteractiveRun) { QMessageBox::information(nullptr, "Information", "Place and move a spherical seed region anywhere in the image by left-clicking and dragging. If the seed region is colored red, tracking is in progress. If the seed region is colored white, tracking is finished.\nPlacing the seed region for the first time in a newly selected dataset might cause a short delay, since the tracker needs to be initialized."); m_FirstInteractiveRun = false; } QApplication::setOverrideCursor(Qt::PointingHandCursor); QApplication::processEvents(); m_InteractivePointSetNode = mitk::DataNode::New(); m_InteractivePointSetNode->SetProperty("color", mitk::ColorProperty::New(1,1,1)); m_InteractivePointSetNode->SetName("InteractiveSeedRegion"); mitk::PointSetShapeProperty::Pointer shape_prop = mitk::PointSetShapeProperty::New(); shape_prop->SetValue(mitk::PointSetShapeProperty::PointSetShape::CIRCLE); m_InteractivePointSetNode->SetProperty("Pointset.2D.shape", shape_prop); GetDataStorage()->Add(m_InteractivePointSetNode); m_SliceChangeListener.RenderWindowPartActivated(this->GetRenderWindowPart()); connect(&m_SliceChangeListener, SIGNAL(SliceChanged()), this, SLOT(OnSliceChanged())); } else { QApplication::restoreOverrideCursor(); QApplication::processEvents(); m_InteractiveNode = nullptr; m_InteractivePointSetNode = nullptr; m_SliceChangeListener.RenderWindowPartActivated(this->GetRenderWindowPart()); disconnect(&m_SliceChangeListener, SIGNAL(SliceChanged()), this, SLOT(OnSliceChanged())); } } void QmitkStreamlineTrackingView::Activated() { } void QmitkStreamlineTrackingView::Deactivated() { } void QmitkStreamlineTrackingView::Visible() { m_Visible = true; } void QmitkStreamlineTrackingView::Hidden() { m_Visible = false; m_Controls->m_InteractiveBox->setChecked(false); ToggleInteractive(); } void QmitkStreamlineTrackingView::OnSliceChanged() { InteractiveSeedChanged(true); } void QmitkStreamlineTrackingView::SetFocus() { } void QmitkStreamlineTrackingView::DeleteTrackingHandler() { if (!m_ThreadIsRunning && m_TrackingHandler != nullptr) { if (m_TrackingPriorHandler != nullptr) { delete m_TrackingPriorHandler; m_TrackingPriorHandler = nullptr; } delete m_TrackingHandler; m_TrackingHandler = nullptr; m_DeleteTrackingHandler = false; m_LastPrior = nullptr; } else if (m_ThreadIsRunning) { m_DeleteTrackingHandler = true; } } void QmitkStreamlineTrackingView::ForestSwitched() { DeleteTrackingHandler(); } void QmitkStreamlineTrackingView::OutputStyleSwitched() { if (m_InteractiveNode.IsNotNull()) GetDataStorage()->Remove(m_InteractiveNode); m_InteractiveNode = nullptr; } void QmitkStreamlineTrackingView::OnSelectionChanged( berry::IWorkbenchPart::Pointer , const QList& nodes ) { std::vector< mitk::DataNode::Pointer > last_nodes = m_InputImageNodes; m_InputImageNodes.clear(); m_AdditionalInputImages.clear(); bool retrack = false; for( auto node : nodes ) { if( node.IsNotNull() && dynamic_cast(node->GetData()) ) { if( dynamic_cast(node->GetData()) || dynamic_cast(node->GetData()) || dynamic_cast(node->GetData()) || dynamic_cast(node->GetData()) || mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(node->GetData()))) { m_InputImageNodes.push_back(node); retrack = true; } else { mitk::Image* img = dynamic_cast(node->GetData()); if (img!=nullptr && img->GetDimension()==3) m_AdditionalInputImages.push_back(dynamic_cast(node->GetData())); } } } // sometimes the OnSelectionChanged event is sent twice and actually no selection has changed for the first event. We need to catch that. if (last_nodes.size() == m_InputImageNodes.size()) { bool same_nodes = true; for (unsigned int i=0; im_TensorImageLabel->setText("select in data-manager"); m_Controls->m_fBox->setEnabled(false); m_Controls->m_fLabel->setEnabled(false); m_Controls->m_gBox->setEnabled(false); m_Controls->m_gLabel->setEnabled(false); m_Controls->m_FaImageSelectionWidget->setEnabled(true); m_Controls->mFaImageLabel->setEnabled(true); m_Controls->m_OdfCutoffBox->setEnabled(false); m_Controls->m_OdfCutoffLabel->setEnabled(false); m_Controls->m_SharpenOdfsBox->setEnabled(false); m_Controls->m_ForestSelectionWidget->setVisible(false); m_Controls->m_ForestLabel->setVisible(false); m_Controls->commandLinkButton->setEnabled(false); m_Controls->m_TrialsPerSeedBox->setEnabled(false); m_Controls->m_TrialsPerSeedLabel->setEnabled(false); m_Controls->m_TargetImageSelectionWidget->setEnabled(false); m_Controls->m_TargetImageLabel->setEnabled(false); if (m_Controls->m_InteractiveBox->isChecked()) { m_Controls->m_InteractiveSeedingFrame->setVisible(true); m_Controls->m_StaticSeedingFrame->setVisible(false); m_Controls->commandLinkButton_2->setVisible(false); m_Controls->commandLinkButton->setVisible(false); } else { m_Controls->m_InteractiveSeedingFrame->setVisible(false); m_Controls->m_StaticSeedingFrame->setVisible(true); m_Controls->commandLinkButton_2->setVisible(m_ThreadIsRunning); m_Controls->commandLinkButton->setVisible(!m_ThreadIsRunning); } if (m_Controls->m_EpConstraintsBox->currentIndex()>0) { m_Controls->m_TargetImageSelectionWidget->setEnabled(true); m_Controls->m_TargetImageLabel->setEnabled(true); } // trials per seed are only important for probabilistic tractography if (m_Controls->m_ModeBox->currentIndex()==1) { m_Controls->m_TrialsPerSeedBox->setEnabled(true); m_Controls->m_TrialsPerSeedLabel->setEnabled(true); } if(!m_InputImageNodes.empty()) { if (m_InputImageNodes.size()>1) m_Controls->m_TensorImageLabel->setText( ( std::to_string(m_InputImageNodes.size()) + " images selected").c_str() ); else m_Controls->m_TensorImageLabel->setText(m_InputImageNodes.at(0)->GetName().c_str()); m_Controls->commandLinkButton->setEnabled(!m_Controls->m_InteractiveBox->isChecked() && !m_ThreadIsRunning); m_Controls->m_ScalarThresholdBox->setEnabled(true); m_Controls->m_FaThresholdLabel->setEnabled(true); if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { m_Controls->m_fBox->setEnabled(true); m_Controls->m_fLabel->setEnabled(true); m_Controls->m_gBox->setEnabled(true); m_Controls->m_gLabel->setEnabled(true); } else if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) || dynamic_cast(m_InputImageNodes.at(0)->GetData())) { m_Controls->m_OdfCutoffBox->setEnabled(true); m_Controls->m_OdfCutoffLabel->setEnabled(true); m_Controls->m_SharpenOdfsBox->setEnabled(true); } else if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(m_InputImageNodes.at(0)->GetData())) ) { m_Controls->m_ForestSelectionWidget->setVisible(true); m_Controls->m_ForestLabel->setVisible(true); m_Controls->m_ScalarThresholdBox->setEnabled(false); m_Controls->m_FaThresholdLabel->setEnabled(false); } } } void QmitkStreamlineTrackingView::StartStopTrackingGui(bool start) { m_ThreadIsRunning = start; if (!m_Controls->m_InteractiveBox->isChecked()) { m_Controls->commandLinkButton_2->setVisible(start); m_Controls->commandLinkButton->setVisible(!start); m_Controls->m_InteractiveBox->setEnabled(!start); m_Controls->m_StatusTextBox->setVisible(start); } } void QmitkStreamlineTrackingView::DoFiberTracking() { if (m_InputImageNodes.empty()) { QMessageBox::information(nullptr, "Information", "Please select an input image in the datamaneger (tensor, ODF, peak or dMRI image)!"); return; } if (m_ThreadIsRunning || !m_Visible) return; if (m_Controls->m_InteractiveBox->isChecked() && m_SeedPoints.empty()) return; StartStopTrackingGui(true); m_Tracker = TrackerType::New(); if( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { if (m_Controls->m_ModeBox->currentIndex()==1) { if (m_InputImageNodes.size()>1) { QMessageBox::information(nullptr, "Information", "Probabilistic tensor tractography is only implemented for single-tensor mode!"); StartStopTrackingGui(false); return; } if (m_TrackingHandler==nullptr) { m_TrackingHandler = new mitk::TrackingHandlerOdf(); typedef itk::TensorImageToOdfImageFilter< float, float > FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetInput( mitk::convert::GetItkTensorFromTensorImage(dynamic_cast(m_InputImageNodes.at(0)->GetData())) ); filter->Update(); dynamic_cast(m_TrackingHandler)->SetOdfImage(filter->GetOutput()); if (m_Controls->m_FaImageSelectionWidget->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer itkImg = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageSelectionWidget->GetSelectedNode()->GetData()), itkImg); dynamic_cast(m_TrackingHandler)->SetGfaImage(itkImg); } } dynamic_cast(m_TrackingHandler)->SetGfaThreshold(static_cast(m_Controls->m_ScalarThresholdBox->value())); dynamic_cast(m_TrackingHandler)->SetOdfThreshold(0); dynamic_cast(m_TrackingHandler)->SetSharpenOdfs(true); dynamic_cast(m_TrackingHandler)->SetIsOdfFromTensor(true); } else { if (m_TrackingHandler==nullptr) { m_TrackingHandler = new mitk::TrackingHandlerTensor(); for (unsigned int i=0; i(m_TrackingHandler)->AddTensorImage(mitk::convert::GetItkTensorFromTensorImage(dynamic_cast(m_InputImageNodes.at(i)->GetData())).GetPointer()); if (m_Controls->m_FaImageSelectionWidget->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer itkImg = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageSelectionWidget->GetSelectedNode()->GetData()), itkImg); dynamic_cast(m_TrackingHandler)->SetFaImage(itkImg); } } dynamic_cast(m_TrackingHandler)->SetFaThreshold(static_cast(m_Controls->m_ScalarThresholdBox->value())); dynamic_cast(m_TrackingHandler)->SetF(static_cast(m_Controls->m_fBox->value())); dynamic_cast(m_TrackingHandler)->SetG(static_cast(m_Controls->m_gBox->value())); } } else if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) || dynamic_cast(m_InputImageNodes.at(0)->GetData())) { if (m_TrackingHandler==nullptr) { m_TrackingHandler = new mitk::TrackingHandlerOdf(); if (dynamic_cast(m_InputImageNodes.at(0)->GetData())) dynamic_cast(m_TrackingHandler)->SetOdfImage(mitk::convert::GetItkOdfFromShImage(dynamic_cast(m_InputImageNodes.at(0)->GetData()))); else dynamic_cast(m_TrackingHandler)->SetOdfImage(mitk::convert::GetItkOdfFromOdfImage(dynamic_cast(m_InputImageNodes.at(0)->GetData()))); if (m_Controls->m_FaImageSelectionWidget->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer itkImg = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageSelectionWidget->GetSelectedNode()->GetData()), itkImg); dynamic_cast(m_TrackingHandler)->SetGfaImage(itkImg); } } dynamic_cast(m_TrackingHandler)->SetGfaThreshold(static_cast(m_Controls->m_ScalarThresholdBox->value())); dynamic_cast(m_TrackingHandler)->SetOdfThreshold(static_cast(m_Controls->m_OdfCutoffBox->value())); dynamic_cast(m_TrackingHandler)->SetSharpenOdfs(m_Controls->m_SharpenOdfsBox->isChecked()); } else if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(m_InputImageNodes.at(0)->GetData())) ) { if ( m_Controls->m_ForestSelectionWidget->GetSelectedNode().IsNull() ) { QMessageBox::information(nullptr, "Information", "Not random forest for machine learning based tractography (raw dMRI tractography) selected. Did you accidentally select the raw diffusion-weighted image in the datamanager?"); StartStopTrackingGui(false); return; } if (m_TrackingHandler==nullptr) { mitk::TractographyForest::Pointer forest = dynamic_cast(m_Controls->m_ForestSelectionWidget->GetSelectedNode()->GetData()); mitk::Image::Pointer dwi = dynamic_cast(m_InputImageNodes.at(0)->GetData()); std::vector< std::vector< ItkFloatImageType::Pointer > > additionalFeatureImages; additionalFeatureImages.push_back(std::vector< ItkFloatImageType::Pointer >()); for (auto img : m_AdditionalInputImages) { ItkFloatImageType::Pointer itkimg = ItkFloatImageType::New(); mitk::CastToItkImage(img, itkimg); additionalFeatureImages.at(0).push_back(itkimg); } bool forest_valid = false; if (forest->GetNumFeatures()>=100) { unsigned int num_previous_directions = static_cast((forest->GetNumFeatures() - (100 + additionalFeatureImages.at(0).size()))/3); m_TrackingHandler = new mitk::TrackingHandlerRandomForest<6, 100>(); dynamic_cast*>(m_TrackingHandler)->AddDwi(dwi); dynamic_cast*>(m_TrackingHandler)->SetAdditionalFeatureImages(additionalFeatureImages); dynamic_cast*>(m_TrackingHandler)->SetForest(forest); dynamic_cast*>(m_TrackingHandler)->SetNumPreviousDirections(num_previous_directions); forest_valid = dynamic_cast*>(m_TrackingHandler)->IsForestValid(); } else { unsigned int num_previous_directions = static_cast((forest->GetNumFeatures() - (28 + additionalFeatureImages.at(0).size()))/3); m_TrackingHandler = new mitk::TrackingHandlerRandomForest<6, 28>(); dynamic_cast*>(m_TrackingHandler)->AddDwi(dwi); dynamic_cast*>(m_TrackingHandler)->SetAdditionalFeatureImages(additionalFeatureImages); dynamic_cast*>(m_TrackingHandler)->SetForest(forest); dynamic_cast*>(m_TrackingHandler)->SetNumPreviousDirections(num_previous_directions); forest_valid = dynamic_cast*>(m_TrackingHandler)->IsForestValid(); } if (!forest_valid) { QMessageBox::information(nullptr, "Information", "Random forest is invalid. The forest signatue does not match the parameters of TrackingHandlerRandomForest."); StartStopTrackingGui(false); return; } } } else { if (m_Controls->m_ModeBox->currentIndex()==1) { QMessageBox::information(nullptr, "Information", "Probabilstic tractography is not implemented for peak images."); StartStopTrackingGui(false); return; } if (m_TrackingHandler==nullptr) { m_TrackingHandler = new mitk::TrackingHandlerPeaks(); dynamic_cast(m_TrackingHandler)->SetPeakImage(mitk::convert::GetItkPeakFromPeakImage(dynamic_cast(m_InputImageNodes.at(0)->GetData()))); } dynamic_cast(m_TrackingHandler)->SetPeakThreshold(static_cast(m_Controls->m_ScalarThresholdBox->value())); } m_TrackingHandler->SetFlipX(m_Controls->m_FlipXBox->isChecked()); m_TrackingHandler->SetFlipY(m_Controls->m_FlipYBox->isChecked()); m_TrackingHandler->SetFlipZ(m_Controls->m_FlipZBox->isChecked()); m_TrackingHandler->SetInterpolate(m_Controls->m_InterpolationBox->isChecked()); switch (m_Controls->m_ModeBox->currentIndex()) { case 0: m_TrackingHandler->SetMode(mitk::TrackingDataHandler::MODE::DETERMINISTIC); break; case 1: m_TrackingHandler->SetMode(mitk::TrackingDataHandler::MODE::PROBABILISTIC); break; default: m_TrackingHandler->SetMode(mitk::TrackingDataHandler::MODE::DETERMINISTIC); } if (m_Controls->m_InteractiveBox->isChecked()) { m_Tracker->SetSeedPoints(m_SeedPoints); } else if (m_Controls->m_SeedImageSelectionWidget->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_SeedImageSelectionWidget->GetSelectedNode()->GetData()), mask); m_Tracker->SetSeedImage(mask); } if (m_Controls->m_MaskImageSelectionWidget->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_MaskImageSelectionWidget->GetSelectedNode()->GetData()), mask); m_Tracker->SetMaskImage(mask); } if (m_Controls->m_StopImageSelectionWidget->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_StopImageSelectionWidget->GetSelectedNode()->GetData()), mask); m_Tracker->SetStoppingRegions(mask); } if (m_Controls->m_TargetImageSelectionWidget->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_TargetImageSelectionWidget->GetSelectedNode()->GetData()), mask); m_Tracker->SetTargetRegions(mask); } if (m_Controls->m_PriorImageSelectionWidget->GetSelectedNode().IsNotNull()) { if (m_LastPrior!=m_Controls->m_PriorImageSelectionWidget->GetSelectedNode() || m_TrackingPriorHandler==nullptr) { typedef mitk::ImageToItk< mitk::TrackingHandlerPeaks::PeakImgType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(dynamic_cast(m_Controls->m_PriorImageSelectionWidget->GetSelectedNode()->GetData())); caster->SetCopyMemFlag(true); caster->Update(); mitk::TrackingHandlerPeaks::PeakImgType::Pointer itkImg = caster->GetOutput(); m_TrackingPriorHandler = new mitk::TrackingHandlerPeaks(); dynamic_cast(m_TrackingPriorHandler)->SetPeakImage(itkImg); dynamic_cast(m_TrackingPriorHandler)->SetPeakThreshold(0.0); m_LastPrior = m_Controls->m_PriorImageSelectionWidget->GetSelectedNode(); } m_TrackingPriorHandler->SetInterpolate(m_Controls->m_InterpolationBox->isChecked()); m_TrackingPriorHandler->SetMode(mitk::TrackingDataHandler::MODE::DETERMINISTIC); m_Tracker->SetTrackingPriorHandler(m_TrackingPriorHandler); m_Tracker->SetTrackingPriorWeight(m_Controls->m_PriorWeightBox->value()); m_Tracker->SetTrackingPriorAsMask(m_Controls->m_PriorAsMaskBox->isChecked()); m_Tracker->SetIntroduceDirectionsFromPrior(m_Controls->m_NewDirectionsFromPriorBox->isChecked()); } else if (m_Controls->m_PriorImageSelectionWidget->GetSelectedNode().IsNull()) m_Tracker->SetTrackingPriorHandler(nullptr); if (m_Controls->m_ExclusionImageSelectionWidget->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_ExclusionImageSelectionWidget->GetSelectedNode()->GetData()), mask); m_Tracker->SetExclusionRegions(mask); } // Endpoint constraints switch (m_Controls->m_EpConstraintsBox->currentIndex()) { case 0: m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::NONE); m_Tracker->SetTargetRegions(nullptr); break; case 1: m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_TARGET); break; case 2: m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_TARGET_LABELDIFF); break; case 3: m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_SEED_AND_TARGET); break; case 4: m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::MIN_ONE_EP_IN_TARGET); break; case 5: m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::ONE_EP_IN_TARGET); break; case 6: m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::NO_EP_IN_TARGET); break; } if (m_Tracker->GetEndpointConstraint()!=itk::StreamlineTrackingFilter::EndpointConstraints::NONE && m_Controls->m_TargetImageSelectionWidget->GetSelectedNode().IsNull()) { QMessageBox::information(nullptr, "Error", "Endpoint constraints are used but no target image is set!"); StartStopTrackingGui(false); return; } else if (m_Tracker->GetEndpointConstraint()==itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_SEED_AND_TARGET && (m_Controls->m_SeedImageSelectionWidget->GetSelectedNode().IsNull()|| m_Controls->m_TargetImageSelectionWidget->GetSelectedNode().IsNull()) ) { QMessageBox::information(nullptr, "Error", "Endpoint constraint EPS_IN_SEED_AND_TARGET is used but no target or no seed image is set!"); StartStopTrackingGui(false); return; } m_Tracker->SetInterpolateMasks(m_Controls->m_MaskInterpolationBox->isChecked()); m_Tracker->SetVerbose(!m_Controls->m_InteractiveBox->isChecked()); m_Tracker->SetSeedsPerVoxel(m_Controls->m_SeedsPerVoxelBox->value()); m_Tracker->SetStepSize(m_Controls->m_StepSizeBox->value()); m_Tracker->SetSamplingDistance(m_Controls->m_SamplingDistanceBox->value()); m_Tracker->SetUseStopVotes(m_Controls->m_StopVotesBox->isChecked()); m_Tracker->SetOnlyForwardSamples(m_Controls->m_FrontalSamplesBox->isChecked()); m_Tracker->SetTrialsPerSeed(m_Controls->m_TrialsPerSeedBox->value()); m_Tracker->SetMaxNumTracts(m_Controls->m_NumFibersBox->value()); m_Tracker->SetNumberOfSamples(m_Controls->m_NumSamplesBox->value()); m_Tracker->SetTrackingHandler(m_TrackingHandler); m_Tracker->SetLoopCheck(m_Controls->m_LoopCheckBox->value()); m_Tracker->SetAngularThreshold(m_Controls->m_AngularThresholdBox->value()); m_Tracker->SetMinTractLength(m_Controls->m_MinTractLengthBox->value()); m_Tracker->SetUseOutputProbabilityMap(m_Controls->m_OutputProbMap->isChecked()); m_Tracker->SetRandom(!m_Controls->m_FixSeedBox->isChecked()); m_ParentNode = m_InputImageNodes.at(0); m_TrackingThread.start(QThread::LowestPriority); } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.h b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.h index f2d3bcb433..96a45af112 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.h @@ -1,154 +1,158 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef QmitkStreamlineTrackingView_h #define QmitkStreamlineTrackingView_h #include #include "ui_QmitkStreamlineTrackingViewControls.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include +#include class QmitkStreamlineTrackingView; class QmitkStreamlineTrackingWorker : public QObject { Q_OBJECT public: QmitkStreamlineTrackingWorker(QmitkStreamlineTrackingView* view); public slots: void run(); private: QmitkStreamlineTrackingView* m_View; }; /*! \brief View for tensor based deterministic streamline fiber tracking. */ class QmitkStreamlineTrackingView : public QmitkAbstractView, public mitk::ILifecycleAwarePart { // this is needed for all Qt objects that should have a Qt meta-object // (everything that derives from QObject and wants to have signal/slots) Q_OBJECT public: static const std::string VIEW_ID; typedef itk::Image< unsigned int, 3 > ItkUintImgType; typedef itk::Image< unsigned char, 3 > ItkUCharImageType; typedef itk::Image< float, 3 > ItkFloatImageType; typedef itk::StreamlineTrackingFilter TrackerType; QmitkStreamlineTrackingView(); virtual ~QmitkStreamlineTrackingView(); virtual void CreateQtPartControl(QWidget *parent) override; /// /// Sets the focus to an internal widget. /// virtual void SetFocus() override; TrackerType::Pointer m_Tracker; QmitkStreamlineTrackingWorker m_TrackingWorker; QThread m_TrackingThread; virtual void Activated() override; virtual void Deactivated() override; virtual void Visible() override; virtual void Hidden() override; protected slots: void DoFiberTracking(); ///< start fiber tracking void UpdateGui(); void ToggleInteractive(); void DeleteTrackingHandler(); void OnParameterChanged(); void InteractiveSeedChanged(bool posChanged=false); void ForestSwitched(); void OutputStyleSwitched(); void AfterThread(); ///< update gui etc. after tracking has finished void BeforeThread(); ///< start timer etc. void TimerUpdate(); void StopTractography(); void OnSliceChanged(); + void SaveParameters(); protected: /// \brief called by QmitkAbstractView when DataManager's selection has changed virtual void OnSelectionChanged(berry::IWorkbenchPart::Pointer part, const QList& nodes) override; Ui::QmitkStreamlineTrackingViewControls* m_Controls; protected slots: private: bool CheckAndStoreLastParams(QObject* obj); void StartStopTrackingGui(bool start); + mitk::StreamlineTractographyParameters GetParametersFromGui(); std::vector< itk::Point > m_SeedPoints; mitk::DataNode::Pointer m_ParentNode; mitk::DataNode::Pointer m_InteractiveNode; mitk::DataNode::Pointer m_InteractivePointSetNode; std::vector< mitk::DataNode::Pointer > m_InputImageNodes; ///< input image nodes std::vector< mitk::Image::ConstPointer > m_AdditionalInputImages; bool m_FirstTensorProbRun; bool m_FirstInteractiveRun; mitk::TrackingDataHandler* m_TrackingHandler; bool m_ThreadIsRunning; QTimer* m_TrackingTimer; bool m_DeleteTrackingHandler; QmitkSliceNavigationListener m_SliceChangeListener; bool m_Visible; mitk::DataNode::Pointer m_LastPrior; mitk::TrackingDataHandler* m_TrackingPriorHandler; std::map< QString, std::string > m_LastTractoParams; + QString m_ParameterFile; }; #endif // _QMITKFIBERTRACKINGVIEW_H_INCLUDED diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingViewControls.ui index e7ec1f72fc..b04abdebc7 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingViewControls.ui @@ -1,1511 +1,1524 @@ QmitkStreamlineTrackingViewControls 0 0 453 859 0 0 QmitkTemplate QCommandLinkButton:disabled { border: none; } QGroupBox { background-color: transparent; } 3 3 0 40 QFrame::NoFrame QFrame::Raised 0 15 0 0 6 15 true 0 0 true QFrame::NoFrame QFrame::Raised 0 0 0 0 Input Image. ODF, tensor and peak images are currently supported. Input Image: Input Image. ODF, tensor, peak, and, in case of ML tractography, raw diffusion-weighted images are currently supported. <html><head/><body><p><span style=" color:#ff0000;">select image in data-manager</span></p></body></html> true Tractography Forest: true Stop tractography and return all fibers reconstructed until now. Stop Tractography false Start Tractography + + + + true + + + Stop tractography and return all fibers reconstructed until now. + + + Save Parameters + + + 0 0 0 0 75 true 0 0 0 421 254 Seeding Specify how, where and how many tractography seed points are placed. QFrame::NoFrame QFrame::Raised 0 0 0 0 QFrame::NoFrame QFrame::Raised 0 0 0 0 Number of seed points equally distributed around selected position. 1 9999999 50 Radius: Seedpoints are equally distributed within a sphere centered at the selected position with the specified radius (in mm). 2 50.000000000000000 0.100000000000000 2.000000000000000 Num. Seeds: true When checked, parameter changes cause instant retracking while in interactive mode. Update on Parameter Change true QFrame::NoFrame QFrame::Raised 0 0 0 0 Try each seed N times until a valid streamline is obtained (only for probabilistic tractography). Minimum fiber length (in mm) 1 999 10 Trials Per Seed: Max. Num. Fibers: Tractography is stopped after the desired number of fibers is reached, even before all seed points are processed (-1 means no limit). -1 999999999 -1 QFrame::NoFrame QFrame::Raised 0 0 0 0 Seeds per Voxel: Seed Image: Number of seed points placed in each voxel. 1 9999999 true Dynamically pick a seed location by click into image. Enable Interactive Tractography Qt::Vertical 20 40 0 0 435 232 ROI Constraints Specify various ROI and mask images to constrain the tractography process. Mask Image: Select which fibers should be accepted or rejected based on the location of their endpoints. QComboBox::AdjustToMinimumContentsLength No Constraints on EP locations Both EPs in Target Image Both EPs in Target Image But Different Label One EP in Seed Image and One EP in Target Image At Least One EP in Target Image Exactly One EP in Target Image No EP in Target Image Endpoint Constraints: Stop ROI Image: Exclusion ROI Image: Target ROI Image: Qt::Vertical 20 40 0 - -132 + 0 421 366 Tractography Parameters Specify the behavior of the tractography at each streamline integration step (step size, deterministic/probabilistic, ...). Angular threshold between two steps (in degree). Default: 90° * step_size -1 90 1 -1 FA/GFA Image: ODF Cutoff: Always produce the same random numbers. Qt::Vertical 20 40 f=1 + g=0 means FACT (depending on the chosen interpolation). f=0 and g=1 means TEND (disable interpolation for this mode!). 2 1.000000000000000 0.100000000000000 0.000000000000000 Minimum tract length in mm. Shorter fibers are discarded. Minimum fiber length (in mm) 1 999.000000000000000 1.000000000000000 20.000000000000000 Sharpen ODFs: Threshold on peak magnitude, FA, GFA, ... 5 1.000000000000000 0.100000000000000 0.100000000000000 g: Loop Check: If you are using dODF images as input, it is advisable to sharpen the ODFs (min-max normalize and raise to the power of 4). This is not necessary for CSD fODFs, since they are naturally much sharper. f=1 + g=0 means FACT (depending on the chosen interpolation). f=0 and g=1 means TEND (disable interpolation for this mode!). 2 1.000000000000000 0.100000000000000 1.000000000000000 Maximum allowed angular SDTEV over 4 voxel lengths. Default: no loop check. -1 180 -1 Min. Tract Length: f parameter of tensor tractography. f=1 + g=0 means FACT (depending on the chosen interpolation). f=0 and g=1 means TEND (disable interpolation for this mode!). f: Angular Threshold: Step size (in voxels) 2 0.010000000000000 10.000000000000000 0.100000000000000 0.500000000000000 Toggle between deterministic and probabilistic tractography. Some modes might not be available for all types of tractography. Deterministic Probabilistic Mode: Step Size: Cutoff: Additional threshold on the ODF magnitude. This is useful in case of CSD fODF tractography. For fODFs a good default value is 0.1, for normalized dODFs, e.g. Q-ball ODFs, this threshold should be very low (0.00025) or 0. 5 1.000000000000000 0.100000000000000 0.000250000000000 Fix Random Seed: 0 0 435 232 Tractography Prior Restrict to Prior: Weight: Peak Image: Weighting factor between prior and data. 1.000000000000000 0.100000000000000 0.500000000000000 Restrict tractography to regions where the prior is valid. true Qt::Vertical 20 40 New Directions from Prior: If unchecked, the prior cannot create directions where there are none in the data. true 0 0 435 232 Neighborhood Sampling Specify if and how information about the current streamline neighborhood should be used. Only neighborhood samples in front of the current streamline position are considered. Use Only Frontal Samples false If checked, the majority of sampling points has to place a stop-vote for the streamline to terminate. If not checked, all sampling positions have to vote for a streamline termination. Use Stop-Votes false QFrame::NoFrame QFrame::Raised 0 0 0 0 Num. Samples: Number of neighborhood samples that are used to determine the next fiber progression direction. 50 Sampling Distance: Sampling distance (in voxels) 2 10.000000000000000 0.100000000000000 0.250000000000000 Qt::Vertical 20 40 0 0 435 232 Data Handling Specify interpolation and direction flips. QFrame::NoFrame QFrame::Raised 0 0 0 0 Trilinearly interpolate the input image used for tractography. Interpolate Tractography Data true Trilinearly interpolate the ROI images used to constrain the tractography. Interpolate ROI Images true QFrame::NoFrame QFrame::Raised 0 0 0 0 QFrame::NoFrame QFrame::Raised 0 0 0 0 Internally flips progression directions. This might be necessary depending on the input data. x Internally flips progression directions. This might be necessary depending on the input data. y Internally flips progression directions. This might be necessary depending on the input data. z Flip directions: Qt::Vertical 20 40 0 0 435 232 Output and Postprocessing Specify the tractography output (streamlines or probability maps) and postprocessing steps. QFrame::NoFrame QFrame::Raised 0 0 0 0 Compress fibers using the specified error constraint. Compress Fibers true Qt::StrongFocus Lossy fiber compression. Recommended for large tractograms. Maximum error in mm. 3 10.000000000000000 0.010000000000000 0.100000000000000 Output map with voxel-wise visitation counts instead of streamlines. Output Probability Map false Qt::Vertical 20 40 QmitkSingleNodeSelectionWidget QWidget
QmitkSingleNodeSelectionWidget.h
1