diff --git a/CMake/BuildConfigurations/DiffusionAll.cmake b/CMake/BuildConfigurations/DiffusionAll.cmake index 26f661e2c9..79807ee352 100644 --- a/CMake/BuildConfigurations/DiffusionAll.cmake +++ b/CMake/BuildConfigurations/DiffusionAll.cmake @@ -1,34 +1,33 @@ message(STATUS "Configuring MITK Diffusion with all Plugins") # Enable non-optional external dependencies set(MITK_USE_Vigra ON CACHE BOOL "MITK Use Vigra Library" FORCE) set(MITK_USE_HDF5 ON CACHE BOOL "MITK Use HDF5 Library" FORCE) -set(MITK_USE_DCMQI ON CACHE BOOL "" FORCE) set(MITK_USE_MatchPoint ON CACHE BOOL "" FORCE) # Disable all apps but MITK Diffusion set(MITK_BUILD_ALL_APPS OFF CACHE BOOL "Build all MITK applications" FORCE) set(MITK_BUILD_APP_CoreApp OFF CACHE BOOL "Build the MITK CoreApp" FORCE) set(MITK_BUILD_APP_Diffusion ON CACHE BOOL "Build MITK Diffusion" FORCE) # Activate Diffusion Mini Apps set(BUILD_DiffusionCoreCmdApps ON CACHE BOOL "Build commandline tools for diffusion" FORCE) set(BUILD_DiffusionFiberProcessingCmdApps ON CACHE BOOL "Build commandline tools for diffusion fiber processing" FORCE) set(BUILD_DiffusionFiberfoxCmdApps ON CACHE BOOL "Build commandline tools for diffusion data simulation (Fiberfox)" FORCE) set(BUILD_DiffusionMiscCmdApps ON CACHE BOOL "Build miscellaneous commandline tools for diffusion" FORCE) set(BUILD_DiffusionQuantificationCmdApps ON CACHE BOOL "Build commandline tools for diffusion quantification (IVIM, ADC, ...)" FORCE) set(BUILD_DiffusionTractographyCmdApps ON CACHE BOOL "Build commandline tools for diffusion fiber tractography" FORCE) set(BUILD_DiffusionTractographyEvaluationCmdApps ON CACHE BOOL "Build commandline tools for diffusion fiber tractography evaluation" FORCE) set(BUILD_DiffusionConnectomicsCmdApps ON CACHE BOOL "Build commandline tools for diffusion connectomics" FORCE) # Build neither all plugins nor examples set(MITK_BUILD_ALL_PLUGINS OFF CACHE BOOL "Build all MITK plugins" FORCE) set(MITK_BUILD_EXAMPLES OFF CACHE BOOL "Build the MITK examples" FORCE) set(BUILD_TESTING ON CACHE BOOL "Build the MITK tests" FORCE) # Activate in-application help generation set(MITK_DOXYGEN_GENERATE_QCH_FILES ON CACHE BOOL "Use doxygen to generate Qt compressed help files for MITK docs" FORCE) set(BLUEBERRY_USE_QT_HELP ON CACHE BOOL "Enable support for integrating bundle documentation into Qt Help" FORCE) # Disable console window set(MITK_SHOW_CONSOLE_WINDOW ON CACHE BOOL "Use this to enable or disable the console window when starting MITK GUI Applications" FORCE) diff --git a/CMake/BuildConfigurations/DiffusionRelease.cmake b/CMake/BuildConfigurations/DiffusionRelease.cmake index d79a8af152..873a68aedf 100644 --- a/CMake/BuildConfigurations/DiffusionRelease.cmake +++ b/CMake/BuildConfigurations/DiffusionRelease.cmake @@ -1,37 +1,36 @@ message(STATUS "Configuring MITK Diffusion Release Build") # Enable non-optional external dependencies set(MITK_USE_Vigra ON CACHE BOOL "MITK Use Vigra Library" FORCE) set(MITK_USE_HDF5 ON CACHE BOOL "MITK Use HDF5 Library" FORCE) -set(MITK_USE_DCMQI ON CACHE BOOL "" FORCE) set(MITK_USE_MatchPoint ON CACHE BOOL "" FORCE) # Disable all apps but MITK Diffusion set(MITK_BUILD_ALL_APPS OFF CACHE BOOL "Build all MITK applications" FORCE) set(MITK_BUILD_APP_CoreApp OFF CACHE BOOL "Build the MITK CoreApp" FORCE) set(MITK_BUILD_APP_Workbench OFF CACHE BOOL "Build the MITK Workbench" FORCE) set(MITK_BUILD_APP_Diffusion ON CACHE BOOL "Build MITK Diffusion" FORCE) # Activate Diffusion Mini Apps set(BUILD_DiffusionCoreCmdApps OFF CACHE BOOL "Build commandline tools for diffusion" FORCE) set(BUILD_DiffusionFiberProcessingCmdApps ON CACHE BOOL "Build commandline tools for diffusion fiber processing" FORCE) set(BUILD_DiffusionFiberfoxCmdApps ON CACHE BOOL "Build commandline tools for diffusion data simulation (Fiberfox)" FORCE) set(BUILD_DiffusionMiscCmdApps OFF CACHE BOOL "Build miscellaneous commandline tools for diffusion" FORCE) set(BUILD_DiffusionQuantificationCmdApps OFF CACHE BOOL "Build commandline tools for diffusion quantification (IVIM, ADC, ...)" FORCE) set(BUILD_DiffusionTractographyCmdApps ON CACHE BOOL "Build commandline tools for diffusion fiber tractography" FORCE) set(BUILD_DiffusionTractographyEvaluationCmdApps OFF CACHE BOOL "Build commandline tools for diffusion fiber tractography evaluation" FORCE) set(BUILD_DiffusionConnectomicsCmdApps ON CACHE BOOL "Build commandline tools for diffusion connectomics" FORCE) # Build neither all plugins nor examples set(MITK_BUILD_ALL_PLUGINS OFF CACHE BOOL "Build all MITK plugins" FORCE) set(MITK_BUILD_EXAMPLES OFF CACHE BOOL "Build the MITK examples" FORCE) set(BUILD_TESTING OFF CACHE BOOL "Build the MITK tests" FORCE) # Activate in-application help generation set(MITK_DOXYGEN_GENERATE_QCH_FILES ON CACHE BOOL "Use doxygen to generate Qt compressed help files for MITK docs" FORCE) set(BLUEBERRY_USE_QT_HELP ON CACHE BOOL "Enable support for integrating bundle documentation into Qt Help" FORCE) # Enable console window set(MITK_SHOW_CONSOLE_WINDOW ON CACHE BOOL "Use this to enable or disable the console window when starting MITK GUI Applications" FORCE) set(CMAKE_BUILD_TYPE Release CACHE STRING "" FORCE) diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetric.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetric.h new file mode 100644 index 0000000000..17fac83341 --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetric.h @@ -0,0 +1,61 @@ +/*=================================================================== + +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 _ClusteringMetric +#define _ClusteringMetric + +#include + +namespace mitk +{ + +/** +* \brief Base class for fiber clustering metrics */ + +class MITKFIBERTRACKING_EXPORT ClusteringMetric +{ + +public: + + ClusteringMetric() + : m_Scale(1.0) + {} + virtual ~ClusteringMetric(){} + + virtual float CalculateDistance(vnl_matrix& s, vnl_matrix& t, bool &flipped) = 0; + + float SetScale() const; + void SetScale(float Scale); + +protected: + + float m_Scale; + +}; + +float ClusteringMetric::SetScale() const +{ + return m_Scale; +} + +void ClusteringMetric::SetScale(float Scale) +{ + m_Scale = Scale; +} + +} + +#endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetricAnatomic.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetricAnatomic.h new file mode 100644 index 0000000000..9e8db560d4 --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetricAnatomic.h @@ -0,0 +1,199 @@ +/*=================================================================== + +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 _ClusteringMetricMdfAnatomic +#define _ClusteringMetricMdfAnatomic + +#include +#include +#include +#include + +namespace mitk +{ + +/** +* \brief Fiber clustering metric based on white matter parcellation histograms along the tracts (Siless et al. https://www.ncbi.nlm.nih.gov/pubmed/29100937) */ + +class MITKFIBERTRACKING_EXPORT ClusteringMetricAnatomic : public ClusteringMetric +{ + +public: + + typedef itk::Image ItkShortImgType; + + ClusteringMetricAnatomic() + : m_Radius(1) + {} + virtual ~ClusteringMetricAnatomic(){} + + float CalculateDistance(vnl_matrix& s, vnl_matrix& t, bool &flipped) + { + vnl_matrix hist1; hist1.set_size( (2*m_Radius+1)*(2*m_Radius+1)*(2*m_Radius+1), m_NumLabels ); hist1.fill(0); + vnl_matrix hist2; hist2.set_size( (2*m_Radius+1)*(2*m_Radius+1)*(2*m_Radius+1), m_NumLabels ); hist2.fill(0); + + + float d_direct = 0; + float d_flipped = 0; + + vnl_vector dists_d; dists_d.set_size(s.cols()); + vnl_vector dists_f; dists_f.set_size(s.cols()); + + for (unsigned int i=0; i p; + p[0] = s[0][i]; + p[1] = s[1][i]; + p[2] = s[2][i]; + GetHistogramAtPoint(p, hist1); + + p[0] = t[0][i]; + p[1] = t[1][i]; + p[2] = t[2][i]; + GetHistogramAtPoint(p, hist2); + + d_direct += (s.get_column(i)-t.get_column(i)).magnitude(); + d_flipped += (s.get_column(i)-t.get_column(s.cols()-i-1)).magnitude(); + } + +// float eudist = 0; + if (d_direct>d_flipped) + { + flipped = true; +// eudist = d_flipped/s.cols(); + } + else + { + flipped = false; +// eudist = d_direct/s.cols(); + } + + float label_intersection = 0; + for (unsigned int c=0; c0) + l1 = true; + if (hist2[r][c]>0) + l2 = true; + if (l1 && l2) + { + label_intersection += 1; + break; + } + } + } + + float similarity = 0; + if (label_intersection>0) + { + hist1.normalize_rows(); + hist2.normalize_rows(); + label_intersection /= m_NumLabels; + + for (unsigned int l=0; l0) + return m_Scale*0.2/similarity; + else + return 9999; + } + + void GetHistogramAtPoint(itk::Point& itkP, vnl_matrix& hist) + { + int parc_idx = 0; + for (auto parc : m_Parcellations) + { + int dir_idx=0; + itk::Index<3> tmp_idx; + itk::Index<3> idx; + parc->TransformPhysicalPointToIndex(itkP, idx); + if (!parc->GetLargestPossibleRegion().IsInside(idx)) + continue; + + short label = parc->GetPixel(idx); + short S0 = label; + hist[dir_idx][m_LabelMaps.at(parc_idx).at(label)] += 1; + + for (int x=-m_Radius; x<=m_Radius; ++x) + for (int y=-m_Radius; y<=m_Radius; ++y) + for (int z=-m_Radius; z<=m_Radius; ++z) + { + if (x==0 && y==0 && z==0) + continue; + + ++dir_idx; + + for (int d=1; d<5; ++d) + { + tmp_idx[0] = idx[0] + x*d; + tmp_idx[1] = idx[1] + y*d; + tmp_idx[2] = idx[2] + z*d; + if (!parc->GetLargestPossibleRegion().IsInside(tmp_idx)) + break; + label = parc->GetPixel(tmp_idx); + + if (label!=S0) + { + hist[dir_idx][m_LabelMaps.at(parc_idx).at(label)] += 1; + break; + } + } + } + ++parc_idx; + } + } + + void SetParcellations(const std::vector &Parcellations) + { + m_Parcellations = Parcellations; + m_NumLabels = 0; + + for (auto parc : m_Parcellations) + { + std::map< short, short > label_map; + itk::ImageRegionConstIterator it(parc, parc->GetLargestPossibleRegion()); + while (!it.IsAtEnd()) + { + if (label_map.count(it.Get())==0) + { + label_map.insert( std::pair< short, short >( it.Get(), m_NumLabels) ); + ++m_NumLabels; + } + ++it; + } + m_LabelMaps.push_back(label_map); + } + } + +protected: + + std::vector< ItkShortImgType::Pointer > m_Parcellations; + short m_NumLabels; + int m_Radius; + std::vector< std::map< short, short > > m_LabelMaps; +}; + +} + +#endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetricEuclideanMax.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetricEuclideanMax.h new file mode 100644 index 0000000000..cd2ffad9f0 --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetricEuclideanMax.h @@ -0,0 +1,71 @@ +/*=================================================================== + +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 _ClusteringMetricMdfMax +#define _ClusteringMetricMdfMax + +#include + +namespace mitk +{ + +/** +* \brief Fiber clustering metric based on the euclidean maximum distance between tracts */ + +class MITKFIBERTRACKING_EXPORT ClusteringMetricEuclideanMax : public ClusteringMetric +{ + +public: + + ClusteringMetricEuclideanMax(){} + virtual ~ClusteringMetricEuclideanMax(){} + + float CalculateDistance(vnl_matrix& s, vnl_matrix& t, bool &flipped) + { + float d_direct = 0; + float d_flipped = 0; + + vnl_vector dists_d; dists_d.set_size(s.cols()); + vnl_vector dists_f; dists_f.set_size(s.cols()); + + for (unsigned int i=0; id_flipped) + { + float d = dists_f.max_value(); + flipped = true; + return m_Scale*d; + } + + float d = dists_d.max_value(); + flipped = false; + return m_Scale*d; + } + +protected: + +}; + +} + +#endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetricEuclideanMean.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetricEuclideanMean.h new file mode 100644 index 0000000000..dd0d6bf624 --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetricEuclideanMean.h @@ -0,0 +1,62 @@ +/*=================================================================== + +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 _ClusteringMetricMdf +#define _ClusteringMetricMdf + +#include + +namespace mitk +{ + +/** +* \brief Fiber clustering metric based on the mean euclidean distance between tracts */ + +class MITKFIBERTRACKING_EXPORT ClusteringMetricEuclideanMean : public ClusteringMetric +{ + +public: + + ClusteringMetricEuclideanMean(){} + virtual ~ClusteringMetricEuclideanMean(){} + + float CalculateDistance(vnl_matrix& s, vnl_matrix& t, bool &flipped) + { + float d_direct = 0; + float d_flipped = 0; + + for (unsigned int i=0; id_flipped) + { + flipped = true; + return m_Scale*d_flipped/s.cols(); + } + flipped = false; + return m_Scale*d_direct/s.cols(); + } + +protected: + +}; + +} + +#endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetricEuclideanStd.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetricEuclideanStd.h new file mode 100644 index 0000000000..7314bd1c88 --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetricEuclideanStd.h @@ -0,0 +1,77 @@ +/*=================================================================== + +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 _ClusteringMetricMdfStd +#define _ClusteringMetricMdfStd + +#include + +namespace mitk +{ + +/** +* \brief Fiber clustering metric based on the euclidean distance between tracts and the corresponding standard deviation of the distances */ + +class MITKFIBERTRACKING_EXPORT ClusteringMetricEuclideanStd : public ClusteringMetric +{ + +public: + + ClusteringMetricEuclideanStd(){} + virtual ~ClusteringMetricEuclideanStd(){} + + float CalculateDistance(vnl_matrix& s, vnl_matrix& t, bool &flipped) + { + float d_direct = 0; + float d_flipped = 0; + + vnl_vector dists_d; dists_d.set_size(s.cols()); + vnl_vector dists_f; dists_f.set_size(s.cols()); + + for (unsigned int i=0; id_flipped) + { + float d = d_flipped/s.cols(); + dists_f -= d; + d += dists_f.magnitude(); + + flipped = true; + return m_Scale*d/2; + } + + float d = d_direct/s.cols(); + dists_d -= d; + d += dists_d.magnitude(); + + flipped = false; + return m_Scale*d/2; + } + +protected: + +}; + +} + +#endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetricScalarMap.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetricScalarMap.h new file mode 100644 index 0000000000..736d705c78 --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetricScalarMap.h @@ -0,0 +1,126 @@ +/*=================================================================== + +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 _ClusteringMetricScalarMap +#define _ClusteringMetricScalarMap + +#include +#include +#include +#include +#include + +namespace mitk +{ + +/** +* \brief Fiber clustering metric based on the scalar image values along a tract */ + +class MITKFIBERTRACKING_EXPORT ClusteringMetricScalarMap : public ClusteringMetric +{ + +public: + + typedef itk::Image ItkFloatImgType; + + ClusteringMetricScalarMap() + {} + virtual ~ClusteringMetricScalarMap(){} + + float CalculateDistance(vnl_matrix& s, vnl_matrix& t, bool &flipped) + { + float d_direct = 0; + float d_flipped = 0; + float map_distance = 0; + + vnl_vector dists_d; dists_d.set_size(s.cols()); + vnl_vector dists_f; dists_f.set_size(s.cols()); + + for (unsigned int i=0; id_flipped) + { + flipped = true; + + for (unsigned int i=0; i p; + p[0] = s[0][i]; + p[1] = s[1][i]; + p[2] = s[2][i]; + vnl_vector vals1 = GetImageValuesAtPoint(p); + + p[0] = t[0][s.cols()-i-1]; + p[1] = t[1][s.cols()-i-1]; + p[2] = t[2][s.cols()-i-1]; + vnl_vector vals2 = GetImageValuesAtPoint(p); + + map_distance += (vals1-vals2).magnitude(); + } + } + else + { + flipped = false; + + for (unsigned int i=0; i p; + p[0] = s[0][i]; + p[1] = s[1][i]; + p[2] = s[2][i]; + vnl_vector vals1 = GetImageValuesAtPoint(p); + + p[0] = t[0][i]; + p[1] = t[1][i]; + p[2] = t[2][i]; + vnl_vector vals2 = GetImageValuesAtPoint(p); + + map_distance += (vals1-vals2).magnitude(); + } + } + + return m_Scale*map_distance; + } + + vnl_vector GetImageValuesAtPoint(itk::Point& itkP) + { + vnl_vector vals; vals.set_size(m_ScalarMaps.size()); + int c = 0; + for (auto map : m_ScalarMaps) + { + vals[c] = mitk::TrackingDataHandler::GetImageValue(itkP, map, true); + ++c; + } + return vals; + } + + void SetImages(const std::vector &Parcellations) + { + m_ScalarMaps = Parcellations; + } + +protected: + + std::vector< ItkFloatImgType::Pointer > m_ScalarMaps; +}; + +} + +#endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingDataHandler.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingDataHandler.h index b7d3749d54..be789cb90d 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingDataHandler.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingDataHandler.h @@ -1,382 +1,401 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef _TrackingDataHandler #define _TrackingDataHandler #include #include #include #include #include #include #include #include namespace mitk { /** * \brief Abstract class for tracking handler. A tracking handler deals with determining the next progression direction of a streamline fiber. There are different handlers for tensor images, peak images, ... */ class MITKFIBERTRACKING_EXPORT TrackingDataHandler { public: enum MODE { DETERMINISTIC, PROBABILISTIC }; TrackingDataHandler(); virtual ~TrackingDataHandler(){} typedef itk::Statistics::MersenneTwisterRandomVariateGenerator ItkRngType; typedef boost::mt19937 BoostRngType; typedef itk::Image ItkUcharImgType; typedef itk::Image ItkShortImgType; typedef itk::Image ItkFloatImgType; typedef itk::Image ItkDoubleImgType; typedef vnl_vector_fixed< float, 3 > TrackingDirectionType; virtual TrackingDirectionType ProposeDirection(const itk::Point& pos, std::deque< TrackingDirectionType >& olddirs, itk::Index<3>& oldIndex) = 0; ///< predicts next progression direction at the given position virtual void InitForTracking() = 0; virtual itk::Vector GetSpacing() = 0; virtual itk::Point GetOrigin() = 0; virtual itk::Matrix GetDirection() = 0; virtual itk::ImageRegion<3> GetLargestPossibleRegion() = 0; virtual bool WorldToIndex(itk::Point& pos, itk::Index<3>& index) = 0; virtual void SetMode(MODE m) = 0; MODE GetMode(){ return m_Mode; } void 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 ) { m_Random = random; if (!random) { m_Rng.seed(0); std::srand(0); m_RngItk->SetSeed(0); } else { m_Rng.seed(); m_RngItk->SetSeed(); std::srand(std::time(0)); } } double GetRandDouble(const double & a, const double & b) { return m_RngItk->GetUniformVariate(a, b); } + static bool IsInsideMask(const itk::Point &pos, ItkUcharImgType::Pointer mask, bool interpolate) + { + if (mask.IsNull()) + return true; + + if (interpolate) + { + if ( mitk::TrackingDataHandler::GetImageValue(pos, mask, true)<0.5) + return false; + } + else + { + if ( mitk::TrackingDataHandler::GetImageValue(pos, mask, false)==0) + return false; + } + + return true; + } + template< class TPixelType > - TPixelType GetImageValue(itk::Point itkP, itk::Image* image, vnl_vector_fixed& interpWeights){ + static TPixelType GetImageValue(const itk::Point& itkP, itk::Image* image, vnl_vector_fixed& interpWeights){ // transform physical point to index coordinates itk::Index<3> idx; itk::ContinuousIndex< float, 3> cIdx; image->TransformPhysicalPointToIndex(itkP, idx); image->TransformPhysicalPointToContinuousIndex(itkP, cIdx); TPixelType pix = 0.0; if ( image->GetLargestPossibleRegion().IsInside(idx) ) { pix = image->GetPixel(idx); } else return pix; float frac_x = cIdx[0] - idx[0]; float frac_y = cIdx[1] - idx[1]; float frac_z = cIdx[2] - idx[2]; if (frac_x<0) { idx[0] -= 1; frac_x += 1; } if (frac_y<0) { idx[1] -= 1; frac_y += 1; } if (frac_z<0) { idx[2] -= 1; frac_z += 1; } frac_x = 1-frac_x; frac_y = 1-frac_y; frac_z = 1-frac_z; // int coordinates inside image? if (idx[0] >= 0 && idx[0] < image->GetLargestPossibleRegion().GetSize(0)-1 && idx[1] >= 0 && idx[1] < image->GetLargestPossibleRegion().GetSize(1)-1 && idx[2] >= 0 && idx[2] < image->GetLargestPossibleRegion().GetSize(2)-1) { // trilinear interpolation interpWeights[0] = ( frac_x)*( frac_y)*( frac_z); interpWeights[1] = (1-frac_x)*( frac_y)*( frac_z); interpWeights[2] = ( frac_x)*(1-frac_y)*( frac_z); interpWeights[3] = ( frac_x)*( frac_y)*(1-frac_z); interpWeights[4] = (1-frac_x)*(1-frac_y)*( frac_z); interpWeights[5] = ( frac_x)*(1-frac_y)*(1-frac_z); interpWeights[6] = (1-frac_x)*( frac_y)*(1-frac_z); interpWeights[7] = (1-frac_x)*(1-frac_y)*(1-frac_z); pix = image->GetPixel(idx) * interpWeights[0]; typename itk::Image::IndexType tmpIdx = idx; tmpIdx[0]++; pix += image->GetPixel(tmpIdx) * interpWeights[1]; tmpIdx = idx; tmpIdx[1]++; pix += image->GetPixel(tmpIdx) * interpWeights[2]; tmpIdx = idx; tmpIdx[2]++; pix += image->GetPixel(tmpIdx) * interpWeights[3]; tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++; pix += image->GetPixel(tmpIdx) * interpWeights[4]; tmpIdx = idx; tmpIdx[1]++; tmpIdx[2]++; pix += image->GetPixel(tmpIdx) * interpWeights[5]; tmpIdx = idx; tmpIdx[2]++; tmpIdx[0]++; pix += image->GetPixel(tmpIdx) * interpWeights[6]; tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++; tmpIdx[2]++; pix += image->GetPixel(tmpIdx) * interpWeights[7]; } if (pix!=pix) { pix = 0; MITK_WARN << "nan values in image!"; } return pix; } template< class TPixelType, class TOutPixelType=TPixelType > - TPixelType GetImageValue(itk::Point itkP, itk::Image* image, bool interpolate){ + static TPixelType GetImageValue(const itk::Point& itkP, itk::Image* image, bool interpolate){ // transform physical point to index coordinates itk::Index<3> idx; itk::ContinuousIndex< float, 3> cIdx; image->TransformPhysicalPointToIndex(itkP, idx); image->TransformPhysicalPointToContinuousIndex(itkP, cIdx); TOutPixelType pix = 0.0; if ( image->GetLargestPossibleRegion().IsInside(idx) ) { pix = (TOutPixelType)image->GetPixel(idx); if (!interpolate) return pix; } else return pix; float frac_x = cIdx[0] - idx[0]; float frac_y = cIdx[1] - idx[1]; float frac_z = cIdx[2] - idx[2]; if (frac_x<0) { idx[0] -= 1; frac_x += 1; } if (frac_y<0) { idx[1] -= 1; frac_y += 1; } if (frac_z<0) { idx[2] -= 1; frac_z += 1; } frac_x = 1-frac_x; frac_y = 1-frac_y; frac_z = 1-frac_z; // int coordinates inside image? if (idx[0] >= 0 && idx[0] < static_cast(image->GetLargestPossibleRegion().GetSize(0) - 1) && idx[1] >= 0 && idx[1] < static_cast(image->GetLargestPossibleRegion().GetSize(1) - 1) && idx[2] >= 0 && idx[2] < static_cast(image->GetLargestPossibleRegion().GetSize(2) - 1)) { // trilinear interpolation vnl_vector_fixed interpWeights; interpWeights[0] = ( frac_x)*( frac_y)*( frac_z); interpWeights[1] = (1-frac_x)*( frac_y)*( frac_z); interpWeights[2] = ( frac_x)*(1-frac_y)*( frac_z); interpWeights[3] = ( frac_x)*( frac_y)*(1-frac_z); interpWeights[4] = (1-frac_x)*(1-frac_y)*( frac_z); interpWeights[5] = ( frac_x)*(1-frac_y)*(1-frac_z); interpWeights[6] = (1-frac_x)*( frac_y)*(1-frac_z); interpWeights[7] = (1-frac_x)*(1-frac_y)*(1-frac_z); pix = image->GetPixel(idx) * interpWeights[0]; typename itk::Image::IndexType tmpIdx = idx; tmpIdx[0]++; pix += (TOutPixelType)image->GetPixel(tmpIdx) * interpWeights[1]; tmpIdx = idx; tmpIdx[1]++; pix += (TOutPixelType)image->GetPixel(tmpIdx) * interpWeights[2]; tmpIdx = idx; tmpIdx[2]++; pix += (TOutPixelType)image->GetPixel(tmpIdx) * interpWeights[3]; tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++; pix += (TOutPixelType)image->GetPixel(tmpIdx) * interpWeights[4]; tmpIdx = idx; tmpIdx[1]++; tmpIdx[2]++; pix += (TOutPixelType)image->GetPixel(tmpIdx) * interpWeights[5]; tmpIdx = idx; tmpIdx[2]++; tmpIdx[0]++; pix += (TOutPixelType)image->GetPixel(tmpIdx) * interpWeights[6]; tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++; tmpIdx[2]++; pix += (TOutPixelType)image->GetPixel(tmpIdx) * interpWeights[7]; } if (pix!=pix) { pix = 0; MITK_WARN << "nan values in image!"; } return pix; } template< class TPixelType, int components > - itk::Vector< TPixelType, components > GetImageValue(itk::Point itkP, itk::Image, 3>* image, bool interpolate){ + static itk::Vector< TPixelType, components > GetImageValue(const itk::Point& itkP, itk::Image, 3>* image, bool interpolate){ // transform physical point to index coordinates itk::Index<3> idx; itk::ContinuousIndex< float, 3> cIdx; image->TransformPhysicalPointToIndex(itkP, idx); image->TransformPhysicalPointToContinuousIndex(itkP, cIdx); itk::Vector< TPixelType, components > pix; pix.Fill(0.0); if ( image->GetLargestPossibleRegion().IsInside(idx) ) { pix = image->GetPixel(idx); if (!interpolate) return pix; } else return pix; float frac_x = cIdx[0] - idx[0]; float frac_y = cIdx[1] - idx[1]; float frac_z = cIdx[2] - idx[2]; if (frac_x<0) { idx[0] -= 1; frac_x += 1; } if (frac_y<0) { idx[1] -= 1; frac_y += 1; } if (frac_z<0) { idx[2] -= 1; frac_z += 1; } frac_x = 1-frac_x; frac_y = 1-frac_y; frac_z = 1-frac_z; // int coordinates inside image? if (idx[0] >= 0 && idx[0] < static_cast(image->GetLargestPossibleRegion().GetSize(0) - 1) && idx[1] >= 0 && idx[1] < static_cast(image->GetLargestPossibleRegion().GetSize(1) - 1) && idx[2] >= 0 && idx[2] < static_cast(image->GetLargestPossibleRegion().GetSize(2) - 1)) { // trilinear interpolation vnl_vector_fixed interpWeights; interpWeights[0] = ( frac_x)*( frac_y)*( frac_z); interpWeights[1] = (1-frac_x)*( frac_y)*( frac_z); interpWeights[2] = ( frac_x)*(1-frac_y)*( frac_z); interpWeights[3] = ( frac_x)*( frac_y)*(1-frac_z); interpWeights[4] = (1-frac_x)*(1-frac_y)*( frac_z); interpWeights[5] = ( frac_x)*(1-frac_y)*(1-frac_z); interpWeights[6] = (1-frac_x)*( frac_y)*(1-frac_z); interpWeights[7] = (1-frac_x)*(1-frac_y)*(1-frac_z); pix = image->GetPixel(idx) * interpWeights[0]; typename itk::Image, 3>::IndexType tmpIdx = idx; tmpIdx[0]++; pix += image->GetPixel(tmpIdx) * interpWeights[1]; tmpIdx = idx; tmpIdx[1]++; pix += image->GetPixel(tmpIdx) * interpWeights[2]; tmpIdx = idx; tmpIdx[2]++; pix += image->GetPixel(tmpIdx) * interpWeights[3]; tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++; pix += image->GetPixel(tmpIdx) * interpWeights[4]; tmpIdx = idx; tmpIdx[1]++; tmpIdx[2]++; pix += image->GetPixel(tmpIdx) * interpWeights[5]; tmpIdx = idx; tmpIdx[2]++; tmpIdx[0]++; pix += image->GetPixel(tmpIdx) * interpWeights[6]; tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++; tmpIdx[2]++; pix += image->GetPixel(tmpIdx) * interpWeights[7]; } if (pix!=pix) { pix.Fill(0.0); mitkThrow() << "nan values in image!"; } return pix; } protected: float m_AngularThreshold; bool m_Interpolate; bool m_FlipX; bool m_FlipY; bool m_FlipZ; MODE m_Mode; BoostRngType m_Rng; ItkRngType::Pointer m_RngItk; bool m_NeedsDataInit; bool m_Random; void DataModified() { m_NeedsDataInit = true; } }; } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp index 1e23175c0a..ece114f169 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp @@ -1,1108 +1,1094 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include "itkStreamlineTrackingFilter.h" #include #include #include #include #include "itkPointShell.h" #include #include #include #include #include #include #define _USE_MATH_DEFINES #include namespace itk { StreamlineTrackingFilter ::StreamlineTrackingFilter() : m_PauseTracking(false) , m_AbortTracking(false) , m_BuildFibersFinished(false) , m_BuildFibersReady(0) , m_FiberPolyData(nullptr) , m_Points(nullptr) , m_Cells(nullptr) , m_StoppingRegions(nullptr) , m_TargetRegions(nullptr) , m_SeedImage(nullptr) , m_MaskImage(nullptr) , m_TissueImage(nullptr) , m_OutputProbabilityMap(nullptr) , m_AngularThresholdDeg(-1) , m_StepSizeVox(-1) , m_SamplingDistanceVox(-1) , m_AngularThreshold(-1) , m_StepSize(0) , m_MaxLength(10000) , m_MinTractLength(20.0) , m_MaxTractLength(400.0) , m_SeedsPerVoxel(1) , m_AvoidStop(true) , m_RandomSampling(false) , m_SamplingDistance(-1) , m_DeflectionMod(1.0) , m_OnlyForwardSamples(true) , m_UseStopVotes(true) , m_NumberOfSamples(30) , m_NumPreviousDirections(1) , m_WmLabel(3) // mrtrix 5ttseg labels , m_GmLabel(1) // mrtrix 5ttseg labels , m_SeedOnlyGm(false) , m_ControlGmEndings(false) , m_MaxNumTracts(-1) , m_Verbose(true) , m_AposterioriCurvCheck(false) , m_DemoMode(false) , m_Random(true) , m_UseOutputProbabilityMap(false) , m_CurrentTracts(0) , m_Progress(0) , m_StopTracking(false) , m_InterpolateMask(false) { this->SetNumberOfRequiredInputs(0); } std::string StreamlineTrackingFilter::GetStatusText() { std::string status = "Seedpoints processed: " + boost::lexical_cast(m_Progress) + "/" + boost::lexical_cast(m_SeedPoints.size()); if (m_SeedPoints.size()>0) status += " (" + boost::lexical_cast(100*m_Progress/m_SeedPoints.size()) + "%)"; if (m_MaxNumTracts>0) status += "\nFibers accepted: " + boost::lexical_cast(m_CurrentTracts) + "/" + boost::lexical_cast(m_MaxNumTracts); else status += "\nFibers accepted: " + boost::lexical_cast(m_CurrentTracts); return status; } void StreamlineTrackingFilter::BeforeTracking() { m_StopTracking = false; m_TrackingHandler->SetRandom(m_Random); m_TrackingHandler->InitForTracking(); m_FiberPolyData = PolyDataType::New(); m_Points = vtkSmartPointer< vtkPoints >::New(); m_Cells = vtkSmartPointer< vtkCellArray >::New(); itk::Vector< double, 3 > imageSpacing = m_TrackingHandler->GetSpacing(); float minSpacing; if(imageSpacing[0]SetAngularThreshold(m_AngularThreshold); if (m_SamplingDistanceVoxGetNumberOfThreads(); i++) { PolyDataType poly = PolyDataType::New(); m_PolyDataContainer.push_back(poly); } if (m_UseOutputProbabilityMap) { m_OutputProbabilityMap = ItkDoubleImgType::New(); m_OutputProbabilityMap->SetSpacing(imageSpacing); m_OutputProbabilityMap->SetOrigin(m_TrackingHandler->GetOrigin()); m_OutputProbabilityMap->SetDirection(m_TrackingHandler->GetDirection()); m_OutputProbabilityMap->SetRegions(m_TrackingHandler->GetLargestPossibleRegion()); m_OutputProbabilityMap->Allocate(); m_OutputProbabilityMap->FillBuffer(0); } if (m_StoppingRegions.IsNull()) { m_StoppingRegions = ItkUcharImgType::New(); m_StoppingRegions->SetSpacing( imageSpacing ); m_StoppingRegions->SetOrigin( m_TrackingHandler->GetOrigin() ); m_StoppingRegions->SetDirection( m_TrackingHandler->GetDirection() ); m_StoppingRegions->SetRegions( m_TrackingHandler->GetLargestPossibleRegion() ); m_StoppingRegions->Allocate(); m_StoppingRegions->FillBuffer(0); } else std::cout << "StreamlineTracking - Using stopping region image" << std::endl; if (m_TargetRegions.IsNull()) { m_TargetRegions = ItkUintImgType::New(); m_TargetRegions->SetSpacing( imageSpacing ); m_TargetRegions->SetOrigin( m_TrackingHandler->GetOrigin() ); m_TargetRegions->SetDirection( m_TrackingHandler->GetDirection() ); m_TargetRegions->SetRegions( m_TrackingHandler->GetLargestPossibleRegion() ); m_TargetRegions->Allocate(); m_TargetRegions->FillBuffer(1); } else std::cout << "StreamlineTracking - Using target region image" << std::endl; if (m_SeedImage.IsNull()) { m_SeedImage = ItkUcharImgType::New(); m_SeedImage->SetSpacing( imageSpacing ); m_SeedImage->SetOrigin( m_TrackingHandler->GetOrigin() ); m_SeedImage->SetDirection( m_TrackingHandler->GetDirection() ); m_SeedImage->SetRegions( m_TrackingHandler->GetLargestPossibleRegion() ); m_SeedImage->Allocate(); m_SeedImage->FillBuffer(1); } else std::cout << "StreamlineTracking - Using seed image" << std::endl; if (m_MaskImage.IsNull()) { // initialize mask image m_MaskImage = ItkUcharImgType::New(); m_MaskImage->SetSpacing( imageSpacing ); m_MaskImage->SetOrigin( m_TrackingHandler->GetOrigin() ); m_MaskImage->SetDirection( m_TrackingHandler->GetDirection() ); m_MaskImage->SetRegions( m_TrackingHandler->GetLargestPossibleRegion() ); m_MaskImage->Allocate(); m_MaskImage->FillBuffer(1); } else std::cout << "StreamlineTracking - Using mask image" << std::endl; if (m_SeedPoints.empty()) GetSeedPointsFromSeedImage(); else m_SeedOnlyGm = false; if (m_TissueImage.IsNull()) { if (m_SeedOnlyGm) { MITK_WARN << "StreamlineTracking - Cannot seed in gray matter. No tissue type image set."; m_SeedOnlyGm = false; } if (m_ControlGmEndings) { MITK_WARN << "StreamlineTracking - Cannot control gray matter endings. No tissue type image set."; m_ControlGmEndings = false; } } else { if (m_ControlGmEndings) m_SeedOnlyGm = true; if (m_ControlGmEndings || m_SeedOnlyGm) std::cout << "StreamlineTracking - Using tissue image" << std::endl; else MITK_WARN << "StreamlineTracking - Tissue image set but no gray matter seeding or fiber endpoint-control enabled" << std::endl; } m_BuildFibersReady = 0; m_BuildFibersFinished = false; m_Tractogram.clear(); m_SamplingPointset = mitk::PointSet::New(); m_AlternativePointset = mitk::PointSet::New(); m_StopVotePointset = mitk::PointSet::New(); m_StartTime = std::chrono::system_clock::now(); if (m_SeedOnlyGm && m_ControlGmEndings) InitGrayMatterEndings(); if (m_DemoMode) omp_set_num_threads(1); if (m_TrackingHandler->GetMode()==mitk::TrackingDataHandler::MODE::DETERMINISTIC) std::cout << "StreamlineTracking - Mode: deterministic" << std::endl; else if(m_TrackingHandler->GetMode()==mitk::TrackingDataHandler::MODE::PROBABILISTIC) std::cout << "StreamlineTracking - Mode: probabilistic" << std::endl; else std::cout << "StreamlineTracking - Mode: ???" << std::endl; std::cout << "StreamlineTracking - Angular threshold: " << m_AngularThreshold << " (" << 180*std::acos( m_AngularThreshold )/M_PI << "°)" << std::endl; std::cout << "StreamlineTracking - Stepsize: " << m_StepSize << "mm (" << m_StepSize/minSpacing << "*vox)" << std::endl; std::cout << "StreamlineTracking - Seeds per voxel: " << m_SeedsPerVoxel << std::endl; std::cout << "StreamlineTracking - Max. tract length: " << m_MaxTractLength << "mm" << std::endl; std::cout << "StreamlineTracking - Min. tract length: " << m_MinTractLength << "mm" << std::endl; std::cout << "StreamlineTracking - Max. num. tracts: " << m_MaxNumTracts << std::endl; std::cout << "StreamlineTracking - Num. neighborhood samples: " << m_NumberOfSamples << std::endl; std::cout << "StreamlineTracking - Max. sampling distance: " << m_SamplingDistance << "mm (" << m_SamplingDistance/minSpacing << "*vox)" << std::endl; std::cout << "StreamlineTracking - Deflection modifier: " << m_DeflectionMod << std::endl; std::cout << "StreamlineTracking - Use stop votes: " << m_UseStopVotes << std::endl; std::cout << "StreamlineTracking - Only frontal samples: " << m_OnlyForwardSamples << std::endl; if (m_DemoMode) { std::cout << "StreamlineTracking - Running in demo mode"; std::cout << "StreamlineTracking - Starting streamline tracking using 1 thread" << std::endl; } else std::cout << "StreamlineTracking - Starting streamline tracking using " << omp_get_max_threads() << " threads" << std::endl; } void StreamlineTrackingFilter::InitGrayMatterEndings() { m_TrackingHandler->SetAngularThreshold(0); m_GmStubs.clear(); if (m_TissueImage.IsNotNull()) { std::cout << "StreamlineTracking - initializing GM endings" << std::endl; ImageRegionConstIterator< ItkUcharImgType > it(m_TissueImage, m_TissueImage->GetLargestPossibleRegion() ); it.GoToBegin(); vnl_vector_fixed d1; d1.fill(0.0); std::deque< vnl_vector_fixed > olddirs; while (olddirs.size() start; m_TissueImage->TransformIndexToPhysicalPoint(s_idx, start); itk::Point wm_p; float max = -1; FiberType fib; for (int x : {-1,0,1}) for (int y : {-1,0,1}) for (int z : {-1,0,1}) { if (x==y && y==z) continue; ItkUcharImgType::IndexType e_idx; e_idx[0] = s_idx[0] + x; e_idx[1] = s_idx[1] + y; e_idx[2] = s_idx[2] + z; if ( !m_TissueImage->GetLargestPossibleRegion().IsInside(e_idx) || m_TissueImage->GetPixel(e_idx)!=m_WmLabel ) continue; itk::ContinuousIndex end; m_TissueImage->TransformIndexToPhysicalPoint(e_idx, end); d1 = m_TrackingHandler->ProposeDirection(end, olddirs, s_idx); if (d1.magnitude()<0.0001) continue; d1.normalize(); vnl_vector_fixed< float, 3 > d2; d2[0] = end[0] - start[0]; d2[1] = end[1] - start[1]; d2[2] = end[2] - start[2]; d2.normalize(); float a = fabs(dot_product(d1,d2)); if (a>max) { max = a; wm_p = end; } } if (max>=0) { fib.push_back(start); fib.push_back(wm_p); m_GmStubs.push_back(fib); } } ++it; } } m_TrackingHandler->SetAngularThreshold(m_AngularThreshold); } void StreamlineTrackingFilter::CalculateNewPosition(itk::Point& pos, vnl_vector_fixed& dir) { pos[0] += dir[0]*m_StepSize; pos[1] += dir[1]*m_StepSize; pos[2] += dir[2]*m_StepSize; } bool StreamlineTrackingFilter ::IsInsideMask(const itk::Point &pos, ItkUcharImgType::Pointer mask) { - if (mask.IsNull()) - return true; - - if (m_InterpolateMask) - { - if (m_TrackingHandler->GetImageValue(pos, mask, m_InterpolateMask)<0.5) - return false; - } - else - { - if (m_TrackingHandler->GetImageValue(pos, mask, m_InterpolateMask)==0) - return false; - } - - return true; + return m_TrackingHandler->IsInsideMask(pos, mask, m_InterpolateMask); } bool StreamlineTrackingFilter ::IsInGm(const itk::Point &pos) { if (m_TissueImage.IsNull()) return true; ItkUcharImgType::IndexType idx; m_TissueImage->TransformPhysicalPointToIndex(pos, idx); if (m_TissueImage->GetLargestPossibleRegion().IsInside(idx) && m_TissueImage->GetPixel(idx) == m_GmLabel) return true; return false; } std::vector< vnl_vector_fixed > StreamlineTrackingFilter::CreateDirections(int NPoints) { std::vector< vnl_vector_fixed > pointshell; if (NPoints<2) return pointshell; std::vector< float > theta; theta.resize(NPoints); std::vector< float > phi; phi.resize(NPoints); float C = sqrt(4*M_PI); phi[0] = 0.0; phi[NPoints-1] = 0.0; for(int i=0; i0 && i d; d[0] = cos(theta[i]) * cos(phi[i]); d[1] = cos(theta[i]) * sin(phi[i]); d[2] = sin(theta[i]); pointshell.push_back(d); } return pointshell; } vnl_vector_fixed StreamlineTrackingFilter::GetNewDirection(itk::Point &pos, std::deque >& olddirs, itk::Index<3> &oldIndex) { if (m_DemoMode) { m_SamplingPointset->Clear(); m_AlternativePointset->Clear(); m_StopVotePointset->Clear(); } vnl_vector_fixed direction; direction.fill(0); if (IsInsideMask(pos, m_MaskImage) && !IsInsideMask(pos, m_StoppingRegions)) direction = m_TrackingHandler->ProposeDirection(pos, olddirs, oldIndex); // get direction proposal at current streamline position else return direction; vnl_vector_fixed olddir = olddirs.back(); std::vector< vnl_vector_fixed > probeVecs = CreateDirections(m_NumberOfSamples); itk::Point sample_pos; int alternatives = 1; int stop_votes = 0; int possible_stop_votes = 0; for (unsigned int i=0; i d; bool is_stop_voter = false; if (m_Random && m_RandomSampling) { d[0] = m_TrackingHandler->GetRandDouble(-0.5, 0.5); d[1] = m_TrackingHandler->GetRandDouble(-0.5, 0.5); d[2] = m_TrackingHandler->GetRandDouble(-0.5, 0.5); d.normalize(); d *= m_TrackingHandler->GetRandDouble(0,m_SamplingDistance); } else { d = probeVecs.at(i); float dot = dot_product(d, olddir); if (m_UseStopVotes && dot>0.7) { is_stop_voter = true; possible_stop_votes++; } else if (m_OnlyForwardSamples && dot<0) continue; d *= m_SamplingDistance; } sample_pos[0] = pos[0] + d[0]; sample_pos[1] = pos[1] + d[1]; sample_pos[2] = pos[2] + d[2]; vnl_vector_fixed tempDir; tempDir.fill(0.0); if (IsInsideMask(sample_pos, m_MaskImage)) tempDir = m_TrackingHandler->ProposeDirection(sample_pos, olddirs, oldIndex); // sample neighborhood if (tempDir.magnitude()>mitk::eps) { direction += tempDir; if(m_DemoMode) m_SamplingPointset->InsertPoint(i, sample_pos); } else if (m_AvoidStop && olddir.magnitude()>0.5) // out of white matter { if (is_stop_voter) stop_votes++; if (m_DemoMode) m_StopVotePointset->InsertPoint(i, sample_pos); float dot = dot_product(d, olddir); if (dot >= 0.0) // in front of plane defined by pos and olddir d = -d + 2*dot*olddir; // reflect else d = -d; // invert // look a bit further into the other direction sample_pos[0] = pos[0] + d[0]; sample_pos[1] = pos[1] + d[1]; sample_pos[2] = pos[2] + d[2]; alternatives++; vnl_vector_fixed tempDir; tempDir.fill(0.0); if (IsInsideMask(sample_pos, m_MaskImage)) tempDir = m_TrackingHandler->ProposeDirection(sample_pos, olddirs, oldIndex); // sample neighborhood if (tempDir.magnitude()>mitk::eps) // are we back in the white matter? { direction += d * m_DeflectionMod; // go into the direction of the white matter direction += tempDir; // go into the direction of the white matter direction at this location if(m_DemoMode) m_AlternativePointset->InsertPoint(alternatives, sample_pos); } else { if (m_DemoMode) m_StopVotePointset->InsertPoint(i, sample_pos); } } else { if (m_DemoMode) m_StopVotePointset->InsertPoint(i, sample_pos); if (is_stop_voter) stop_votes++; } } if (direction.magnitude()>0.001 && (possible_stop_votes==0 || (float)stop_votes/possible_stop_votes<0.5) ) direction.normalize(); else direction.fill(0); return direction; } float StreamlineTrackingFilter::FollowStreamline(itk::Point pos, vnl_vector_fixed dir, FiberType* fib, float tractLength, bool front) { vnl_vector_fixed zero_dir; zero_dir.fill(0.0); std::deque< vnl_vector_fixed > last_dirs; for (unsigned int i=0; i oldIndex; m_TrackingHandler->WorldToIndex(pos, oldIndex); // get new position CalculateNewPosition(pos, dir); // is new position inside of image and mask if (m_AbortTracking) // if not end streamline { return tractLength; } else // if yes, add new point to streamline { tractLength += m_StepSize; if (front) fib->push_front(pos); else fib->push_back(pos); if (m_AposterioriCurvCheck) { int curv = CheckCurvature(fib, front); if (curv>0) { tractLength -= m_StepSize*curv; while (curv>0) { if (front) fib->pop_front(); else fib->pop_back(); curv--; } return tractLength; } } if (tractLength>m_MaxTractLength) return tractLength; } if (m_DemoMode && !m_UseOutputProbabilityMap) // CHECK: warum sind die samplingpunkte der streamline in der visualisierung immer einen schritt voras? { #pragma omp critical { m_BuildFibersReady++; m_Tractogram.push_back(*fib); BuildFibers(true); m_Stop = true; while (m_Stop){ } } } dir.normalize(); last_dirs.push_back(dir); if (last_dirs.size()>m_NumPreviousDirections) last_dirs.pop_front(); dir = GetNewDirection(pos, last_dirs, oldIndex); while (m_PauseTracking){} if (dir.magnitude()<0.0001) return tractLength; } return tractLength; } int StreamlineTrackingFilter::CheckCurvature(FiberType* fib, bool front) { float m_Distance = 5; if (fib->size()<3) return 0; float dist = 0; std::vector< vnl_vector_fixed< float, 3 > > vectors; vnl_vector_fixed< float, 3 > meanV; meanV.fill(0); float dev = 0; if (front) { int c=0; while(distsize()-1) { itk::Point p1 = fib->at(c); itk::Point p2 = fib->at(c+1); vnl_vector_fixed< float, 3 > v; v[0] = p2[0]-p1[0]; v[1] = p2[1]-p1[1]; v[2] = p2[2]-p1[2]; dist += v.magnitude(); v.normalize(); vectors.push_back(v); if (c==0) meanV += v; c++; } } else { int c=fib->size()-1; while(dist0) { itk::Point p1 = fib->at(c); itk::Point p2 = fib->at(c-1); vnl_vector_fixed< float, 3 > v; v[0] = p2[0]-p1[0]; v[1] = p2[1]-p1[1]; v[2] = p2[2]-p1[2]; dist += v.magnitude(); v.normalize(); vectors.push_back(v); if (c==(int)fib->size()-1) meanV += v; c--; } } meanV.normalize(); for (unsigned int c=0; c1.0) angle = 1.0; if (angle<-1.0) angle = -1.0; dev += acos(angle)*180/M_PI; } if (vectors.size()>0) dev /= vectors.size(); if (dev<30) return 0; else return vectors.size(); } void StreamlineTrackingFilter::GetSeedPointsFromSeedImage() { MITK_INFO << "StreamlineTracking - Calculating seed points."; m_SeedPoints.clear(); if (!m_ControlGmEndings) { typedef ImageRegionConstIterator< ItkUcharImgType > MaskIteratorType; MaskIteratorType sit(m_SeedImage, m_SeedImage->GetLargestPossibleRegion()); sit.GoToBegin(); while (!sit.IsAtEnd()) { if (sit.Value()>0) { ItkUcharImgType::IndexType index = sit.GetIndex(); itk::ContinuousIndex start; start[0] = index[0]; start[1] = index[1]; start[2] = index[2]; itk::Point worldPos; m_SeedImage->TransformContinuousIndexToPhysicalPoint(start, worldPos); if (IsInsideMask(worldPos, m_MaskImage) && (!m_SeedOnlyGm || IsInGm(worldPos))) { m_SeedPoints.push_back(worldPos); for (int s = 1; s < m_SeedsPerVoxel; s++) { start[0] = index[0] + m_TrackingHandler->GetRandDouble(-0.5, 0.5); start[1] = index[1] + m_TrackingHandler->GetRandDouble(-0.5, 0.5); start[2] = index[2] + m_TrackingHandler->GetRandDouble(-0.5, 0.5); itk::Point worldPos; m_SeedImage->TransformContinuousIndexToPhysicalPoint(start, worldPos); m_SeedPoints.push_back(worldPos); } } } ++sit; } } else { for (auto s : m_GmStubs) m_SeedPoints.push_back(s[1]); } } void StreamlineTrackingFilter::GenerateData() { this->BeforeTracking(); if (m_Random) std::random_shuffle(m_SeedPoints.begin(), m_SeedPoints.end()); m_CurrentTracts = 0; int num_seeds = m_SeedPoints.size(); itk::Index<3> zeroIndex; zeroIndex.Fill(0); m_Progress = 0; int i = 0; int print_interval = num_seeds/100; if (print_interval<100) m_Verbose=false; #pragma omp parallel while (i=num_seeds || m_StopTracking) continue; else if (m_Verbose && i%print_interval==0) #pragma omp critical { m_Progress += print_interval; std::cout << " \r"; if (m_MaxNumTracts>0) std::cout << "Tried: " << m_Progress << "/" << num_seeds << " | Accepted: " << m_CurrentTracts << "/" << m_MaxNumTracts << '\r'; else std::cout << "Tried: " << m_Progress << "/" << num_seeds << " | Accepted: " << m_CurrentTracts << '\r'; cout.flush(); } const itk::Point worldPos = m_SeedPoints.at(temp_i); FiberType fib; float tractLength = 0; unsigned int counter = 0; // get starting direction vnl_vector_fixed dir; dir.fill(0.0); std::deque< vnl_vector_fixed > olddirs; while (olddirs.size() gm_start_dir; if (m_ControlGmEndings) { gm_start_dir[0] = m_GmStubs[temp_i][1][0] - m_GmStubs[temp_i][0][0]; gm_start_dir[1] = m_GmStubs[temp_i][1][1] - m_GmStubs[temp_i][0][1]; gm_start_dir[2] = m_GmStubs[temp_i][1][2] - m_GmStubs[temp_i][0][2]; gm_start_dir.normalize(); olddirs.pop_back(); olddirs.push_back(gm_start_dir); } if (IsInsideMask(worldPos, m_MaskImage)) dir = m_TrackingHandler->ProposeDirection(worldPos, olddirs, zeroIndex); if (dir.magnitude()>0.0001) { if (m_ControlGmEndings) { float a = dot_product(gm_start_dir, dir); if (a<0) dir = -dir; } // forward tracking tractLength = FollowStreamline(worldPos, dir, &fib, 0, false); fib.push_front(worldPos); if (m_ControlGmEndings) { fib.push_front(m_GmStubs[temp_i][0]); CheckFiberForGmEnding(&fib); } else { // backward tracking (only if we don't explicitely start in the GM) tractLength = FollowStreamline(worldPos, -dir, &fib, tractLength, true); if (m_ControlGmEndings) { CheckFiberForGmEnding(&fib); std::reverse(fib.begin(),fib.end()); CheckFiberForGmEnding(&fib); } } counter = fib.size(); if (tractLength>=m_MinTractLength && counter>=2) { ItkUintImgType::IndexType idx_begin, idx_end; m_TargetRegions->TransformPhysicalPointToIndex(fib.front(), idx_begin); m_TargetRegions->TransformPhysicalPointToIndex(fib.back(), idx_end); #pragma omp critical if ( !m_TargetRegions->GetLargestPossibleRegion().IsInside(idx_end) || !m_TargetRegions->GetLargestPossibleRegion().IsInside(idx_begin) || (m_TargetRegions->GetPixel(idx_begin)>0 && m_TargetRegions->GetPixel(idx_end)==m_TargetRegions->GetPixel(idx_begin)) ) { if (!m_StopTracking) { if (!m_UseOutputProbabilityMap) m_Tractogram.push_back(fib); else FiberToProbmap(&fib); m_CurrentTracts++; } if (m_MaxNumTracts > 0 && m_CurrentTracts>=static_cast(m_MaxNumTracts)) { if (!m_StopTracking) { std::cout << " \r"; MITK_INFO << "Reconstructed maximum number of tracts (" << m_CurrentTracts << "). Stopping tractography."; } m_StopTracking = true; } } } } } this->AfterTracking(); } void StreamlineTrackingFilter::CheckFiberForGmEnding(FiberType* fib) { if (m_TissueImage.IsNull()) return; // first check if the current fibe rendpoint is located inside of the white matter // if not, remove last fiber point and repeat bool in_wm = false; while (!in_wm && fib->size()>2) { ItkUcharImgType::IndexType idx; m_TissueImage->TransformPhysicalPointToIndex(fib->back(), idx); if (m_TissueImage->GetPixel(idx)==m_WmLabel) in_wm = true; else fib->pop_back(); } if (fib->size()<3 || !in_wm) { fib->clear(); return; } // get fiber direction at end point vnl_vector_fixed< float, 3 > d1; d1[0] = fib->back()[0] - fib->at(fib->size()-2)[0]; d1[1] = fib->back()[1] - fib->at(fib->size()-2)[1]; d1[2] = fib->back()[2] - fib->at(fib->size()-2)[2]; d1.normalize(); // find closest gray matter voxel ItkUcharImgType::IndexType s_idx; m_TissueImage->TransformPhysicalPointToIndex(fib->back(), s_idx); itk::Point gm_endp; float max = -1; for (int x : {-1,0,1}) for (int y : {-1,0,1}) for (int z : {-1,0,1}) { if (x==y && y==z) continue; ItkUcharImgType::IndexType e_idx; e_idx[0] = s_idx[0] + x; e_idx[1] = s_idx[1] + y; e_idx[2] = s_idx[2] + z; if ( !m_TissueImage->GetLargestPossibleRegion().IsInside(e_idx) || m_TissueImage->GetPixel(e_idx)!=m_GmLabel ) continue; itk::ContinuousIndex end; m_TissueImage->TransformIndexToPhysicalPoint(e_idx, end); vnl_vector_fixed< float, 3 > d2; d2[0] = end[0] - fib->back()[0]; d2[1] = end[1] - fib->back()[1]; d2[2] = end[2] - fib->back()[2]; d2.normalize(); float a = dot_product(d1,d2); if (a>max) { max = a; gm_endp = end; } } if (max>=0) fib->push_back(gm_endp); else // no gray matter enpoint found -> delete fiber fib->clear(); } void StreamlineTrackingFilter::FiberToProbmap(FiberType* fib) { ItkDoubleImgType::IndexType last_idx; last_idx.Fill(0); for (auto p : *fib) { ItkDoubleImgType::IndexType idx; m_OutputProbabilityMap->TransformPhysicalPointToIndex(p, idx); if (idx != last_idx) { if (m_OutputProbabilityMap->GetLargestPossibleRegion().IsInside(idx)) m_OutputProbabilityMap->SetPixel(idx, m_OutputProbabilityMap->GetPixel(idx)+1); last_idx = idx; } } } void StreamlineTrackingFilter::BuildFibers(bool check) { if (m_BuildFibersReady::New(); vtkSmartPointer vNewLines = vtkSmartPointer::New(); vtkSmartPointer vNewPoints = vtkSmartPointer::New(); for (unsigned int i=0; i container = vtkSmartPointer::New(); FiberType fib = m_Tractogram.at(i); for (FiberType::iterator it = fib.begin(); it!=fib.end(); ++it) { vtkIdType id = vNewPoints->InsertNextPoint((*it).GetDataPointer()); container->GetPointIds()->InsertNextId(id); } vNewLines->InsertNextCell(container); } if (check) for (int i=0; iSetPoints(vNewPoints); m_FiberPolyData->SetLines(vNewLines); m_BuildFibersFinished = true; } void StreamlineTrackingFilter::AfterTracking() { if (m_Verbose) std::cout << " \r"; if (!m_UseOutputProbabilityMap) { MITK_INFO << "Reconstructed " << m_Tractogram.size() << " fibers."; MITK_INFO << "Generating polydata "; BuildFibers(false); } else { itk::RescaleIntensityImageFilter< ItkDoubleImgType, ItkDoubleImgType >::Pointer filter = itk::RescaleIntensityImageFilter< ItkDoubleImgType, ItkDoubleImgType >::New(); filter->SetInput(m_OutputProbabilityMap); filter->SetOutputMaximum(1.0); filter->SetOutputMinimum(0.0); filter->Update(); m_OutputProbabilityMap = filter->GetOutput(); } MITK_INFO << "done"; m_EndTime = std::chrono::system_clock::now(); std::chrono::hours hh = std::chrono::duration_cast(m_EndTime - m_StartTime); std::chrono::minutes mm = std::chrono::duration_cast(m_EndTime - m_StartTime); std::chrono::seconds ss = std::chrono::duration_cast(m_EndTime - m_StartTime); mm %= 60; ss %= 60; MITK_INFO << "Tracking took " << hh.count() << "h, " << mm.count() << "m and " << ss.count() << "s"; m_SeedPoints.clear(); } 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()) { algo_code_value = "sup181_ee04"; algo_code_meaning = "FACT"; } else if (m_TrackingHandler->GetMode()==mitk::TrackingDataHandler::DETERMINISTIC) { algo_code_value = "sup181_ee01"; algo_code_meaning = "Deterministic"; } else if (m_TrackingHandler->GetMode()==mitk::TrackingDataHandler::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/itkTractClusteringFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractClusteringFilter.cpp index ebf2399ada..41e656680e 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractClusteringFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractClusteringFilter.cpp @@ -1,594 +1,486 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "itkTractClusteringFilter.h" #define _USE_MATH_DEFINES #include #include +#include namespace itk{ TractClusteringFilter::TractClusteringFilter() : m_NumPoints(12) , m_InCentroids(nullptr) , m_MinClusterSize(1) , m_MaxClusters(0) - , m_Metric(Metric::MDF) - , m_ScalarMap(nullptr) - , m_Scale(0) , m_MergeDuplicateThreshold(-1) , m_DoResampling(true) , m_FilterMask(nullptr) , m_OverlapThreshold(0.0) { } TractClusteringFilter::~TractClusteringFilter() { + for (auto m : m_Metrics) + delete m; +} +void TractClusteringFilter::SetMetrics(const std::vector &Metrics) +{ + m_Metrics = Metrics; } std::vector TractClusteringFilter::GetOutClusters() const { return m_OutClusters; } std::vector TractClusteringFilter::GetOutCentroids() const { return m_OutCentroids; } -void TractClusteringFilter::SetMetric(const Metric &Metric) -{ - m_Metric = Metric; -} - std::vector TractClusteringFilter::GetOutTractograms() const { return m_OutTractograms; } void TractClusteringFilter::SetDistances(const std::vector &Distances) { m_Distances = Distances; } float TractClusteringFilter::CalcOverlap(vnl_matrix& t) { float overlap = 0; if (m_FilterMask.IsNotNull()) { for (unsigned int i=0; i p = t.get_column(i); itk::Point point; point[0] = p[0]; point[1] = p[1]; point[2] = p[2]; itk::Index<3> idx; m_FilterMask->TransformPhysicalPointToIndex(point, idx); if (m_FilterMask->GetLargestPossibleRegion().IsInside(idx) && m_FilterMask->GetPixel(idx)>0) overlap += 1; } overlap /= m_NumPoints; } else return 1.0; return overlap; } -float TractClusteringFilter::CalcMDF(vnl_matrix& s, vnl_matrix& t, bool& flipped) -{ - float d_direct = 0; - float d_flipped = 0; - - for (unsigned int i=0; id_flipped) - { - flipped = true; - return d_flipped/m_NumPoints; - } - flipped = false; - return d_direct/m_NumPoints; -} - -float TractClusteringFilter::CalcMDF_STD(vnl_matrix& s, vnl_matrix& t, bool& flipped) -{ - float d_direct = 0; - float d_flipped = 0; - - vnl_vector dists_d; dists_d.set_size(m_NumPoints); - vnl_vector dists_f; dists_f.set_size(m_NumPoints); - - for (unsigned int i=0; id_flipped) - { - float d = d_flipped/m_NumPoints; - dists_f -= d; - d += dists_f.magnitude(); - - flipped = true; - return d/2; - } - - float d = d_direct/m_NumPoints; - dists_d -= d; - d += dists_d.magnitude(); - - flipped = false; - return d/2; -} - -float TractClusteringFilter::CalcMAX_MDF(vnl_matrix& s, vnl_matrix& t, bool& flipped) -{ - float d_direct = 0; - float d_flipped = 0; - - vnl_vector dists_d; dists_d.set_size(m_NumPoints); - vnl_vector dists_f; dists_f.set_size(m_NumPoints); - - for (unsigned int i=0; id_flipped) - { - float d = dists_f.max_value(); - flipped = true; - return d; - } - - float d = dists_d.max_value(); - flipped = false; - return d; -} - std::vector > TractClusteringFilter::ResampleFibers(mitk::FiberBundle::Pointer tractogram) { mitk::FiberBundle::Pointer temp_fib = tractogram->GetDeepCopy(); if (m_DoResampling) temp_fib->ResampleToNumPoints(m_NumPoints); std::vector< vnl_matrix > out_fib; for (int i=0; iGetFiberPolyData()->GetNumberOfCells(); i++) { vtkCell* cell = temp_fib->GetFiberPolyData()->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vnl_matrix streamline; - if (m_ScalarMap.IsNull()) - streamline.set_size(3, m_NumPoints); - else - streamline.set_size(4, m_NumPoints); + streamline.set_size(3, m_NumPoints); streamline.fill(0.0); for (int j=0; jGetPoint(j, cand); - if (m_ScalarMap.IsNull()) - { - vnl_vector_fixed< float, 3 > candV; - candV[0]=cand[0]; candV[1]=cand[1]; candV[2]=cand[2]; - streamline.set_column(j, candV); - } - else - { - vnl_vector_fixed< float, 4 > candV; - candV[0]=cand[0]; candV[1]=cand[1]; candV[2]=cand[2]; candV[3]=0; - itk::Point wp; wp[0]=cand[0]; wp[1]=cand[1]; wp[2]=cand[2]; - - itk::Index<3> idx; - m_ScalarMap->TransformPhysicalPointToIndex(wp, idx); - if (m_ScalarMap->GetLargestPossibleRegion().IsInside(idx)) - candV[3]=m_ScalarMap->GetPixel(idx)*m_Scale; - - streamline.set_column(j, candV); - } + vnl_vector_fixed< float, 3 > candV; + candV[0]=cand[0]; candV[1]=cand[1]; candV[2]=cand[2]; + streamline.set_column(j, candV); } out_fib.push_back(streamline); } return out_fib; } std::vector< TractClusteringFilter::Cluster > TractClusteringFilter::ClusterStep(std::vector< long > f_indices, std::vector distances) { float dist_thres = distances.back(); distances.pop_back(); std::vector< Cluster > C; int N = f_indices.size(); Cluster c1; c1.I.push_back(f_indices.at(0)); c1.h = T[f_indices.at(0)]; c1.n = 1; C.push_back(c1); if (f_indices.size()==1) return C; for (int i=1; i t = T.at(f_indices.at(i)); int min_cluster_index = -1; float min_cluster_distance = 99999; bool flip = false; for (unsigned int k=0; k v = C.at(k).h / C.at(k).n; bool f = false; float d = 0; - if (m_Metric==Metric::MDF) - d = CalcMDF(t, v, f); - else if (m_Metric==Metric::MDF_STD) - d = CalcMDF_STD(t, v, f); - else if (m_Metric==Metric::MAX_MDF) - d = CalcMAX_MDF(t, v, f); + for (auto m : m_Metrics) + d += m->CalculateDistance(t, v, f); + d /= m_Metrics.size(); if (d=0 && min_cluster_distance outC; #pragma omp parallel for for (int c=0; c<(int)C.size(); c++) { std::vector< Cluster > tempC = ClusterStep(C.at(c).I, distances); AppendCluster(outC, tempC); } return outC; } else return C; } void TractClusteringFilter::AppendCluster(std::vector< Cluster >& a, std::vector< Cluster >&b) { for (auto c : b) a.push_back(c); } void TractClusteringFilter::MergeDuplicateClusters(std::vector< TractClusteringFilter::Cluster >& clusters) { if (m_MergeDuplicateThreshold<0) - m_MergeDuplicateThreshold = m_Distances.at(0); + m_MergeDuplicateThreshold = m_Distances.at(0)/2; bool found = true; MITK_INFO << "Merging duplicate clusters with distance threshold " << m_MergeDuplicateThreshold; int start = 0; while (found && m_MergeDuplicateThreshold>mitk::eps) { std::cout << " \r"; std::cout << "Number of clusters: " << clusters.size() << '\r'; cout.flush(); found = false; for (int k1=start; k1<(int)clusters.size(); ++k1) { Cluster c1 = clusters.at(k1); vnl_matrix t = c1.h / c1.n; std::vector< int > merge_indices; std::vector< bool > flip_indices; #pragma omp parallel for for (int k2=0; k2<(int)clusters.size(); ++k2) { if (k1!=k2) { Cluster c2 = clusters.at(k2); vnl_matrix v = c2.h / c2.n; bool f = false; - float d = 0; - // if (m_Metric==Metric::MDF) - // d = CalcMDF(t, v, f); - // else if (m_Metric==Metric::MDF_STD) - // d = CalcMDF_STD(t, v, f); - // else if (m_Metric==Metric::MAX_MDF) - d = CalcMAX_MDF(t, v, f); +// d = m_Metric->CalculateDistance(t, v, f); // alwayse use MDF? + float d = 0; + for (auto m : m_Metrics) + d += m->CalculateDistance(t, v, f); + d /= m_Metrics.size(); #pragma omp critical if (d TractClusteringFilter::AddToKnownClusters(std::vector< long > f_indices, std::vector >& centroids) { float dist_thres = m_Distances.at(0); int N = f_indices.size(); std::vector< Cluster > C; vnl_matrix zero_h; zero_h.set_size(T.at(0).rows(), T.at(0).cols()); zero_h.fill(0.0); Cluster no_fit; no_fit.h = zero_h; for (unsigned int i=0; i t = T.at(f_indices.at(i)); int min_cluster_index = -1; float min_cluster_distance = 99999; bool flip = false; if (CalcOverlap(t)>=m_OverlapThreshold) { int c_idx = 0; for (vnl_matrix centroid : centroids) { bool f = false; float d = 0; - - if (m_Metric==Metric::MDF) - d = CalcMDF(t, centroid, f); - else if (m_Metric==Metric::MDF_STD) - d = CalcMDF_STD(t, centroid, f); - else if (m_Metric==Metric::MAX_MDF) - d = CalcMAX_MDF(t, centroid, f); + for (auto m : m_Metrics) + d += m->CalculateDistance(t, centroid, f); + d /= m_Metrics.size(); if (d=0 && min_cluster_distance f_indices; for (unsigned int i=0; i clusters; if (m_InCentroids.IsNull()) { MITK_INFO << "Clustering fibers"; clusters = ClusterStep(f_indices, m_Distances); MITK_INFO << "Number of clusters: " << clusters.size(); MergeDuplicateClusters(clusters); std::sort(clusters.begin(),clusters.end()); } else { std::vector > centroids = ResampleFibers(m_InCentroids); if (centroids.empty()) { MITK_INFO << "No fibers in centroid tractogram!"; return; } MITK_INFO << "Clustering with input centroids"; clusters = AddToKnownClusters(f_indices, centroids); no_match = clusters.back(); clusters.pop_back(); MITK_INFO << "Number of clusters: " << clusters.size(); MergeDuplicateClusters(clusters); } MITK_INFO << "Clustering finished"; int max = clusters.size()-1; if (m_MaxClusters>0 && clusters.size()-1>m_MaxClusters) max = m_MaxClusters; for (int i=clusters.size()-1; i>=0; --i) { Cluster c = clusters.at(i); if ( c.n>=(int)m_MinClusterSize && !(m_MaxClusters>0 && clusters.size()-i>=m_MaxClusters) ) { m_OutClusters.push_back(c); vtkSmartPointer weights = vtkSmartPointer::New(); vtkSmartPointer pTmp = m_Tractogram->GeneratePolyDataByIds(c.I, weights); mitk::FiberBundle::Pointer fib = mitk::FiberBundle::New(pTmp); if (max>0) fib->SetFiberWeights((float)i/max); m_OutTractograms.push_back(fib); // create centroid vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); vtkSmartPointer polyData = vtkSmartPointer::New(); vtkSmartPointer container = vtkSmartPointer::New(); vnl_matrix centroid_points = c.h / c.n; for (unsigned int j=0; jInsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); polyData->SetPoints(vtkNewPoints); polyData->SetLines(vtkNewCells); mitk::FiberBundle::Pointer centroid = mitk::FiberBundle::New(polyData); centroid->SetFiberColors(255, 255, 255); m_OutCentroids.push_back(centroid); } } MITK_INFO << "Final number of clusters: " << m_OutTractograms.size(); int w = 0; for (auto fib : m_OutTractograms) { if (m_OutTractograms.size()>1) { fib->SetFiberWeights((float)w/(m_OutTractograms.size()-1)); m_OutCentroids.at(w)->SetFiberWeights((float)w/(m_OutTractograms.size()-1)); } else { fib->SetFiberWeights(1); m_OutCentroids.at(w)->SetFiberWeights(1); } fib->ColorFibersByFiberWeights(false, false); ++w; } if (no_match.n>0) { vtkSmartPointer weights = vtkSmartPointer::New(); vtkSmartPointer pTmp = m_Tractogram->GeneratePolyDataByIds(no_match.I, weights); mitk::FiberBundle::Pointer fib = mitk::FiberBundle::New(pTmp); fib->SetFiberColors(0, 0, 0); m_OutTractograms.push_back(fib); } } } diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractClusteringFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractClusteringFilter.h index 4de44967dd..dd4bb041ed 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractClusteringFilter.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractClusteringFilter.h @@ -1,151 +1,136 @@ /*=================================================================== 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 itkTractClusteringFilter_h #define itkTractClusteringFilter_h // MITK #include #include #include +#include // ITK #include // VTK #include #include #include #include #include namespace itk{ /** * \brief */ class TractClusteringFilter : public ProcessObject { public: struct Cluster { Cluster() : n(0), f_id(-1) {} vnl_matrix h; std::vector< long > I; int n; int f_id; bool operator <(Cluster const& b) const { return this->n < b.n; } }; - enum Metric - { - MDF, - MDF_STD, - MAX_MDF - }; - typedef TractClusteringFilter Self; typedef ProcessObject Superclass; typedef SmartPointer< Self > Pointer; typedef SmartPointer< const Self > ConstPointer; typedef itk::Image< float, 3 > FloatImageType; typedef itk::Image< unsigned char, 3 > UcharImageType; itkFactorylessNewMacro(Self) itkCloneMacro(Self) itkTypeMacro( TractClusteringFilter, ProcessObject ) itkSetMacro(NumPoints, unsigned int) ///< Fibers are resampled to the specified number of points. If scalar maps are used, a larger number of points is recommended. itkGetMacro(NumPoints, unsigned int) ///< Fibers are resampled to the specified number of points. If scalar maps are used, a larger number of points is recommended. itkSetMacro(MinClusterSize, unsigned int) ///< Clusters with too few fibers are discarded itkGetMacro(MinClusterSize, unsigned int) ///< Clusters with too few fibers are discarded itkSetMacro(MaxClusters, unsigned int) ///< Only the N largest clusters are kept itkGetMacro(MaxClusters, unsigned int) ///< Only the N largest clusters are kept - itkSetMacro(Scale, float) ///< Scaling factor for scalar map. By default, the map is scaled by a factor corresponding to the cluster distance - itkGetMacro(Scale, float) ///< Scaling factor for scalar map. By default, the map is scaled by a factor corresponding to the cluster distance itkSetMacro(MergeDuplicateThreshold, float) ///< Clusters with centroids very close to each other are merged. Set to 0 to avoid merging and to -1 to use the original cluster size. itkGetMacro(MergeDuplicateThreshold, float) ///< Clusters with centroids very close to each other are merged. Set to 0 to avoid merging and to -1 to use the original cluster size. itkSetMacro(DoResampling, bool) ///< Resample fibers to equal number of points. This is mandatory, but can be performed outside of the filter if desired. itkGetMacro(DoResampling, bool) ///< Resample fibers to equal number of points. This is mandatory, but can be performed outside of the filter if desired. itkSetMacro(OverlapThreshold, float) ///< Overlap threshold used in conjunction with the filter mask when clustering around known centroids. itkGetMacro(OverlapThreshold, float) ///< Overlap threshold used in conjunction with the filter mask when clustering around known centroids. itkSetMacro(Tractogram, mitk::FiberBundle::Pointer) ///< The streamlines to be clustered - itkSetMacro(ScalarMap, FloatImageType::Pointer) ///< Scalar values along the streamlines are included into the distance calculation. If this feature is used, it is recommended to use a larger number of fiber sampling points. itkSetMacro(InCentroids, mitk::FiberBundle::Pointer) ///< If a tractogram containing known tract centroids is set, the input fibers are assigned to the closest centroid. If no centroid is found within the specified smallest clustering distance, the fiber is assigned to the no-fit cluster. itkSetMacro(FilterMask, UcharImageType::Pointer) ///< If fibers are clustered around the nearest input centroids (see SetInCentroids), the complete input tractogram can additionally be pre-filtered with this binary mask and a given overlap threshold (see SetOverlapThreshold). virtual void Update() override{ this->GenerateData(); } void SetDistances(const std::vector &Distances); ///< Set clustering distances that are traversed recoursively. The distances have to be sorted in an ascending manner. The actual cluster size is determined by the smallest entry in the distance-list (index 0). std::vector GetOutTractograms() const; - - void SetMetric(const Metric &Metric); - std::vector GetOutCentroids() const; - std::vector GetOutClusters() const; + void SetMetrics(const std::vector &Metrics); + protected: void GenerateData() override; std::vector< vnl_matrix > ResampleFibers(FiberBundle::Pointer tractogram); - float CalcMDF(vnl_matrix& s, vnl_matrix& t, bool &flipped); - float CalcMDF_STD(vnl_matrix& s, vnl_matrix& t, bool &flipped); - float CalcMAX_MDF(vnl_matrix& s, vnl_matrix& t, bool &flipped); float CalcOverlap(vnl_matrix& t); std::vector< Cluster > ClusterStep(std::vector< long > f_indices, std::vector< float > distances); + void MergeDuplicateClusters(std::vector< TractClusteringFilter::Cluster >& clusters); std::vector< Cluster > AddToKnownClusters(std::vector< long > f_indices, std::vector > ¢roids); void AppendCluster(std::vector< Cluster >& a, std::vector< Cluster >&b); TractClusteringFilter(); virtual ~TractClusteringFilter(); unsigned int m_NumPoints; std::vector< float > m_Distances; mitk::FiberBundle::Pointer m_Tractogram; mitk::FiberBundle::Pointer m_InCentroids; std::vector< mitk::FiberBundle::Pointer > m_OutTractograms; std::vector< mitk::FiberBundle::Pointer > m_OutCentroids; std::vector > T; unsigned int m_MinClusterSize; unsigned int m_MaxClusters; - Metric m_Metric; - FloatImageType::Pointer m_ScalarMap; - float m_Scale; float m_MergeDuplicateThreshold; std::vector< Cluster > m_OutClusters; bool m_DoResampling; UcharImageType::Pointer m_FilterMask; float m_OverlapThreshold; + std::vector< mitk::ClusteringMetric* > m_Metrics; }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkTractClusteringFilter.cpp" #endif #endif diff --git a/Modules/DiffusionImaging/FiberTracking/CMakeLists.txt b/Modules/DiffusionImaging/FiberTracking/CMakeLists.txt index f4c88f461a..3d2ef66df8 100644 --- a/Modules/DiffusionImaging/FiberTracking/CMakeLists.txt +++ b/Modules/DiffusionImaging/FiberTracking/CMakeLists.txt @@ -1,54 +1,54 @@ set(_module_deps MitkDiffusionCore MitkGraphAlgorithms MitkCLVigraRandomForest) mitk_check_module_dependencies( MODULES ${_module_deps} MISSING_DEPENDENCIES_VAR _missing_deps ) if(NOT _missing_deps) set(lut_url http://mitk.org/download/data/FibertrackingLUT.tar.gz) set(lut_tarball ${CMAKE_CURRENT_BINARY_DIR}/FibertrackingLUT.tar.gz) #message("Downloading FiberTracking LUT ${lut_url}...") file(DOWNLOAD ${lut_url} ${lut_tarball} EXPECTED_MD5 38ecb6d4a826c9ebb0f4965eb9aeee44 TIMEOUT 60 STATUS status SHOW_PROGRESS ) list(GET status 0 status_code) list(GET status 1 status_msg) if(NOT status_code EQUAL 0) message(SEND_ERROR "${status_msg} (error code ${status_code})") endif() file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/Resources) #message("Unpacking FiberTracking LUT tarball...") execute_process(COMMAND ${CMAKE_COMMAND} -E tar xzf ../FibertrackingLUT.tar.gz WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/Resources RESULT_VARIABLE result ERROR_VARIABLE err_msg) if(result) message(SEND_ERROR "Unpacking FibertrackingLUT.tar.gz failed: ${err_msg}") endif() endif() MITK_CREATE_MODULE( SUBPROJECTS MITK-DTI - INCLUDE_DIRS Fiberfox Fiberfox/SignalModels Fiberfox/Sequences Algorithms Algorithms/TrackingHandlers Algorithms/GibbsTracking Algorithms/StochasticTracking IODataStructures IODataStructures/FiberBundle IODataStructures/PlanarFigureComposite Rendering ${CMAKE_CURRENT_BINARY_DIR} + INCLUDE_DIRS Fiberfox Fiberfox/SignalModels Fiberfox/Sequences Algorithms Algorithms/TrackingHandlers Algorithms/ClusteringMetrics Algorithms/GibbsTracking Algorithms/StochasticTracking IODataStructures IODataStructures/FiberBundle IODataStructures/PlanarFigureComposite Rendering ${CMAKE_CURRENT_BINARY_DIR} DEPENDS ${_module_deps} PACKAGE_DEPENDS PUBLIC ITK|ITKFFT ITK|ITKDiffusionTensorImage Vigra HDF5 ) if(MODULE_IS_ENABLED) add_subdirectory(cmdapps/Fiberfox) add_subdirectory(cmdapps/FiberProcessing) add_subdirectory(cmdapps/Tractography) add_subdirectory(cmdapps/TractographyEvaluation) add_subdirectory(cmdapps/Misc) add_subdirectory(Testing) endif() diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberClustering.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberClustering.cpp index fe6b191d6f..30e2b5e2cc 100644 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberClustering.cpp +++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberClustering.cpp @@ -1,190 +1,225 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include +#include +#include +#include +#include +#include mitk::FiberBundle::Pointer LoadFib(std::string filename) { std::vector fibInfile = mitk::IOUtil::Load(filename); if( fibInfile.empty() ) std::cout << "File " << filename << " could not be read!"; mitk::BaseData::Pointer baseData = fibInfile.at(0); return dynamic_cast(baseData.GetPointer()); } /*! \brief Spatially cluster fibers */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Fiber Clustering"); parser.setCategory("Fiber Processing"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkCommandLineParser::InputFile, "Input:", "input fiber bundle (.fib, .trk, .tck)", us::Any(), false); parser.addArgument("", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output root", us::Any(), false); parser.addArgument("cluster_size", "", mitkCommandLineParser::Int, "Cluster size:", "", 10); parser.addArgument("fiber_points", "", mitkCommandLineParser::Int, "Fiber points:", "", 12); parser.addArgument("min_fibers", "", mitkCommandLineParser::Int, "Min. fibers per cluster:", "", 1); parser.addArgument("max_clusters", "", mitkCommandLineParser::Int, "Max. clusters:", ""); parser.addArgument("merge_clusters", "", mitkCommandLineParser::Float, "Merge clusters:", "", -1.0); parser.addArgument("output_centroids", "", mitkCommandLineParser::Bool, "Output centroids:", ""); - parser.addArgument("metric", "", mitkCommandLineParser::String, "Metric:", ""); + parser.addArgument("metrics", "", mitkCommandLineParser::StringList, "Metrics:", "EU_MEAN, EU_STD, EU_MAX, ANAT, MAP"); parser.addArgument("input_centroids", "", mitkCommandLineParser::String, "Input centroids:", ""); parser.addArgument("scalar_map", "", mitkCommandLineParser::String, "Scalar map:", ""); - parser.addArgument("map_scale", "", mitkCommandLineParser::Float, "Map scale:", "", 0.0); + parser.addArgument("parcellation", "", mitkCommandLineParser::String, "Parcellation:", ""); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string inFileName = us::any_cast(parsedArgs["i"]); std::string out_root = us::any_cast(parsedArgs["o"]); int cluster_size = 10; if (parsedArgs.count("cluster_size")) cluster_size = us::any_cast(parsedArgs["cluster_size"]); int fiber_points = 12; if (parsedArgs.count("fiber_points")) fiber_points = us::any_cast(parsedArgs["fiber_points"]); int min_fibers = 1; if (parsedArgs.count("min_fibers")) min_fibers = us::any_cast(parsedArgs["min_fibers"]); int max_clusters = 0; if (parsedArgs.count("max_clusters")) max_clusters = us::any_cast(parsedArgs["max_clusters"]); - float map_scale = 0.0; - if (parsedArgs.count("map_scale")) - map_scale = us::any_cast(parsedArgs["map_scale"]); - float merge_clusters = -1.0; if (parsedArgs.count("merge_clusters")) merge_clusters = us::any_cast(parsedArgs["merge_clusters"]); bool output_centroids = false; if (parsedArgs.count("output_centroids")) output_centroids = us::any_cast(parsedArgs["output_centroids"]); - std::string metric = "MDF_STD"; - if (parsedArgs.count("metric")) - metric = us::any_cast(parsedArgs["metric"]); + std::vector< std::string > metric_strings = {"EUCL"}; + if (parsedArgs.count("metrics")) + metric_strings = us::any_cast(parsedArgs["metrics"]); std::string input_centroids = ""; if (parsedArgs.count("input_centroids")) input_centroids = us::any_cast(parsedArgs["input_centroids"]); std::string scalar_map = ""; if (parsedArgs.count("scalar_map")) scalar_map = us::any_cast(parsedArgs["scalar_map"]); + std::string parcellation = ""; + if (parsedArgs.count("parcellation")) + parcellation = us::any_cast(parsedArgs["parcellation"]); + try { typedef itk::Image< float, 3 > FloatImageType; + typedef itk::Image< short, 3 > ShortImageType; mitk::FiberBundle::Pointer fib = LoadFib(inFileName); float max_d = 0; int i=1; std::vector< float > distances; while (max_d < fib->GetGeometry()->GetDiagonalLength()/2) { distances.push_back(cluster_size*i); max_d = cluster_size*i; ++i; } itk::TractClusteringFilter::Pointer clusterer = itk::TractClusteringFilter::New(); clusterer->SetDistances(distances); clusterer->SetTractogram(fib); if (input_centroids!="") { mitk::FiberBundle::Pointer in_centroids = LoadFib(input_centroids); clusterer->SetInCentroids(in_centroids); } - if (metric=="MDF") - clusterer->SetMetric(itk::TractClusteringFilter::Metric::MDF); - else if (metric=="MDF_STD") - clusterer->SetMetric(itk::TractClusteringFilter::Metric::MDF_STD); - else if (metric=="MAX_MDF") - clusterer->SetMetric(itk::TractClusteringFilter::Metric::MAX_MDF); + std::vector< mitk::ClusteringMetric* > metrics; - if (scalar_map!="") + for (auto m : metric_strings) { - mitk::Image::Pointer mitk_map = dynamic_cast(mitk::IOUtil::Load(scalar_map)[0].GetPointer()); - if (mitk_map->GetDimension()==3) + MITK_INFO << "Metric: " << m; + if (m=="EU_MEAN") + metrics.push_back({new mitk::ClusteringMetricEuclideanMean()}); + else if (m=="EU_STD") + metrics.push_back({new mitk::ClusteringMetricEuclideanStd()}); + else if (m=="EU_MAX") + metrics.push_back({new mitk::ClusteringMetricEuclideanMax()}); + else if (m=="MAP" && scalar_map!="") + { + mitk::Image::Pointer mitk_map = dynamic_cast(mitk::IOUtil::Load(scalar_map)[0].GetPointer()); + if (mitk_map->GetDimension()==3) + { + FloatImageType::Pointer itk_map = FloatImageType::New(); + mitk::CastToItkImage(mitk_map, itk_map); + + mitk::ClusteringMetricScalarMap* metric = new mitk::ClusteringMetricScalarMap(); + metric->SetImages({itk_map}); + metric->SetScale(distances.at(0)); + metrics.push_back(metric); + } + } + else if (m=="ANAT" && parcellation!="") { - FloatImageType::Pointer itk_map = FloatImageType::New(); - mitk::CastToItkImage(mitk_map, itk_map); - clusterer->SetScalarMap(itk_map); + mitk::Image::Pointer mitk_map = dynamic_cast(mitk::IOUtil::Load(parcellation)[0].GetPointer()); + if (mitk_map->GetDimension()==3) + { + ShortImageType::Pointer itk_map = ShortImageType::New(); + mitk::CastToItkImage(mitk_map, itk_map); + + mitk::ClusteringMetricAnatomic* metric = new mitk::ClusteringMetricAnatomic(); + metric->SetParcellations({itk_map}); + metrics.push_back(metric); + } } } + + if (metrics.empty()) + { + MITK_INFO << "No metric selected!"; + return EXIT_FAILURE; + } + + clusterer->SetMetrics(metrics); clusterer->SetMergeDuplicateThreshold(merge_clusters); - clusterer->SetScale(map_scale); clusterer->SetNumPoints(fiber_points); clusterer->SetMaxClusters(max_clusters); clusterer->SetMinClusterSize(min_fibers); clusterer->Update(); std::vector tracts = clusterer->GetOutTractograms(); std::vector centroids = clusterer->GetOutCentroids(); MITK_INFO << "Saving clusters"; std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect unsigned int c = 0; for (auto f : tracts) { mitk::IOUtil::Save(f, out_root + "Cluster_" + boost::lexical_cast(c) + ".fib"); if (output_centroids) mitk::IOUtil::Save(centroids.at(c), out_root + "Centroid_" + boost::lexical_cast(c) + ".fib"); ++c; } std::cout.rdbuf (old); // <-- restore } catch (itk::ExceptionObject e) { std::cout << e; return EXIT_FAILURE; } catch (std::exception e) { std::cout << e.what(); return EXIT_FAILURE; } catch (...) { std::cout << "ERROR!?!"; return EXIT_FAILURE; } return EXIT_SUCCESS; } diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/ExtractSimilarTracts.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/ExtractSimilarTracts.cpp index 74caed4caa..dbbe833487 100644 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/ExtractSimilarTracts.cpp +++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/ExtractSimilarTracts.cpp @@ -1,211 +1,216 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include +#include +#include +#include +#include +#include typedef itksys::SystemTools ist; typedef itk::Image ItkUcharImgType; mitk::FiberBundle::Pointer LoadFib(std::string filename) { std::vector fibInfile = mitk::IOUtil::Load(filename); if( fibInfile.empty() ) std::cout << "File " << filename << " could not be read!"; mitk::BaseData::Pointer baseData = fibInfile.at(0); return dynamic_cast(baseData.GetPointer()); } ItkUcharImgType::Pointer LoadItkMaskImage(const std::string& filename) { mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::Load(filename)[0].GetPointer()); ItkUcharImgType::Pointer itkMask = ItkUcharImgType::New(); mitk::CastToItkImage(img, itkMask); return itkMask; } /*! \brief Spatially cluster fibers */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Extract Similar Tracts"); parser.setCategory("Fiber Tracking Evaluation"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkCommandLineParser::InputFile, "Input:", "input fiber bundle (.fib, .trk, .tck)", us::Any(), false); parser.addArgument("ref_tracts", "", mitkCommandLineParser::StringList, "Ref. Tracts:", "reference tracts (.fib, .trk, .tck)", us::Any(), false); parser.addArgument("ref_masks", "", mitkCommandLineParser::StringList, "Ref. Masks:", "reference bundle masks", us::Any()); parser.addArgument("", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output root", us::Any(), false); parser.addArgument("distance", "", mitkCommandLineParser::Int, "Distance:", "", 10); parser.addArgument("metric", "", mitkCommandLineParser::String, "Metric:", ""); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string in_fib = us::any_cast(parsedArgs["i"]); std::string out_root = us::any_cast(parsedArgs["o"]); mitkCommandLineParser::StringContainerType ref_bundle_files = us::any_cast(parsedArgs["ref_tracts"]); mitkCommandLineParser::StringContainerType ref_mask_files; if (parsedArgs.count("ref_masks")) ref_mask_files = us::any_cast(parsedArgs["ref_masks"]); if (ref_mask_files.size()>0 && ref_mask_files.size()!=ref_bundle_files.size()) { MITK_INFO << "If reference masks are used, there has to be one mask per reference tract."; return EXIT_FAILURE; } int distance = 10; if (parsedArgs.count("distance")) distance = us::any_cast(parsedArgs["distance"]); - std::string metric = "MDF_STD"; + std::string metric = "EU_MEAN"; if (parsedArgs.count("metric")) metric = us::any_cast(parsedArgs["metric"]); try { mitk::FiberBundle::Pointer fib = LoadFib(in_fib); fib->ResampleToNumPoints(12); std::vector< mitk::FiberBundle::Pointer > ref_fibs; std::vector< ItkUcharImgType::Pointer > ref_masks; for (std::size_t i=0; i distances; distances.push_back(distance); mitk::FiberBundle::Pointer anchor_tractogram = mitk::FiberBundle::New(nullptr); unsigned int c = 0; for (auto ref_fib : ref_fibs) { MITK_INFO << "Extracting " << ist::GetFilenameName(ref_bundle_files.at(c)); std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect try { itk::TractClusteringFilter::Pointer segmenter = itk::TractClusteringFilter::New(); // calculate centroids from reference bundle { itk::TractClusteringFilter::Pointer clusterer = itk::TractClusteringFilter::New(); clusterer->SetDistances({10,20,30}); clusterer->SetTractogram(ref_fib); - clusterer->SetMetric(itk::TractClusteringFilter::Metric::MDF_STD); + clusterer->SetMetrics({new mitk::ClusteringMetricEuclideanStd()}); clusterer->Update(); std::vector tracts = clusterer->GetOutCentroids(); ref_fib = mitk::FiberBundle::New(nullptr); ref_fib = ref_fib->AddBundles(tracts); mitk::IOUtil::Save(ref_fib, out_root + "centroids_" + ist::GetFilenameName(ref_bundle_files.at(c))); segmenter->SetInCentroids(ref_fib); } // segment tract segmenter->SetFilterMask(ref_masks.at(c)); segmenter->SetOverlapThreshold(0.8); segmenter->SetDistances(distances); segmenter->SetTractogram(fib); segmenter->SetDoResampling(false); - if (metric=="MDF") - segmenter->SetMetric(itk::TractClusteringFilter::Metric::MDF); - else if (metric=="MDF_STD") - segmenter->SetMetric(itk::TractClusteringFilter::Metric::MDF_STD); - else if (metric=="MAX_MDF") - segmenter->SetMetric(itk::TractClusteringFilter::Metric::MAX_MDF); + if (metric=="EU_MEAN") + segmenter->SetMetrics({new mitk::ClusteringMetricEuclideanMean()}); + else if (metric=="EU_STD") + segmenter->SetMetrics({new mitk::ClusteringMetricEuclideanStd()}); + else if (metric=="EU_MAX") + segmenter->SetMetrics({new mitk::ClusteringMetricEuclideanMax()}); segmenter->Update(); std::vector clusters = segmenter->GetOutTractograms(); if (clusters.size()>0) { fib = clusters.back(); clusters.pop_back(); mitk::FiberBundle::Pointer result = mitk::FiberBundle::New(nullptr); result = result->AddBundles(clusters); anchor_tractogram = anchor_tractogram->AddBundle(result); mitk::IOUtil::Save(result, out_root + "anchor_" + ist::GetFilenameName(ref_bundle_files.at(c))); } } catch(itk::ExceptionObject& excpt) { MITK_INFO << "Exception while processing " << ist::GetFilenameName(ref_bundle_files.at(c)); MITK_INFO << excpt.GetDescription(); } catch(std::exception& excpt) { MITK_INFO << "Exception while processing " << ist::GetFilenameName(ref_bundle_files.at(c)); MITK_INFO << excpt.what(); } std::cout.rdbuf (old); // <-- restore if (fib->GetNumFibers()==0) break; ++c; } MITK_INFO << "Streamlines in anchor tractogram: " << anchor_tractogram->GetNumFibers(); mitk::IOUtil::Save(anchor_tractogram, out_root + "anchor_tractogram.trk"); MITK_INFO << "Streamlines remaining in candidate tractogram: " << fib->GetNumFibers(); mitk::IOUtil::Save(fib, out_root + "candidate_tractogram.trk"); } catch (itk::ExceptionObject e) { std::cout << e; return EXIT_FAILURE; } catch (std::exception e) { std::cout << e.what(); return EXIT_FAILURE; } catch (...) { std::cout << "ERROR!?!"; return EXIT_FAILURE; } return EXIT_SUCCESS; } diff --git a/Modules/DiffusionImaging/FiberTracking/files.cmake b/Modules/DiffusionImaging/FiberTracking/files.cmake index 5b0d93c0a0..6eba166f25 100644 --- a/Modules/DiffusionImaging/FiberTracking/files.cmake +++ b/Modules/DiffusionImaging/FiberTracking/files.cmake @@ -1,93 +1,101 @@ set(CPP_FILES mitkFiberTrackingModuleActivator.cpp ## IO datastructures IODataStructures/FiberBundle/mitkFiberBundle.cpp IODataStructures/FiberBundle/mitkTrackvis.cpp IODataStructures/PlanarFigureComposite/mitkPlanarFigureComposite.cpp IODataStructures/mitkTractographyForest.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 # 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 # 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 + # Fiberfox Fiberfox/itkFibersFromPlanarFiguresFilter.h Fiberfox/itkTractsToDWIImageFilter.h Fiberfox/itkKspaceImageFilter.h Fiberfox/itkDftImageFilter.h Fiberfox/itkFieldmapGeneratorFilter.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/mitkCartesianReadout.h ) set(RESOURCE_FILES # Binary directory resources FiberTrackingLUTBaryCoords.bin FiberTrackingLUTIndices.bin ) diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberClusteringView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberClusteringView.cpp index 1b04352336..1d7a4d2b05 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberClusteringView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberClusteringView.cpp @@ -1,175 +1,204 @@ /*=================================================================== 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 "QmitkFiberClusteringView.h" #include +#include +#include #include #include #include #include +#include +#include +#include +#include +#include +#include const std::string QmitkFiberClusteringView::VIEW_ID = "org.mitk.views.fiberclustering"; using namespace mitk; QmitkFiberClusteringView::QmitkFiberClusteringView() : QmitkAbstractView() , m_Controls( nullptr ) { } // Destructor QmitkFiberClusteringView::~QmitkFiberClusteringView() { } void QmitkFiberClusteringView::CreateQtPartControl( QWidget *parent ) { // build up qt view, unless already done if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkFiberClusteringViewControls; m_Controls->setupUi( parent ); connect( m_Controls->m_StartButton, SIGNAL(clicked()), this, SLOT(StartClustering()) ); connect( m_Controls->m_TractBox, SIGNAL(currentIndexChanged(int)), this, SLOT(FiberSelectionChanged()) ); mitk::TNodePredicateDataType::Pointer isFib = mitk::TNodePredicateDataType::New(); m_Controls->m_TractBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_TractBox->SetPredicate(isFib); m_Controls->m_InCentroidsBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_InCentroidsBox->SetPredicate(isFib); m_Controls->m_InCentroidsBox->SetZeroEntryText("--"); - mitk::TNodePredicateDataType::Pointer isImage = mitk::TNodePredicateDataType::New(); + + mitk::TNodePredicateDataType::Pointer imageP = mitk::TNodePredicateDataType::New(); + mitk::NodePredicateDimension::Pointer dimP = mitk::NodePredicateDimension::New(3); + m_Controls->m_MapBox->SetDataStorage(this->GetDataStorage()); - m_Controls->m_MapBox->SetPredicate(isImage); + m_Controls->m_MapBox->SetPredicate(mitk::NodePredicateAnd::New(imageP, dimP)); m_Controls->m_MapBox->SetZeroEntryText("--"); + m_Controls->m_ParcellationBox->SetDataStorage(this->GetDataStorage()); + m_Controls->m_ParcellationBox->SetPredicate(mitk::NodePredicateAnd::New(imageP, dimP)); + m_Controls->m_ParcellationBox->SetZeroEntryText("--"); + FiberSelectionChanged(); } } void QmitkFiberClusteringView::SetFocus() { FiberSelectionChanged(); } void QmitkFiberClusteringView::FiberSelectionChanged() { if (m_Controls->m_TractBox->GetSelectedNode().IsNull()) m_Controls->m_StartButton->setEnabled(false); else m_Controls->m_StartButton->setEnabled(true); } void QmitkFiberClusteringView::StartClustering() { if (m_Controls->m_TractBox->GetSelectedNode().IsNull()) return; mitk::DataNode::Pointer node = m_Controls->m_TractBox->GetSelectedNode(); float clusterSize = m_Controls->m_ClusterSizeBox->value(); mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); float max_d = 0; int i=1; std::vector< float > distances; while (max_d < fib->GetGeometry()->GetDiagonalLength()/2) { distances.push_back(clusterSize*i); max_d = clusterSize*i; ++i; } itk::TractClusteringFilter::Pointer clusterer = itk::TractClusteringFilter::New(); clusterer->SetDistances(distances); clusterer->SetTractogram(fib); if (m_Controls->m_InCentroidsBox->GetSelectedNode().IsNotNull()) { mitk::FiberBundle::Pointer in_centroids = dynamic_cast(m_Controls->m_InCentroidsBox->GetSelectedNode()->GetData()); clusterer->SetInCentroids(in_centroids); } - switch (m_Controls->m_MetricBox->currentIndex()) + + std::vector< mitk::ClusteringMetric* > metrics; + + if (m_Controls->m_MetricBox1->isChecked()) + metrics.push_back(new mitk::ClusteringMetricEuclideanMean()); + if (m_Controls->m_MetricBox2->isChecked()) + metrics.push_back(new mitk::ClusteringMetricEuclideanStd()); + if (m_Controls->m_MetricBox3->isChecked()) + metrics.push_back(new mitk::ClusteringMetricEuclideanMax()); + if (m_Controls->m_ParcellationBox->GetSelectedNode().IsNotNull() && m_Controls->m_MetricBox4->isChecked()) { - case 0: - clusterer->SetMetric(itk::TractClusteringFilter::Metric::MDF); - break; - case 1: - clusterer->SetMetric(itk::TractClusteringFilter::Metric::MDF_STD); - break; - case 2: - clusterer->SetMetric(itk::TractClusteringFilter::Metric::MAX_MDF); - break; - } + mitk::Image::Pointer mitk_map = dynamic_cast(m_Controls->m_ParcellationBox->GetSelectedNode()->GetData()); - if (m_Controls->m_MapBox->GetSelectedNode().IsNotNull()) + ShortImageType::Pointer itk_map = ShortImageType::New(); + mitk::CastToItkImage(mitk_map, itk_map); + mitk::ClusteringMetricAnatomic* metric = new mitk::ClusteringMetricAnatomic(); + metric->SetParcellations({itk_map}); + metrics.push_back(metric); + } + if (m_Controls->m_MapBox->GetSelectedNode().IsNotNull() && m_Controls->m_MetricBox5->isChecked()) { mitk::Image::Pointer mitk_map = dynamic_cast(m_Controls->m_MapBox->GetSelectedNode()->GetData()); - if (mitk_map->GetDimension()==3) - { - FloatImageType::Pointer itk_map = FloatImageType::New(); - mitk::CastToItkImage(mitk_map, itk_map); - clusterer->SetScalarMap(itk_map); - } + + FloatImageType::Pointer itk_map = FloatImageType::New(); + mitk::CastToItkImage(mitk_map, itk_map); + mitk::ClusteringMetricScalarMap* metric = new mitk::ClusteringMetricScalarMap(); + metric->SetImages({itk_map}); + metric->SetScale(distances.at(0)); + metrics.push_back(metric); } + + if (metrics.empty()) + { + QMessageBox::warning(nullptr, "Warning", "No metric selected!"); + return; + } + + clusterer->SetMetrics(metrics); clusterer->SetMergeDuplicateThreshold(m_Controls->m_MergeDuplicatesBox->value()); - clusterer->SetScale(m_Controls->m_MapScaleBox->value()); clusterer->SetNumPoints(m_Controls->m_FiberPointsBox->value()); clusterer->SetMaxClusters(m_Controls->m_MaxClustersBox->value()); clusterer->SetMinClusterSize(m_Controls->m_MinFibersBox->value()); clusterer->Update(); std::vector tracts = clusterer->GetOutTractograms(); std::vector centroids = clusterer->GetOutCentroids(); unsigned int c = 0; for (auto f : tracts) { mitk::DataNode::Pointer new_node = mitk::DataNode::New(); new_node->SetData(f); new_node->SetName("Cluster_" + boost::lexical_cast(c)); if (m_Controls->m_InCentroidsBox->GetSelectedNode().IsNotNull()) this->GetDataStorage()->Add(new_node, m_Controls->m_InCentroidsBox->GetSelectedNode()); else this->GetDataStorage()->Add(new_node, node); node->SetVisibility(false); if (m_Controls->m_CentroidsBox->isChecked()) { mitk::DataNode::Pointer new_node2 = mitk::DataNode::New(); new_node2->SetData(centroids.at(c)); new_node2->SetName("Centroid_" + boost::lexical_cast(c)); this->GetDataStorage()->Add(new_node2, new_node); } ++c; } } void QmitkFiberClusteringView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList& ) { } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberClusteringView.h b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberClusteringView.h index 049e4df15b..585ff5f55e 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberClusteringView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberClusteringView.h @@ -1,67 +1,68 @@ /*=================================================================== 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 QmitkFiberClusteringView_h #define QmitkFiberClusteringView_h #include #include "ui_QmitkFiberClusteringViewControls.h" #include /*! \brief Cluster fibers by shape and image properties */ class QmitkFiberClusteringView : public QmitkAbstractView { // this is needed for all Qt objects that should have a Qt meta-object // (everything that derives from QObject and wants to have signal/slots) Q_OBJECT public: static const std::string VIEW_ID; QmitkFiberClusteringView(); virtual ~QmitkFiberClusteringView(); virtual void CreateQtPartControl(QWidget *parent) override; typedef itk::Image< float, 3 > FloatImageType; + typedef itk::Image< short, 3 > ShortImageType; /// /// Sets the focus to an internal widget. /// virtual void SetFocus() override; protected slots: void StartClustering(); void FiberSelectionChanged(); protected: /// \brief called by QmitkAbstractView when DataManager's selection has changed virtual void OnSelectionChanged(berry::IWorkbenchPart::Pointer part, const QList& nodes) override; Ui::QmitkFiberClusteringViewControls* m_Controls; }; #endif // _QMITKFIBERTRACKINGVIEW_H_INCLUDED diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberClusteringViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberClusteringViewControls.ui index 722d0a0a65..b7440970f3 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberClusteringViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberClusteringViewControls.ui @@ -1,313 +1,359 @@ QmitkFiberClusteringViewControls 0 0 484 574 Form - + + + + Qt::Vertical + + + + 20 + 40 + + + + + QFrame::NoFrame QFrame::Raised 0 0 0 0 6 - - + + + + If set, the input tractogram is clustered around the input centroids and no new clusters are created. + + + + + - Metric: + Output Centroids: - + false 0 0 200 16777215 11 Start - + + + + Cluster Size: + + + + + + + Fiber Points: + + + + Only output the N largest clusters. Zero means no limit. 99999999 - 0 + 10 - + Fibers are resampled to the desired number of points for clustering. Smaller is faster but less accurate. 2 9999999 12 - - - - 1 + + + + Input Centroids: - - - MDF - - - - - MDF_STD - - - - - MAX_MDF - - - - - - If set, distance metrics include distance between scalar map values along the tract. + + + + - - + + - Tractogram: + Min. Fibers per Cluster: - + - Merge duplicate clusters withthe specified distance threshold. If threshold is < 0, the threshold is set to the specified cluster size. + Merge duplicate clusters withthe specified distance threshold. If threshold is < 0, the threshold is set to half of the specified cluster size. -1.000000000000000 99999.000000000000000 -1.000000000000000 - - - - - - - Cluster size: - - - - - - - Map scale: - - - - - + + - Scalar map: + Max. Clusters: - - - Min. fibers per cluster: - - - - - - - Output centroids: - - - - - + - Fiber points: + Merge Duplicate Clusters: - - + + - Mergde duplicate clusters: + Tractogram: - + Only output clusters with ate least the specified number of fibers. 1 9999999 - 1 - - - - - - - Max. clusters: - - - - - - - + 50 - + Cluster size in mm. 1 9999999 - 10 + 20 - - - - Scaling factor for the scalar map (0=auto). + + + + + + + + + + Metrics + + + + 0 + + + 0 + + + 0 + + + 0 + + + + + Euclidean - - 99999.000000000000000 + + true + + + + + + + Euclidean Maximum - + - Input centroids: + Euclidean with STDEV - - - - If set, the input tractogram is clustered around the input centroids and no new clusters are created. + + + + QFrame::NoFrame + + + QFrame::Raised + + + 0 + + + 0 + + + 0 + + + 0 + + + 6 + + + + + Distance is based on the selected parcellation. + + + Anatomical + + + + + + + Distance is based on the selected parcellation. + + + + + + + Distance is based on the scalar map values along the tract. + + + Scalar Map + + + + + + + Distance is based on the scalar map values along the tract. + + + + - - - - Qt::Vertical - - - - 20 - 40 - - - - QmitkDataStorageComboBox QComboBox
QmitkDataStorageComboBox.h
 QmitkDataStorageComboBoxWithSelectNone QComboBox
