diff --git a/Modules/FiberTracking/Algorithms/itkTractParcellationFilter.cpp b/Modules/FiberTracking/Algorithms/itkTractParcellationFilter.cpp
index 0b11005..0513246 100644
--- a/Modules/FiberTracking/Algorithms/itkTractParcellationFilter.cpp
+++ b/Modules/FiberTracking/Algorithms/itkTractParcellationFilter.cpp
@@ -1,468 +1,563 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center.
 
 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 "itkTractParcellationFilter.h"
 
 // VTK
 #include <vtkPolyLine.h>
 #include <vtkCellArray.h>
 #include <vtkCellData.h>
 #include <boost/progress.hpp>
 #include <mitkDiffusionFunctionCollection.h>
 #include <itkTractDensityImageFilter.h>
 #include <itkImageRegionIterator.h>
 #include <mitkTractClusteringFilter.h>
 #include <mitkClusteringMetricEuclideanStd.h>
 #include <mitkClusteringMetricEuclideanMean.h>
 #include <mitkClusteringMetricLength.h>
 #include <itkImageRegionConstIterator.h>
 #include <itkTractsToVectorImageFilter.h>
 #include <mitkTractometry.h>
 
 namespace itk{
 
 template< class OutImageType, class InputImageType >
 TractParcellationFilter< OutImageType, InputImageType >::TractParcellationFilter()
   : m_UpsamplingFactor(1)
   , m_NumParcels(0)
   , m_NumCentroids(0)
   , m_StartClusterSize(5)
   , m_InputImage(nullptr)
 {
   this->SetNumberOfRequiredOutputs(2);
 }
 
 template< class OutImageType, class InputImageType >
 TractParcellationFilter< OutImageType, InputImageType >::~TractParcellationFilter()
 {
 }
 
 template< class OutImageType, class InputImageType >
 bool TractParcellationFilter< OutImageType, InputImageType >::Flip(vtkSmartPointer< vtkPolyData > polydata1, int i, vtkSmartPointer< vtkPolyData > ref_poly, int ref_i)
 {
   double d_direct = 0;
   double d_flipped = 0;
 
   vtkCell* cell1 = polydata1->GetCell(ref_i);
   if (ref_poly!=nullptr)
     cell1 = ref_poly->GetCell(ref_i);
   auto numPoints1 = cell1->GetNumberOfPoints();
   vtkPoints* points1 = cell1->GetPoints();
 
   std::vector<itk::Point<double, 3>> ref_points;
   for (int j=0; j<numPoints1; ++j)
   {
     double* p1 = points1->GetPoint(j);
     itk::Point<double, 3> itk_p;
     itk_p[0] = p1[0];
     itk_p[1] = p1[1];
     itk_p[2] = p1[2];
     ref_points.push_back(itk_p);
   }
 
   vtkCell* cell2 = polydata1->GetCell(i);
   vtkPoints* points2 = cell2->GetPoints();
 
   for (int j=0; j<numPoints1; ++j)
   {
     auto p1 = ref_points.at(j);
 
     double* p2 = points2->GetPoint(j);
     d_direct += (p1[0]-p2[0])*(p1[0]-p2[0]) + (p1[1]-p2[1])*(p1[1]-p2[1]) + (p1[2]-p2[2])*(p1[2]-p2[2]);
 
     double* p3 = points2->GetPoint(numPoints1-j-1);
     d_flipped += (p1[0]-p3[0])*(p1[0]-p3[0]) + (p1[1]-p3[1])*(p1[1]-p3[1]) + (p1[2]-p3[2])*(p1[2]-p3[2]);
   }
 
   if (d_direct>d_flipped)
     return true;
   return false;
 }
 
 template< class OutImageType, class InputImageType >
 mitk::FiberBundle::Pointer TractParcellationFilter< OutImageType, InputImageType >::GetWorkingFib()
 {
   mitk::FiberBundle::Pointer fib_static_resampled = m_InputTract->GetDeepCopy();
   mitk::Tractometry::ResampleIfNecessary(fib_static_resampled, m_NumParcels);
 
   // clustering
   std::vector< mitk::ClusteringMetric* > metrics;
   metrics.push_back({new mitk::ClusteringMetricEuclideanStd()});
 //  metrics.push_back({new mitk::ClusteringMetricEuclideanMean()});
 //  metrics.push_back({new mitk::ClusteringMetricLength()});
 
   if (m_NumCentroids>0)
   {
     std::vector<mitk::FiberBundle::Pointer> centroids;
     std::shared_ptr< mitk::TractClusteringFilter > clusterer = std::make_shared<mitk::TractClusteringFilter>();
     int c=0;
     while (c<30 && (centroids.empty() || centroids.size()>static_cast<unsigned long>(m_NumCentroids)))
     {
       float cluster_size = m_StartClusterSize + m_StartClusterSize*0.2*c;
       float max_d = 0;
       int i=1;
       std::vector< float > distances;
       while (max_d < m_InputTract->GetGeometry()->GetDiagonalLength()/2)
       {
         distances.push_back(cluster_size*i);
         max_d = cluster_size*i;
         ++i;
       }
 
       clusterer->SetDistances(distances);
       clusterer->SetTractogram(fib_static_resampled);
       clusterer->SetMetrics(metrics);
       clusterer->SetMergeDuplicateThreshold(cluster_size);
       clusterer->SetDoResampling(false);
       clusterer->SetNumPoints(m_NumParcels);
       //  clusterer->SetMaxClusters(m_Controls->m_MaxCentroids->value());
       clusterer->SetMinClusterSize(1);
       clusterer->Update();
       centroids = clusterer->GetOutCentroids();
       ++c;
     }
 
     mitk::FiberBundle::Pointer centroid_bundle = mitk::FiberBundle::New();
     centroid_bundle = centroid_bundle->AddBundles(centroids);
     return centroid_bundle;
   }
 
   return fib_static_resampled;
 }
 
 template< class OutImageType, class InputImageType >
 std::vector< typename itk::Image<unsigned char, 3>::Pointer > TractParcellationFilter< OutImageType, InputImageType >::GetBinarySplit(typename OutImageType::Pointer inImage)
 {
   std::vector< typename itk::Image<unsigned char, 3>::Pointer > binary_maps;
 
   for (unsigned int i=0; i<m_NumParcels; ++i)
   {
     typename itk::Image<unsigned char, 3>::Pointer parcel_image = itk::Image<unsigned char, 3>::New();
     parcel_image->SetSpacing( inImage->GetSpacing() );
     parcel_image->SetOrigin( inImage->GetOrigin() );
     parcel_image->SetDirection( inImage->GetDirection() );
     parcel_image->SetRegions( inImage->GetLargestPossibleRegion() );
     parcel_image->Allocate();
     parcel_image->FillBuffer(0);
 
     binary_maps.push_back(parcel_image);
   }
 
   itk::ImageRegionIterator< itk::Image<unsigned char, 3> > p_it(inImage, inImage->GetLargestPossibleRegion());
   while(!p_it.IsAtEnd())
   {
     if (p_it.Get()>0)
     {
       binary_maps.at(p_it.Get()-1)->SetPixel(p_it.GetIndex(), 1);
     }
     ++p_it;
   }
 
   return binary_maps;
 }
 
 template< class OutImageType, class InputImageType >
 typename OutImageType::Pointer TractParcellationFilter< OutImageType, InputImageType >::PostprocessParcellation(typename OutImageType::Pointer inImage)
 {
   itk::ImageRegionConstIterator< OutImageType > in_it(inImage, inImage->GetLargestPossibleRegion());
 
   typename OutImageType::Pointer outImage = OutImageType::New();
   outImage->SetSpacing( inImage->GetSpacing() );
   outImage->SetOrigin( inImage->GetOrigin() );
   outImage->SetDirection( inImage->GetDirection() );
   outImage->SetRegions( inImage->GetLargestPossibleRegion() );
   outImage->Allocate();
   outImage->FillBuffer(0);
   itk::ImageRegionIterator< OutImageType > out_it(outImage, outImage->GetLargestPossibleRegion());
 
   MITK_INFO << "Postprocessing parcellation";
   while( !in_it.IsAtEnd() )
   {
     if (in_it.Get()>0)
     {
       std::vector< OutPixelType > vote; vote.resize(m_NumParcels + 1, 0);
 
       typename OutImageType::SizeType regionSize;
       regionSize[0] = 3;
       regionSize[1] = 3;
       regionSize[2] = 3;
 
       typename OutImageType::IndexType regionIndex = in_it.GetIndex();
       regionIndex[0] -= 1;
       regionIndex[1] -= 1;
       regionIndex[2] -= 1;
 
       typename OutImageType::RegionType region;
       region.SetSize(regionSize);
       region.SetIndex(regionIndex);
 
       itk::ImageRegionConstIterator<OutImageType> rit(inImage, region);
       while( !rit.IsAtEnd() )
       {
         if (rit.GetIndex()!=regionIndex)
           vote[rit.Get()]++;
         ++rit;
       }
 
       OutPixelType max = 0;
       unsigned int max_parcel = -1;
       for (unsigned int i=1; i<vote.size(); ++i)
       {
         if (vote[i]>max)
         {
           max = vote[i];
           max_parcel = i;
         }
       }
 
       out_it.Set(max_parcel);
     }
     ++out_it;
     ++in_it;
   }
   MITK_INFO << "DONE";
   return outImage;
 }
 
 template< class OutImageType, class InputImageType >
-void TractParcellationFilter< OutImageType, InputImageType >::StaticResampleParcelVoting(typename OutImageType::Pointer outImage)
+void TractParcellationFilter< OutImageType, InputImageType >::CentroidBasedParcellation(typename OutImageType::Pointer outImage)
 {
   typename itk::TractDensityImageFilter< OutImageType >::Pointer generator = itk::TractDensityImageFilter< OutImageType >::New();
   generator->SetFiberBundle(m_InputTract);
   generator->SetMode(TDI_MODE::BINARY);
   generator->SetInputImage(outImage);
   generator->SetUseImageGeometry(true);
   generator->Update();
   auto tdi = generator->GetOutput();
 
   if (m_NumParcels < 3)
-    m_NumParcels = mitk::Tractometry::EstimateNumSamplingPoints(tdi, m_InputTract, 3);
+    m_NumParcels = mitk::Tractometry::EstimateNumSamplingPoints(tdi, m_InputTract, 5);
+
+  m_WorkingTract = GetWorkingFib();
+  vtkSmartPointer< vtkPolyData > polydata = m_WorkingTract->GetFiberPolyData();
+
+  itk::ImageRegionIterator< OutImageType > it(outImage, outImage->GetLargestPossibleRegion());
+  itk::ImageRegionIterator< OutImageType > it_tdi(tdi, tdi->GetLargestPossibleRegion());
+
+  vtkSmartPointer< vtkPolyData > reference_polydata = nullptr;
+  if (m_ReferenceTract.IsNotNull())
+    reference_polydata = m_ReferenceTract->GetFiberPolyData();
+
+  unsigned long num_vox = 0;
+  while( !it_tdi.IsAtEnd() )
+  {
+    if (it_tdi.Get()>0)
+      ++num_vox;
+    ++it_tdi;
+  }
+  it_tdi.GoToBegin();
+
+  MITK_INFO << "Parcellating tract (centroid-based)";
+  boost::progress_display disp(num_vox);
+  while( !it.IsAtEnd() )
+  {
+    if (it_tdi.Get()>0)
+    {
+      int final_seg_id = -1;
+      float min_d = 99999999;
+
+      itk::Image< float, 4 >::IndexType idx4;
+      idx4[0] = it_tdi.GetIndex()[0];
+      idx4[1] = it_tdi.GetIndex()[1];
+      idx4[2] = it_tdi.GetIndex()[2];
+
+      typename OutImageType::PointType image_point;
+      tdi->TransformIndexToPhysicalPoint(it.GetIndex(), image_point);
+
+      for (unsigned int i=0; i<m_WorkingTract->GetNumFibers(); ++i)
+      {
+        vtkCell* cell = polydata->GetCell(i);
+        auto numPoints = cell->GetNumberOfPoints();
+        vtkPoints* points = cell->GetPoints();
+
+        bool flip = Flip(polydata, i, reference_polydata);
+
+        for (int j=0; j<numPoints; j++)
+        {
+          itk::Point<float, 3> p;
+          int segment_id = -1;
+          if (flip)
+          {
+            segment_id = numPoints - j - 1;
+            p = mitk::imv::GetItkPoint(points->GetPoint(segment_id));
+          }
+          else
+          {
+            p = mitk::imv::GetItkPoint(points->GetPoint(j));
+            segment_id = j;
+          }
+
+          float d = std::fabs( (p[0]-image_point[0]) ) + std::fabs( (p[1]-image_point[1]) ) + std::fabs( (p[2]-image_point[2]) );
+
+          if (d<min_d)
+          {
+            min_d = d;
+            final_seg_id = j;
+          }
+        }
+      }
+
+      it.Set(final_seg_id + 1);
+      ++disp;
+    }
+    ++it;
+    ++it_tdi;
+  }
+  MITK_INFO << "DONE";
+}
+
+template< class OutImageType, class InputImageType >
+void TractParcellationFilter< OutImageType, InputImageType >::VotingBasedParcellation(typename OutImageType::Pointer outImage)
+{
+  typename itk::TractDensityImageFilter< OutImageType >::Pointer generator = itk::TractDensityImageFilter< OutImageType >::New();
+  generator->SetFiberBundle(m_InputTract);
+  generator->SetMode(TDI_MODE::BINARY);
+  generator->SetInputImage(outImage);
+  generator->SetUseImageGeometry(true);
+  generator->Update();
+  auto tdi = generator->GetOutput();
+
+  if (m_NumParcels < 3)
+    m_NumParcels = mitk::Tractometry::EstimateNumSamplingPoints(tdi, m_InputTract, 5);
 
   itk::TractsToVectorImageFilter<float>::Pointer fOdfFilter = itk::TractsToVectorImageFilter<float>::New();
   fOdfFilter->SetMaskImage(tdi);
   fOdfFilter->SetFiberBundle(m_InputTract);
   fOdfFilter->SetNormalizationMethod(itk::TractsToVectorImageFilter<float>::NormalizationMethods::SINGLE_VEC_NORM);
   fOdfFilter->SetMaxNumDirections(1);
   fOdfFilter->SetOnlyUseMaskGeometry(true);
   fOdfFilter->Update();
   itk::Image< float, 4 >::Pointer dir_image = fOdfFilter->GetDirectionImage();
 
   m_WorkingTract = GetWorkingFib();
   vtkSmartPointer< vtkPolyData > polydata = m_WorkingTract->GetFiberPolyData();
 
   float maxd = m_InputTract->GetMeanFiberLength()/(0.5*m_NumParcels);
   itk::ImageRegionIterator< OutImageType > it(outImage, outImage->GetLargestPossibleRegion());
   itk::ImageRegionIterator< OutImageType > it_tdi(tdi, tdi->GetLargestPossibleRegion());
 
   vtkSmartPointer< vtkPolyData > reference_polydata = nullptr;
   if (m_ReferenceTract.IsNotNull())
     reference_polydata = m_ReferenceTract->GetFiberPolyData();
 
+  // count voxels in mask
   unsigned long num_vox = 0;
   while( !it_tdi.IsAtEnd() )
   {
     if (it_tdi.Get()>0)
       ++num_vox;
     ++it_tdi;
   }
   it_tdi.GoToBegin();
 
-  MITK_INFO << "Parcellating tract";
+  //
+  std::vector< std::vector< itk::Point<float, 3 > > > streamlines;
+  for (unsigned int i=0; i<m_WorkingTract->GetNumFibers(); ++i)
+  {
+    vtkCell* cell = polydata->GetCell(i);
+    auto numPoints = cell->GetNumberOfPoints();
+    vtkPoints* points = cell->GetPoints();
+
+    bool flip = Flip(polydata, i, reference_polydata);
+
+    std::vector< itk::Point<float, 3 > > streamline;
+
+    for (int j=0; j<numPoints; ++j)
+    {
+      itk::Point<float, 3> p;
+      if (flip)
+        p = mitk::imv::GetItkPoint(points->GetPoint(numPoints - j - 1));
+      else
+        p = mitk::imv::GetItkPoint(points->GetPoint(j));
+      streamline.push_back(p);
+    }
+    streamlines.push_back(streamline);
+  }
+
+  MITK_INFO << "Parcellating tract (parcel-voting)";
   boost::progress_display disp(num_vox);
   while( !it.IsAtEnd() )
   {
     if (it_tdi.Get()>0)
     {
       int final_seg_id = -1;
-      int mult = 1;
+      int mult = 2;
 
       itk::Image< float, 4 >::IndexType idx4;
       idx4[0] = it_tdi.GetIndex()[0];
       idx4[1] = it_tdi.GetIndex()[1];
       idx4[2] = it_tdi.GetIndex()[2];
 
       vnl_vector_fixed<float, 3> ref_dir;
       idx4[3] = 0;
       ref_dir[0] = dir_image->GetPixel(idx4);
       idx4[3] = 1;
       ref_dir[1] = dir_image->GetPixel(idx4);
       idx4[3] = 2;
       ref_dir[2] = dir_image->GetPixel(idx4);
 
       if (ref_dir.magnitude()>0.01)
       {
         ref_dir.normalize();
 
         while(final_seg_id<0 && mult<5)
         {
           std::vector<float> seg_vote; seg_vote.resize(m_NumParcels, 0);
           typename OutImageType::PointType image_point;
           tdi->TransformIndexToPhysicalPoint(it.GetIndex(), image_point);
 
-          for (unsigned int i=0; i<m_WorkingTract->GetNumFibers(); ++i)
+#pragma omp parallel for
+          for (int sid = 0; sid<static_cast<int>(streamlines.size()); ++sid)
           {
-            vtkCell* cell = polydata->GetCell(i);
-            auto numPoints = cell->GetNumberOfPoints();
-            vtkPoints* points = cell->GetPoints();
-
-            bool flip = Flip(polydata, i, reference_polydata);
 
             float local_d = 99999999;
             int local_closest_seg = -1;
-  //          float weight = 1.0;
 
-            for (int j=0; j<numPoints; j++)
+            auto streamline = streamlines[sid];
+            for (unsigned int segment_id=0; segment_id<streamline.size(); ++segment_id)
             {
-              itk::Point<float, 3> p;
-              int segment_id = -1;
-              if (flip)
-              {
-                segment_id = numPoints - j - 1;
-                p = mitk::imv::GetItkPoint(points->GetPoint(segment_id));
-              }
-              else
-              {
-                p = mitk::imv::GetItkPoint(points->GetPoint(j));
-                segment_id = j;
-              }
-
+              itk::Point<float, 3> p = streamline[segment_id];
               float d = std::fabs( (p[0]-image_point[0]) ) + std::fabs( (p[1]-image_point[1]) ) + std::fabs( (p[2]-image_point[2]) );
 
-
               itk::Point<float, 3> p2;
-              if (segment_id<numPoints-1)
-              {
-                p2 = mitk::imv::GetItkPoint(points->GetPoint(segment_id+1));
-              }
+              if (segment_id<streamline.size()-1)
+                p2 = streamline.at(segment_id+1);
               else
-              {
-                p2 = mitk::imv::GetItkPoint(points->GetPoint(segment_id-1));
-              }
+                p2 = streamline.at(segment_id-1);
 
               vnl_vector_fixed<float, 3> dir;
               dir[0] = p[0]-p2[0];
               dir[1] = p[1]-p2[1];
               dir[2] = p[2]-p2[2];
 
               if (dir.magnitude()<0.0000001)
                 continue;
               dir.normalize();
 
               float a = std::fabs(dot_product(dir, ref_dir));
               if (a<0.0000001)
                 a += 0.0000001;
               d += (1.0/a - 1.0) * maxd;
 
               if (d<local_d)
               {
                 local_d = d;
-                local_closest_seg = j;
-  //              typename OutImageType::IndexType p_idx;
-  //              typename OutImageType::PointType mitk_p;
-  //              mitk_p[0] = p[0];
-  //              mitk_p[1] = p[2];
-  //              mitk_p[2] = p[1];
-  //              tdi->TransformPhysicalPointToIndex(mitk_p, p_idx);
-  //              weight = tdi->GetPixel(p_idx);
+                local_closest_seg = segment_id;
               }
             }
 
+#pragma omp critical
             if (local_d<maxd*mult)
+            {
               seg_vote[local_closest_seg] += 1.0/(local_d);
+            }
           }
 
           float max_count = 0;
           for (unsigned int i=0; i<seg_vote.size(); ++i)
           {
             if (seg_vote.at(i)>max_count)
             {
               final_seg_id = i;
               max_count = seg_vote.at(i);
             }
           }
 
           if (final_seg_id>=0)
             it.Set(final_seg_id + 1);
 
           ++mult;
         }
       }
 
       ++disp;
     }
     ++it;
     ++it_tdi;
   }
   MITK_INFO << "DONE";
 }
 
 template< class OutImageType, class InputImageType >
 void TractParcellationFilter< OutImageType, InputImageType >::GenerateData()
 {
   // generate upsampled image
   mitk::BaseGeometry::Pointer geometry = m_InputTract->GetGeometry();
 
   // calculate new image parameters
   itk::Vector<double, 3> newSpacing;
   mitk::Point3D newOrigin;
   itk::Matrix<double, 3, 3> newDirection;
   ImageRegion<3> upsampledRegion;
   if (!m_InputImage.IsNull())
   {
     newSpacing = m_InputImage->GetSpacing()/m_UpsamplingFactor;
     upsampledRegion = m_InputImage->GetLargestPossibleRegion();
     newOrigin = m_InputImage->GetOrigin();
     typename OutImageType::RegionType::SizeType size = upsampledRegion.GetSize();
     size[0] *= m_UpsamplingFactor;
     size[1] *= m_UpsamplingFactor;
     size[2] *= m_UpsamplingFactor;
     upsampledRegion.SetSize(size);
     newDirection = m_InputImage->GetDirection();
   }
   else
   {
     newSpacing = geometry->GetSpacing()/m_UpsamplingFactor;
     newOrigin = geometry->GetOrigin();
     mitk::Geometry3D::BoundsArrayType bounds = geometry->GetBounds();
     newOrigin[0] += bounds.GetElement(0);
     newOrigin[1] += bounds.GetElement(2);
     newOrigin[2] += bounds.GetElement(4);
 
     for (int i=0; i<3; i++)
       for (int j=0; j<3; j++)
         newDirection[j][i] = geometry->GetMatrixColumn(i)[j];
     upsampledRegion.SetSize(0, geometry->GetExtent(0)*m_UpsamplingFactor);
     upsampledRegion.SetSize(1, geometry->GetExtent(1)*m_UpsamplingFactor);
     upsampledRegion.SetSize(2, geometry->GetExtent(2)*m_UpsamplingFactor);
   }
 
   // apply new image parameters
   typename OutImageType::Pointer outImage = OutImageType::New();
   outImage->SetSpacing( newSpacing );
   outImage->SetOrigin( newOrigin );
   outImage->SetDirection( newDirection );
   outImage->SetRegions( upsampledRegion );
   outImage->Allocate();
   outImage->FillBuffer(0);
   this->SetNthOutput(0, outImage);
 
-  StaticResampleParcelVoting(outImage);
+  if (m_NumCentroids > 0)
+    CentroidBasedParcellation(outImage);
+  else
+    VotingBasedParcellation(outImage);
 
   auto outImage_pp = PostprocessParcellation(outImage);
   outImage_pp = PostprocessParcellation(outImage_pp);
   outImage_pp = PostprocessParcellation(outImage_pp);
   this->SetNthOutput(1, outImage_pp);
 }
 }
