diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractClusteringFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractClusteringFilter.cpp
index b7d9c8e655..ebf2399ada 100644
--- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractClusteringFilter.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractClusteringFilter.cpp
@@ -1,542 +1,594 @@
 /*===================================================================
 
 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 <math.h>
 #include <boost/progress.hpp>
 
 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()
 {
 
 }
 
 std::vector<TractClusteringFilter::Cluster> TractClusteringFilter::GetOutClusters() const
 {
   return m_OutClusters;
 }
 
 std::vector<mitk::FiberBundle::Pointer> TractClusteringFilter::GetOutCentroids() const
 {
   return m_OutCentroids;
 }
 
 void TractClusteringFilter::SetMetric(const Metric &Metric)
 {
   m_Metric = Metric;
 }
 
 std::vector<mitk::FiberBundle::Pointer> TractClusteringFilter::GetOutTractograms() const
 {
   return m_OutTractograms;
 }
 
 void TractClusteringFilter::SetDistances(const std::vector<float> &Distances)
 {
   m_Distances = Distances;
 }
 
+float TractClusteringFilter::CalcOverlap(vnl_matrix<float>& t)
+{
+  float overlap = 0;
+  if (m_FilterMask.IsNotNull())
+  {
+
+    for (unsigned int i=0; i<m_NumPoints; ++i)
+    {
+      vnl_vector<float> p = t.get_column(i);
+      itk::Point<float,3> 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<float>& s, vnl_matrix<float>& t, bool& flipped)
 {
   float d_direct = 0;
   float d_flipped = 0;
 
   for (unsigned int i=0; i<m_NumPoints; ++i)
   {
     d_direct += (s.get_column(i)-t.get_column(i)).magnitude();
     d_flipped += (s.get_column(i)-t.get_column(m_NumPoints-i-1)).magnitude();
   }
 
   if (d_direct>d_flipped)
   {
     flipped = true;
     return d_flipped/m_NumPoints;
   }
   flipped = false;
   return d_direct/m_NumPoints;
 }
 
 float TractClusteringFilter::CalcMDF_STD(vnl_matrix<float>& s, vnl_matrix<float>& t, bool& flipped)
 {
   float d_direct = 0;
   float d_flipped = 0;
 
   vnl_vector<float> dists_d; dists_d.set_size(m_NumPoints);
   vnl_vector<float> dists_f; dists_f.set_size(m_NumPoints);
 
   for (unsigned int i=0; i<m_NumPoints; ++i)
   {
     dists_d[i] = (s.get_column(i)-t.get_column(i)).magnitude();
     d_direct += dists_d[i];
 
     dists_f[i] = (s.get_column(i)-t.get_column(m_NumPoints-i-1)).magnitude();
     d_flipped += dists_f[i];
   }
 
   if (d_direct>d_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<float>& s, vnl_matrix<float>& t, bool& flipped)
 {
   float d_direct = 0;
   float d_flipped = 0;
 
   vnl_vector<float> dists_d; dists_d.set_size(m_NumPoints);
   vnl_vector<float> dists_f; dists_f.set_size(m_NumPoints);
 
   for (unsigned int i=0; i<m_NumPoints; ++i)
   {
     dists_d[i] = (s.get_column(i)-t.get_column(i)).magnitude();
     d_direct += dists_d[i];
 
     dists_f[i] = (s.get_column(i)-t.get_column(m_NumPoints-i-1)).magnitude();
     d_flipped += dists_f[i];
   }
 
   if (d_direct>d_flipped)
   {
     float d = dists_f.max_value();
     flipped = true;
     return d;
   }
 
   float d = dists_d.max_value();
   flipped = false;
   return d;
 }
 
 std::vector<vnl_matrix<float> > TractClusteringFilter::ResampleFibers(mitk::FiberBundle::Pointer tractogram)
 {
   mitk::FiberBundle::Pointer temp_fib = tractogram->GetDeepCopy();
-  temp_fib->ResampleToNumPoints(m_NumPoints);
+  if (m_DoResampling)
+    temp_fib->ResampleToNumPoints(m_NumPoints);
 
   std::vector< vnl_matrix<float> > out_fib;
 
   for (int i=0; i<temp_fib->GetFiberPolyData()->GetNumberOfCells(); i++)
   {
     vtkCell* cell = temp_fib->GetFiberPolyData()->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     vnl_matrix<float> streamline;
     if (m_ScalarMap.IsNull())
       streamline.set_size(3, m_NumPoints);
     else
       streamline.set_size(4, m_NumPoints);
     streamline.fill(0.0);
 
     for (int j=0; j<numPoints; j++)
     {
       double cand[3];
       points->GetPoint(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<float,3> wp; wp[0]=cand[0]; wp[1]=cand[1]; wp[2]=cand[2];
 
         itk::Index<3> idx;
         m_ScalarMap->TransformPhysicalPointToIndex(wp, idx);
         if (m_ScalarMap->GetLargestPossibleRegion().IsInside(idx))
           candV[3]=m_ScalarMap->GetPixel(idx)*m_Scale;
 
         streamline.set_column(j, candV);
       }
     }
 
     out_fib.push_back(streamline);
   }
 
   return out_fib;
 }
 
 std::vector< TractClusteringFilter::Cluster > TractClusteringFilter::ClusterStep(std::vector< long > f_indices, std::vector<float> 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<N; ++i)
   {
     vnl_matrix<float> 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<C.size(); ++k)
     {
       vnl_matrix<float> v = C.at(k).h / C.at(k).n;
       bool f = false;
       float d = 0;
       if (m_Metric==Metric::MDF)
         d = CalcMDF(t, v, f);
       else if (m_Metric==Metric::MDF_STD)
         d = CalcMDF_STD(t, v, f);
       else if (m_Metric==Metric::MAX_MDF)
         d = CalcMAX_MDF(t, v, f);
 
       if (d<min_cluster_distance)
       {
         min_cluster_distance = d;
         min_cluster_index = k;
         flip = f;
       }
     }
 
     if (min_cluster_index>=0 && min_cluster_distance<dist_thres)
     {
       C[min_cluster_index].I.push_back(f_indices.at(i));
       if (!flip)
         C[min_cluster_index].h += t;
       else
         C[min_cluster_index].h += t.fliplr();
       C[min_cluster_index].n += 1;
     }
     else
     {
       Cluster c;
       c.I.push_back(f_indices.at(i));
       c.h = t;
       c.n = 1;
       C.push_back(c);
     }
   }
 
   if (!distances.empty())
   {
     std::vector< Cluster > outC;
 #pragma omp parallel for
     for (int c=0; c<(int)C.size(); c++)
     {
 
       std::vector< Cluster > tempC = ClusterStep(C.at(c).I, distances);
       AppendCluster(outC, tempC);
     }
     return outC;
   }
   else
     return C;
 }
 
 void TractClusteringFilter::AppendCluster(std::vector< Cluster >& a, std::vector< Cluster >&b)
 {
   for (auto c : b)
     a.push_back(c);
 }
 
 void TractClusteringFilter::MergeDuplicateClusters(std::vector< TractClusteringFilter::Cluster >& clusters)
 {
   if (m_MergeDuplicateThreshold<0)
     m_MergeDuplicateThreshold = m_Distances.at(0);
   bool found = true;
 
   MITK_INFO << "Merging duplicate clusters with distance threshold " << m_MergeDuplicateThreshold;
+  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=0; k1<(int)clusters.size(); ++k1)
+    for (int k1=start; k1<(int)clusters.size(); ++k1)
     {
       Cluster c1 = clusters.at(k1);
       vnl_matrix<float> 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<float> 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);
+          //          if (m_Metric==Metric::MDF)
+          //            d = CalcMDF(t, v, f);
+          //          else if (m_Metric==Metric::MDF_STD)
+          //            d = CalcMDF_STD(t, v, f);
+          //          else if (m_Metric==Metric::MAX_MDF)
+          d = CalcMAX_MDF(t, v, f);
 
 #pragma omp critical
           if (d<m_MergeDuplicateThreshold)
           {
             merge_indices.push_back(k2);
             flip_indices.push_back(f);
           }
         }
       }
 
       for (unsigned int i=0; i<merge_indices.size(); ++i)
       {
         Cluster c2 = clusters.at(merge_indices.at(i));
         for (int i=0; i<c2.n; ++i)
         {
           clusters[k1].I.push_back(c2.I.at(i));
           clusters[k1].n += 1;
         }
         if (!flip_indices.at(i))
           clusters[k1].h += c2.h;
         else
           clusters[k1].h += c2.h.fliplr();
       }
 
       std::sort(merge_indices.begin(), merge_indices.end());
       for (unsigned int i=0; i<merge_indices.size(); ++i)
       {
-        clusters.erase (clusters.begin()+merge_indices.at(i));
+        clusters.erase (clusters.begin()+merge_indices.at(i)-i);
         found = true;
       }
       if (found)
         break;
     }
+    ++start;
   }
   MITK_INFO << "\nNumber of clusters after merging duplicates: " << clusters.size();
 }
 
 std::vector<TractClusteringFilter::Cluster> TractClusteringFilter::AddToKnownClusters(std::vector< long > f_indices, std::vector<vnl_matrix<float> >& centroids)
 {
   float dist_thres = m_Distances.at(0);
   int N = f_indices.size();
 
   std::vector< Cluster > C;
   vnl_matrix<float> 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<centroids.size(); ++i)
   {
     Cluster c;
     c.h.set_size(T.at(0).rows(), T.at(0).cols()); c.h.fill(0.0);
     c.f_id = i;
     C.push_back(c);
   }
 
+#pragma omp parallel for
   for (int i=0; i<N; ++i)
   {
     vnl_matrix<float> t = T.at(f_indices.at(i));
 
     int min_cluster_index = -1;
     float min_cluster_distance = 99999;
     bool flip = false;
 
-    int c_idx = 0;
-    for (vnl_matrix<float> centroid : centroids)
+    if (CalcOverlap(t)>=m_OverlapThreshold)
     {
-      bool f = false;
-      float d = 0;
+      int c_idx = 0;
+      for (vnl_matrix<float> centroid : centroids)
+      {
+        bool f = false;
+        float d = 0;
 
-      if (m_Metric==Metric::MDF)
-        d = CalcMDF(t, centroid, f);
-      else if (m_Metric==Metric::MDF_STD)
-        d = CalcMDF_STD(t, centroid, f);
-      else if (m_Metric==Metric::MAX_MDF)
-        d = CalcMAX_MDF(t, centroid, f);
+        if (m_Metric==Metric::MDF)
+          d = CalcMDF(t, centroid, f);
+        else if (m_Metric==Metric::MDF_STD)
+          d = CalcMDF_STD(t, centroid, f);
+        else if (m_Metric==Metric::MAX_MDF)
+          d = CalcMAX_MDF(t, centroid, f);
 
-      if (d<min_cluster_distance)
-      {
-        min_cluster_distance = d;
-        min_cluster_index = c_idx;
-        flip = f;
+        if (d<min_cluster_distance)
+        {
+          min_cluster_distance = d;
+          min_cluster_index = c_idx;
+          flip = f;
+        }
+        ++c_idx;
       }
-      ++c_idx;
     }
 
     if (min_cluster_index>=0 && min_cluster_distance<dist_thres)
     {
-      C[min_cluster_index].I.push_back(f_indices.at(i));
-      if (!flip)
-        C[min_cluster_index].h += t;
-      else
-        C[min_cluster_index].h += t.fliplr();
-      C[min_cluster_index].n += 1;
+#pragma omp critical
+      {
+        C[min_cluster_index].I.push_back(f_indices.at(i));
+        if (!flip)
+          C[min_cluster_index].h += t;
+        else
+          C[min_cluster_index].h += t.fliplr();
+        C[min_cluster_index].n += 1;
+      }
     }
     else
     {
-      no_fit.I.push_back(f_indices.at(i));
-      no_fit.n++;
+#pragma omp critical
+      {
+        no_fit.I.push_back(f_indices.at(i));
+        no_fit.n++;
+      }
     }
   }
   C.push_back(no_fit);
   return C;
 }
 
 void TractClusteringFilter::GenerateData()
 {
   m_OutTractograms.clear();
   m_OutCentroids.clear();
   m_OutClusters.clear();
 
   if (m_Scale==0)
     m_Scale = m_Distances.at(0);
-  MITK_INFO << "Scalar map scaling: " << m_Scale;
 
   T = ResampleFibers(m_Tractogram);
+  if (T.empty())
+  {
+    MITK_INFO << "No fibers in tractogram!";
+    return;
+  }
 
   std::vector< long > f_indices;
   for (unsigned int i=0; i<T.size(); ++i)
     f_indices.push_back(i);
   //  std::random_shuffle(f_indices.begin(), f_indices.end());
 
   Cluster no_match;
   std::vector< Cluster > clusters;
   if (m_InCentroids.IsNull())
   {
     MITK_INFO << "Clustering fibers";
     clusters = ClusterStep(f_indices, m_Distances);
     MITK_INFO << "Number of clusters: " << clusters.size();
     MergeDuplicateClusters(clusters);
     std::sort(clusters.begin(),clusters.end());
   }
   else
   {
-    MITK_INFO << "Clustering with input centroids";
     std::vector<vnl_matrix<float> > 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<vtkFloatArray> weights = vtkSmartPointer<vtkFloatArray>::New();
       vtkSmartPointer<vtkPolyData> 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<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
       vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
       vtkSmartPointer<vtkPolyData> polyData = vtkSmartPointer<vtkPolyData>::New();
       vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
       vnl_matrix<float> centroid_points = c.h / c.n;
       for (unsigned int j=0; j<centroid_points.cols(); j++)
       {
         double p[3];
         p[0] = centroid_points.get(0,j);
         p[1] = centroid_points.get(1,j);
         p[2] = centroid_points.get(2,j);
 
         vtkIdType id = vtkNewPoints->InsertNextPoint(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<vtkFloatArray> weights = vtkSmartPointer<vtkFloatArray>::New();
     vtkSmartPointer<vtkPolyData> 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 032c24ad5e..4de44967dd 100644
--- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractClusteringFilter.h
+++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractClusteringFilter.h
@@ -1,141 +1,151 @@
 /*===================================================================
 
 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 <mitkPlanarEllipse.h>
 #include <mitkFiberBundle.h>
 #include <mitkFiberfoxParameters.h>
 
 // ITK
 #include <itkProcessObject.h>
 
 // VTK
 #include <vtkSmartPointer.h>
 #include <vtkPolyData.h>
 #include <vtkCellArray.h>
 #include <vtkPoints.h>
 #include <vtkPolyLine.h>
 
 namespace itk{
 
 /**
 * \brief    */
 
 class TractClusteringFilter : public ProcessObject
 {
 public:
 
   struct Cluster
   {
     Cluster() : n(0), f_id(-1) {}
 
     vnl_matrix<float> 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)
-  itkGetMacro(NumPoints, unsigned int)
-  itkSetMacro(MinClusterSize, unsigned int)
-  itkGetMacro(MinClusterSize, unsigned int)
-  itkSetMacro(MaxClusters, unsigned int)
-  itkGetMacro(MaxClusters, unsigned int)
-  itkSetMacro(Scale, float)
-  itkGetMacro(Scale, float)
-  itkSetMacro(MergeDuplicateThreshold, float)
-  itkGetMacro(MergeDuplicateThreshold, float)
-
-  itkSetMacro(Tractogram, mitk::FiberBundle::Pointer)
-  itkSetMacro(InCentroids, mitk::FiberBundle::Pointer)
-  itkSetMacro(ScalarMap, FloatImageType::Pointer)
+  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<float> &Distances);
+  void SetDistances(const std::vector<float> &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<mitk::FiberBundle::Pointer> GetOutTractograms() const;
 
   void SetMetric(const Metric &Metric);
 
   std::vector<mitk::FiberBundle::Pointer> GetOutCentroids() const;
 
   std::vector<Cluster> GetOutClusters() const;
 
 protected:
 
   void GenerateData() override;
   std::vector< vnl_matrix<float> > ResampleFibers(FiberBundle::Pointer tractogram);
   float CalcMDF(vnl_matrix<float>& s, vnl_matrix<float>& t, bool &flipped);
   float CalcMDF_STD(vnl_matrix<float>& s, vnl_matrix<float>& t, bool &flipped);
   float CalcMAX_MDF(vnl_matrix<float>& s, vnl_matrix<float>& t, bool &flipped);
+  float CalcOverlap(vnl_matrix<float>& 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<vnl_matrix<float> > &centroids);
   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<vnl_matrix<float> >             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;
 };
 }
 
 #ifndef ITK_MANUAL_INSTANTIATION
 #include "itkTractClusteringFilter.cpp"
 #endif
 
 #endif
diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/CMakeLists.txt b/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/CMakeLists.txt
index 49e27af547..0c57ddd8c9 100755
--- a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/CMakeLists.txt
+++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/CMakeLists.txt
@@ -1,45 +1,46 @@
 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
     )
 
     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/FiberExtractionRoi.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberExtractionRoi.cpp
new file mode 100755
index 0000000000..7511ba5d2e
--- /dev/null
+++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberExtractionRoi.cpp
@@ -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.
+
+===================================================================*/
+
+#include <metaCommand.h>
+#include "mitkCommandLineParser.h"
+#include <usAny.h>
+#include <mitkIOUtil.h>
+#include <boost/lexical_cast.hpp>
+#include <mitkCoreObjectFactory.h>
+#include <mitkPlanarFigure.h>
+#include <mitkPlanarFigureComposite.h>
+#include <mitkFiberBundle.h>
+#include <mitkImageCast.h>
+#include <mitkImageToItk.h>
+#include <mitkStandaloneDataStorage.h>
+#include <itksys/SystemTools.hxx>
+
+#define _USE_MATH_DEFINES
+#include <math.h>
+
+typedef itksys::SystemTools ist;
+typedef itk::Image<unsigned char, 3>    ItkUcharImgType;
+
+ItkUcharImgType::Pointer LoadItkMaskImage(const std::string& filename)
+{
+  mitk::Image::Pointer img = dynamic_cast<mitk::Image*>(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.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<std::string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
+  if (parsedArgs.size()==0)
+    return EXIT_FAILURE;
+
+  std::string inFib = us::any_cast<std::string>(parsedArgs["input"]);
+  std::string outFib = us::any_cast<std::string>(parsedArgs["out"]);
+  mitkCommandLineParser::StringContainerType roi_files = us::any_cast<mitkCommandLineParser::StringContainerType>(parsedArgs["rois"]);
+
+  bool both_ends = false;
+  if (parsedArgs.count("both_ends"))
+    both_ends = us::any_cast<bool>(parsedArgs["both_ends"]);
+
+  bool invert = false;
+  if (parsedArgs.count("invert"))
+    invert = us::any_cast<bool>(parsedArgs["invert"]);
+
+  float overlap_fraction = -1;
+  if (parsedArgs.count("overlap_fraction"))
+    overlap_fraction = us::any_cast<float>(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::FiberBundle*>(mitk::IOUtil::Load(inFib)[0].GetPointer());
+
+    for (std::size_t i=0; i<roi_files.size(); ++i)
+    {
+      try
+      {
+        ItkUcharImgType::Pointer mask = LoadItkMaskImage(roi_files.at(i));
+        mitk::FiberBundle::Pointer result = inputTractogram->ExtractFiberSubset(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/Misc/CMakeLists.txt b/Modules/DiffusionImaging/FiberTracking/cmdapps/Misc/CMakeLists.txt
index e759910e7f..8342f17ffc 100644
--- a/Modules/DiffusionImaging/FiberTracking/cmdapps/Misc/CMakeLists.txt
+++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/Misc/CMakeLists.txt
@@ -1,36 +1,37 @@
 option(BUILD_DiffusionMiscCmdApps "Build miscellaneous commandline tools for diffusion" OFF)
 
 if(BUILD_DiffusionMiscCmdApps 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( diffusionMisccmdapps
     PeakExtraction^^MitkFiberTracking
     FileFormatConverter^^MitkFiberTracking
+    FlipPeaks^^MitkFiberTracking
     )
 
     foreach(diffusionMisccmdapp ${diffusionMisccmdapps})
       # extract cmd app name and dependencies
       string(REPLACE "^^" "\\;" cmdapp_info ${diffusionMisccmdapp})
       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()
     
 endif()
 
diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/Misc/FlipPeaks.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/Misc/FlipPeaks.cpp
new file mode 100644
index 0000000000..66d4857af8
--- /dev/null
+++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/Misc/FlipPeaks.cpp
@@ -0,0 +1,102 @@
+/*===================================================================
+
+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 <mitkImageCast.h>
+#include <mitkImage.h>
+#include <mitkIOUtil.h>
+#include "mitkCommandLineParser.h"
+#include <mitkPeakImage.h>
+#include <itkFlipPeaksFilter.h>
+#include <mitkImageToItk.h>
+
+/*!
+\brief Copies transformation matrix of one image to another
+*/
+int main(int argc, char* argv[])
+{
+  mitkCommandLineParser parser;
+
+  parser.setTitle("Flip Peaks");
+  parser.setCategory("Preprocessing Tools");
+  parser.setDescription("Flips the peaks of the input peak image along the given dimensions.");
+  parser.setContributor("MIC");
+
+  parser.setArgumentPrefix("--", "-");
+  parser.addArgument("", "i", mitkCommandLineParser::InputFile, "Input", "input image", us::Any(), false);
+  parser.addArgument("", "o", mitkCommandLineParser::OutputFile, "Output", "output image", us::Any(), false);
+  parser.addArgument("", "x", mitkCommandLineParser::Bool, "Flip x", "flip along x-axis");
+  parser.addArgument("", "y", mitkCommandLineParser::Bool, "Flip y", "flip along y-axis");
+  parser.addArgument("", "z", mitkCommandLineParser::Bool, "Flip z", "flip along z-axis");
+
+  std::map<std::string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
+  if (parsedArgs.size()==0)
+    return EXIT_FAILURE;
+
+  std::string imageName = us::any_cast<std::string>(parsedArgs["i"]);
+  std::string outImage = us::any_cast<std::string>(parsedArgs["o"]);
+
+  bool x = false;
+  if (parsedArgs.count("x"))
+    x = us::any_cast<bool>(parsedArgs["x"]);
+
+  bool y = false;
+  if (parsedArgs.count("y"))
+    y = us::any_cast<bool>(parsedArgs["y"]);
+
+  bool z = false;
+  if (parsedArgs.count("z"))
+    z = us::any_cast<bool>(parsedArgs["z"]);
+
+  try
+  {
+    mitk::PeakImage::Pointer image = dynamic_cast<mitk::PeakImage*>(mitk::IOUtil::Load(imageName)[0].GetPointer());
+
+    typedef mitk::ImageToItk< mitk::PeakImage::ItkPeakImageType > CasterType;
+    CasterType::Pointer caster = CasterType::New();
+    caster->SetInput(image);
+    caster->Update();
+    mitk::PeakImage::ItkPeakImageType::Pointer itkImg = caster->GetOutput();
+
+    itk::FlipPeaksFilter< float >::Pointer flipper = itk::FlipPeaksFilter< float >::New();
+    flipper->SetInput(itkImg);
+    flipper->SetFlipX(x);
+    flipper->SetFlipY(y);
+    flipper->SetFlipZ(z);
+    flipper->Update();
+
+    mitk::Image::Pointer resultImage = dynamic_cast<mitk::Image*>(mitk::PeakImage::New().GetPointer());
+    mitk::CastToMitkImage(flipper->GetOutput(), resultImage);
+    resultImage->SetVolume(flipper->GetOutput()->GetBufferPointer());
+
+    mitk::IOUtil::Save(resultImage, outImage);
+  }
+  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/AnchorBasedScoring.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/AnchorBasedScoring.cpp
index fb08c52ef5..3f9208ccc5 100755
--- a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/AnchorBasedScoring.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/AnchorBasedScoring.cpp
@@ -1,443 +1,443 @@
 /*===================================================================
 
 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 <mitkBaseData.h>
 #include <mitkImageCast.h>
 #include <mitkImageToItk.h>
 #include <metaCommand.h>
 #include <mitkCommandLineParser.h>
 #include <usAny.h>
 #include <mitkIOUtil.h>
 #include <itksys/SystemTools.hxx>
 #include <itkDirectory.h>
 #include <mitkFiberBundle.h>
 #include <mitkPreferenceListReaderOptionsFunctor.h>
 #include <itkImageFileWriter.h>
 #include <mitkPeakImage.h>
 #include <itkFitFibersToImageFilter.h>
 #include <boost/lexical_cast.hpp>
 #include <itkTractDensityImageFilter.h>
 
 typedef itksys::SystemTools ist;
 typedef itk::Point<float, 4> PointType4;
 typedef itk::Image< float, 4 >  PeakImgType;
 typedef itk::Image< unsigned char, 3 >  ItkUcharImageType;
 
 std::vector< std::string > get_file_list(const std::string& path, std::vector< std::string > extensions={".fib", ".trk"})
 {
   std::vector< std::string > file_list;
   itk::Directory::Pointer dir = itk::Directory::New();
 
   if (dir->Load(path.c_str()))
   {
     int n = dir->GetNumberOfFiles();
     for (int r = 0; r < n; r++)
     {
       const char *filename = dir->GetFile(r);
       std::string ext = ist::GetFilenameExtension(filename);
       for (auto e : extensions)
       {
         if (ext==e)
         {
           file_list.push_back(path + '/' + filename);
           break;
         }
       }
     }
   }
   return file_list;
 }
 
 std::vector< mitk::FiberBundle::Pointer > CombineTractograms(std::vector< mitk::FiberBundle::Pointer > reference, std::vector< mitk::FiberBundle::Pointer > candidates, int skip=-1)
 {
   std::vector< mitk::FiberBundle::Pointer > fib;
   for (auto f : reference)
     fib.push_back(f);
 
   int c = 0;
   for (auto f : candidates)
   {
     if (c!=skip)
       fib.push_back(f);
     ++c;
   }
 
   return fib;
 }
 
 /*!
 \brief Fits the tractogram to the input peak image by assigning a weight to each fiber (similar to https://doi.org/10.1016/j.neuroimage.2015.06.092).
 */
 int main(int argc, char* argv[])
 {
   mitkCommandLineParser parser;
 
   parser.setTitle("Anchor Based Scoring");
   parser.setCategory("Fiber Tracking Evaluation");
   parser.setDescription("");
   parser.setContributor("MIC");
 
   parser.setArgumentPrefix("--", "-");
   parser.addArgument("", "i1", mitkCommandLineParser::InputFile, "Anchor tractogram:", "anchor tracts in one tractogram file", us::Any(), false);
   parser.addArgument("", "i2", mitkCommandLineParser::InputFile, "Input peaks:", "input peak image", us::Any(), false);
   parser.addArgument("", "i3", mitkCommandLineParser::InputFile, "Candidate folder:", "folder containing all candidate tracts (separate files)", us::Any(), false);
   parser.addArgument("", "o", mitkCommandLineParser::OutputDirectory, "Output folder:", "output folder", us::Any(), false);
 
   parser.addArgument("reference_mask_folder", "", mitkCommandLineParser::String, "Reference Mask Folder:", "reference tract masks for accuracy evaluation");
   parser.addArgument("mask", "", mitkCommandLineParser::InputFile, "Mask image:", "scoring is only performed inside the mask image");
   parser.addArgument("greedy_add", "", mitkCommandLineParser::Bool, "Greedy:", "if enabled, the candidate tracts are not jointly fitted to the residual image but one after the other employing a greedy scheme", false);
   parser.addArgument("lambda", "", mitkCommandLineParser::Float, "Lambda:", "modifier for regularization", 0.1);
   parser.addArgument("dont_filter_outliers", "", mitkCommandLineParser::Bool, "Don't filter outliers:", "don't perform second optimization run with an upper weight bound based on the first weight estimation (95% quantile)", false);
 
   std::map<std::string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
   if (parsedArgs.size()==0)
     return EXIT_FAILURE;
 
   std::string anchors_file = us::any_cast<std::string>(parsedArgs["i1"]);
   std::string peak_file_name = us::any_cast<std::string>(parsedArgs["i2"]);
   std::string candidate_folder = us::any_cast<std::string>(parsedArgs["i3"]);
   std::string out_folder = us::any_cast<std::string>(parsedArgs["o"]);
 
   bool greedy_add = false;
   if (parsedArgs.count("greedy_add"))
     greedy_add = us::any_cast<bool>(parsedArgs["greedy_add"]);
 
   float lambda = 0.1;
   if (parsedArgs.count("lambda"))
     lambda = us::any_cast<float>(parsedArgs["lambda"]);
 
   bool filter_outliers = true;
   if (parsedArgs.count("dont_filter_outliers"))
     filter_outliers = !us::any_cast<bool>(parsedArgs["dont_filter_outliers"]);
 
   std::string mask_file = "";
   if (parsedArgs.count("mask"))
     mask_file = us::any_cast<std::string>(parsedArgs["mask"]);
 
   std::string reference_mask_folder = "";
   if (parsedArgs.count("reference_mask_folder"))
     reference_mask_folder = us::any_cast<std::string>(parsedArgs["reference_mask_folder"]);
 
   try
   {
     itk::TimeProbe clock;
     clock.Start();
 
     MITK_INFO << "Creating directory structure";
     if (ist::PathExists(out_folder))
       ist::RemoveADirectory(out_folder);
     ist::MakeDirectory(out_folder);
 
     MITK_INFO << "Loading data";
 
     std::streambuf *old = cout.rdbuf(); // <-- save
     std::stringstream ss;
     std::cout.rdbuf (ss.rdbuf());       // <-- redirect
 
     ofstream logfile;
     logfile.open (out_folder + "log.txt");
 
     itk::ImageFileWriter< PeakImgType >::Pointer peak_image_writer = itk::ImageFileWriter< PeakImgType >::New();
     mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Peak Image", "Fiberbundles"}, {});
     mitk::Image::Pointer inputImage = dynamic_cast<mitk::PeakImage*>(mitk::IOUtil::Load(peak_file_name, &functor)[0].GetPointer());
 
     float minSpacing = 1;
     if(inputImage->GetGeometry()->GetSpacing()[0]<inputImage->GetGeometry()->GetSpacing()[1] && inputImage->GetGeometry()->GetSpacing()[0]<inputImage->GetGeometry()->GetSpacing()[2])
       minSpacing = inputImage->GetGeometry()->GetSpacing()[0];
     else if (inputImage->GetGeometry()->GetSpacing()[1] < inputImage->GetGeometry()->GetSpacing()[2])
       minSpacing = inputImage->GetGeometry()->GetSpacing()[1];
     else
       minSpacing = inputImage->GetGeometry()->GetSpacing()[2];
 
     // Load mask file. Fit is only performed inside the mask
     itk::FitFibersToImageFilter::UcharImgType::Pointer mask = nullptr;
     if (mask_file.compare("")!=0)
     {
       mitk::Image::Pointer mitk_mask = dynamic_cast<mitk::Image*>(mitk::IOUtil::Load(mask_file)[0].GetPointer());
       mitk::CastToItkImage(mitk_mask, mask);
     }
 
     // Load masks covering the true positives for evaluation purposes
     std::vector< std::string > reference_mask_filenames;
     std::vector< itk::FitFibersToImageFilter::UcharImgType::Pointer > reference_masks;
     if (reference_mask_folder.compare("")!=0)
     {
       reference_mask_filenames = get_file_list(reference_mask_folder, {".nii", ".nii.gz", ".nrrd"});
       for (auto filename : reference_mask_filenames)
       {
         itk::FitFibersToImageFilter::UcharImgType::Pointer ref_mask = nullptr;
         mitk::Image::Pointer ref_mitk_mask = dynamic_cast<mitk::Image*>(mitk::IOUtil::Load(filename)[0].GetPointer());
         mitk::CastToItkImage(ref_mitk_mask, ref_mask);
         reference_masks.push_back(ref_mask);
       }
     }
 
     // Load peak image
     typedef mitk::ImageToItk< PeakImgType > CasterType;
     CasterType::Pointer caster = CasterType::New();
     caster->SetInput(inputImage);
     caster->Update();
     PeakImgType::Pointer peak_image = caster->GetOutput();
 
     // Load all candidate tracts
     std::vector< mitk::FiberBundle::Pointer > input_candidates;
     std::vector< std::string > candidate_tract_files = get_file_list(candidate_folder);
     for (std::string f : candidate_tract_files)
     {
       mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(mitk::IOUtil::Load(f)[0].GetPointer());
       if (fib.IsNull())
         continue;
       if (fib->GetNumFibers()<=0)
         return 0;
       fib->ResampleLinear(minSpacing/10.0);
       input_candidates.push_back(fib);
     }
     std::cout.rdbuf (old);              // <-- restore
 
     double rmse = 0.0;
     int iteration = 0;
     std::string name = "NOANCHOR";
     // Load reference tractogram consisting of all known tracts
     std::vector< mitk::FiberBundle::Pointer > input_reference;
     mitk::FiberBundle::Pointer anchor_tractogram = dynamic_cast<mitk::FiberBundle*>(mitk::IOUtil::Load(anchors_file)[0].GetPointer());
     if ( !(anchor_tractogram.IsNull() || anchor_tractogram->GetNumFibers()==0) )
     {
       std::streambuf *old = cout.rdbuf(); // <-- save
       std::stringstream ss;
       std::cout.rdbuf (ss.rdbuf());       // <-- redirect
       anchor_tractogram->ResampleLinear(minSpacing/10.0);
       std::cout.rdbuf (old);              // <-- restore
       input_reference.push_back(anchor_tractogram);
 
       // Fit known tracts to peak image to obtain underexplained image
       MITK_INFO << "Fit anchor tracts";
       itk::FitFibersToImageFilter::Pointer fitter = itk::FitFibersToImageFilter::New();
       fitter->SetTractograms(input_reference);
       fitter->SetLambda(lambda);
       fitter->SetFilterOutliers(filter_outliers);
       fitter->SetPeakImage(peak_image);
       fitter->SetVerbose(true);
       fitter->SetResampleFibers(false);
       fitter->SetMaskImage(mask);
       fitter->Update();
       rmse = fitter->GetRMSE();
 
       name = ist::GetFilenameWithoutExtension(anchors_file);
       mitk::FiberBundle::Pointer anchor_tracts = fitter->GetTractograms().at(0);
       anchor_tracts->SetFiberColors(255,255,255);
       mitk::IOUtil::Save(anchor_tracts, out_folder + "0_" + name + ".fib");
 
       peak_image = fitter->GetUnderexplainedImage();
       peak_image_writer->SetInput(peak_image);
       peak_image_writer->SetFileName(out_folder + boost::lexical_cast<std::string>(iteration) + "_" + name + ".nrrd");
       peak_image_writer->Update();
     }
 
     if (!greedy_add)
     {
       MITK_INFO << "Fit candidate tracts";
       itk::FitFibersToImageFilter::Pointer fitter = itk::FitFibersToImageFilter::New();
       fitter->SetLambda(lambda);
       fitter->SetFilterOutliers(filter_outliers);
       fitter->SetVerbose(true);
       fitter->SetPeakImage(peak_image);
       fitter->SetResampleFibers(false);
       fitter->SetMaskImage(mask);
       fitter->SetTractograms(input_candidates);
       fitter->Update();
       vnl_vector<double> rms_diff = fitter->GetRmsDiffPerBundle();
 
       vnl_vector<double> log_rms_diff = rms_diff-rms_diff.min_value() + 1;
       log_rms_diff = log_rms_diff.apply(std::log);
       log_rms_diff /= log_rms_diff.max_value();
       int c = 0;
       for (auto fib : input_candidates)
       {
         fib->SetFiberWeights( log_rms_diff[c] );
         fib->ColorFibersByFiberWeights(false, true);
 
         std::string bundle_name = ist::GetFilenameWithoutExtension(candidate_tract_files.at(c));
 
         std::streambuf *old = cout.rdbuf(); // <-- save
         std::stringstream ss;
         std::cout.rdbuf (ss.rdbuf());       // <-- redirect
-        mitk::IOUtil::Save(fib, out_folder + boost::lexical_cast<std::string>((int)(10000*rms_diff[c])) + "_" + bundle_name + ".fib");
+        mitk::IOUtil::Save(fib, out_folder + boost::lexical_cast<std::string>((int)(100000*rms_diff[c])) + "_" + bundle_name + ".fib");
 
         float best_overlap = 0;
         int best_overlap_index = -1;
         int m_idx = 0;
         for (auto ref_mask : reference_masks)
         {
           float overlap = fib->GetOverlap(ref_mask, false);
           if (overlap>best_overlap)
           {
             best_overlap = overlap;
             best_overlap_index = m_idx;
           }
           ++m_idx;
         }
 
         unsigned int num_voxels = 0;
         {
           itk::TractDensityImageFilter< ItkUcharImageType >::Pointer masks_filter = itk::TractDensityImageFilter< ItkUcharImageType >::New();
           masks_filter->SetInputImage(mask);
           masks_filter->SetBinaryOutput(true);
           masks_filter->SetFiberBundle(fib);
           masks_filter->SetUseImageGeometry(true);
           masks_filter->Update();
           num_voxels = masks_filter->GetNumCoveredVoxels();
         }
 
         std::cout.rdbuf (old);              // <-- restore
 
         logfile << "RMS_DIFF: " << setprecision(5) << rms_diff[c] << " " << bundle_name << " " << num_voxels << "\n";
         if (best_overlap_index>=0)
           logfile << "Best_overlap: " << setprecision(5) << best_overlap << " " << ist::GetFilenameWithoutExtension(reference_mask_filenames.at(best_overlap_index)) << "\n";
         else
           logfile << "No_overlap\n";
 
         ++c;
       }
 
       mitk::FiberBundle::Pointer out_fib = mitk::FiberBundle::New();
       out_fib = out_fib->AddBundles(input_candidates);
       out_fib->ColorFibersByFiberWeights(false, true);
       mitk::IOUtil::Save(out_fib, out_folder + "AllCandidates.fib");
     }
     else
     {
       MITK_INFO << "RMSE: " << setprecision(5) << rmse;
       //    fitter->SetPeakImage(peak_image);
 
       // Iteratively add candidate bundles in a greedy manner
       while (!input_candidates.empty())
       {
         double next_rmse = rmse;
         double num_peaks = 0;
         mitk::FiberBundle::Pointer best_candidate = nullptr;
         PeakImgType::Pointer best_candidate_peak_image = nullptr;
 
         for (int i=0; i<(int)input_candidates.size(); ++i)
         {
           // WHY NECESSARY AGAIN??
           itk::FitFibersToImageFilter::Pointer fitter = itk::FitFibersToImageFilter::New();
           fitter->SetLambda(lambda);
           fitter->SetFilterOutliers(filter_outliers);
           fitter->SetVerbose(false);
           fitter->SetPeakImage(peak_image);
           fitter->SetResampleFibers(false);
           fitter->SetMaskImage(mask);
           // ******************************
           fitter->SetTractograms({input_candidates.at(i)});
 
           std::streambuf *old = cout.rdbuf(); // <-- save
           std::stringstream ss;
           std::cout.rdbuf (ss.rdbuf());       // <-- redirect
           fitter->Update();
           std::cout.rdbuf (old);              // <-- restore
 
           double candidate_rmse = fitter->GetRMSE();
 
           if (candidate_rmse<next_rmse)
           {
             next_rmse = candidate_rmse;
             num_peaks = fitter->GetNumCoveredDirections();
             best_candidate = fitter->GetTractograms().at(0);
             best_candidate_peak_image = fitter->GetUnderexplainedImage();
           }
         }
 
         if (best_candidate.IsNull())
           break;
 
         //      fitter->SetPeakImage(peak_image);
         peak_image = best_candidate_peak_image;
 
         int i=0;
         std::vector< mitk::FiberBundle::Pointer > remaining_candidates;
         std::vector< std::string > remaining_candidate_files;
         for (auto fib : input_candidates)
         {
           if (fib!=best_candidate)
           {
             remaining_candidates.push_back(fib);
             remaining_candidate_files.push_back(candidate_tract_files.at(i));
           }
           else
             name = ist::GetFilenameWithoutExtension(candidate_tract_files.at(i));
           ++i;
         }
         input_candidates = remaining_candidates;
         candidate_tract_files = remaining_candidate_files;
 
         iteration++;
         std::streambuf *old = cout.rdbuf(); // <-- save
         std::stringstream ss;
         std::cout.rdbuf (ss.rdbuf());       // <-- redirect
         // Save winning candidate
         mitk::IOUtil::Save(best_candidate, out_folder + boost::lexical_cast<std::string>(iteration) + "_" + name + ".fib");
         peak_image_writer->SetInput(peak_image);
         peak_image_writer->SetFileName(out_folder + boost::lexical_cast<std::string>(iteration) + "_" + name + ".nrrd");
         peak_image_writer->Update();
 
         // Calculate best overlap with reference masks for evaluation purposes
         float best_overlap = 0;
         int best_overlap_index = -1;
         i = 0;
         for (auto ref_mask : reference_masks)
         {
           float overlap = best_candidate->GetOverlap(ref_mask, false);
           if (overlap>best_overlap)
           {
             best_overlap = overlap;
             best_overlap_index = i;
           }
           ++i;
         }
         std::cout.rdbuf (old);              // <-- restore
 
         logfile << "RMSE: " << setprecision(5) << rmse << " " << name << " " << num_peaks << "\n";
         if (best_overlap_index>=0)
           logfile << "Best_overlap: " << setprecision(5) << best_overlap << " " << ist::GetFilenameWithoutExtension(reference_mask_filenames.at(best_overlap_index)) << "\n";
         else
           logfile << "No_overlap\n";
       }
     }
 
     clock.Stop();
     int h = clock.GetTotal()/3600;
     int m = ((int)clock.GetTotal()%3600)/60;
     int s = (int)clock.GetTotal()%60;
     MITK_INFO << "Plausibility estimation took " << h << "h, " << m << "m and " << s << "s";
     logfile.close();
   }
   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/CMakeLists.txt b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/CMakeLists.txt
index a72860b7ab..85af60ba39 100755
--- a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/CMakeLists.txt
+++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/CMakeLists.txt
@@ -1,44 +1,45 @@
 option(BUILD_DiffusionTractographyEvaluationCmdApps "Build commandline tools for diffusion fiber tractography evaluation" OFF)
 
 if(BUILD_DiffusionTractographyEvaluationCmdApps 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( diffusionTractographyEvaluationcmdapps
     PeaksAngularError^^MitkFiberTracking
     TractometerMetrics^^MitkFiberTracking
     MergeOverlappingTracts^^MitkFiberTracking
-    ExtractAnchorTracts^^MitkFiberTracking
+    ExtractOverlappingTracts^^MitkFiberTracking
+    ExtractSimilarTracts^^MitkFiberTracking
     AnchorBasedScoring^^MitkFiberTracking
     EvaluateLiFE^^MitkFiberTracking
     # LocalDirectionalFiberPlausibility^^MitkFiberTracking # HAS TO USE NEW PEAK IMAGE FORMAT
     )
 
     foreach(diffusionTractographyEvaluationcmdapp ${diffusionTractographyEvaluationcmdapps})
       # extract cmd app name and dependencies
       string(REPLACE "^^" "\\;" cmdapp_info ${diffusionTractographyEvaluationcmdapp})
       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/TractographyEvaluation/ExtractAnchorTracts.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/ExtractOverlappingTracts.cpp
similarity index 99%
rename from Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/ExtractAnchorTracts.cpp
rename to Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/ExtractOverlappingTracts.cpp
index a38418adeb..70ecc3ac90 100755
--- a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/ExtractAnchorTracts.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/ExtractOverlappingTracts.cpp
@@ -1,268 +1,268 @@
 /*===================================================================
 
 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 <mitkBaseData.h>
 #include <mitkImageCast.h>
 #include <mitkImageToItk.h>
 #include <metaCommand.h>
 #include <mitkCommandLineParser.h>
 #include <usAny.h>
 #include <mitkIOUtil.h>
 #include <boost/lexical_cast.hpp>
 #include <itksys/SystemTools.hxx>
 #include <itkDirectory.h>
 #include <mitkFiberBundle.h>
 #include <vtkTransformPolyDataFilter.h>
 #include <fstream>
 #include <chrono>
 #include <boost/progress.hpp>
 
 typedef itksys::SystemTools ist;
 typedef itk::Image<unsigned char, 3>    ItkUcharImgType;
 typedef std::tuple< ItkUcharImgType::Pointer, std::string > MaskType;
 
 void CreateFolderStructure(const std::string& path)
 {
   if (ist::PathExists(path))
     ist::RemoveADirectory(path);
 
   ist::MakeDirectory(path);
   ist::MakeDirectory(path + "/anchor_tracts/");
   ist::MakeDirectory(path + "/candidate_tracts/");
   ist::MakeDirectory(path + "/implausible_tracts/");
   ist::MakeDirectory(path + "/skipped_masks/");
 }
 
 ItkUcharImgType::Pointer LoadItkMaskImage(const std::string& filename)
 {
   mitk::Image::Pointer img = dynamic_cast<mitk::Image*>(mitk::IOUtil::Load(filename)[0].GetPointer());
   ItkUcharImgType::Pointer itkMask = ItkUcharImgType::New();
   mitk::CastToItkImage(img, itkMask);
   return itkMask;
 }
 
 std::vector< MaskType > get_file_list(const std::string& path, float anchor_fraction, const std::string& skipped_path)
 {
   if (anchor_fraction<0)
     anchor_fraction = 0;
   else if (anchor_fraction>1.0)
     anchor_fraction = 1.0;
 
   std::chrono::milliseconds ms = std::chrono::duration_cast< std::chrono::milliseconds >(std::chrono::system_clock::now().time_since_epoch());
   std::srand(ms.count());
 
   std::vector< MaskType > mask_list;
 
   itk::Directory::Pointer dir = itk::Directory::New();
 
   int skipped = 0;
   if (dir->Load(path.c_str()))
   {
     int n = dir->GetNumberOfFiles();
     int num_images = 0;
 
     std::vector< int > im_indices;
     for (int r = 0; r < n; r++)
     {
       const char *filename = dir->GetFile(r);
       std::string ext = ist::GetFilenameExtension(filename);
       if (ext==".nii" || ext==".nii.gz" || ext==".nrrd")
       {
         ++num_images;
         im_indices.push_back(r);
       }
     }
 
     int skipping_num = num_images * (1.0 - anchor_fraction);
     std::random_shuffle(im_indices.begin(), im_indices.end());
     MITK_INFO << "Skipping " << skipping_num << " images";
     MITK_INFO << "Number of anchors: " << num_images-skipping_num;
 
     int c = -1;
     for (int r : im_indices)
     {
       c++;
       const char *filename = dir->GetFile(r);
 
       if (c<skipping_num && (skipping_num==num_images || ist::GetFilenameWithoutExtension(filename)!="CC"))
       {
         MITK_INFO << "Skipping " << ist::GetFilenameWithoutExtension(filename);
         ist::CopyAFile(path + '/' + filename, skipped_path + ist::GetFilenameName(filename));
         skipped++;
         continue;
       }
 
       MITK_INFO << "Loading " << ist::GetFilenameWithoutExtension(filename);
       std::streambuf *old = cout.rdbuf(); // <-- save
       std::stringstream ss;
       std::cout.rdbuf (ss.rdbuf());       // <-- redirect
       MaskType m(LoadItkMaskImage(path + '/' + filename), ist::GetFilenameName(filename));
       mask_list.push_back(m);
       std::cout.rdbuf (old);              // <-- restore
     }
   }
 
   return mask_list;
 }
 
 mitk::FiberBundle::Pointer TransformToMRtrixSpace(mitk::FiberBundle::Pointer fib)
 {
   mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New();
   vtkSmartPointer< vtkMatrix4x4 > matrix = vtkSmartPointer< vtkMatrix4x4 >::New();
   matrix->Identity();
   matrix->SetElement(0,0,-matrix->GetElement(0,0));
   matrix->SetElement(1,1,-matrix->GetElement(1,1));
   geometry->SetIndexToWorldTransformByVtkMatrix(matrix);
 
   vtkSmartPointer<vtkTransformPolyDataFilter> transformFilter = vtkSmartPointer<vtkTransformPolyDataFilter>::New();
   transformFilter->SetInputData(fib->GetFiberPolyData());
   transformFilter->SetTransform(geometry->GetVtkTransform());
   transformFilter->Update();
   mitk::FiberBundle::Pointer transformed_fib = mitk::FiberBundle::New(transformFilter->GetOutput());
   return transformed_fib;
 }
 
 /*!
 \brief
 */
 int main(int argc, char* argv[])
 {
   mitkCommandLineParser parser;
 
-  parser.setTitle("Extract Anchor Tracts");
+  parser.setTitle("Extract Overlapping Tracts");
   parser.setCategory("Fiber Tracking Evaluation");
   parser.setDescription("");
   parser.setContributor("MIC");
 
   parser.setArgumentPrefix("--", "-");
   parser.addArgument("input", "i", mitkCommandLineParser::InputFile, "Input:", "input tractogram (.fib, vtk ascii file format)", us::Any(), false);
   parser.addArgument("out", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output folder", us::Any(), false);
   parser.addArgument("reference_mask_folder", "m", mitkCommandLineParser::String, "Reference Mask Folder:", "reference masks of known bundles", us::Any(), false);
   parser.addArgument("gray_matter_mask", "gm", mitkCommandLineParser::String, "GM mask:", "remove fibers not ending in the gray matter");
   parser.addArgument("anchor_fraction", "", mitkCommandLineParser::Float, "Anchor fraction:", "Fraction of tracts used as anchors", 0.5);
   parser.addArgument("overlap", "", mitkCommandLineParser::Float, "Overlap threshold:", "Overlap threshold used to identify the anchor tracts", 0.8);
   parser.addArgument("subsample", "", mitkCommandLineParser::Float, "Subsampling factor:", "Only use specified fraction of input fibers for the analysis", 1.0);
 
   std::map<std::string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
   if (parsedArgs.size()==0)
     return EXIT_FAILURE;
 
   std::string fibFile = us::any_cast<std::string>(parsedArgs["input"]);
   std::string reference_mask_folder = us::any_cast<std::string>(parsedArgs["reference_mask_folder"]);
   std::string out_folder = us::any_cast<std::string>(parsedArgs["out"]);
 
   std::string gray_matter_mask = "";
   if (parsedArgs.count("gray_matter_mask"))
     gray_matter_mask = us::any_cast<std::string>(parsedArgs["gray_matter_mask"]);
 
   float anchor_fraction = 0.5;
   if (parsedArgs.count("anchor_fraction"))
     anchor_fraction = us::any_cast<float>(parsedArgs["anchor_fraction"]);
 
   float overlap = 0.8;
   if (parsedArgs.count("overlap"))
     overlap = us::any_cast<float>(parsedArgs["overlap"]);
 
   float subsample = 1.0;
   if (parsedArgs.count("subsample"))
     subsample = us::any_cast<float>(parsedArgs["subsample"]);
 
   try
   {
     CreateFolderStructure(out_folder);
     std::vector< MaskType > known_tract_masks = get_file_list(reference_mask_folder, anchor_fraction, out_folder + "/skipped_masks/");
     mitk::FiberBundle::Pointer inputTractogram = dynamic_cast<mitk::FiberBundle*>(mitk::IOUtil::Load(fibFile)[0].GetPointer());
 
     MITK_INFO << "Removing fibers not ending inside of GM";
     if (gray_matter_mask.compare("")!=0)
     {
       std::streambuf *old = cout.rdbuf(); // <-- save
       std::stringstream ss;
       std::cout.rdbuf (ss.rdbuf());       // <-- redirect
       ItkUcharImgType::Pointer gm_image = LoadItkMaskImage(gray_matter_mask);
       std::cout.rdbuf (old);              // <-- restore
 
       mitk::FiberBundle::Pointer not_gm_fibers = inputTractogram->ExtractFiberSubset(gm_image, false, true, true);
 
       old = cout.rdbuf(); // <-- save
       std::cout.rdbuf (ss.rdbuf());       // <-- redirect
       mitk::IOUtil::Save(not_gm_fibers, out_folder + "/implausible_tracts/no_gm_endings.trk");
       inputTractogram = inputTractogram->ExtractFiberSubset(gm_image, false, false, true);
       std::cout.rdbuf (old);              // <-- restore
     }
 
     if (subsample<1.0)
       inputTractogram = inputTractogram->SubsampleFibers(subsample);
 
 
     // resample fibers
     float minSpacing = 1;
     if (!known_tract_masks.empty())
     {
       if(std::get<0>(known_tract_masks.at(0))->GetSpacing()[0]<std::get<0>(known_tract_masks.at(0))->GetSpacing()[1] && std::get<0>(known_tract_masks.at(0))->GetSpacing()[0]<std::get<0>(known_tract_masks.at(0))->GetSpacing()[2])
         minSpacing = std::get<0>(known_tract_masks.at(0))->GetSpacing()[0];
       else if (std::get<0>(known_tract_masks.at(0))->GetSpacing()[1] < std::get<0>(known_tract_masks.at(0))->GetSpacing()[2])
         minSpacing = std::get<0>(known_tract_masks.at(0))->GetSpacing()[1];
       else
         minSpacing = std::get<0>(known_tract_masks.at(0))->GetSpacing()[2];
     }
     inputTractogram->ResampleLinear(minSpacing/5);
 
     mitk::FiberBundle::Pointer all_anchor_tracts = mitk::FiberBundle::New(nullptr);
     if (!known_tract_masks.empty())
     {
       MITK_INFO << "Find known tracts via overlap match";
       boost::progress_display disp(known_tract_masks.size());
       for ( MaskType mask : known_tract_masks )
       {
         ++disp;
         std::streambuf *old = cout.rdbuf(); // <-- save
         std::stringstream ss;
         std::cout.rdbuf (ss.rdbuf());       // <-- redirect
 
         ItkUcharImgType::Pointer mask_image = std::get<0>(mask);
         std::string mask_name = std::get<1>(mask);
         mitk::FiberBundle::Pointer known_tract = inputTractogram->ExtractFiberSubset(mask_image, true, false, false, overlap, false);
         mitk::IOUtil::Save(known_tract, out_folder + "/anchor_tracts/" + mask_name + ".trk");
 
         all_anchor_tracts = all_anchor_tracts->AddBundle(known_tract);
 
         std::cout.rdbuf (old);              // <-- restore
       }
     }
 
     mitk::IOUtil::Save(all_anchor_tracts, out_folder + "/anchor_tracts/anchor_tractogram.trk");
     mitk::FiberBundle::Pointer remaining_tracts = inputTractogram->SubtractBundle(all_anchor_tracts);
     mitk::IOUtil::Save(remaining_tracts, out_folder + "/candidate_tracts/candidate_tractogram.trk");
     mitk::IOUtil::Save(inputTractogram, out_folder + "/filtered_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/cmdapps/TractographyEvaluation/ExtractSimilarTracts.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/ExtractSimilarTracts.cpp
new file mode 100644
index 0000000000..74caed4caa
--- /dev/null
+++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/ExtractSimilarTracts.cpp
@@ -0,0 +1,211 @@
+/*===================================================================
+
+The Medical Imaging Interaction Toolkit (MITK)
+
+Copyright (c) German Cancer Research Center,
+Division of Medical and Biological Informatics.
+All rights reserved.
+
+This software is distributed WITHOUT ANY WARRANTY; without
+even the implied warranty of MERCHANTABILITY or FITNESS FOR
+A PARTICULAR PURPOSE.
+
+See LICENSE.txt or http://www.mitk.org for details.
+
+===================================================================*/
+
+#include <mitkFiberBundle.h>
+#include <mitkCommandLineParser.h>
+#include <boost/lexical_cast.hpp>
+#include <mitkIOUtil.h>
+#include <itkTractClusteringFilter.h>
+
+typedef itksys::SystemTools ist;
+typedef itk::Image<unsigned char, 3>    ItkUcharImgType;
+
+mitk::FiberBundle::Pointer LoadFib(std::string filename)
+{
+  std::vector<mitk::BaseData::Pointer> 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<mitk::FiberBundle*>(baseData.GetPointer());
+}
+
+ItkUcharImgType::Pointer LoadItkMaskImage(const std::string& filename)
+{
+  mitk::Image::Pointer img = dynamic_cast<mitk::Image*>(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<std::string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
+  if (parsedArgs.size()==0)
+    return EXIT_FAILURE;
+
+  std::string in_fib = us::any_cast<std::string>(parsedArgs["i"]);
+  std::string out_root = us::any_cast<std::string>(parsedArgs["o"]);
+  mitkCommandLineParser::StringContainerType ref_bundle_files = us::any_cast<mitkCommandLineParser::StringContainerType>(parsedArgs["ref_tracts"]);
+  mitkCommandLineParser::StringContainerType ref_mask_files;
+  if (parsedArgs.count("ref_masks"))
+    ref_mask_files = us::any_cast<mitkCommandLineParser::StringContainerType>(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<int>(parsedArgs["distance"]);
+
+  std::string metric = "MDF_STD";
+  if (parsedArgs.count("metric"))
+    metric = us::any_cast<std::string>(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<ref_bundle_files.size(); ++i)
+    {
+      MITK_INFO << "Loading " << ist::GetFilenameName(ref_bundle_files.at(i));
+      if (i<ref_mask_files.size())
+        MITK_INFO << "Loading " << ist::GetFilenameName(ref_mask_files.at(i));
+
+      std::streambuf *old = cout.rdbuf(); // <-- save
+      std::stringstream ss;
+      std::cout.rdbuf (ss.rdbuf());       // <-- redirect
+      try
+      {
+        ref_fibs.push_back(LoadFib(ref_bundle_files.at(i)));
+        if (i<ref_mask_files.size())
+          ref_masks.push_back(LoadItkMaskImage(ref_mask_files.at(i)));
+        else
+          ref_masks.push_back(nullptr);
+        std::cout.rdbuf (old);              // <-- restore
+      }
+      catch(...){
+        std::cout.rdbuf (old);              // <-- restore
+        std::cout << "could not load: " << ref_bundle_files.at(i);
+        return EXIT_FAILURE;
+      }
+    }
+
+    std::vector< float > 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->Update();
+          std::vector<mitk::FiberBundle::Pointer> 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);
+        segmenter->Update();
+
+        std::vector<mitk::FiberBundle::Pointer> 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;
+}