QmitkDataStorageComboBoxWithSelectNone.h
 diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberFitViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberFitViewControls.ui index b7ded5b31e..7905b595a8 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberFitViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberFitViewControls.ui @@ -1,179 +1,179 @@ QmitkFiberFitViewControls 0 0 484 574 Form QFrame::NoFrame QFrame::Raised 0 0 0 0 6 Select a peak or raw diffusion-weighted image. true - Input image: + Input Image: - Suppress outliers: + Suppress Outliers: Modifier for regularization. 999999.000000000000000 0.100000000000000 0.100000000000000 false 0 0 200 16777215 11 Start Tractogram: λ: - Output residuals: + Output Residuals: false Qt::Vertical 20 40 QmitkDataStorageComboBox QComboBox
QmitkDataStorageComboBox.h
 diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkTractometryView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkTractometryView.cpp index 05c6ec0191..90d3a9ed04 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkTractometryView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkTractometryView.cpp @@ -1,277 +1,287 @@ /*=================================================================== 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 "QmitkTractometryView.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include const std::string QmitkTractometryView::VIEW_ID = "org.mitk.views.tractometry"; using namespace mitk; QmitkTractometryView::QmitkTractometryView() : QmitkAbstractView() , m_Controls( nullptr ) { } // Destructor QmitkTractometryView::~QmitkTractometryView() { } void QmitkTractometryView::CreateQtPartControl( QWidget *parent ) { // build up qt view, unless already done if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkTractometryViewControls; m_Controls->setupUi( parent ); mitk::TNodePredicateDataType::Pointer imageP = mitk::TNodePredicateDataType::New(); mitk::NodePredicateDimension::Pointer dimP = mitk::NodePredicateDimension::New(3); m_Controls->m_ImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_ImageBox->SetPredicate(mitk::NodePredicateAnd::New(imageP, dimP)); this->m_Controls->m_ChartWidget->SetXAxisLabel("Tract position"); this->m_Controls->m_ChartWidget->SetYAxisLabel("Image Value"); } } void QmitkTractometryView::OnPageSuccessfullyLoaded() { -// berry::IPreferencesService* prefService = berry::WorkbenchPlugin::GetDefault()->GetPreferencesService(); -// berry::IPreferences::Pointer m_StylePref = prefService->GetSystemPreferences()->Node(berry::QtPreferences::QT_STYLES_NODE); - -// QString styleName = m_StylePref->Get(berry::QtPreferences::QT_STYLE_NAME, ""); - -// if (styleName == ":/org.blueberry.ui.qt/darkstyle.qss") -// { -// this->m_Controls->m_ChartWidget->SetTheme(QmitkChartWidget::ChartStyle::darkstyle); -// } -// else -// { -// this->m_Controls->m_ChartWidget->SetTheme(QmitkChartWidget::ChartStyle::lightstyle); -// } + berry::IPreferencesService* prefService = berry::WorkbenchPlugin::GetDefault()->GetPreferencesService(); + berry::IPreferences::Pointer m_StylePref = prefService->GetSystemPreferences()->Node(berry::QtPreferences::QT_STYLES_NODE); + + QString styleName = m_StylePref->Get(berry::QtPreferences::QT_STYLE_NAME, ""); + + if (styleName == ":/org.blueberry.ui.qt/darkstyle.qss") + { + this->m_Controls->m_ChartWidget->SetTheme(QmitkChartWidget::ChartStyle::darkstyle); + } + else + { + this->m_Controls->m_ChartWidget->SetTheme(QmitkChartWidget::ChartStyle::lightstyle); + } } void QmitkTractometryView::SetFocus() { } -bool QmitkTractometryView::Flip(vtkSmartPointer< vtkPolyData > polydata, int i) +bool QmitkTractometryView::Flip(vtkSmartPointer< vtkPolyData > polydata1, int i, vtkSmartPointer< vtkPolyData > ref_poly) { float d_direct = 0; float d_flipped = 0; - vtkCell* cell1 = polydata->GetCell(0); + vtkCell* cell1 = polydata1->GetCell(0); + if (ref_poly!=nullptr) + cell1 = ref_poly->GetCell(0); int numPoints1 = cell1->GetNumberOfPoints(); vtkPoints* points1 = cell1->GetPoints(); - vtkCell* cell2 = polydata->GetCell(i); + vtkCell* cell2 = polydata1->GetCell(i); vtkPoints* points2 = cell2->GetPoints(); for (int j=0; jGetPoint(j); double* p2 = points2->GetPoint(j); d_direct = (p1[0]-p2[0])*(p1[0]-p2[0]) + (p1[1]-p2[1])*(p1[1]-p2[1]) + (p1[2]-p2[2])*(p1[2]-p2[2]); double* p3 = points2->GetPoint(numPoints1-j-1); d_flipped = (p1[0]-p3[0])*(p1[0]-p3[0]) + (p1[1]-p3[1])*(p1[1]-p3[1]) + (p1[2]-p3[2])*(p1[2]-p3[2]); } if (d_direct>d_flipped) return true; return false; } template void QmitkTractometryView::ImageValuesAlongTract(const mitk::PixelType, mitk::Image::Pointer image, mitk::FiberBundle::Pointer fib, std::vector > &data) { int num_points = 100; mitk::ImagePixelReadAccessor readimage(image, image->GetVolumeData(0)); mitk::FiberBundle::Pointer working_fib = fib->GetDeepCopy(); working_fib->ResampleToNumPoints(num_points); vtkSmartPointer< vtkPolyData > polydata = working_fib->GetFiberPolyData(); std::vector > all_values; std::vector< double > mean_values; for (int i=0; iGetNumFibers(); ++i) { vtkCell* cell = polydata->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); std::vector< double > fib_vals; bool flip = false; if (i>0) flip = Flip(polydata, i); + else if (m_ReferencePolyData!=nullptr) + flip = Flip(polydata, 0, m_ReferencePolyData); for (int j=0; jGetPoint(numPoints - j - 1); else p = points->GetPoint(j); Point3D px; px[0] = p[0]; px[1] = p[1]; px[2] = p[2]; double pixelValue = readimage.GetPixelByWorldCoordinates(px); fib_vals.push_back(pixelValue); mean += pixelValue; if (pixelValuemax) max = pixelValue; mean_values.at(j) += pixelValue; } all_values.push_back(fib_vals); } + if (m_ReferencePolyData==nullptr) + m_ReferencePolyData = polydata; + std::vector< double > std_values1; std::vector< double > std_values2; for (int i=0; iGetNumFibers(); double stdev = 0; for (unsigned int j=0; jGetNumberOfPoints(); } std::string QmitkTractometryView::RGBToHexString(double *rgb) { std::ostringstream os; for (int i = 0; i < 3; ++i) { os << std::setw(2) << std::setfill('0') << std::hex << static_cast(rgb[i] * 255); } return os.str(); } void QmitkTractometryView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList& nodes) { if(m_Controls->m_ImageBox->GetSelectedNode().IsNull()) return; + m_ReferencePolyData = nullptr; mitk::Image::Pointer image = dynamic_cast(m_Controls->m_ImageBox->GetSelectedNode()->GetData()); vtkSmartPointer lookupTable = vtkSmartPointer::New(); - lookupTable->SetTableRange(0.0, 0.8); + lookupTable->SetTableRange(0.0, 1.0); lookupTable->Build(); int num_tracts = 0; for (auto node: nodes) if ( dynamic_cast(node->GetData()) ) num_tracts++; int c = 1; this->m_Controls->m_ChartWidget->Clear(); for (auto node: nodes) { if ( dynamic_cast(node->GetData()) ) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); std::vector< std::vector< double > > data; mitkPixelTypeMultiplex3( ImageValuesAlongTract, image->GetPixelType(), image, fib, data ); m_Controls->m_ChartWidget->AddData1D(data.at(0), node->GetName() + " Mean", QmitkChartWidget::ChartType::line); if (m_Controls->m_StDevBox->isChecked()) { this->m_Controls->m_ChartWidget->AddData1D(data.at(1), node->GetName() + " +STDEV", QmitkChartWidget::ChartType::line); this->m_Controls->m_ChartWidget->AddData1D(data.at(2), node->GetName() + " -STDEV", QmitkChartWidget::ChartType::line); } double color[3]; if (num_tracts>1) { float scalar_color = ( (float)c/num_tracts - 1.0/num_tracts )/(1.0-1.0/num_tracts); lookupTable->GetColor(1.0 - scalar_color, color); } else lookupTable->GetColor(0, color); this->m_Controls->m_ChartWidget->SetColor(node->GetName() + " Mean", RGBToHexString(color)); if (m_Controls->m_StDevBox->isChecked()) { color[0] *= 0.8; color[1] *= 0.8; color[2] *= 0.8; this->m_Controls->m_ChartWidget->SetColor(node->GetName() + " +STDEV", RGBToHexString(color)); this->m_Controls->m_ChartWidget->SetColor(node->GetName() + " -STDEV", RGBToHexString(color)); } this->m_Controls->m_ChartWidget->Show(true); this->m_Controls->m_ChartWidget->SetShowDataPoints(false); ++c; } } + + MITK_INFO << "OnSelectionChanged DONE"; } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkTractometryView.h b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkTractometryView.h index 00f1c1a986..a13f7ae875 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkTractometryView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkTractometryView.h @@ -1,69 +1,72 @@ /*=================================================================== 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 QmitkTractometryView_h #define QmitkTractometryView_h #include #include "ui_QmitkTractometryViewControls.h" #include #include #include #include /*! \brief Weight fibers by linearly fitting them to the image data. */ class QmitkTractometryView : public QmitkAbstractView { // this is needed for all Qt objects that should have a Qt meta-object // (everything that derives from QObject and wants to have signal/slots) Q_OBJECT public: static const std::string VIEW_ID; QmitkTractometryView(); virtual ~QmitkTractometryView(); virtual void CreateQtPartControl(QWidget *parent) override; template void ImageValuesAlongTract(const mitk::PixelType, mitk::Image::Pointer image, mitk::FiberBundle::Pointer fib, std::vector< std::vector< double > >& data); virtual void SetFocus() override; void OnPageSuccessfullyLoaded(); + protected slots: protected: /// \brief called by QmitkAbstractView when DataManager's selection has changed virtual void OnSelectionChanged(berry::IWorkbenchPart::Pointer part, const QList& nodes) override; Ui::QmitkTractometryViewControls* m_Controls; - bool Flip(vtkSmartPointer< vtkPolyData > polydata, int i); + bool Flip(vtkSmartPointer< vtkPolyData > polydata1, int i, vtkSmartPointer ref_poly=nullptr); std::string RGBToHexString(double *rgb); + + vtkSmartPointer< vtkPolyData > m_ReferencePolyData; }; #endif // _QMITKFIBERTRACKINGVIEW_H_INCLUDED diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkTractometryViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkTractometryViewControls.ui index f97c61b283..924023a42e 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkTractometryViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkTractometryViewControls.ui @@ -1,107 +1,107 @@ QmitkTractometryViewControls 0 0 484 574 Form QFrame::NoFrame QFrame::Raised 0 0 0 0 6 - Select a peak or raw diffusion-weighted image. + - Input image: + Input Image: 0 100 Qt::Vertical 20 40 Show STDEV true QmitkChartWidget QWidget
QmitkChartWidget.h
 QmitkDataStorageComboBox QComboBox
QmitkDataStorageComboBox.h