diff --git a/Modules/FiberTracking/Algorithms/itkTractParcellationFilter.h b/Modules/FiberTracking/Algorithms/itkTractParcellationFilter.h
index 15f3875..3b0d258 100644
--- a/Modules/FiberTracking/Algorithms/itkTractParcellationFilter.h
+++ b/Modules/FiberTracking/Algorithms/itkTractParcellationFilter.h
@@ -1,96 +1,97 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center.
 
 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 __itkTractParcellationFilter_h__
 #define __itkTractParcellationFilter_h__
 
 #include <itkImageSource.h>
 #include <itkImage.h>
 #include <mitkFiberBundle.h>
 
 namespace itk{
 
 /**
 * \brief Generates image where the pixel values are set according to the position along a fiber bundle.   */
 
 template< class OutImageType=itk::Image<unsigned char, 3>, class InputImageType=OutImageType >
 class TractParcellationFilter : public ImageSource< OutImageType >
 {
 
 public:
   typedef TractParcellationFilter Self;
   typedef ProcessObject Superclass;
   typedef SmartPointer< Self > Pointer;
   typedef SmartPointer< const Self > ConstPointer;
 
   typedef typename OutImageType::PixelType OutPixelType;
 
   itkFactorylessNewMacro(Self)
   itkCloneMacro(Self)
   itkTypeMacro( TractParcellationFilter, ImageSource )
 
   itkSetMacro( UpsamplingFactor, float)
   itkGetMacro( UpsamplingFactor, float)
 
   itkSetMacro( NumParcels, unsigned int )
   itkGetMacro( NumParcels, unsigned int )
 
   itkSetMacro( NumCentroids, unsigned int )
   itkGetMacro( NumCentroids, unsigned int )
 
   itkSetMacro( StartClusterSize, float )
   itkGetMacro( StartClusterSize, float )
 
   itkSetMacro( InputTract, mitk::FiberBundle::Pointer)
   itkSetMacro( ReferenceTract, mitk::FiberBundle::Pointer)
   itkGetMacro( WorkingTract, mitk::FiberBundle::Pointer)
 
   itkSetMacro( InputImage, typename InputImageType::Pointer)
 
 
   std::vector< typename itk::Image<unsigned char, 3>::Pointer > GetBinarySplit(typename OutImageType::Pointer inImage);
 
   void GenerateData() override;
 
 protected:
 
   TractParcellationFilter();
   virtual ~TractParcellationFilter();
 
   bool Flip(vtkSmartPointer< vtkPolyData > polydata1, int i, vtkSmartPointer< vtkPolyData > ref_poly=nullptr, int ref_i=0);
 
   mitk::FiberBundle::Pointer GetWorkingFib();
 
   typename OutImageType::Pointer PostprocessParcellation(typename OutImageType::Pointer outImage);
 
-  void StaticResampleParcelVoting(typename OutImageType::Pointer outImage);
+  void VotingBasedParcellation(typename OutImageType::Pointer outImage);
+  void CentroidBasedParcellation(typename OutImageType::Pointer outImage);
 
   mitk::FiberBundle::Pointer        m_ReferenceTract;
   mitk::FiberBundle::Pointer        m_InputTract;          ///< input fiber bundle
   mitk::FiberBundle::Pointer        m_WorkingTract;
   float                             m_UpsamplingFactor;     ///< use higher resolution for ouput image
   unsigned int                      m_NumParcels;
   unsigned int                      m_NumCentroids;
   float                             m_StartClusterSize;
   typename InputImageType::Pointer  m_InputImage;
 };
 
 }
 
 #ifndef ITK_MANUAL_INSTANTIATION
 #include "itkTractParcellationFilter.cpp"
 #endif
 
 #endif // __itkTractParcellationFilter_h__
diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkTractometryView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkTractometryView.cpp
index 01fba29..a6a01a0 100644
--- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkTractometryView.cpp
+++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkTractometryView.cpp
@@ -1,450 +1,450 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center.
 
 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 <berryISelectionService.h>
 #include <berryIWorkbenchWindow.h>
 
 #include "QmitkTractometryView.h"
 
 #include <mitkNodePredicateDataType.h>
 #include <mitkNodePredicateDimension.h>
 #include <mitkNodePredicateAnd.h>
 #include <mitkImageCast.h>
 #include <mitkImageAccessByItk.h>
 #include <mitkImage.h>
 #include <mitkFiberBundle.h>
 #include <mitkDiffusionPropertyHelper.h>
 #include <mitkITKImageImport.h>
 #include <mitkPixelTypeMultiplex.h>
 #include <mitkImagePixelReadAccessor.h>
 #include <berryIPreferences.h>
 #include <berryWorkbenchPlugin.h>
 #include <berryQtPreferences.h>
 #include <vtkLookupTable.h>
 #include <QClipboard>
 #include <mitkLexicalCast.h>
 #include <QmitkChartWidget.h>
 #include <mitkLookupTable.h>
 #include <mitkTractClusteringFilter.h>
 #include <mitkClusteringMetricEuclideanStd.h>
 #include <itkTractDensityImageFilter.h>
 #include <itkTractParcellationFilter.h>
 #include <mitkLookupTableProperty.h>
 #include <mitkLevelWindowProperty.h>
 #include <itkMaskedStatisticsImageFilter.h>
 #include <mitkDiffusionFunctionCollection.h>
 #include <mitkTractometry.h>
 
 
 const std::string QmitkTractometryView::VIEW_ID = "org.mitk.views.tractometry";
 using namespace mitk;
 
 QmitkTractometryView::QmitkTractometryView()
   : QmitkAbstractView()
   , m_Controls( nullptr )
   , m_Visible(false)
 {
 
 }
 
 // Destructor
 QmitkTractometryView::~QmitkTractometryView()
 {
 
 }
 
 void QmitkTractometryView::CreateQtPartControl( QWidget *parent )
 {
   // build up qt view, unless already done
   if ( !m_Controls )
   {
     // create GUI widgets from the Qt Designer's .ui file
     m_Controls = new Ui::QmitkTractometryViewControls;
     m_Controls->setupUi( parent );
 
     connect( m_Controls->m_MethodBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui()) );
     connect( m_Controls->m_StartButton, SIGNAL(clicked()), this, SLOT(StartTractometry()) );
 
     mitk::TNodePredicateDataType<mitk::Image>::Pointer imageP = mitk::TNodePredicateDataType<mitk::Image>::New();
     mitk::NodePredicateDimension::Pointer dimP = mitk::NodePredicateDimension::New(3);
 
     m_Controls->m_ImageBox->SetDataStorage(this->GetDataStorage());
     m_Controls->m_ImageBox->SetPredicate(mitk::NodePredicateAnd::New(imageP, dimP));
 
     m_Controls->m_ChartWidget->SetXAxisLabel("Tract position");
     m_Controls->m_ChartWidget->SetYAxisLabel("Image Value");
 
     m_Controls->m_ChartWidget->SetTheme(QmitkChartWidget::ColorTheme::darkstyle);
   }
 }
 
 void QmitkTractometryView::SetFocus()
 {
 }
 
 void QmitkTractometryView::UpdateGui()
 {
   berry::IWorkbenchPart::Pointer nullPart;
   OnSelectionChanged(nullPart, QList<mitk::DataNode::Pointer>(m_CurrentSelection));
 }
 
 void QmitkTractometryView::StaticResamplingTractometry(mitk::Image::Pointer image, mitk::DataNode::Pointer node, std::vector<std::vector<double> > &data, std::string& clipboard_string)
 {
   itk::Image<float, 3>::Pointer itkImage = itk::Image<float, 3>::New();
   CastToItkImage(image, itkImage);
 
   mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
 
   unsigned int num_points = m_NumSamplingPoints;
   mitk::FiberBundle::Pointer working_fib = fib->GetDeepCopy();
 
   vnl_matrix<float> output = mitk::Tractometry::StaticResamplingTractometry(itkImage, working_fib, num_points, m_ReferenceFib);
 
   std::vector< double > std_values1;
   std::vector< double > std_values2;
   std::vector< double > mean_values;
 
   for (unsigned int i=0; i<output.rows(); ++i)
   {
     auto row = output.get_row(i);
     float mean = row.mean();
     double stdev = 0;
 
     for (unsigned int j=0; j<row.size(); ++j)
     {
       double diff = mean - row.get(j);
       diff *= diff;
       stdev += diff;
     }
     stdev /= row.size();
     stdev = std::sqrt(stdev);
 
     clipboard_string += boost::lexical_cast<std::string>(mean);
     clipboard_string += " ";
     clipboard_string += boost::lexical_cast<std::string>(stdev);
     clipboard_string += "\n";
 
     mean_values.push_back(mean);
     std_values1.push_back(mean + stdev);
     std_values2.push_back(mean - stdev);
   }
   clipboard_string += "\n";
 
   data.push_back(mean_values);
   data.push_back(std_values1);
   data.push_back(std_values2);
 
   if (m_Controls->m_ShowBinned->isChecked())
   {
     mitk::DataNode::Pointer new_node = mitk::DataNode::New();
     new_node->SetData(working_fib);
     new_node->SetName("binned_static");
     new_node->SetVisibility(true);
     node->SetVisibility(false);
     GetDataStorage()->Add(new_node, node);
   }
 }
 
 void QmitkTractometryView::NearestCentroidPointTractometry(mitk::Image::Pointer image, mitk::DataNode::Pointer node, std::vector< std::vector< double > >& data, std::string& clipboard_string)
 {
   mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
 
   unsigned int num_points = m_NumSamplingPoints;
 
   mitk::FiberBundle::Pointer working_fib = fib->GetDeepCopy();
   working_fib->ResampleSpline(1.0);
 
   itk::Image<float, 3>::Pointer itkImage = itk::Image<float, 3>::New();
   CastToItkImage(image, itkImage);
 
   auto output = mitk::Tractometry::NearestCentroidPointTractometry(itkImage, working_fib, num_points, m_Controls->m_MaxCentroids->value(), m_Controls->m_ClusterSize->value(), m_ReferenceFib);
 
   std::vector< double > std_values1;
   std::vector< double > std_values2;
   std::vector< double > mean_values;
 
   for (auto row : output)
   {
     float mean = row.mean();
     double stdev = 0;
 
     for (unsigned int j=0; j<row.size(); ++j)
     {
       double diff = mean - row.get(j);
       diff *= diff;
       stdev += diff;
     }
     stdev /= row.size();
     stdev = std::sqrt(stdev);
 
     clipboard_string += boost::lexical_cast<std::string>(mean);
     clipboard_string += " ";
     clipboard_string += boost::lexical_cast<std::string>(stdev);
     clipboard_string += "\n";
 
     mean_values.push_back(mean);
     std_values1.push_back(mean + stdev);
     std_values2.push_back(mean - stdev);
   }
   clipboard_string += "\n";
 
   data.push_back(mean_values);
   data.push_back(std_values1);
   data.push_back(std_values2);
 
   if (m_Controls->m_ShowBinned->isChecked())
   {
     mitk::DataNode::Pointer new_node = mitk::DataNode::New();
     new_node->SetData(working_fib);
     new_node->SetName("binned_centroid");
     new_node->SetVisibility(true);
     node->SetVisibility(false);
     GetDataStorage()->Add(new_node, node);
   }
 }
 
 void QmitkTractometryView::Activated()
 {
 
 }
 
 void QmitkTractometryView::Deactivated()
 {
 
 }
 
 void QmitkTractometryView::Visible()
 {
   m_Visible = true;
   QList<mitk::DataNode::Pointer> selection = GetDataManagerSelection();
   berry::IWorkbenchPart::Pointer nullPart;
   OnSelectionChanged(nullPart, selection);
 }
 
 void QmitkTractometryView::Hidden()
 {
   m_Visible = false;
 }
 
 std::string QmitkTractometryView::RGBToHexString(double *rgb)
 {
   std::ostringstream os;
   for (int i = 0; i < 3; ++i)
     {
     os << std::setw(2) << std::setfill('0') << std::hex
        << static_cast<int>(rgb[i] * 255);
     }
   return os.str();
 }
 
 void QmitkTractometryView::AlongTractRadiomicsPreprocessing(mitk::Image::Pointer image, mitk::DataNode::Pointer node, std::vector< std::vector< double > >& data, std::string& clipboard_string)
 {
   mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
 
   // calculate mask
   typedef itk::Image<unsigned char, 3> ParcellationImageType;
   ParcellationImageType::Pointer itkImage = ParcellationImageType::New();
   CastToItkImage(image, itkImage);
 
   itk::TractParcellationFilter< >::Pointer parcellator = itk::TractParcellationFilter< >::New();
   parcellator->SetInputImage(itkImage);
   parcellator->SetNumParcels(m_NumSamplingPoints);
   parcellator->SetInputTract(fib);
   parcellator->SetNumCentroids(m_Controls->m_MaxCentroids->value());
   parcellator->SetStartClusterSize(m_Controls->m_ClusterSize->value());
   parcellator->Update();
   ParcellationImageType::Pointer out_image = parcellator->GetOutput(0);
   ParcellationImageType::Pointer out_image_pp = parcellator->GetOutput(1);
   auto binary_masks = parcellator->GetBinarySplit(out_image_pp);
 
   mitk::Image::Pointer seg_img = mitk::Image::New();
   seg_img->InitializeByItk(out_image.GetPointer());
   seg_img->SetVolume(out_image->GetBufferPointer());
 
   mitk::Image::Pointer seg_img_pp = mitk::Image::New();
   seg_img_pp->InitializeByItk(out_image_pp.GetPointer());
   seg_img_pp->SetVolume(out_image_pp->GetBufferPointer());
 
   std::vector< double > std_values1;
   std::vector< double > std_values2;
   std::vector< double > mean_values;
   for (auto mask : binary_masks)
   {
     itk::Image<float, 3>::Pointer data_image = itk::Image<float, 3>::New();
     CastToItkImage(image, data_image);
     itk::MaskedStatisticsImageFilter<itk::Image<float, 3>>::Pointer statisticsImageFilter = itk::MaskedStatisticsImageFilter<itk::Image<float, 3>>::New();
     statisticsImageFilter->SetInput(data_image);
     statisticsImageFilter->SetMask(mask);
     statisticsImageFilter->Update();
     double mean = statisticsImageFilter->GetMean();
     double stdev = std::sqrt(statisticsImageFilter->GetVariance());
 
     std_values1.push_back(mean + stdev);
     std_values2.push_back(mean - stdev);
     mean_values.push_back(mean);
 
     clipboard_string += boost::lexical_cast<std::string>(mean);
     clipboard_string += " ";
     clipboard_string += boost::lexical_cast<std::string>(stdev);
     clipboard_string += "\n";
   }
   clipboard_string += "\n";
 
   data.push_back(mean_values);
   data.push_back(std_values1);
   data.push_back(std_values2);
 
   mitk::LookupTable::Pointer lut = mitk::LookupTable::New();
   lut->SetType( mitk::LookupTable::MULTILABEL );
   mitk::LookupTableProperty::Pointer lut_prop = mitk::LookupTableProperty::New();
   lut_prop->SetLookupTable( lut );
 
   mitk::LevelWindow lw;
   lw.SetRangeMinMax(0, parcellator->GetNumParcels());
 
 //  mitk::DataNode::Pointer new_node = mitk::DataNode::New();
 //  new_node->SetData(seg_img);
 //  new_node->SetName("tract parcellation");
 //  new_node->SetVisibility(true);
 //  new_node->SetProperty("LookupTable", lut_prop );
 //  new_node->SetProperty( "levelwindow", mitk::LevelWindowProperty::New( lw ) );
 //  node->SetVisibility(false);
 //  GetDataStorage()->Add(new_node, node);
 
   mitk::DataNode::Pointer new_node2 = mitk::DataNode::New();
   new_node2->SetData(seg_img_pp);
   new_node2->SetName("tract parcellation pp");
   new_node2->SetVisibility(false);
   new_node2->SetProperty("LookupTable", lut_prop );
   new_node2->SetProperty( "levelwindow", mitk::LevelWindowProperty::New( lw ) );
   GetDataStorage()->Add(new_node2, node);
 
   mitk::DataNode::Pointer new_node3 = mitk::DataNode::New();
   auto working_fib = fib->GetDeepCopy();
   working_fib->ColorFibersByScalarMap(seg_img, false, false, mitk::LookupTable::LookupTableType::MULTILABEL, 0.9);
   new_node3->SetData(working_fib);
-  new_node3->SetName("centroids");
+  new_node3->SetName("bin-colored fib");
   GetDataStorage()->Add(new_node3, node);
 }
 
 void QmitkTractometryView::StartTractometry()
 {
   m_ReferenceFib = dynamic_cast<mitk::FiberBundle*>(m_CurrentSelection.at(0)->GetData())->GetDeepCopy();
 
   mitk::Image::Pointer image = dynamic_cast<mitk::Image*>(m_Controls->m_ImageBox->GetSelectedNode()->GetData());
 
   MITK_INFO << "Resanmpling reference fibers";
   if (m_Controls->m_SamplingPointsBox->value()<3)
   {
     typedef itk::Image<unsigned char, 3> ParcellationImageType;
     ParcellationImageType::Pointer itkImage = ParcellationImageType::New();
     CastToItkImage(image, itkImage);
 
-    m_NumSamplingPoints = mitk::Tractometry::EstimateNumSamplingPoints(itkImage, m_ReferenceFib, 3);
+    m_NumSamplingPoints = mitk::Tractometry::EstimateNumSamplingPoints(itkImage, m_ReferenceFib, 5);
   }
   else
     m_NumSamplingPoints = m_Controls->m_SamplingPointsBox->value();
   m_ReferenceFib->ResampleToNumPoints(m_NumSamplingPoints);
 
   double color[3] = {0,0,0};
   mitk::LookupTable::Pointer lookupTable = mitk::LookupTable::New();
   lookupTable->SetType(mitk::LookupTable::MULTILABEL);
 
   this->m_Controls->m_ChartWidget->Clear();
   std::string clipboardString = "";
   int c = 1;
   for (auto node : m_CurrentSelection)
   {
     clipboardString += node->GetName() + "\n";
     clipboardString += "mean stdev\n";
 
     std::vector< std::vector< double > > data;
     switch (m_Controls->m_MethodBox->currentIndex())
     {
     case 0:
     {
       StaticResamplingTractometry( image, node, data, clipboardString );
       break;
     }
     case 1:
     {
       NearestCentroidPointTractometry( image, node, data, clipboardString );
       break;
     }
     case 2:
     {
       AlongTractRadiomicsPreprocessing(image, node, data, clipboardString);
       break;
     }
     default:
     {
       StaticResamplingTractometry( image, node, data, clipboardString );
     }
     }
 
     m_Controls->m_ChartWidget->AddData1D(data.at(0), node->GetName() + " Mean", QmitkChartWidget::ChartType::line);
     m_Controls->m_ChartWidget->SetLineStyle(node->GetName() + " Mean", QmitkChartWidget::LineStyle::solid);
     if (m_Controls->m_StDevBox->isChecked())
     {
       m_Controls->m_ChartWidget->AddData1D(data.at(1), node->GetName() + " +STDEV", QmitkChartWidget::ChartType::line);
       m_Controls->m_ChartWidget->AddData1D(data.at(2), node->GetName() + " -STDEV", QmitkChartWidget::ChartType::line);
       m_Controls->m_ChartWidget->SetLineStyle(node->GetName() + " +STDEV", QmitkChartWidget::LineStyle::dashed);
       m_Controls->m_ChartWidget->SetLineStyle(node->GetName() + " -STDEV", QmitkChartWidget::LineStyle::dashed);
     }
 
     lookupTable->GetTableValue(c, color);
     this->m_Controls->m_ChartWidget->SetColor(node->GetName() + " Mean", RGBToHexString(color));
 
     if (m_Controls->m_StDevBox->isChecked())
     {
       color[0] *= 0.5;
       color[1] *= 0.5;
       color[2] *= 0.5;
       this->m_Controls->m_ChartWidget->SetColor(node->GetName() + " +STDEV", RGBToHexString(color));
       this->m_Controls->m_ChartWidget->SetColor(node->GetName() + " -STDEV", RGBToHexString(color));
     }
 
     this->m_Controls->m_ChartWidget->Show(true);
     this->m_Controls->m_ChartWidget->SetShowDataPoints(false);
     ++c;
   }
 
   QApplication::clipboard()->setText(clipboardString.c_str(), QClipboard::Clipboard);
 }
 
 void QmitkTractometryView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList<mitk::DataNode::Pointer>& nodes)
 {
   m_Controls->m_StartButton->setEnabled(false);
 
   if (!m_Visible)
     return;
 
   if (m_Controls->m_MethodBox->currentIndex()==0)
     m_Controls->m_ClusterFrame->setVisible(false);
   else
     m_Controls->m_ClusterFrame->setVisible(true);
 
   m_CurrentSelection.clear();
   if(m_Controls->m_ImageBox->GetSelectedNode().IsNull())
     return;
 
   for (auto node: nodes)
     if ( dynamic_cast<mitk::FiberBundle*>(node->GetData()) )
       m_CurrentSelection.push_back(node);
   if (!m_CurrentSelection.empty())
     m_Controls->m_StartButton->setEnabled(true);
 }