diff --git a/Modules/DiffusionCmdApps/FiberProcessing/CMakeLists.txt b/Modules/DiffusionCmdApps/FiberProcessing/CMakeLists.txt
index 2acf97a..a15b7a2 100644
--- a/Modules/DiffusionCmdApps/FiberProcessing/CMakeLists.txt
+++ b/Modules/DiffusionCmdApps/FiberProcessing/CMakeLists.txt
@@ -1,50 +1,51 @@
 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^^
     Sift2WeightCopy^^
     FiberExtraction^^
     FiberExtractionRoi^^
     FiberProcessing^^
     FitFibersToImage^^
     FiberDirectionExtraction^^
     FiberJoin^^
     FiberClustering^^
     TractDensityFilter^^
     PrintFiberStatistics^^
     DistanceToSegmentation^^
     TractParcellation^^
+    TractometryComparison^^
     )
 
     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 MitkDiffusionCmdApps MitkMriSimulation MitkFiberTracking ${dependencies_list}
         PACKAGE_DEPENDS
       )
     endforeach()
 
   if(EXECUTABLE_IS_ENABLED)
     MITK_INSTALL_TARGETS(EXECUTABLES ${EXECUTABLE_TARGET})
   endif()
 
   endif()
diff --git a/Modules/DiffusionCmdApps/FiberProcessing/TractometryComparison.cpp b/Modules/DiffusionCmdApps/FiberProcessing/TractometryComparison.cpp
new file mode 100644
index 0000000..8b980a8
--- /dev/null
+++ b/Modules/DiffusionCmdApps/FiberProcessing/TractometryComparison.cpp
@@ -0,0 +1,133 @@
+/*===================================================================
+
+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 <mitkIOUtil.h>
+#include <metaCommand.h>
+#include "mitkDiffusionCommandLineParser.h"
+#include <usAny.h>
+#include <mitkIOUtil.h>
+#include <mitkLexicalCast.h>
+#include <mitkCoreObjectFactory.h>
+#include <mitkPlanarFigure.h>
+#include <mitkPlanarFigureComposite.h>
+#include <mitkFiberBundle.h>
+#include <mitkDiffusionDataIOHelper.h>
+#include <mitkTractometry.h>
+#include <itksys/SystemTools.hxx>
+#include <itkTractParcellationFilter.h>
+
+#define _USE_MATH_DEFINES
+#include <math.h>
+
+/*!
+\brief Join multiple tractograms
+*/
+int main(int argc, char* argv[])
+{
+  mitkDiffusionCommandLineParser parser;
+
+  parser.setTitle("Fiber Statistics");
+  parser.setCategory("Fiber Tracking and Processing Methods");
+  parser.setContributor("MIC");
+  parser.setDescription("Output statistics about fiber lengths, weights, etc. to a file.");
+
+  parser.setArgumentPrefix("--", "-");
+  parser.addArgument("tracts", "", mitkDiffusionCommandLineParser::StringList, "Input tracts:", "input tracts", us::Any(), false);
+  parser.addArgument("images", "", mitkDiffusionCommandLineParser::StringList, "Input images:", "input images", us::Any(), false);
+  parser.addArgument("out_file", "", mitkDiffusionCommandLineParser::String, "Output File:", "output file", us::Any(), false);
+
+  std::map<std::string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
+  if (parsedArgs.size()==0)
+    return EXIT_FAILURE;
+
+  mitkDiffusionCommandLineParser::StringContainerType inFibs = us::any_cast<mitkDiffusionCommandLineParser::StringContainerType>(parsedArgs["tracts"]);
+  mitkDiffusionCommandLineParser::StringContainerType inImages = us::any_cast<mitkDiffusionCommandLineParser::StringContainerType>(parsedArgs["images"]);
+  std::string outfile = us::any_cast<std::string>(parsedArgs["out_file"]);
+  int num_parcels = 20;
+
+  try
+  {
+
+    std::vector< std::string > fib_names, img_names;
+    auto input_tracts = mitk::DiffusionDataIOHelper::load_fibs(inFibs, &fib_names);
+    auto input_images = mitk::DiffusionDataIOHelper::load_itk_images<itk::Image<float, 3>>(inImages, &img_names);
+
+    bool add_header = !itksys::SystemTools::FileExists(outfile, true);
+
+
+    std::ofstream statistics_file;
+    statistics_file.open (outfile, std::ios_base::app);
+    if (add_header)
+      statistics_file << "method;tract_file;image_file;parcel;value" << std::endl;
+
+    auto ref = input_tracts.at(0)->GetDeepCopy();
+    for (unsigned int i=0; i<inFibs.size(); ++i)
+    {
+      auto fib = input_tracts.at(i);
+      fib->ResampleSpline(1.0);
+
+      {
+        auto output = mitk::Tractometry::StaticResamplingTractometry(input_images.at(i), fib->GetDeepCopy(), num_parcels, ref);
+        for (unsigned int r=0; r<output.rows(); ++r)
+          for (unsigned int c=0; c<output.cols(); ++c)
+            statistics_file << "StaticResamplingTractometry;" << fib_names.at(i) << ";" << img_names.at(i) << ";" << boost::lexical_cast<std::string>(r) << ";" << boost::lexical_cast<std::string>(output.get(r,c)) << std::endl;
+      }
+
+      {
+        auto output = mitk::Tractometry::NearestCentroidPointTractometry(input_images.at(i), fib->GetDeepCopy(), num_parcels, 1, 99, ref);
+        for (unsigned int r=0; r<output.size(); ++r)
+          for (unsigned int c=0; c<output.at(r).size(); ++c)
+            statistics_file << "NearestCentroidPointTractometry;" << fib_names.at(i) << ";" << img_names.at(i) << ";" << boost::lexical_cast<std::string>(r) << ";" << boost::lexical_cast<std::string>(output.at(r).get(c)) << std::endl;
+      }
+      {
+        itk::TractParcellationFilter< itk::Image<unsigned char, 3>, itk::Image<float, 3> >::Pointer parcellator = itk::TractParcellationFilter< itk::Image<unsigned char, 3>, itk::Image<float, 3> >::New();
+        parcellator->SetInputImage(input_images.at(i));
+        parcellator->SetNumParcels(num_parcels);
+        parcellator->SetInputTract(fib->GetDeepCopy());
+        parcellator->SetReferenceTract(ref);
+        parcellator->Update();
+        itk::Image<unsigned char, 3>::Pointer out_image_pp = parcellator->GetOutput(1);
+
+        itk::ImageRegionConstIterator< itk::Image<float, 3> > it(input_images.at(i), input_images.at(i)->GetLargestPossibleRegion());
+        itk::ImageRegionConstIterator< itk::Image<unsigned char, 3> > it2(out_image_pp, out_image_pp->GetLargestPossibleRegion());
+        while (!it.IsAtEnd()) {
+          if (it2.Get()>0)
+            statistics_file << "TractParcellationFilter;" << fib_names.at(i) << ";" << img_names.at(i) << ";" << boost::lexical_cast<std::string>(it2.Get()-1) << ";" << boost::lexical_cast<std::string>(it.Get()) << std::endl;
+          ++it;
+          ++it2;
+        }
+      }
+    }
+
+    statistics_file.close();
+  }
+  catch (const itk::ExceptionObject& e)
+  {
+    std::cout << e.what();
+    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/FiberTracking/Algorithms/itkTractParcellationFilter.cpp b/Modules/FiberTracking/Algorithms/itkTractParcellationFilter.cpp
index 2f50812..d042415 100644
--- a/Modules/FiberTracking/Algorithms/itkTractParcellationFilter.cpp
+++ b/Modules/FiberTracking/Algorithms/itkTractParcellationFilter.cpp
@@ -1,469 +1,469 @@
 /*===================================================================
 
 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>
 
 namespace itk{
 
-template< class OutImageType >
-TractParcellationFilter< OutImageType >::TractParcellationFilter()
+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 >
-TractParcellationFilter< OutImageType >::~TractParcellationFilter()
+template< class OutImageType, class InputImageType >
+TractParcellationFilter< OutImageType, InputImageType >::~TractParcellationFilter()
 {
 }
 
-template< class OutImageType >
-bool TractParcellationFilter< OutImageType >::Flip(vtkSmartPointer< vtkPolyData > polydata1, int i, vtkSmartPointer< vtkPolyData > ref_poly, int ref_i)
+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 >
-mitk::FiberBundle::Pointer TractParcellationFilter< OutImageType >::GetWorkingFib()
+template< class OutImageType, class InputImageType >
+mitk::FiberBundle::Pointer TractParcellationFilter< OutImageType, InputImageType >::GetWorkingFib()
 {
   mitk::FiberBundle::Pointer fib_static_resampled = m_InputTract->GetDeepCopy();
   fib_static_resampled->ResampleToNumPoints(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 >
-std::vector< typename itk::Image<unsigned char, 3>::Pointer > TractParcellationFilter< OutImageType >::GetBinarySplit(typename OutImageType::Pointer inImage)
+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 >
-typename OutImageType::Pointer TractParcellationFilter< OutImageType >::PostprocessParcellation(typename OutImageType::Pointer inImage)
+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 >
-void TractParcellationFilter< OutImageType >::StaticResampleParcelVoting(typename OutImageType::Pointer outImage)
+template< class OutImageType, class InputImageType >
+void TractParcellationFilter< OutImageType, InputImageType >::StaticResampleParcelVoting(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)
   {
     auto spacing = outImage->GetSpacing();
     float f = (spacing[0] + spacing[1] + spacing[2])/3;
     f /= generator->GetAverageSegmentLength();
     MITK_INFO << generator->GetAverageSegmentLength();
     MITK_INFO << generator->GetAverageNumTraversedVoxels();
     MITK_INFO << f;
     m_NumParcels = std::ceil(generator->GetAverageNumTraversedVoxels()/(3.0*f));
 
     MITK_INFO << "Automatically setting number of parcels to " << m_NumParcels;
   }
 
   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->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();
 
   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";
   boost::progress_display disp(num_vox);
   while( !it.IsAtEnd() )
   {
     if (it_tdi.Get()>0)
     {
       int final_seg_id = -1;
       int mult = 1;
 
       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);
 
       while(final_seg_id<0)
       {
         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)
         {
           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++)
           {
             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]) );
 
 
             itk::Point<float, 3> p2;
             if (segment_id<numPoints-1)
             {
               p2 = mitk::imv::GetItkPoint(points->GetPoint(segment_id+1));
             }
             else
             {
               p2 = mitk::imv::GetItkPoint(points->GetPoint(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];
             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);
             }
           }
 
           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 >
-void TractParcellationFilter< OutImageType >::GenerateData()
+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);
 
   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 9bfba01..15f3875 100644
--- a/Modules/FiberTracking/Algorithms/itkTractParcellationFilter.h
+++ b/Modules/FiberTracking/Algorithms/itkTractParcellationFilter.h
@@ -1,96 +1,96 @@
 /*===================================================================
 
 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> >
+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 OutImageType::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);
 
   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 OutImageType::Pointer    m_InputImage;
+  typename InputImageType::Pointer  m_InputImage;
 };
 
 }
 
 #ifndef ITK_MANUAL_INSTANTIATION
 #include "itkTractParcellationFilter.cpp"
 #endif
 
 #endif // __itkTractParcellationFilter_h__