diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetric.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetric.h index 17fac83341..7faaf59f6c 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetric.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetric.h @@ -1,61 +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; + float GetScale() const; void SetScale(float Scale); protected: float m_Scale; }; -float ClusteringMetric::SetScale() const +float ClusteringMetric::GetScale() 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 index 9e8db560d4..cb8be80db4 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetricAnatomic.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetricAnatomic.h @@ -1,199 +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 +#ifndef _ClusteringMetricAnatomic +#define _ClusteringMetricAnatomic #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 index cd2ffad9f0..24f1b52e7c 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetricEuclideanMax.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetricEuclideanMax.h @@ -1,71 +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 +#ifndef _ClusteringMetricEuMax +#define _ClusteringMetricEuMax #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 index dd0d6bf624..9cddeccfff 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetricEuclideanMean.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetricEuclideanMean.h @@ -1,62 +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 +#ifndef _ClusteringMetricEuMean +#define _ClusteringMetricEuMean #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 index 7314bd1c88..1f4ca84b3e 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetricEuclideanStd.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetricEuclideanStd.h @@ -1,77 +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 +#ifndef _ClusteringMetricEuStd +#define _ClusteringMetricEuStd #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/mitkClusteringMetricInnerAngles.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetricInnerAngles.h new file mode 100644 index 0000000000..90206f1163 --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetricInnerAngles.h @@ -0,0 +1,125 @@ +/*=================================================================== + +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 _ClusteringMetricInnerAngles +#define _ClusteringMetricInnerAngles + +#include + +namespace mitk +{ + +/** +* \brief Fiber clustering metric based on the angles between certain parts of the streamline */ + +class MITKFIBERTRACKING_EXPORT ClusteringMetricInnerAngles : public ClusteringMetric +{ + +public: + + ClusteringMetricInnerAngles(){} + virtual ~ClusteringMetricInnerAngles(){} + + float CalculateDistance(vnl_matrix& s, vnl_matrix& t, bool &flipped) + { + int p1 = 0; + int p2 = s.cols()/4; + int p3 = s.cols()/2; + int p4 = 3*s.cols()/4; + int p5 = s.cols()-1; + + float a1_s = 0; + float a2_s = 0; + float a3_s = 0; + float a1_t = 0; + float a2_t = 0; + float a3_t = 0; + + { + vnl_vector v1 = s.get_column(p1)-s.get_column(p2); v1.normalize(); + vnl_vector v2 = s.get_column(p3)-s.get_column(p2); v2.normalize(); + a1_s = dot_product(v1,v2); + a1_s = std::acos( a1_s ) * 180.0/M_PI; + } + + { + vnl_vector v1 = s.get_column(p1)-s.get_column(p3); v1.normalize(); + vnl_vector v2 = s.get_column(p5)-s.get_column(p3); v2.normalize(); + a2_s = dot_product(v1,v2); + a2_s = std::acos( a2_s ) * 180.0/M_PI; + } + + { + vnl_vector v1 = s.get_column(p3)-s.get_column(p4); v1.normalize(); + vnl_vector v2 = s.get_column(p5)-s.get_column(p4); v2.normalize(); + a3_s = dot_product(v1,v2); + a3_s = std::acos( a3_s ) * 180.0/M_PI; + } + + // + { + vnl_vector v1 = t.get_column(p1)-t.get_column(p2); v1.normalize(); + vnl_vector v2 = t.get_column(p3)-t.get_column(p2); v2.normalize(); + a1_t = dot_product(v1,v2); + a1_t = std::acos( a1_t ) * 180.0/M_PI; + } + + { + vnl_vector v1 = t.get_column(p1)-t.get_column(p3); v1.normalize(); + vnl_vector v2 = t.get_column(p5)-t.get_column(p3); v2.normalize(); + a2_t = dot_product(v1,v2); + a2_t = std::acos( a2_t ) * 180.0/M_PI; + } + + { + vnl_vector v1 = t.get_column(p3)-t.get_column(p4); v1.normalize(); + vnl_vector v2 = t.get_column(p5)-t.get_column(p4); v2.normalize(); + a3_t = dot_product(v1,v2); + a3_t = std::acos( a3_t ) * 180.0/M_PI; + } + + float d_direct = 0; + float d_flipped = 0; + + int inc = s.cols()/4; + for (unsigned int i=0; id_flipped) + { + flipped = true; + float d1 = std::fabs(a1_s-a3_t); + float d2 = std::fabs(a2_s-a2_t); + float d3 = std::fabs(a3_s-a1_t); + return m_Scale * std::max( d1, std::max(d2,d3) ); + } + flipped = false; + float d1 = std::fabs(a1_s-a1_t); + float d2 = std::fabs(a2_s-a2_t); + float d3 = std::fabs(a3_s-a3_t); + return m_Scale * std::max( d1, std::max(d2,d3) ); + } + +protected: + +}; + +} + +#endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetricScalarMap.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetricScalarMap.h index 736d705c78..ff8108be40 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetricScalarMap.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/ClusteringMetrics/mitkClusteringMetricScalarMap.h @@ -1,126 +1,127 @@ /*=================================================================== 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/itkTractClusteringFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractClusteringFilter.cpp index 41e656680e..fbade99dae 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractClusteringFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractClusteringFilter.cpp @@ -1,486 +1,497 @@ /*=================================================================== 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_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; } 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; } 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; streamline.set_size(3, m_NumPoints); streamline.fill(0.0); for (int j=0; jGetPoint(j, cand); 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; 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)/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; // 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; 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); + while (clusters.size()>5000) + { + MITK_INFO << "Number of clusters: " << clusters.size(); + MITK_INFO << "Increasing cluster size"; + for (unsigned int i=0; i > 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) ) + 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/cmdapps/FiberProcessing/CMakeLists.txt b/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/CMakeLists.txt index 0c57ddd8c9..642b53de57 100755 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/CMakeLists.txt +++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/CMakeLists.txt @@ -1,46 +1,47 @@ option(BUILD_DiffusionFiberProcessingCmdApps "Build commandline tools for diffusion fiber processing" OFF) if(BUILD_DiffusionFiberProcessingCmdApps OR MITK_BUILD_ALL_APPS) # needed include directories include_directories( ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR} ) # list of diffusion cmdapps # if an app requires additional dependencies # they are added after a "^^" and separated by "_" set( diffusionFiberProcessingcmdapps TractDensity^^MitkFiberTracking Sift2WeightCopy^^MitkFiberTracking FiberExtraction^^MitkFiberTracking FiberExtractionRoi^^MitkFiberTracking FiberProcessing^^MitkFiberTracking FitFibersToImage^^MitkFiberTracking FiberDirectionExtraction^^MitkFiberTracking FiberJoin^^MitkFiberTracking FiberClustering^^MitkFiberTracking + GetOverlappingTracts^^MitkFiberTracking ) foreach(diffusionFiberProcessingcmdapp ${diffusionFiberProcessingcmdapps}) # extract cmd app name and dependencies string(REPLACE "^^" "\\;" cmdapp_info ${diffusionFiberProcessingcmdapp}) set(cmdapp_info_list ${cmdapp_info}) list(GET cmdapp_info_list 0 appname) list(GET cmdapp_info_list 1 raw_dependencies) string(REPLACE "_" "\\;" dependencies "${raw_dependencies}") set(dependencies_list ${dependencies}) mitkFunctionCreateCommandLineApp( NAME ${appname} DEPENDS MitkCore MitkDiffusionCore ${dependencies_list} PACKAGE_DEPENDS ITK ) endforeach() if(EXECUTABLE_IS_ENABLED) MITK_INSTALL_TARGETS(EXECUTABLES ${EXECUTABLE_TARGET}) endif() endif() diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberClustering.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberClustering.cpp index 30e2b5e2cc..f645349ce5 100644 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberClustering.cpp +++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberClustering.cpp @@ -1,225 +1,228 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include #include #include +#include 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("metrics", "", mitkCommandLineParser::StringList, "Metrics:", "EU_MEAN, EU_STD, EU_MAX, ANAT, MAP"); + parser.addArgument("metrics", "", mitkCommandLineParser::StringList, "Metrics:", "EU_MEAN, EU_STD, EU_MAX, ANAT, MAP, ANGLES"); parser.addArgument("input_centroids", "", mitkCommandLineParser::String, "Input centroids:", ""); parser.addArgument("scalar_map", "", mitkCommandLineParser::String, "Scalar map:", ""); 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 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::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); } std::vector< mitk::ClusteringMetric* > metrics; for (auto m : metric_strings) { 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=="ANGLES") + metrics.push_back({new mitk::ClusteringMetricInnerAngles()}); 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!="") { 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->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/FiberProcessing/FiberExtractionRoi.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberExtractionRoi.cpp index 7511ba5d2e..cff0bed217 100755 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberExtractionRoi.cpp +++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberExtractionRoi.cpp @@ -1,126 +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. ===================================================================*/ #include #include "mitkCommandLineParser.h" #include #include #include #include #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include typedef itksys::SystemTools ist; typedef itk::Image ItkUcharImgType; 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 Extract fibers from a tractogram using binary image ROIs */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; - parser.setTitle("Fiber Extraction"); + parser.setTitle("Fiber Extraction With ROI Image"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setContributor("MIC"); parser.setDescription("Extract fibers from a tractogram using binary image ROIs"); parser.setArgumentPrefix("--", "-"); parser.addArgument("input", "i", mitkCommandLineParser::String, "Input:", "input tractogram (.fib/.trk/.tck/.dcm)", us::Any(), false); parser.addArgument("out", "o", mitkCommandLineParser::String, "Output:", "output tractogram", us::Any(), false); parser.addArgument("rois", "", mitkCommandLineParser::StringList, "ROI images:", "Images with binary ROIs", us::Any(), false); parser.addArgument("both_ends", "", mitkCommandLineParser::Bool, "Both ends:", "Fibers are extracted if both endpoints are located in the ROI.", false, false); parser.addArgument("overlap_fraction", "", mitkCommandLineParser::Float, "Overlap fraction:", "Extract by overlap, not by endpoints. Extract fibers that overlap to at least the provided factor (0-1) with the ROI.", -1, false); parser.addArgument("invert", "", mitkCommandLineParser::Bool, "Invert:", "invert mask image", false, false); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string inFib = us::any_cast(parsedArgs["input"]); std::string outFib = us::any_cast(parsedArgs["out"]); mitkCommandLineParser::StringContainerType roi_files = us::any_cast(parsedArgs["rois"]); bool both_ends = false; if (parsedArgs.count("both_ends")) both_ends = us::any_cast(parsedArgs["both_ends"]); bool invert = false; if (parsedArgs.count("invert")) invert = us::any_cast(parsedArgs["invert"]); float overlap_fraction = -1; if (parsedArgs.count("overlap_fraction")) overlap_fraction = us::any_cast(parsedArgs["overlap_fraction"]); bool any_point = false; if (overlap_fraction>=0) any_point = true; try { // load fiber bundle mitk::FiberBundle::Pointer inputTractogram = dynamic_cast(mitk::IOUtil::Load(inFib)[0].GetPointer()); for (std::size_t i=0; iExtractFiberSubset(mask, any_point, invert, both_ends, overlap_fraction); mitk::IOUtil::Save(result, outFib); } catch(...) { std::cout << "could not load: " << roi_files.at(i); return EXIT_FAILURE; } } } 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/FiberProcessing/GetOverlappingTracts.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/GetOverlappingTracts.cpp new file mode 100755 index 0000000000..4db8ad3a87 --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/GetOverlappingTracts.cpp @@ -0,0 +1,135 @@ +/*=================================================================== + +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 "mitkCommandLineParser.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#define _USE_MATH_DEFINES +#include + +typedef itksys::SystemTools ist; +typedef itk::Image ItkUcharImgType; + +/*! +\brief Extract fibers from a tractogram using binary image ROIs +*/ +int main(int argc, char* argv[]) +{ + mitkCommandLineParser parser; + + parser.setTitle("Move Overlapping Tracts"); + parser.setCategory("Fiber Tracking and Processing Methods"); + parser.setContributor("MIC"); + parser.setDescription("Move tracts that overlap with one of the reference tracts"); + + parser.setArgumentPrefix("--", "-"); + parser.addArgument("input", "i", mitkCommandLineParser::StringList, "Input:", "input tractograms (.fib/.trk/.tck/.dcm)", us::Any(), false); + parser.addArgument("reference", "r", mitkCommandLineParser::StringList, "Reference:", "reference tractograms (.fib/.trk/.tck/.dcm)", us::Any(), false); + parser.addArgument("out", "o", mitkCommandLineParser::OutputDirectory, "Output Folder:", "move input tracts that do/don't overlap here", us::Any(), false); + + parser.addArgument("overlap_fraction", "", mitkCommandLineParser::Float, "Overlap fraction:", "", 0.9); + parser.addArgument("move_overlapping", "", mitkCommandLineParser::Bool, ":", ""); + + std::map parsedArgs = parser.parseArguments(argc, argv); + if (parsedArgs.size()==0) + return EXIT_FAILURE; + + mitkCommandLineParser::StringContainerType input = us::any_cast(parsedArgs["input"]); + mitkCommandLineParser::StringContainerType reference = us::any_cast(parsedArgs["reference"]); + + std::string out_folder = us::any_cast(parsedArgs["out"]); + + bool move_overlapping = false; + if (parsedArgs.count("move_overlapping")) + move_overlapping = us::any_cast(parsedArgs["move_overlapping"]); + + float overlap_threshold = 0.9; + if (parsedArgs.count("overlap_fraction")) + overlap_threshold = us::any_cast(parsedArgs["overlap_fraction"]); + + try + { + itk::TractDensityImageFilter< ItkUcharImgType >::Pointer filter = itk::TractDensityImageFilter< ItkUcharImgType >::New(); + filter->SetDoFiberResampling(true); + filter->SetUpsamplingFactor(0.25); + filter->SetBinaryOutput(true); + std::vector< ItkUcharImgType::Pointer > masks; + for (auto f : reference) + { + mitk::FiberBundle::Pointer fib = dynamic_cast(mitk::IOUtil::Load(f)[0].GetPointer()); + filter->SetFiberBundle(fib); + + std::streambuf *old = cout.rdbuf(); // <-- save + std::stringstream ss; + std::cout.rdbuf (ss.rdbuf()); // <-- redirect + filter->Update(); + masks.push_back(filter->GetOutput()); + std::cout.rdbuf (old); // <-- restore + } + + boost::progress_display disp(input.size()); + for (auto f : input) + { + ++disp; + std::streambuf *old = cout.rdbuf(); // <-- save + std::stringstream ss; + std::cout.rdbuf (ss.rdbuf()); // <-- redirect + bool is_overlapping = false; + mitk::FiberBundle::Pointer fib = dynamic_cast(mitk::IOUtil::Load(f)[0].GetPointer()); + fib->ResampleLinear(2); + for (auto m : masks) + { + float overlap = fib->GetOverlap(m, false); + + if (overlap>=overlap_threshold) + { + is_overlapping = true; + break; + } + } + if (move_overlapping && is_overlapping) + ist::CopyAFile(f, out_folder + ist::GetFilenameName(f)); + else if (!move_overlapping && !is_overlapping) + ist::CopyAFile(f, out_folder + ist::GetFilenameName(f)); + 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/files.cmake b/Modules/DiffusionImaging/FiberTracking/files.cmake index 6eba166f25..1fe24d7019 100644 --- a/Modules/DiffusionImaging/FiberTracking/files.cmake +++ b/Modules/DiffusionImaging/FiberTracking/files.cmake @@ -1,101 +1,102 @@ 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 + Algorithms/ClusteringMetrics/mitkClusteringMetricInnerAngles.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 1d7a4d2b05..5f5c02575f 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,204 +1,224 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ // Blueberry #include #include #include "QmitkFiberClusteringView.h" #include #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 imageP = mitk::TNodePredicateDataType::New(); mitk::NodePredicateDimension::Pointer dimP = mitk::NodePredicateDimension::New(3); m_Controls->m_MapBox->SetDataStorage(this->GetDataStorage()); 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("--"); + + m_Controls->m_MetricBox6->setVisible(false); + m_Controls->m_MetricWeight6->setVisible(false); 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); } std::vector< mitk::ClusteringMetric* > metrics; - if (m_Controls->m_MetricBox1->isChecked()) - metrics.push_back(new mitk::ClusteringMetricEuclideanMean()); + { + mitk::ClusteringMetricEuclideanMean* metric = new mitk::ClusteringMetricEuclideanMean(); + metric->SetScale(m_Controls->m_MetricWeight1->value()); + metrics.push_back(metric); + } if (m_Controls->m_MetricBox2->isChecked()) - metrics.push_back(new mitk::ClusteringMetricEuclideanStd()); + { + mitk::ClusteringMetricEuclideanStd* metric = new mitk::ClusteringMetricEuclideanStd(); + metric->SetScale(m_Controls->m_MetricWeight2->value()); + metrics.push_back(metric); + } if (m_Controls->m_MetricBox3->isChecked()) - metrics.push_back(new mitk::ClusteringMetricEuclideanMax()); + { + mitk::ClusteringMetricEuclideanMax* metric = new mitk::ClusteringMetricEuclideanMax(); + metric->SetScale(m_Controls->m_MetricWeight3->value()); + metrics.push_back(metric); + } + if (m_Controls->m_MetricBox6->isChecked()) + { + mitk::ClusteringMetricInnerAngles* metric = new mitk::ClusteringMetricInnerAngles(); + metric->SetScale(m_Controls->m_MetricWeight6->value()); + metrics.push_back(metric); + } if (m_Controls->m_ParcellationBox->GetSelectedNode().IsNotNull() && m_Controls->m_MetricBox4->isChecked()) { mitk::Image::Pointer mitk_map = dynamic_cast(m_Controls->m_ParcellationBox->GetSelectedNode()->GetData()); ShortImageType::Pointer itk_map = ShortImageType::New(); mitk::CastToItkImage(mitk_map, itk_map); mitk::ClusteringMetricAnatomic* metric = new mitk::ClusteringMetricAnatomic(); metric->SetParcellations({itk_map}); + metric->SetScale(m_Controls->m_MetricWeight4->value()); 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()); 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)); + metric->SetScale(m_Controls->m_MetricWeight5->value()); 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->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/QmitkFiberClusteringViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberClusteringViewControls.ui index b7440970f3..48f6321379 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,359 +1,486 @@ 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. Output Centroids: false 0 0 200 16777215 11 Start Cluster Size: Fiber Points: Only output the N largest clusters. Zero means no limit. 99999999 10 Fibers are resampled to the desired number of points for clustering. Smaller is faster but less accurate. 2 9999999 12 Input Centroids: Min. Fibers per Cluster: 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 Max. Clusters: Merge Duplicate Clusters: Tractogram: Only output clusters with ate least the specified number of fibers. 1 9999999 50 Cluster size in mm. 1 9999999 20 Metrics - + 0 0 0 0 + + + Weighting factor for metric values. + + + 1 + + + 999.000000000000000 + + + 1.000000000000000 + + + + Euclidean true - + + + + Euclidean with STDEV + + + + Euclidean Maximum - - + + - Euclidean with STDEV + Inner Angles - + QFrame::NoFrame QFrame::Raised 0 0 0 0 6 - + Distance is based on the selected parcellation. Anatomical - - + + Distance is based on the selected parcellation. - + + + + + + + Weighting factor for metric values. + + + 1 + + + 999.000000000000000 + + + 1.000000000000000 + + + + + + + Weighting factor for metric values. + + + 1 + + + 999.000000000000000 + + + 1.000000000000000 + + + + + + + Weighting factor for metric values. + + + 1 + + + 999.000000000000000 + + + 1.000000000000000 + + + + + + + Weighting factor for metric values. + + + 1 + + + 999.000000000000000 + + + 1.000000000000000 + + + + + + + Weighting factor for metric values. + + + 1 + + + 999.000000000000000 + + + 30.000000000000000 + + + + + + + QFrame::NoFrame + + + QFrame::Raised + + + + 0 + + + 0 + + + 0 + + + 0 + + Distance is based on the scalar map values along the tract. Scalar Map - - + + Distance is based on the scalar map values along the tract. QmitkDataStorageComboBox QComboBox
QmitkDataStorageComboBox.h
 QmitkDataStorageComboBoxWithSelectNone QComboBox
QmitkDataStorageComboBoxWithSelectNone.h