diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.cpp
index 422bd35376..0e143d3a09 100644
--- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.cpp
@@ -1,467 +1,459 @@
 #include "itkFitFibersToImageFilter.h"
 
 #include <boost/progress.hpp>
 
 namespace itk{
 
 FitFibersToImageFilter::FitFibersToImageFilter()
   : m_PeakImage(nullptr)
   , m_MaskImage(nullptr)
   , m_FitIndividualFibers(true)
   , m_GradientTolerance(1e-5)
   , m_Lambda(0.1)
   , m_MaxIterations(20)
   , m_FiberSampling(10)
   , m_Coverage(0)
   , m_Overshoot(0)
+  , m_RMSE(0.0)
   , m_FilterOutliers(true)
   , m_MeanWeight(1.0)
   , m_MedianWeight(1.0)
   , m_MinWeight(1.0)
   , m_MaxWeight(1.0)
   , m_Verbose(true)
   , m_DeepCopy(true)
   , m_ResampleFibers(true)
   , m_NumUnknowns(0)
   , m_NumResiduals(0)
   , m_NumCoveredDirections(0)
 {
   this->SetNumberOfRequiredOutputs(3);
 }
 
 FitFibersToImageFilter::~FitFibersToImageFilter()
 {
 
 }
 
 void FitFibersToImageFilter::GenerateData()
 {
   int sz_x = m_PeakImage->GetLargestPossibleRegion().GetSize(0);
   int sz_y = m_PeakImage->GetLargestPossibleRegion().GetSize(1);
   int sz_z = m_PeakImage->GetLargestPossibleRegion().GetSize(2);
   int sz_peaks = m_PeakImage->GetLargestPossibleRegion().GetSize(3)/3 + 1; // +1 for zero - peak
   int num_voxels = sz_x*sz_y*sz_z;
 
   float minSpacing = 1;
   if(m_PeakImage->GetSpacing()[0]<m_PeakImage->GetSpacing()[1] && m_PeakImage->GetSpacing()[0]<m_PeakImage->GetSpacing()[2])
     minSpacing = m_PeakImage->GetSpacing()[0];
   else if (m_PeakImage->GetSpacing()[1] < m_PeakImage->GetSpacing()[2])
     minSpacing = m_PeakImage->GetSpacing()[1];
   else
     minSpacing = m_PeakImage->GetSpacing()[2];
 
   m_NumUnknowns = m_Tractograms.size();
   if (m_FitIndividualFibers)
   {
     m_NumUnknowns = 0;
     for (unsigned int bundle=0; bundle<m_Tractograms.size(); bundle++)
       m_NumUnknowns += m_Tractograms.at(bundle)->GetNumFibers();
   }
 
   if (m_ResampleFibers)
     for (unsigned int bundle=0; bundle<m_Tractograms.size(); bundle++)
     {
       std::streambuf *old = cout.rdbuf(); // <-- save
       std::stringstream ss;
       std::cout.rdbuf (ss.rdbuf());
       if (m_DeepCopy)
         m_Tractograms[bundle] = m_Tractograms.at(bundle)->GetDeepCopy();
       m_Tractograms.at(bundle)->ResampleLinear(minSpacing/m_FiberSampling);
       std::cout.rdbuf (old);
     }
 
   m_NumResiduals = num_voxels * sz_peaks;
 
   MITK_INFO << "Num. unknowns: " << m_NumUnknowns;
   MITK_INFO << "Num. residuals: " << m_NumResiduals;
   MITK_INFO << "Creating system ...";
 
   vnl_sparse_matrix<double> A;
   vnl_vector<double> b;
   A.set_size(m_NumResiduals, m_NumUnknowns);
   b.set_size(m_NumResiduals); b.fill(0.0);
 
   double TD = 0;
   double FD = 0;
   m_NumCoveredDirections = 0;
   unsigned int fiber_count = 0;
 
   for (unsigned int bundle=0; bundle<m_Tractograms.size(); bundle++)
   {
     vtkSmartPointer<vtkPolyData> polydata = m_Tractograms.at(bundle)->GetFiberPolyData();
 
     for (int i=0; i<m_Tractograms.at(bundle)->GetNumFibers(); ++i)
     {
       vtkCell* cell = polydata->GetCell(i);
       int numPoints = cell->GetNumberOfPoints();
       vtkPoints* points = cell->GetPoints();
 
       if (numPoints<2)
         MITK_INFO << "FIBER WITH ONLY ONE POINT ENCOUNTERED!";
 
       for (int j=0; j<numPoints-1; ++j)
       {
         double* p1 = points->GetPoint(j);
         PointType4 p;
         p[0]=p1[0];
         p[1]=p1[1];
         p[2]=p1[2];
         p[3]=0;
 
         itk::Index<4> idx4;
         m_PeakImage->TransformPhysicalPointToIndex(p, idx4);
         itk::Index<3> idx3; idx3[0] = idx4[0]; idx3[1] = idx4[1]; idx3[2] = idx4[2];
         if (!m_PeakImage->GetLargestPossibleRegion().IsInside(idx4) || (m_MaskImage.IsNotNull() && m_MaskImage->GetPixel(idx3)==0))
           continue;
 
         double* p2 = points->GetPoint(j+1);
         vnl_vector_fixed<float,3> fiber_dir;
         fiber_dir[0] = p[0]-p2[0];
         fiber_dir[1] = p[1]-p2[1];
         fiber_dir[2] = p[2]-p2[2];
         fiber_dir.normalize();
 
         double w = 1;
         int peak_id = sz_peaks-1;
         vnl_vector_fixed<float,3> odf_peak = GetClosestPeak(idx4, m_PeakImage, fiber_dir, peak_id, w);
         float peak_mag = odf_peak.magnitude();
 
         int x = idx4[0];
         int y = idx4[1];
         int z = idx4[2];
 
         unsigned int linear_index = x + sz_x*y + sz_x*sz_y*z + sz_x*sz_y*sz_z*peak_id;
 
         if (b[linear_index] == 0 && peak_id<3)
         {
           m_NumCoveredDirections++;
           FD += peak_mag;
         }
         TD += w;
 
         if (m_FitIndividualFibers)
         {
           b[linear_index] = peak_mag;
           A.put(linear_index, fiber_count, A.get(linear_index, fiber_count) + w);
         }
         else
         {
           b[linear_index] = peak_mag;
           A.put(linear_index, bundle, A.get(linear_index, bundle) + w);
         }
       }
 
       ++fiber_count;
     }
   }
 
   TD /= (m_NumCoveredDirections*fiber_count);
   FD /= m_NumCoveredDirections;
   A /= TD;
   b *= 100.0/FD;  // times 100 because we want to avoid too small values for computational reasons
 
   double init_lambda = fiber_count;  // initialization for lambda estimation
 
   itk::TimeProbe clock;
   clock.Start();
 
   VnlCostFunction cost(m_NumUnknowns);
   cost.SetProblem(A, b, init_lambda);
   m_Weights.set_size(m_NumUnknowns); // m_Weights.fill( TD/100.0 * FD/2.0 );
   m_Weights.fill( 0.0 );
   vnl_lbfgsb minimizer(cost);
   vnl_vector<double> l; l.set_size(m_NumUnknowns); l.fill(0);
   vnl_vector<long> bound_selection; bound_selection.set_size(m_NumUnknowns); bound_selection.fill(1);
   minimizer.set_bound_selection(bound_selection);
   minimizer.set_lower_bound(l);
   minimizer.set_projected_gradient_tolerance(m_GradientTolerance);
 
   MITK_INFO << "Estimating regularization";
   minimizer.set_trace(false);
   minimizer.set_max_function_evals(2);
   minimizer.minimize(m_Weights);
   vnl_vector<double> dx; dx.set_size(m_NumUnknowns); dx.fill(0.0);
   cost.grad_regu_localMSE(m_Weights, dx);
   double r = dx.magnitude()/m_Weights.magnitude();
   cost.m_Lambda *= m_Lambda*55.0/r;
   MITK_INFO << r << " - " << m_Lambda*55.0/r;
   if (cost.m_Lambda>10e7)
   {
     MITK_INFO << "Regularization estimation failed. Using default value.";
     cost.m_Lambda = fiber_count;
   }
   MITK_INFO << "Using regularization factor of " << cost.m_Lambda << " (λ: " << m_Lambda << ")";
 
   MITK_INFO << "Fitting fibers";
   minimizer.set_trace(m_Verbose);
 
   minimizer.set_max_function_evals(m_MaxIterations);
   minimizer.minimize(m_Weights);
 
   std::vector< double > weights;
   if (m_FilterOutliers)
   {
     for (auto w : m_Weights)
       weights.push_back(w);
     std::sort(weights.begin(), weights.end());
     MITK_INFO << "Setting upper weight bound to " << weights.at(m_NumUnknowns*0.99);
     vnl_vector<double> u; u.set_size(m_NumUnknowns); u.fill(weights.at(m_NumUnknowns*0.99));
     minimizer.set_upper_bound(u);
     bound_selection.fill(2);
     minimizer.set_bound_selection(bound_selection);
     minimizer.minimize(m_Weights);
     weights.clear();
   }
 
   for (auto w : m_Weights)
     weights.push_back(w);
   std::sort(weights.begin(), weights.end());
 
   m_MeanWeight = m_Weights.mean();
   m_MedianWeight = weights.at(m_NumUnknowns*0.5);
   m_MinWeight = weights.at(0);
   m_MaxWeight = weights.at(m_NumUnknowns-1);
 
   MITK_INFO << "*************************";
   MITK_INFO << "Weight statistics";
   MITK_INFO << "Mean: " << m_MeanWeight;
   MITK_INFO << "Median: " << m_MedianWeight;
   MITK_INFO << "75% quantile: " << weights.at(m_NumUnknowns*0.75);
   MITK_INFO << "95% quantile: " << weights.at(m_NumUnknowns*0.95);
   MITK_INFO << "99% quantile: " << weights.at(m_NumUnknowns*0.99);
   MITK_INFO << "Min: " << m_MinWeight;
   MITK_INFO << "Max: " << m_MaxWeight;
   MITK_INFO << "*************************";
   MITK_INFO << "NumEvals: " << minimizer.get_num_evaluations();
   MITK_INFO << "NumIterations: " << minimizer.get_num_iterations();
   MITK_INFO << "Residual cost: " << minimizer.get_end_error();
-  double rms = cost.S->get_rms_error(m_Weights);
-  MITK_INFO << "Final RMS: " << rms;
+  m_RMSE = cost.S->get_rms_error(m_Weights);
+  MITK_INFO << "Final RMS: " << m_RMSE;
 
   clock.Stop();
   int h = clock.GetTotal()/3600;
   int m = ((int)clock.GetTotal()%3600)/60;
   int s = (int)clock.GetTotal()%60;
   MITK_INFO << "Optimization took " << h << "h, " << m << "m and " << s << "s";
 
   MITK_INFO << "Weighting fibers";
   m_RmsDiffPerFiber.set_size(m_Weights.size());
   m_RmsDiffPerBundle.set_size(m_Tractograms.size());
   std::streambuf *old = cout.rdbuf(); // <-- save
   std::stringstream ss;
   std::cout.rdbuf (ss.rdbuf());
   if (m_FitIndividualFibers)
   {
     unsigned int fiber_count = 0;
     for (unsigned int bundle=0; bundle<m_Tractograms.size(); bundle++)
     {
-//      double mean_d = 0;
       vnl_vector< double > temp_weights;
       temp_weights.set_size(m_Weights.size());
       temp_weights.copy_in(m_Weights.data_block());
 
       for (int i=0; i<m_Tractograms.at(bundle)->GetNumFibers(); i++)
       {
         m_Tractograms.at(bundle)->SetFiberWeight(i, m_Weights[fiber_count]);
         temp_weights[fiber_count] = 0;
-//        double w = m_Weights[fiber_count];
-//        m_Weights[fiber_count] = 0;
-//        double d_rms = cost.S->get_rms_error(m_Weights) - rms;
-//        m_RmsDiffPerFiber[fiber_count] = d_rms;
-//        mean_d += d_rms;
-//        m_Weights[fiber_count] = w;
 
         ++fiber_count;
       }
-      double d_rms = cost.S->get_rms_error(temp_weights) - rms;
-      MITK_INFO << d_rms;
-//      mean_d /= m_Tractograms.at(bundle)->GetNumFibers();
+      double d_rms = cost.S->get_rms_error(temp_weights) - m_RMSE;
       m_RmsDiffPerBundle[bundle] = d_rms;
       m_Tractograms.at(bundle)->Compress(0.1);
       m_Tractograms.at(bundle)->ColorFibersByFiberWeights(false, true);
     }
   }
   else
   {
     for (unsigned int i=0; i<m_Tractograms.size(); ++i)
     {
       m_Tractograms.at(i)->SetFiberWeights(m_Weights[i]);
       m_Tractograms.at(i)->Compress(0.1);
       m_Tractograms.at(i)->ColorFibersByFiberWeights(false, true);
     }
   }
   std::cout.rdbuf (old);
 
   // transform back for peak image creation
   A *= FD/100.0;
   b *= FD/100.0;
 
   MITK_INFO << "Generating output images ...";
 
   itk::ImageDuplicator< PeakImgType >::Pointer duplicator = itk::ImageDuplicator< PeakImgType >::New();
   duplicator->SetInputImage(m_PeakImage);
   duplicator->Update();
   m_UnderexplainedImage = duplicator->GetOutput();
   m_UnderexplainedImage->FillBuffer(0.0);
 
   duplicator->SetInputImage(m_UnderexplainedImage);
   duplicator->Update();
   m_OverexplainedImage = duplicator->GetOutput();
   m_OverexplainedImage->FillBuffer(0.0);
 
   duplicator->SetInputImage(m_OverexplainedImage);
   duplicator->Update();
   m_ResidualImage = duplicator->GetOutput();
   m_ResidualImage->FillBuffer(0.0);
 
   duplicator->SetInputImage(m_ResidualImage);
   duplicator->Update();
   m_FittedImage = duplicator->GetOutput();
   m_FittedImage->FillBuffer(0.0);
 
   vnl_vector<double> fitted_b; fitted_b.set_size(b.size());
   cost.S->multiply(m_Weights, fitted_b);
 
   for (unsigned int r=0; r<b.size(); r++)
   {
     itk::Index<4> idx4;
     unsigned int linear_index = r;
     idx4[0] = linear_index % sz_x; linear_index /= sz_x;
     idx4[1] = linear_index % sz_y; linear_index /= sz_y;
     idx4[2] = linear_index % sz_z; linear_index /= sz_z;
     int peak_id = linear_index % sz_peaks;
 
     if (peak_id<sz_peaks-1)
     {
       vnl_vector_fixed<float,3> peak_dir;
 
       idx4[3] = peak_id*3;
       peak_dir[0] = m_PeakImage->GetPixel(idx4);
       idx4[3] += 1;
       peak_dir[1] = m_PeakImage->GetPixel(idx4);
       idx4[3] += 1;
       peak_dir[2] = m_PeakImage->GetPixel(idx4);
 
       peak_dir.normalize();
       peak_dir *= fitted_b[r];
 
       idx4[3] = peak_id*3;
       m_FittedImage->SetPixel(idx4, peak_dir[0]);
 
       idx4[3] += 1;
       m_FittedImage->SetPixel(idx4, peak_dir[1]);
 
       idx4[3] += 1;
       m_FittedImage->SetPixel(idx4, peak_dir[2]);
     }
   }
 
   FD = 0;
   m_Coverage = 0;
   m_Overshoot = 0;
 
   itk::Index<4> idx4;
   for (idx4[0]=0; idx4[0]<sz_x; ++idx4[0])
     for (idx4[1]=0; idx4[1]<sz_y; ++idx4[1])
       for (idx4[2]=0; idx4[2]<sz_z; ++idx4[2])
       {
         itk::Index<3> idx3; idx3[0] = idx4[0]; idx3[1] = idx4[1]; idx3[2] = idx4[2];
         if (m_MaskImage.IsNotNull() && m_MaskImage->GetPixel(idx3)==0)
           continue;
 
         vnl_vector_fixed<float,3> peak_dir;
         vnl_vector_fixed<float,3> fitted_dir;
         vnl_vector_fixed<float,3> overshoot_dir;
         for (idx4[3]=0; idx4[3]<(itk::IndexValueType)m_PeakImage->GetLargestPossibleRegion().GetSize(3); ++idx4[3])
         {
           peak_dir[idx4[3]%3] = m_PeakImage->GetPixel(idx4);
           fitted_dir[idx4[3]%3] = m_FittedImage->GetPixel(idx4);
           m_ResidualImage->SetPixel(idx4, m_PeakImage->GetPixel(idx4) - m_FittedImage->GetPixel(idx4));
 
           if (idx4[3]%3==2)
           {
             FD += peak_dir.magnitude();
 
             itk::Index<4> tidx= idx4;
             if (peak_dir.magnitude()>fitted_dir.magnitude())
             {
               m_Coverage += fitted_dir.magnitude();
               m_UnderexplainedImage->SetPixel(tidx, peak_dir[2]-fitted_dir[2]); tidx[3]--;
               m_UnderexplainedImage->SetPixel(tidx, peak_dir[1]-fitted_dir[1]); tidx[3]--;
               m_UnderexplainedImage->SetPixel(tidx, peak_dir[0]-fitted_dir[0]);
             }
             else
             {
               overshoot_dir[0] = fitted_dir[0]-peak_dir[0];
               overshoot_dir[1] = fitted_dir[1]-peak_dir[1];
               overshoot_dir[2] = fitted_dir[2]-peak_dir[2];
               m_Coverage += peak_dir.magnitude();
               m_Overshoot += overshoot_dir.magnitude();
               m_OverexplainedImage->SetPixel(tidx, overshoot_dir[2]); tidx[3]--;
               m_OverexplainedImage->SetPixel(tidx, overshoot_dir[1]); tidx[3]--;
               m_OverexplainedImage->SetPixel(tidx, overshoot_dir[0]);
             }
           }
         }
       }
 
   m_Coverage = m_Coverage/FD;
   m_Overshoot = m_Overshoot/FD;
 
   MITK_INFO << std::fixed << "Coverage: " << setprecision(2) << 100.0*m_Coverage << "%";
   MITK_INFO << std::fixed << "Overshoot: " << setprecision(2) << 100.0*m_Overshoot << "%";
 }
 
 vnl_vector_fixed<float,3> FitFibersToImageFilter::GetClosestPeak(itk::Index<4> idx, PeakImgType::Pointer peak_image , vnl_vector_fixed<float,3> fiber_dir, int& id, double& w )
 {
   int m_NumDirs = peak_image->GetLargestPossibleRegion().GetSize()[3]/3;
   vnl_vector_fixed<float,3> out_dir; out_dir.fill(0);
   float angle = 0.9;
 
   for (int i=0; i<m_NumDirs; i++)
   {
     vnl_vector_fixed<float,3> dir;
     idx[3] = i*3;
     dir[0] = peak_image->GetPixel(idx);
     idx[3] += 1;
     dir[1] = peak_image->GetPixel(idx);
     idx[3] += 1;
     dir[2] = peak_image->GetPixel(idx);
 
     float mag = dir.magnitude();
     if (mag<mitk::eps)
       continue;
 
     dir.normalize();
 
     float a = dot_product(dir, fiber_dir);
     if (fabs(a)>angle)
     {
       angle = fabs(a);
       w = angle;
       if (a<0)
         out_dir = -dir;
       else
         out_dir = dir;
       out_dir *= mag;
       id = i;
     }
   }
 
   return out_dir;
 }
 
 std::vector<mitk::FiberBundle::Pointer> FitFibersToImageFilter::GetTractograms() const
 {
   return m_Tractograms;
 }
 
 void FitFibersToImageFilter::SetTractograms(const std::vector<mitk::FiberBundle::Pointer> &tractograms)
 {
   m_Tractograms = tractograms;
 }
 
 }
 
 
 
diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.h
index 7a8a9be830..5907d9e874 100644
--- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.h
+++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.h
@@ -1,269 +1,271 @@
 #ifndef __itkFitFibersToImageFilter_h__
 #define __itkFitFibersToImageFilter_h__
 
 // MITK
 #include <mitkFiberBundle.h>
 #include <itkImageSource.h>
 #include <mitkPeakImage.h>
 #include <vnl/algo/vnl_lbfgsb.h>
 #include <vnl/vnl_sparse_matrix.h>
 #include <vnl/vnl_sparse_matrix_linear_system.h>
 #include <itkImageDuplicator.h>
 #include <itkTimeProbe.h>
 #include <itkMersenneTwisterRandomVariateGenerator.h>
 
 namespace itk{
 
 /**
 * \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).  */
 
 class FitFibersToImageFilter : public ImageSource< mitk::PeakImage::ItkPeakImageType >
 {
 
 public:
 
   typedef FitFibersToImageFilter Self;
   typedef ProcessObject Superclass;
   typedef SmartPointer< Self > Pointer;
   typedef SmartPointer< const Self > ConstPointer;
 
   typedef itk::Point<float, 4> PointType4;
   typedef mitk::PeakImage::ItkPeakImageType       PeakImgType;
   typedef itk::Image<unsigned char, 3>            UcharImgType;
 
   itkFactorylessNewMacro(Self)
   itkCloneMacro(Self)
   itkTypeMacro( FitFibersToImageFilter, ImageSource )
 
   itkSetMacro( PeakImage, PeakImgType::Pointer)
   itkGetMacro( PeakImage, PeakImgType::Pointer)
   itkSetMacro( MaskImage, UcharImgType::Pointer)
   itkGetMacro( MaskImage, UcharImgType::Pointer)
   itkSetMacro( FitIndividualFibers, bool)
   itkGetMacro( FitIndividualFibers, bool)
   itkSetMacro( GradientTolerance, double)
   itkGetMacro( GradientTolerance, double)
   itkSetMacro( Lambda, double)
   itkGetMacro( Lambda, double)
   itkSetMacro( MaxIterations, int)
   itkGetMacro( MaxIterations, int)
   itkSetMacro( FiberSampling, float)
   itkGetMacro( FiberSampling, float)
   itkSetMacro( FilterOutliers, bool)
   itkGetMacro( FilterOutliers, bool)
   itkSetMacro( Verbose, bool)
   itkGetMacro( Verbose, bool)
   itkSetMacro( DeepCopy, bool)
   itkGetMacro( DeepCopy, bool)
   itkSetMacro( ResampleFibers, bool)
   itkGetMacro( ResampleFibers, bool)
 
   itkGetMacro( Weights, vnl_vector<double>)
   itkGetMacro( RmsDiffPerBundle, vnl_vector<double>)
   itkGetMacro( RmsDiffPerFiber, vnl_vector<double>)
   itkGetMacro( FittedImage, PeakImgType::Pointer)
   itkGetMacro( ResidualImage, PeakImgType::Pointer)
   itkGetMacro( OverexplainedImage, PeakImgType::Pointer)
   itkGetMacro( UnderexplainedImage, PeakImgType::Pointer)
   itkGetMacro( Coverage, double)
   itkGetMacro( Overshoot, double)
+  itkGetMacro( RMSE, double)
   itkGetMacro( MeanWeight, double)
   itkGetMacro( MedianWeight, double)
   itkGetMacro( MinWeight, double)
   itkGetMacro( MaxWeight, double)
   itkGetMacro( NumUnknowns, unsigned int)
   itkGetMacro( NumResiduals, unsigned int)
   itkGetMacro( NumCoveredDirections, unsigned int)
 
   void SetTractograms(const std::vector<mitk::FiberBundle::Pointer> &tractograms);
 
   void GenerateData() override;
 
   std::vector<mitk::FiberBundle::Pointer> GetTractograms() const;
 
 protected:
 
   FitFibersToImageFilter();
   virtual ~FitFibersToImageFilter();
 
   vnl_vector_fixed<float,3> GetClosestPeak(itk::Index<4> idx, PeakImgType::Pointer m_PeakImage , vnl_vector_fixed<float,3> fiber_dir, int& id, double& w );
 
   std::vector< mitk::FiberBundle::Pointer >   m_Tractograms;
   PeakImgType::Pointer                        m_PeakImage;
   UcharImgType::Pointer                       m_MaskImage;
   bool                                        m_FitIndividualFibers;
   double                                      m_GradientTolerance;
   double                                      m_Lambda;
   int                                         m_MaxIterations;
   float                                       m_FiberSampling;
   double                                      m_Coverage;
   double                                      m_Overshoot;
+  double                                      m_RMSE;
   bool                                        m_FilterOutliers;
   double                                      m_MeanWeight;
   double                                      m_MedianWeight;
   double                                      m_MinWeight;
   double                                      m_MaxWeight;
   bool                                        m_Verbose;
   bool                                        m_DeepCopy;
   bool                                        m_ResampleFibers;
   unsigned int                                m_NumUnknowns;
   unsigned int                                m_NumResiduals;
   unsigned int                                m_NumCoveredDirections;
 
   // output
   vnl_vector<double>                          m_RmsDiffPerBundle;
   vnl_vector<double>                          m_Weights;
   vnl_vector<double>                          m_RmsDiffPerFiber;
   PeakImgType::Pointer                        m_UnderexplainedImage;
   PeakImgType::Pointer                        m_OverexplainedImage;
   PeakImgType::Pointer                        m_ResidualImage;
   PeakImgType::Pointer                        m_FittedImage;
 };
 
 }
 
 
 class VnlCostFunction : public vnl_cost_function
 {
 public:
 
   vnl_sparse_matrix_linear_system< double >* S;
   vnl_sparse_matrix< double > m_A;
   vnl_sparse_matrix< double > m_A_Ones; // matrix indicating active weights with 1
   vnl_vector< double > m_b;
   double m_Lambda;  // regularization factor
 
   vnl_vector<double> row_sums;  // number of active weights per row
   vnl_vector<double> local_weight_means;  // mean weight of each row
 
   void SetProblem(vnl_sparse_matrix< double >& A, vnl_vector<double>& b, double lambda)
   {
     S = new vnl_sparse_matrix_linear_system<double>(A, b);
     m_A = A;
     m_b = b;
     m_Lambda = lambda;
 
     m_A_Ones.set_size(m_A.rows(), m_A.cols());
     m_A.reset();
     while (m_A.next())
       m_A_Ones.put(m_A.getrow(), m_A.getcolumn(), 1);
 
     unsigned int N = m_b.size();
     vnl_vector<double> ones; ones.set_size(dim); ones.fill(1.0);
     row_sums.set_size(N);
     m_A_Ones.mult(ones, row_sums);
     local_weight_means.set_size(N);
   }
 
   VnlCostFunction(const int NumVars) : vnl_cost_function(NumVars)
   {
   }
 
   void regu_MSE(vnl_vector<double> const &x, double& cost)
   {
     double mean = x.mean();
     vnl_vector<double> tx = x-mean;
     cost += m_Lambda*1e8*tx.squared_magnitude()/x.size();
   }
 
   void regu_MSM(vnl_vector<double> const &x, double& cost)
   {
     cost += m_Lambda*1e8*x.squared_magnitude()/x.size();
   }
 
   void regu_localMSE(vnl_vector<double> const &x, double& cost)
   {
     m_A_Ones.mult(x, local_weight_means);
     local_weight_means = element_quotient(local_weight_means, row_sums);
 
     m_A_Ones.reset();
     double regu = 0;
     while (m_A_Ones.next())
     {
       double d = 0;
       if (x[m_A_Ones.getcolumn()]>local_weight_means[m_A_Ones.getrow()])
         d = std::exp(x[m_A_Ones.getcolumn()]) - std::exp(local_weight_means[m_A_Ones.getrow()]);
       else
         d = x[m_A_Ones.getcolumn()] - local_weight_means[m_A_Ones.getrow()];
       regu += d*d;
     }
     cost += m_Lambda*regu/dim;
   }
 
   void grad_regu_MSE(vnl_vector<double> const &x, vnl_vector<double> &dx)
   {
     double mean = x.mean();
     vnl_vector<double> tx = x-mean;
 
     vnl_vector<double> tx2(dim, 0.0);
     vnl_vector<double> h(dim, 1.0);
     for (int c=0; c<dim; c++)
     {
       h[c] = dim-1;
       tx2[c] += dot_product(h,tx);
       h[c] = 1;
     }
     dx += tx2*m_Lambda*1e8*2.0/(dim*dim);
 
   }
 
   void grad_regu_MSM(vnl_vector<double> const &x, vnl_vector<double> &dx)
   {
     dx += m_Lambda*1e8*2.0*x/dim;
   }
 
   void grad_regu_localMSE(vnl_vector<double> const &x, vnl_vector<double> &dx)
   {
     m_A_Ones.mult(x, local_weight_means);
     local_weight_means = element_quotient(local_weight_means, row_sums);
 
     vnl_vector<double> exp_x = x.apply(std::exp);
     vnl_vector<double> exp_means = local_weight_means.apply(std::exp);
 
     vnl_vector<double> tdx(dim, 0);
     m_A_Ones.reset();
     while (m_A_Ones.next())
     {
       int c = m_A_Ones.getcolumn();
       int r = m_A_Ones.getrow();
 
       if (x[c]>local_weight_means[r])
         tdx[c] += exp_x[c] * ( exp_x[c] - exp_means[r] );
       else
         tdx[c] += x[c] - local_weight_means[r];
     }
     dx += tdx*2.0*m_Lambda/dim;
   }
 
 
   double f(vnl_vector<double> const &x)
   {
     double cost = S->get_rms_error(x);
     cost *= cost;
 
     regu_localMSE(x, cost);
 
     return cost;
   }
 
   void gradf(vnl_vector<double> const &x, vnl_vector<double> &dx)
   {
     dx.fill(0.0);
     unsigned int N = m_b.size();
 
     vnl_vector<double> d; d.set_size(N);
     S->multiply(x,d);
     d -= m_b;
 
     S->transpose_multiply(d, dx);
     dx *= 2.0/N;
 
     grad_regu_localMSE(x,dx);
   }
 };
 
 #ifndef ITK_MANUAL_INSTANTIATION
 #include "itkFitFibersToImageFilter.cpp"
 #endif
 
 #endif // __itkFitFibersToImageFilter_h__
diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp
index d2e8e5b73a..fedc54a446 100755
--- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp
@@ -1,2584 +1,2588 @@
 /*===================================================================
 
 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.
 
 ===================================================================*/
 
 
 #define _USE_MATH_DEFINES
 #include "mitkFiberBundle.h"
 
 #include <mitkPlanarCircle.h>
 #include <mitkPlanarPolygon.h>
 #include <mitkPlanarFigureComposite.h>
 #include "mitkImagePixelReadAccessor.h"
 #include <mitkPixelTypeMultiplex.h>
 
 #include <vtkPointData.h>
 #include <vtkDataArray.h>
 #include <vtkUnsignedCharArray.h>
 #include <vtkPolyLine.h>
 #include <vtkCellArray.h>
 #include <vtkCellData.h>
 #include <vtkIdFilter.h>
 #include <vtkClipPolyData.h>
 #include <vtkPlane.h>
 #include <vtkDoubleArray.h>
 #include <vtkKochanekSpline.h>
 #include <vtkParametricFunctionSource.h>
 #include <vtkParametricSpline.h>
 #include <vtkPolygon.h>
 #include <vtkCleanPolyData.h>
 #include <cmath>
 #include <boost/progress.hpp>
 #include <vtkTransformPolyDataFilter.h>
 #include <mitkTransferFunction.h>
 #include <vtkLookupTable.h>
 #include <mitkLookupTable.h>
 #include <vtkCardinalSpline.h>
 #include <vtkAppendPolyData.h>
 
 const char* mitk::FiberBundle::FIBER_ID_ARRAY = "Fiber_IDs";
 
 using namespace std;
 
 mitk::FiberBundle::FiberBundle( vtkPolyData* fiberPolyData )
   : m_NumFibers(0)
 {
   m_FiberWeights = vtkSmartPointer<vtkFloatArray>::New();
   m_FiberWeights->SetName("FIBER_WEIGHTS");
 
   m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
   if (fiberPolyData != nullptr)
     m_FiberPolyData = fiberPolyData;
   else
   {
     this->m_FiberPolyData->SetPoints(vtkSmartPointer<vtkPoints>::New());
     this->m_FiberPolyData->SetLines(vtkSmartPointer<vtkCellArray>::New());
   }
 
   this->UpdateFiberGeometry();
   this->GenerateFiberIds();
   this->ColorFibersByOrientation();
 }
 
 mitk::FiberBundle::~FiberBundle()
 {
 
 }
 
 mitk::FiberBundle::Pointer mitk::FiberBundle::GetDeepCopy()
 {
   mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(m_FiberPolyData);
   newFib->SetFiberColors(this->m_FiberColors);
   newFib->SetFiberWeights(this->m_FiberWeights);
   return newFib;
 }
 
 vtkSmartPointer<vtkPolyData> mitk::FiberBundle::GeneratePolyDataByIds(std::vector<long> fiberIds, vtkSmartPointer<vtkFloatArray> weights)
 {
   vtkSmartPointer<vtkPolyData> newFiberPolyData = vtkSmartPointer<vtkPolyData>::New();
   vtkSmartPointer<vtkCellArray> newLineSet = vtkSmartPointer<vtkCellArray>::New();
   vtkSmartPointer<vtkPoints> newPointSet = vtkSmartPointer<vtkPoints>::New();
   weights->SetNumberOfValues(fiberIds.size());
 
   int counter = 0;
   auto finIt = fiberIds.begin();
   while ( finIt != fiberIds.end() )
   {
     if (*finIt < 0 || *finIt>GetNumFibers()){
       MITK_INFO << "FiberID can not be negative or >NumFibers!!! check id Extraction!" << *finIt;
       break;
     }
 
     vtkSmartPointer<vtkCell> fiber = m_FiberIdDataSet->GetCell(*finIt);//->DeepCopy(fiber);
     vtkSmartPointer<vtkPoints> fibPoints = fiber->GetPoints();
     vtkSmartPointer<vtkPolyLine> newFiber = vtkSmartPointer<vtkPolyLine>::New();
     newFiber->GetPointIds()->SetNumberOfIds( fibPoints->GetNumberOfPoints() );
 
     for(int i=0; i<fibPoints->GetNumberOfPoints(); i++)
     {
       newFiber->GetPointIds()->SetId(i, newPointSet->GetNumberOfPoints());
       newPointSet->InsertNextPoint(fibPoints->GetPoint(i)[0], fibPoints->GetPoint(i)[1], fibPoints->GetPoint(i)[2]);
     }
 
     weights->InsertValue(counter, this->GetFiberWeight(*finIt));
     newLineSet->InsertNextCell(newFiber);
     ++finIt;
     ++counter;
   }
 
   newFiberPolyData->SetPoints(newPointSet);
   newFiberPolyData->SetLines(newLineSet);
   return newFiberPolyData;
 }
 
 // merge two fiber bundles
 mitk::FiberBundle::Pointer mitk::FiberBundle::AddBundles(std::vector< mitk::FiberBundle::Pointer > fibs)
 {
   vtkSmartPointer<vtkPolyData> vNewPolyData = vtkSmartPointer<vtkPolyData>::New();
   vtkSmartPointer<vtkCellArray> vNewLines = vtkSmartPointer<vtkCellArray>::New();
   vtkSmartPointer<vtkPoints> vNewPoints = vtkSmartPointer<vtkPoints>::New();
 
   // add current fiber bundle
   vtkSmartPointer<vtkFloatArray> weights = vtkSmartPointer<vtkFloatArray>::New();
 
   int num_weights = this->GetNumFibers();
   for (auto fib : fibs)
     num_weights += fib->GetNumFibers();
   weights->SetNumberOfValues(num_weights);
 
   unsigned int counter = 0;
   for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     for (int j=0; j<numPoints; j++)
     {
       double p[3];
       points->GetPoint(j, p);
 
       vtkIdType id = vNewPoints->InsertNextPoint(p);
       container->GetPointIds()->InsertNextId(id);
     }
     weights->InsertValue(counter, this->GetFiberWeight(i));
     vNewLines->InsertNextCell(container);
     counter++;
   }
 
   for (auto fib : fibs)
   {
     // add new fiber bundle
     for (int i=0; i<fib->GetFiberPolyData()->GetNumberOfCells(); i++)
     {
       vtkCell* cell = fib->GetFiberPolyData()->GetCell(i);
       int numPoints = cell->GetNumberOfPoints();
       vtkPoints* points = cell->GetPoints();
 
       vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
       for (int j=0; j<numPoints; j++)
       {
         double p[3];
         points->GetPoint(j, p);
 
         vtkIdType id = vNewPoints->InsertNextPoint(p);
         container->GetPointIds()->InsertNextId(id);
       }
       weights->InsertValue(counter, fib->GetFiberWeight(i));
       vNewLines->InsertNextCell(container);
       counter++;
     }
   }
 
   // initialize PolyData
   vNewPolyData->SetPoints(vNewPoints);
   vNewPolyData->SetLines(vNewLines);
 
   // initialize fiber bundle
   mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(vNewPolyData);
   newFib->SetFiberWeights(weights);
   return newFib;
 }
 
 // merge two fiber bundles
 mitk::FiberBundle::Pointer mitk::FiberBundle::AddBundle(mitk::FiberBundle* fib)
 {
   if (fib==nullptr)
     return this->GetDeepCopy();
 
   MITK_INFO << "Adding fibers";
 
   vtkSmartPointer<vtkPolyData> vNewPolyData = vtkSmartPointer<vtkPolyData>::New();
   vtkSmartPointer<vtkCellArray> vNewLines = vtkSmartPointer<vtkCellArray>::New();
   vtkSmartPointer<vtkPoints> vNewPoints = vtkSmartPointer<vtkPoints>::New();
 
   // add current fiber bundle
   vtkSmartPointer<vtkFloatArray> weights = vtkSmartPointer<vtkFloatArray>::New();
   weights->SetNumberOfValues(this->GetNumFibers()+fib->GetNumFibers());
 
   unsigned int counter = 0;
   for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     for (int j=0; j<numPoints; j++)
     {
       double p[3];
       points->GetPoint(j, p);
 
       vtkIdType id = vNewPoints->InsertNextPoint(p);
       container->GetPointIds()->InsertNextId(id);
     }
     weights->InsertValue(counter, this->GetFiberWeight(i));
     vNewLines->InsertNextCell(container);
     counter++;
   }
 
   // add new fiber bundle
   for (int i=0; i<fib->GetFiberPolyData()->GetNumberOfCells(); i++)
   {
     vtkCell* cell = fib->GetFiberPolyData()->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     for (int j=0; j<numPoints; j++)
     {
       double p[3];
       points->GetPoint(j, p);
 
       vtkIdType id = vNewPoints->InsertNextPoint(p);
       container->GetPointIds()->InsertNextId(id);
     }
     weights->InsertValue(counter, fib->GetFiberWeight(i));
     vNewLines->InsertNextCell(container);
     counter++;
   }
 
   // initialize PolyData
   vNewPolyData->SetPoints(vNewPoints);
   vNewPolyData->SetLines(vNewLines);
 
   // initialize fiber bundle
   mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(vNewPolyData);
   newFib->SetFiberWeights(weights);
   return newFib;
 }
 
 // Only retain fibers with a weight larger than the specified threshold
 mitk::FiberBundle::Pointer mitk::FiberBundle::FilterByWeights(float weight_thr, bool invert)
 {
   vtkSmartPointer<vtkPolyData> vNewPolyData = vtkSmartPointer<vtkPolyData>::New();
   vtkSmartPointer<vtkCellArray> vNewLines = vtkSmartPointer<vtkCellArray>::New();
   vtkSmartPointer<vtkPoints> vNewPoints = vtkSmartPointer<vtkPoints>::New();
 
   std::vector<float> weights;
 
   for (int i=0; i<this->GetNumFibers(); i++)
   {
     if ( (invert && this->GetFiberWeight(i)>weight_thr) || (!invert && this->GetFiberWeight(i)<=weight_thr))
       continue;
 
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     for (int j=0; j<numPoints; j++)
     {
       double p[3];
       points->GetPoint(j, p);
 
       vtkIdType id = vNewPoints->InsertNextPoint(p);
       container->GetPointIds()->InsertNextId(id);
     }
     vNewLines->InsertNextCell(container);
     weights.push_back(this->GetFiberWeight(i));
   }
 
   // initialize PolyData
   vNewPolyData->SetPoints(vNewPoints);
   vNewPolyData->SetLines(vNewLines);
 
   // initialize fiber bundle
   mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(vNewPolyData);
   for (unsigned int i=0; i<weights.size(); ++i)
     newFib->SetFiberWeight(i, weights.at(i));
   return newFib;
 }
 
 // Only retain a subsample of the fibers
 mitk::FiberBundle::Pointer mitk::FiberBundle::SubsampleFibers(float factor)
 {
   vtkSmartPointer<vtkPolyData> vNewPolyData = vtkSmartPointer<vtkPolyData>::New();
   vtkSmartPointer<vtkCellArray> vNewLines = vtkSmartPointer<vtkCellArray>::New();
   vtkSmartPointer<vtkPoints> vNewPoints = vtkSmartPointer<vtkPoints>::New();
 
   int new_num_fibs = this->GetNumFibers()*factor;
   MITK_INFO << "Subsampling fibers with factor " << factor << "(" << new_num_fibs << "/" << this->GetNumFibers() << ")";
 
   // add current fiber bundle
   vtkSmartPointer<vtkFloatArray> weights = vtkSmartPointer<vtkFloatArray>::New();
   weights->SetNumberOfValues(new_num_fibs);
 
   std::vector< int > ids;
   for (int i=0; i<this->GetNumFibers(); i++)
     ids.push_back(i);
   std::random_shuffle(ids.begin(), ids.end());
 
   unsigned int counter = 0;
   for (int i=0; i<new_num_fibs; i++)
   {
     vtkCell* cell = m_FiberPolyData->GetCell(ids.at(i));
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     for (int j=0; j<numPoints; j++)
     {
       double p[3];
       points->GetPoint(j, p);
 
       vtkIdType id = vNewPoints->InsertNextPoint(p);
       container->GetPointIds()->InsertNextId(id);
     }
     weights->InsertValue(counter, this->GetFiberWeight(ids.at(i)));
     vNewLines->InsertNextCell(container);
     counter++;
   }
 
   // initialize PolyData
   vNewPolyData->SetPoints(vNewPoints);
   vNewPolyData->SetLines(vNewLines);
 
   // initialize fiber bundle
   mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(vNewPolyData);
   newFib->SetFiberWeights(weights);
   return newFib;
 }
 
 // subtract two fiber bundles
 mitk::FiberBundle::Pointer mitk::FiberBundle::SubtractBundle(mitk::FiberBundle* fib)
 {
   if (fib==nullptr)
     return this->GetDeepCopy();
 
   MITK_INFO << "Subtracting fibers";
   vtkSmartPointer<vtkPolyData> vNewPolyData = vtkSmartPointer<vtkPolyData>::New();
   vtkSmartPointer<vtkCellArray> vNewLines = vtkSmartPointer<vtkCellArray>::New();
   vtkSmartPointer<vtkPoints> vNewPoints = vtkSmartPointer<vtkPoints>::New();
 
   std::vector< std::vector< itk::Point<float, 3> > > points1;
   for( int i=0; i<m_NumFibers; i++ )
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     if (points==nullptr || numPoints<=0)
       continue;
 
     itk::Point<float, 3> start = GetItkPoint(points->GetPoint(0));
     itk::Point<float, 3> end = GetItkPoint(points->GetPoint(numPoints-1));
 
     points1.push_back( {start, end} );
   }
 
   std::vector< std::vector< itk::Point<float, 3> > > points2;
   for( int i=0; i<fib->GetNumFibers(); i++ )
   {
     vtkCell* cell = fib->GetFiberPolyData()->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     if (points==nullptr || numPoints<=0)
       continue;
 
     itk::Point<float, 3> start = GetItkPoint(points->GetPoint(0));
     itk::Point<float, 3> end = GetItkPoint(points->GetPoint(numPoints-1));
 
     points2.push_back( {start, end} );
   }
 
 //  int progress = 0;
   std::vector< int > ids;
 #pragma omp parallel for
   for (int i=0; i<(int)points1.size(); i++)
   {
 //#pragma omp critical
 //    {
 //      progress++;
 //      std::cout << (int)(100*(float)progress/points1.size()) << "%" << '\r';
 //      cout.flush();
 //    }
 
     bool match = false;
     for (unsigned int j=0; j<points2.size(); j++)
     {
       auto v1 = points1.at(i);
       auto v2 = points2.at(j);
 
       float dist=0;
       for (unsigned int c=0; c<v1.size(); c++)
       {
         float d = v1[c][0]-v2[c][0];
         dist += d*d;
 
         d = v1[c][1]-v2[c][1];
         dist += d*d;
 
         d = v1[c][2]-v2[c][2];
         dist += d*d;
       }
       dist /= v1.size();
 
       if (dist<mitk::eps)
       {
         match = true;
         break;
       }
 
       dist=0;
       for (unsigned int c=0; c<v1.size(); c++)
       {
         float d = v1[v1.size()-1-c][0]-v2[c][0];
         dist += d*d;
 
         d = v1[v1.size()-1-c][1]-v2[c][1];
         dist += d*d;
 
         d = v1[v1.size()-1-c][2]-v2[c][2];
         dist += d*d;
       }
       dist /= v1.size();
 
       if (dist<mitk::eps)
       {
         match = true;
         break;
       }
     }
 
 #pragma omp critical
     if (!match)
       ids.push_back(i);
   }
 
   for( int i : ids )
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     if (points==nullptr || numPoints<=0)
       continue;
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     for( int j=0; j<numPoints; j++)
     {
       vtkIdType id = vNewPoints->InsertNextPoint(points->GetPoint(j));
       container->GetPointIds()->InsertNextId(id);
     }
     vNewLines->InsertNextCell(container);
   }
   if(vNewLines->GetNumberOfCells()==0)
     return mitk::FiberBundle::New();
   // initialize PolyData
   vNewPolyData->SetPoints(vNewPoints);
   vNewPolyData->SetLines(vNewLines);
 
   // initialize fiber bundle
   return mitk::FiberBundle::New(vNewPolyData);
 }
 
 itk::Point<float, 3> mitk::FiberBundle::GetItkPoint(double point[3])
 {
   itk::Point<float, 3> itkPoint;
   itkPoint[0] = point[0];
   itkPoint[1] = point[1];
   itkPoint[2] = point[2];
   return itkPoint;
 }
 
 /*
  * set PolyData (additional flag to recompute fiber geometry, default = true)
  */
 void mitk::FiberBundle::SetFiberPolyData(vtkSmartPointer<vtkPolyData> fiberPD, bool updateGeometry)
 {
   if (fiberPD == nullptr)
     this->m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
   else
     m_FiberPolyData->DeepCopy(fiberPD);
 
   m_NumFibers = m_FiberPolyData->GetNumberOfLines();
 
   if (updateGeometry)
     UpdateFiberGeometry();
   GenerateFiberIds();
   ColorFibersByOrientation();
 }
 
 /*
  * return vtkPolyData
  */
 vtkSmartPointer<vtkPolyData> mitk::FiberBundle::GetFiberPolyData() const
 {
   return m_FiberPolyData;
 }
 
 void mitk::FiberBundle::ColorFibersByOrientation()
 {
   //===== FOR WRITING A TEST ========================
   //  colorT size == tupelComponents * tupelElements
   //  compare color results
   //  to cover this code 100% also PolyData needed, where colorarray already exists
   //  + one fiber with exactly 1 point
   //  + one fiber with 0 points
   //=================================================
 
   vtkPoints* extrPoints = nullptr;
   extrPoints = m_FiberPolyData->GetPoints();
   int numOfPoints = 0;
   if (extrPoints!=nullptr)
     numOfPoints = extrPoints->GetNumberOfPoints();
 
   //colors and alpha value for each single point, RGBA = 4 components
   unsigned char rgba[4] = {0,0,0,0};
   int componentSize = 4;
   m_FiberColors = vtkSmartPointer<vtkUnsignedCharArray>::New();
   m_FiberColors->Allocate(numOfPoints * componentSize);
   m_FiberColors->SetNumberOfComponents(componentSize);
   m_FiberColors->SetName("FIBER_COLORS");
 
   int numOfFibers = m_FiberPolyData->GetNumberOfLines();
   if (numOfFibers < 1)
     return;
 
   /* extract single fibers of fiberBundle */
   vtkCellArray* fiberList = m_FiberPolyData->GetLines();
   fiberList->InitTraversal();
   for (int fi=0; fi<numOfFibers; ++fi) {
 
     vtkIdType* idList; // contains the point id's of the line
     vtkIdType pointsPerFiber; // number of points for current line
     fiberList->GetNextCell(pointsPerFiber, idList);
 
     /* single fiber checkpoints: is number of points valid */
     if (pointsPerFiber > 1)
     {
       /* operate on points of single fiber */
       for (int i=0; i <pointsPerFiber; ++i)
       {
         /* process all points elastV[0]ept starting and endpoint for calculating color value take current point, previous point and next point */
         if (i<pointsPerFiber-1 && i > 0)
         {
           /* The color value of the current point is influenced by the previous point and next point. */
           vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]);
           vnl_vector_fixed< double, 3 > nextPntvtk(extrPoints->GetPoint(idList[i+1])[0], extrPoints->GetPoint(idList[i+1])[1], extrPoints->GetPoint(idList[i+1])[2]);
           vnl_vector_fixed< double, 3 > prevPntvtk(extrPoints->GetPoint(idList[i-1])[0], extrPoints->GetPoint(idList[i-1])[1], extrPoints->GetPoint(idList[i-1])[2]);
 
           vnl_vector_fixed< double, 3 > diff1;
           diff1 = currentPntvtk - nextPntvtk;
 
           vnl_vector_fixed< double, 3 > diff2;
           diff2 = currentPntvtk - prevPntvtk;
 
           vnl_vector_fixed< double, 3 > diff;
           diff = (diff1 - diff2) / 2.0;
           diff.normalize();
 
           rgba[0] = (unsigned char) (255.0 * std::fabs(diff[0]));
           rgba[1] = (unsigned char) (255.0 * std::fabs(diff[1]));
           rgba[2] = (unsigned char) (255.0 * std::fabs(diff[2]));
           rgba[3] = (unsigned char) (255.0);
         }
         else if (i==0)
         {
           /* First point has no previous point, therefore only diff1 is taken */
 
           vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]);
           vnl_vector_fixed< double, 3 > nextPntvtk(extrPoints->GetPoint(idList[i+1])[0], extrPoints->GetPoint(idList[i+1])[1], extrPoints->GetPoint(idList[i+1])[2]);
 
           vnl_vector_fixed< double, 3 > diff1;
           diff1 = currentPntvtk - nextPntvtk;
           diff1.normalize();
 
           rgba[0] = (unsigned char) (255.0 * std::fabs(diff1[0]));
           rgba[1] = (unsigned char) (255.0 * std::fabs(diff1[1]));
           rgba[2] = (unsigned char) (255.0 * std::fabs(diff1[2]));
           rgba[3] = (unsigned char) (255.0);
         }
         else if (i==pointsPerFiber-1)
         {
           /* Last point has no next point, therefore only diff2 is taken */
           vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]);
           vnl_vector_fixed< double, 3 > prevPntvtk(extrPoints->GetPoint(idList[i-1])[0], extrPoints->GetPoint(idList[i-1])[1], extrPoints->GetPoint(idList[i-1])[2]);
 
           vnl_vector_fixed< double, 3 > diff2;
           diff2 = currentPntvtk - prevPntvtk;
           diff2.normalize();
 
           rgba[0] = (unsigned char) (255.0 * std::fabs(diff2[0]));
           rgba[1] = (unsigned char) (255.0 * std::fabs(diff2[1]));
           rgba[2] = (unsigned char) (255.0 * std::fabs(diff2[2]));
           rgba[3] = (unsigned char) (255.0);
         }
         m_FiberColors->InsertTypedTuple(idList[i], rgba);
       }
     }
     else if (pointsPerFiber == 1)
     {
       /* a single point does not define a fiber (use vertex mechanisms instead */
       continue;
     }
     else
     {
       MITK_DEBUG << "Fiber with 0 points detected... please check your tractography algorithm!" ;
       continue;
     }
   }
   m_UpdateTime3D.Modified();
   m_UpdateTime2D.Modified();
 }
 
 void mitk::FiberBundle::ColorFibersByCurvature(bool, bool normalize)
 {
   double window = 5;
 
   //colors and alpha value for each single point, RGBA = 4 components
   unsigned char rgba[4] = {0,0,0,0};
   int componentSize = 4;
   m_FiberColors = vtkSmartPointer<vtkUnsignedCharArray>::New();
   m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * componentSize);
   m_FiberColors->SetNumberOfComponents(componentSize);
   m_FiberColors->SetName("FIBER_COLORS");
 
   mitk::LookupTable::Pointer mitkLookup = mitk::LookupTable::New();
   vtkSmartPointer<vtkLookupTable> lookupTable = vtkSmartPointer<vtkLookupTable>::New();
   lookupTable->SetTableRange(0.0, 0.8);
   lookupTable->Build();
   mitkLookup->SetVtkLookupTable(lookupTable);
   mitkLookup->SetType(mitk::LookupTable::JET);
 
   vector< double > values;
   double min = 1;
   double max = 0;
   MITK_INFO << "Coloring fibers by curvature";
   boost::progress_display disp(m_FiberPolyData->GetNumberOfCells());
   for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
   {
     ++disp;
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     // calculate curvatures
     for (int j=0; j<numPoints; j++)
     {
       double dist = 0;
       int c = j;
       std::vector< vnl_vector_fixed< float, 3 > > vectors;
       vnl_vector_fixed< float, 3 > meanV; meanV.fill(0.0);
       while(dist<window/2 && c>1)
       {
         double p1[3];
         points->GetPoint(c-1, p1);
         double p2[3];
         points->GetPoint(c, p2);
 
         vnl_vector_fixed< float, 3 > v;
         v[0] = p2[0]-p1[0];
         v[1] = p2[1]-p1[1];
         v[2] = p2[2]-p1[2];
         dist += v.magnitude();
         v.normalize();
         vectors.push_back(v);
         if (c==j)
           meanV += v;
         c--;
       }
       c = j;
       dist = 0;
       while(dist<window/2 && c<numPoints-1)
       {
         double p1[3];
         points->GetPoint(c, p1);
         double p2[3];
         points->GetPoint(c+1, p2);
 
         vnl_vector_fixed< float, 3 > v;
         v[0] = p2[0]-p1[0];
         v[1] = p2[1]-p1[1];
         v[2] = p2[2]-p1[2];
         dist += v.magnitude();
         v.normalize();
         vectors.push_back(v);
         if (c==j)
           meanV += v;
         c++;
       }
       meanV.normalize();
 
       double dev = 0;
       for (unsigned int c=0; c<vectors.size(); c++)
       {
         double angle = dot_product(meanV, vectors.at(c));
         if (angle>1.0)
           angle = 1.0;
         if (angle<-1.0)
           angle = -1.0;
         dev += acos(angle)*180/M_PI;
       }
       if (vectors.size()>0)
         dev /= vectors.size();
 
       dev = 1.0-dev/180.0;
       values.push_back(dev);
       if (dev<min)
         min = dev;
       if (dev>max)
         max = dev;
     }
   }
   unsigned int count = 0;
   for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     for (int j=0; j<numPoints; j++)
     {
       double color[3];
       double dev = values.at(count);
       if (normalize)
         dev = (dev-min)/(max-min);
       else if (dev>1)
         dev = 1;
       lookupTable->GetColor(dev, color);
 
       rgba[0] = (unsigned char) (255.0 * color[0]);
       rgba[1] = (unsigned char) (255.0 * color[1]);
       rgba[2] = (unsigned char) (255.0 * color[2]);
       rgba[3] = (unsigned char) (255.0);
       m_FiberColors->InsertTypedTuple(cell->GetPointId(j), rgba);
       count++;
     }
   }
   m_UpdateTime3D.Modified();
   m_UpdateTime2D.Modified();
 }
 
 void mitk::FiberBundle::SetFiberOpacity(vtkDoubleArray* FAValArray)
 {
   for(long i=0; i<m_FiberColors->GetNumberOfTuples(); i++)
   {
     double faValue = FAValArray->GetValue(i);
     faValue = faValue * 255.0;
     m_FiberColors->SetComponent(i,3, (unsigned char) faValue );
   }
   m_UpdateTime3D.Modified();
   m_UpdateTime2D.Modified();
 }
 
 void mitk::FiberBundle::ResetFiberOpacity()
 {
   for(long i=0; i<m_FiberColors->GetNumberOfTuples(); i++)
     m_FiberColors->SetComponent(i,3, 255.0 );
   m_UpdateTime3D.Modified();
   m_UpdateTime2D.Modified();
 }
 
 void mitk::FiberBundle::ColorFibersByScalarMap(mitk::Image::Pointer FAimage, bool opacity, bool normalize)
 {
   mitkPixelTypeMultiplex3( ColorFibersByScalarMap, FAimage->GetPixelType(), FAimage, opacity, normalize );
   m_UpdateTime3D.Modified();
   m_UpdateTime2D.Modified();
 }
 
 template <typename TPixel>
 void mitk::FiberBundle::ColorFibersByScalarMap(const mitk::PixelType, mitk::Image::Pointer image, bool opacity, bool normalize)
 {
   m_FiberColors = vtkSmartPointer<vtkUnsignedCharArray>::New();
   m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * 4);
   m_FiberColors->SetNumberOfComponents(4);
   m_FiberColors->SetName("FIBER_COLORS");
 
   mitk::ImagePixelReadAccessor<TPixel,3> readimage(image, image->GetVolumeData(0));
 
   unsigned char rgba[4] = {0,0,0,0};
   vtkPoints* pointSet = m_FiberPolyData->GetPoints();
 
   mitk::LookupTable::Pointer mitkLookup = mitk::LookupTable::New();
   vtkSmartPointer<vtkLookupTable> lookupTable = vtkSmartPointer<vtkLookupTable>::New();
   lookupTable->SetTableRange(0.0, 0.8);
   lookupTable->Build();
   mitkLookup->SetVtkLookupTable(lookupTable);
   mitkLookup->SetType(mitk::LookupTable::JET);
 
   double min = 9999999;
   double max = -9999999;
   for(long i=0; i<m_FiberPolyData->GetNumberOfPoints(); ++i)
   {
     Point3D px;
     px[0] = pointSet->GetPoint(i)[0];
     px[1] = pointSet->GetPoint(i)[1];
     px[2] = pointSet->GetPoint(i)[2];
     double pixelValue = readimage.GetPixelByWorldCoordinates(px);
     if (pixelValue>max)
       max = pixelValue;
     if (pixelValue<min)
       min = pixelValue;
   }
 
   for(long i=0; i<m_FiberPolyData->GetNumberOfPoints(); ++i)
   {
     Point3D px;
     px[0] = pointSet->GetPoint(i)[0];
     px[1] = pointSet->GetPoint(i)[1];
     px[2] = pointSet->GetPoint(i)[2];
     double pixelValue = readimage.GetPixelByWorldCoordinates(px);
 
     if (normalize)
       pixelValue = (pixelValue-min)/(max-min);
     else if (pixelValue>1)
       pixelValue = 1;
 
     double color[3];
     lookupTable->GetColor(1-pixelValue, color);
 
     rgba[0] = (unsigned char) (255.0 * color[0]);
     rgba[1] = (unsigned char) (255.0 * color[1]);
     rgba[2] = (unsigned char) (255.0 * color[2]);
     if (opacity)
       rgba[3] = (unsigned char) (255.0 * pixelValue);
     else
       rgba[3] = (unsigned char) (255.0);
     m_FiberColors->InsertTypedTuple(i, rgba);
   }
   m_UpdateTime3D.Modified();
   m_UpdateTime2D.Modified();
 }
 
 
 void mitk::FiberBundle::ColorFibersByFiberWeights(bool opacity, bool normalize)
 {
   m_FiberColors = vtkSmartPointer<vtkUnsignedCharArray>::New();
   m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * 4);
   m_FiberColors->SetNumberOfComponents(4);
   m_FiberColors->SetName("FIBER_COLORS");
 
   mitk::LookupTable::Pointer mitkLookup = mitk::LookupTable::New();
   vtkSmartPointer<vtkLookupTable> lookupTable = vtkSmartPointer<vtkLookupTable>::New();
   lookupTable->SetTableRange(0.0, 0.8);
   lookupTable->Build();
   mitkLookup->SetVtkLookupTable(lookupTable);
   mitkLookup->SetType(mitk::LookupTable::JET);
 
   unsigned char rgba[4] = {0,0,0,0};
   unsigned int counter = 0;
 
   float max = -999999;
   float min = 999999;
   for (int i=0; i<m_NumFibers; i++)
   {
     double weight = this->GetFiberWeight(i);
     if (weight>max)
       max = weight;
     if (weight<min)
       min = weight;
   }
   if (fabs(max-min)<0.00001)
   {
     max = 1;
     min = 0;
   }
 
   for (int i=0; i<m_NumFibers; i++)
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     double weight = this->GetFiberWeight(i);
 
     for (int j=0; j<numPoints; j++)
     {
       float v = weight;
       if (normalize)
         v = (v-min)/(max-min);
       else if (v>1)
         v = 1;
       double color[3];
       lookupTable->GetColor(1-v, color);
 
       rgba[0] = (unsigned char) (255.0 * color[0]);
       rgba[1] = (unsigned char) (255.0 * color[1]);
       rgba[2] = (unsigned char) (255.0 * color[2]);
       if (opacity)
         rgba[3] = (unsigned char) (255.0 * v);
       else
         rgba[3] = (unsigned char) (255.0);
 
       m_FiberColors->InsertTypedTuple(counter, rgba);
       counter++;
     }
   }
 
   m_UpdateTime3D.Modified();
   m_UpdateTime2D.Modified();
 }
 
 void mitk::FiberBundle::SetFiberColors(float r, float g, float b, float alpha)
 {
   m_FiberColors = vtkSmartPointer<vtkUnsignedCharArray>::New();
   m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * 4);
   m_FiberColors->SetNumberOfComponents(4);
   m_FiberColors->SetName("FIBER_COLORS");
 
   unsigned char rgba[4] = {0,0,0,0};
   for(long i=0; i<m_FiberPolyData->GetNumberOfPoints(); ++i)
   {
     rgba[0] = (unsigned char) r;
     rgba[1] = (unsigned char) g;
     rgba[2] = (unsigned char) b;
     rgba[3] = (unsigned char) alpha;
     m_FiberColors->InsertTypedTuple(i, rgba);
   }
   m_UpdateTime3D.Modified();
   m_UpdateTime2D.Modified();
 }
 
 void mitk::FiberBundle::GenerateFiberIds()
 {
   if (m_FiberPolyData == nullptr)
     return;
 
   vtkSmartPointer<vtkIdFilter> idFiberFilter = vtkSmartPointer<vtkIdFilter>::New();
   idFiberFilter->SetInputData(m_FiberPolyData);
   idFiberFilter->CellIdsOn();
   //  idFiberFilter->PointIdsOn(); // point id's are not needed
   idFiberFilter->SetIdsArrayName(FIBER_ID_ARRAY);
   idFiberFilter->FieldDataOn();
   idFiberFilter->Update();
 
   m_FiberIdDataSet = idFiberFilter->GetOutput();
 
 }
 
 float mitk::FiberBundle::GetOverlap(ItkUcharImgType* mask, bool do_resampling)
 {
   vtkSmartPointer<vtkPolyData> PolyData = m_FiberPolyData;
   mitk::FiberBundle::Pointer fibCopy = this;
   if (do_resampling)
   {
     float minSpacing = 1;
     if(mask->GetSpacing()[0]<mask->GetSpacing()[1] && mask->GetSpacing()[0]<mask->GetSpacing()[2])
       minSpacing = mask->GetSpacing()[0];
     else if (mask->GetSpacing()[1] < mask->GetSpacing()[2])
       minSpacing = mask->GetSpacing()[1];
     else
       minSpacing = mask->GetSpacing()[2];
 
     fibCopy = this->GetDeepCopy();
     fibCopy->ResampleLinear(minSpacing/5);
     PolyData = fibCopy->GetFiberPolyData();
   }
 
   MITK_INFO << "Calculating overlap";
   int inside = 0;
   int outside = 0;
   boost::progress_display disp(m_NumFibers);
   for (int i=0; i<m_NumFibers; i++)
   {
     ++disp;
     vtkCell* cell = PolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     for (int j=0; j<numPoints; j++)
     {
       double* p = points->GetPoint(j);
       itk::Point<float, 3> itkP;
       itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2];
       itk::Index<3> idx;
       mask->TransformPhysicalPointToIndex(itkP, idx);
 
       if ( mask->GetLargestPossibleRegion().IsInside(idx) && mask->GetPixel(idx) != 0 )
         inside++;
       else
         outside++;
     }
   }
 
   if (inside+outside==0)
     outside = 1;
   return (float)inside/(inside+outside);
 }
 
 mitk::FiberBundle::Pointer mitk::FiberBundle::ExtractFiberSubset(ItkUcharImgType* mask, bool anyPoint, bool invert, bool bothEnds, float fraction, bool do_resampling)
 {
+  if (m_NumFibers==0 || mask==nullptr)
+    return mitk::FiberBundle::New(nullptr);
+
   vtkSmartPointer<vtkPolyData> PolyData = m_FiberPolyData;
   mitk::FiberBundle::Pointer fibCopy = this;
   if (anyPoint && do_resampling)
   {
     float minSpacing = 1;
     if(mask->GetSpacing()[0]<mask->GetSpacing()[1] && mask->GetSpacing()[0]<mask->GetSpacing()[2])
       minSpacing = mask->GetSpacing()[0];
     else if (mask->GetSpacing()[1] < mask->GetSpacing()[2])
       minSpacing = mask->GetSpacing()[1];
     else
       minSpacing = mask->GetSpacing()[2];
 
     fibCopy = this->GetDeepCopy();
     fibCopy->ResampleLinear(minSpacing/5);
     PolyData = fibCopy->GetFiberPolyData();
   }
   vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
   vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
 
   std::vector< float > new_weights;
 
   MITK_INFO << "Extracting fibers with mask image";
   boost::progress_display disp(m_NumFibers);
   for (int i=0; i<m_NumFibers; i++)
   {
     ++disp;
 
     vtkCell* cell = PolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     vtkCell* cellOriginal = m_FiberPolyData->GetCell(i);
     int numPointsOriginal = cellOriginal->GetNumberOfPoints();
     vtkPoints* pointsOriginal = cellOriginal->GetPoints();
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
 
     if (numPoints>1 && numPointsOriginal)
     {
       if (anyPoint)
       {
         int inside = 0;
         int outside = 0;
         if (!invert)
         {
           for (int j=0; j<numPoints; j++)
           {
             double* p = points->GetPoint(j);
 
             itk::Point<float, 3> itkP;
             itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2];
             itk::Index<3> idx;
             mask->TransformPhysicalPointToIndex(itkP, idx);
 
             if ( mask->GetLargestPossibleRegion().IsInside(idx) && mask->GetPixel(idx) != 0 )
             {
               inside++;
               if (fraction==0)
                 break;
             }
             else
               outside++;
           }
 
           float current_fraction = 0.0;
           if (inside+outside>0)
             current_fraction = (float)inside/(inside+outside);
 
           if (current_fraction>fraction)
           {
             for (int k=0; k<numPoints; k++)
             {
               double* p = points->GetPoint(k);
               vtkIdType id = vtkNewPoints->InsertNextPoint(p);
               container->GetPointIds()->InsertNextId(id);
             }
           }
         }
         else
         {
           bool includeFiber = true;
           for (int j=0; j<numPoints; j++)
           {
             double* p = points->GetPoint(j);
 
             itk::Point<float, 3> itkP;
             itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2];
             itk::Index<3> idx;
             mask->TransformPhysicalPointToIndex(itkP, idx);
 
             if ( mask->GetPixel(idx) != 0 && mask->GetLargestPossibleRegion().IsInside(idx) )
             {
               inside++;
               includeFiber = false;
               break;
             }
             else
               outside++;
           }
           if (includeFiber)
           {
 
             for (int k=0; k<numPoints; k++)
             {
               double* p = points->GetPoint(k);
               vtkIdType id = vtkNewPoints->InsertNextPoint(p);
               container->GetPointIds()->InsertNextId(id);
             }
           }
         }
       }
       else
       {
         double* start = pointsOriginal->GetPoint(0);
         itk::Point<float, 3> itkStart;
         itkStart[0] = start[0]; itkStart[1] = start[1]; itkStart[2] = start[2];
         itk::Index<3> idxStart;
         mask->TransformPhysicalPointToIndex(itkStart, idxStart);
 
         double* end = pointsOriginal->GetPoint(numPointsOriginal-1);
         itk::Point<float, 3> itkEnd;
         itkEnd[0] = end[0]; itkEnd[1] = end[1]; itkEnd[2] = end[2];
         itk::Index<3> idxEnd;
         mask->TransformPhysicalPointToIndex(itkEnd, idxEnd);
 
         if (invert)
         {
           if (bothEnds)
           {
             if ( mask->GetPixel(idxStart) == 0 && mask->GetPixel(idxEnd) == 0 )
             {
               for (int j=0; j<numPointsOriginal; j++)
               {
                 double* p = pointsOriginal->GetPoint(j);
                 vtkIdType id = vtkNewPoints->InsertNextPoint(p);
                 container->GetPointIds()->InsertNextId(id);
               }
             }
           }
           else if ( mask->GetPixel(idxStart) == 0 || mask->GetPixel(idxEnd) == 0 )
           {
             for (int j=0; j<numPointsOriginal; j++)
             {
               double* p = pointsOriginal->GetPoint(j);
               vtkIdType id = vtkNewPoints->InsertNextPoint(p);
               container->GetPointIds()->InsertNextId(id);
             }
           }
         }
         else
         {
           if (bothEnds)
           {
             if ( mask->GetPixel(idxStart) != 0 && mask->GetPixel(idxEnd) != 0 && mask->GetLargestPossibleRegion().IsInside(idxStart) && mask->GetLargestPossibleRegion().IsInside(idxEnd) )
             {
               for (int j=0; j<numPointsOriginal; j++)
               {
                 double* p = pointsOriginal->GetPoint(j);
                 vtkIdType id = vtkNewPoints->InsertNextPoint(p);
                 container->GetPointIds()->InsertNextId(id);
               }
             }
           }
           else if ( (mask->GetPixel(idxStart) != 0 && mask->GetLargestPossibleRegion().IsInside(idxStart)) || (mask->GetPixel(idxEnd) != 0 && mask->GetLargestPossibleRegion().IsInside(idxEnd)) )
           {
             for (int j=0; j<numPointsOriginal; j++)
             {
               double* p = pointsOriginal->GetPoint(j);
               vtkIdType id = vtkNewPoints->InsertNextPoint(p);
               container->GetPointIds()->InsertNextId(id);
             }
           }
         }
       }
     }
 
     if (container->GetNumberOfPoints()>0)
     {
       new_weights.push_back(fibCopy->GetFiberWeight(i));
       vtkNewCells->InsertNextCell(container);
     }
   }
 
   if (vtkNewCells->GetNumberOfCells()<=0)
-    return nullptr;
+    return mitk::FiberBundle::New(nullptr);
 
   vtkSmartPointer<vtkPolyData> newPolyData = vtkSmartPointer<vtkPolyData>::New();
   newPolyData->SetPoints(vtkNewPoints);
   newPolyData->SetLines(vtkNewCells);
+
   mitk::FiberBundle::Pointer newfib = mitk::FiberBundle::New(newPolyData);
   for (unsigned int i=0; i<new_weights.size(); i++)
     newfib->SetFiberWeight(i, new_weights.at(i));
   return newfib;
 }
 
 mitk::FiberBundle::Pointer mitk::FiberBundle::RemoveFibersOutside(ItkUcharImgType* mask, bool invert)
 {
   float minSpacing = 1;
   if(mask->GetSpacing()[0]<mask->GetSpacing()[1] && mask->GetSpacing()[0]<mask->GetSpacing()[2])
     minSpacing = mask->GetSpacing()[0];
   else if (mask->GetSpacing()[1] < mask->GetSpacing()[2])
     minSpacing = mask->GetSpacing()[1];
   else
     minSpacing = mask->GetSpacing()[2];
 
   mitk::FiberBundle::Pointer fibCopy = this->GetDeepCopy();
   fibCopy->ResampleLinear(minSpacing/10);
   vtkSmartPointer<vtkPolyData> PolyData =fibCopy->GetFiberPolyData();
 
   vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
   vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
 
   MITK_INFO << "Cutting fibers";
   boost::progress_display disp(m_NumFibers);
   for (int i=0; i<m_NumFibers; i++)
   {
     ++disp;
 
     vtkCell* cell = PolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     if (numPoints>1)
     {
       int newNumPoints = 0;
       for (int j=0; j<numPoints; j++)
       {
         double* p = points->GetPoint(j);
 
         itk::Point<float, 3> itkP;
         itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2];
         itk::Index<3> idx;
         mask->TransformPhysicalPointToIndex(itkP, idx);
 
         if ( mask->GetPixel(idx) != 0 && mask->GetLargestPossibleRegion().IsInside(idx) && !invert )
         {
           vtkIdType id = vtkNewPoints->InsertNextPoint(p);
           container->GetPointIds()->InsertNextId(id);
           newNumPoints++;
         }
         else if ( (mask->GetPixel(idx) == 0 || !mask->GetLargestPossibleRegion().IsInside(idx)) && invert )
         {
           vtkIdType id = vtkNewPoints->InsertNextPoint(p);
           container->GetPointIds()->InsertNextId(id);
           newNumPoints++;
         }
         else if (newNumPoints>0)
         {
           vtkNewCells->InsertNextCell(container);
 
           newNumPoints = 0;
           container = vtkSmartPointer<vtkPolyLine>::New();
         }
       }
 
       if (newNumPoints>0)
         vtkNewCells->InsertNextCell(container);
     }
 
   }
 
   if (vtkNewCells->GetNumberOfCells()<=0)
     return nullptr;
 
   vtkSmartPointer<vtkPolyData> newPolyData = vtkSmartPointer<vtkPolyData>::New();
   newPolyData->SetPoints(vtkNewPoints);
   newPolyData->SetLines(vtkNewCells);
   mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(newPolyData);
   newFib->Compress(0.1);
   return newFib;
 }
 
 mitk::FiberBundle::Pointer mitk::FiberBundle::ExtractFiberSubset(DataNode* roi, DataStorage* storage)
 {
   if (roi==nullptr || !(dynamic_cast<PlanarFigure*>(roi->GetData()) || dynamic_cast<PlanarFigureComposite*>(roi->GetData())) )
     return nullptr;
 
   std::vector<long> tmp = ExtractFiberIdSubset(roi, storage);
 
   if (tmp.size()<=0)
     return mitk::FiberBundle::New();
   vtkSmartPointer<vtkFloatArray> weights = vtkSmartPointer<vtkFloatArray>::New();
   vtkSmartPointer<vtkPolyData> pTmp = GeneratePolyDataByIds(tmp, weights);
   mitk::FiberBundle::Pointer fib = mitk::FiberBundle::New(pTmp);
   fib->SetFiberWeights(weights);
   return fib;
 }
 
 std::vector<long> mitk::FiberBundle::ExtractFiberIdSubset(DataNode *roi, DataStorage* storage)
 {
   std::vector<long> result;
   if (roi==nullptr || roi->GetData()==nullptr)
     return result;
 
   mitk::PlanarFigureComposite::Pointer pfc = dynamic_cast<mitk::PlanarFigureComposite*>(roi->GetData());
   if (!pfc.IsNull()) // handle composite
   {
     DataStorage::SetOfObjects::ConstPointer children = storage->GetDerivations(roi);
     if (children->size()==0)
       return result;
 
     switch (pfc->getOperationType())
     {
     case 0: // AND
     {
       MITK_INFO << "AND";
       result = this->ExtractFiberIdSubset(children->ElementAt(0), storage);
       std::vector<long>::iterator it;
       for (unsigned int i=1; i<children->Size(); ++i)
       {
         std::vector<long> inRoi = this->ExtractFiberIdSubset(children->ElementAt(i), storage);
 
         std::vector<long> rest(std::min(result.size(),inRoi.size()));
         it = std::set_intersection(result.begin(), result.end(), inRoi.begin(), inRoi.end(), rest.begin() );
         rest.resize( it - rest.begin() );
         result = rest;
       }
       break;
     }
     case 1: // OR
     {
       MITK_INFO << "OR";
       result = ExtractFiberIdSubset(children->ElementAt(0), storage);
       std::vector<long>::iterator it;
       for (unsigned int i=1; i<children->Size(); ++i)
       {
         it = result.end();
         std::vector<long> inRoi = ExtractFiberIdSubset(children->ElementAt(i), storage);
         result.insert(it, inRoi.begin(), inRoi.end());
       }
 
       // remove duplicates
       sort(result.begin(), result.end());
       it = unique(result.begin(), result.end());
       result.resize( it - result.begin() );
       break;
     }
     case 2: // NOT
     {
       MITK_INFO << "NOT";
       for(long i=0; i<this->GetNumFibers(); i++)
         result.push_back(i);
 
       std::vector<long>::iterator it;
       for (unsigned int i=0; i<children->Size(); ++i)
       {
         std::vector<long> inRoi = ExtractFiberIdSubset(children->ElementAt(i), storage);
 
         std::vector<long> rest(result.size()-inRoi.size());
         it = std::set_difference(result.begin(), result.end(), inRoi.begin(), inRoi.end(), rest.begin() );
         rest.resize( it - rest.begin() );
         result = rest;
       }
       break;
     }
     }
   }
   else if ( dynamic_cast<mitk::PlanarFigure*>(roi->GetData()) )  // actual extraction
   {
     if ( dynamic_cast<mitk::PlanarPolygon*>(roi->GetData()) )
     {
       mitk::PlanarFigure::Pointer planarPoly = dynamic_cast<mitk::PlanarFigure*>(roi->GetData());
 
       //create vtkPolygon using controlpoints from planarFigure polygon
       vtkSmartPointer<vtkPolygon> polygonVtk = vtkSmartPointer<vtkPolygon>::New();
       for (unsigned int i=0; i<planarPoly->GetNumberOfControlPoints(); ++i)
       {
         itk::Point<double,3> p = planarPoly->GetWorldControlPoint(i);
         vtkIdType id = polygonVtk->GetPoints()->InsertNextPoint(p[0], p[1], p[2] );
         polygonVtk->GetPointIds()->InsertNextId(id);
       }
 
       MITK_INFO << "Extracting with polygon";
       boost::progress_display disp(m_NumFibers);
       for (int i=0; i<m_NumFibers; i++)
       {
         ++disp ;
         vtkCell* cell = m_FiberPolyData->GetCell(i);
         int numPoints = cell->GetNumberOfPoints();
         vtkPoints* points = cell->GetPoints();
 
         for (int j=0; j<numPoints-1; j++)
         {
           // Inputs
           double p1[3] = {0,0,0};
           points->GetPoint(j, p1);
           double p2[3] = {0,0,0};
           points->GetPoint(j+1, p2);
           double tolerance = 0.001;
 
           // Outputs
           double t = 0; // Parametric coordinate of intersection (0 (corresponding to p1) to 1 (corresponding to p2))
           double x[3] = {0,0,0}; // The coordinate of the intersection
           double pcoords[3] = {0,0,0};
           int subId = 0;
 
           int iD = polygonVtk->IntersectWithLine(p1, p2, tolerance, t, x, pcoords, subId);
           if (iD!=0)
           {
             result.push_back(i);
             break;
           }
         }
       }
     }
     else if ( dynamic_cast<mitk::PlanarCircle*>(roi->GetData()) )
     {
       mitk::PlanarFigure::Pointer planarFigure = dynamic_cast<mitk::PlanarFigure*>(roi->GetData());
       Vector3D planeNormal = planarFigure->GetPlaneGeometry()->GetNormal();
       planeNormal.Normalize();
 
       //calculate circle radius
       mitk::Point3D V1w = planarFigure->GetWorldControlPoint(0); //centerPoint
       mitk::Point3D V2w  = planarFigure->GetWorldControlPoint(1); //radiusPoint
 
       double radius = V1w.EuclideanDistanceTo(V2w);
       radius *= radius;
 
       MITK_INFO << "Extracting with circle";
       boost::progress_display disp(m_NumFibers);
       for (int i=0; i<m_NumFibers; i++)
       {
         ++disp ;
         vtkCell* cell = m_FiberPolyData->GetCell(i);
         int numPoints = cell->GetNumberOfPoints();
         vtkPoints* points = cell->GetPoints();
 
         for (int j=0; j<numPoints-1; j++)
         {
           // Inputs
           double p1[3] = {0,0,0};
           points->GetPoint(j, p1);
           double p2[3] = {0,0,0};
           points->GetPoint(j+1, p2);
 
           // Outputs
           double t = 0; // Parametric coordinate of intersection (0 (corresponding to p1) to 1 (corresponding to p2))
           double x[3] = {0,0,0}; // The coordinate of the intersection
 
           int iD = vtkPlane::IntersectWithLine(p1,p2,planeNormal.GetDataPointer(),V1w.GetDataPointer(),t,x);
 
           if (iD!=0)
           {
             double dist = (x[0]-V1w[0])*(x[0]-V1w[0])+(x[1]-V1w[1])*(x[1]-V1w[1])+(x[2]-V1w[2])*(x[2]-V1w[2]);
             if( dist <= radius)
             {
               result.push_back(i);
               break;
             }
           }
         }
       }
     }
     return result;
   }
 
   return result;
 }
 
 void mitk::FiberBundle::UpdateFiberGeometry()
 {
   vtkSmartPointer<vtkCleanPolyData> cleaner = vtkSmartPointer<vtkCleanPolyData>::New();
   cleaner->SetInputData(m_FiberPolyData);
   cleaner->PointMergingOff();
   cleaner->Update();
   m_FiberPolyData = cleaner->GetOutput();
 
   m_FiberLengths.clear();
   m_MeanFiberLength = 0;
   m_MedianFiberLength = 0;
   m_LengthStDev = 0;
   m_NumFibers = m_FiberPolyData->GetNumberOfCells();
 
   if (m_FiberColors==nullptr || m_FiberColors->GetNumberOfTuples()!=m_FiberPolyData->GetNumberOfPoints())
     this->ColorFibersByOrientation();
 
   if (m_FiberWeights->GetNumberOfValues()!=m_NumFibers)
   {
     m_FiberWeights = vtkSmartPointer<vtkFloatArray>::New();
     m_FiberWeights->SetName("FIBER_WEIGHTS");
     m_FiberWeights->SetNumberOfValues(m_NumFibers);
     this->SetFiberWeights(1);
   }
 
   if (m_NumFibers<=0) // no fibers present; apply default geometry
   {
     m_MinFiberLength = 0;
     m_MaxFiberLength = 0;
     mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New();
     geometry->SetImageGeometry(false);
     float b[] = {0, 1, 0, 1, 0, 1};
     geometry->SetFloatBounds(b);
     SetGeometry(geometry);
     return;
   }
   double b[6];
   m_FiberPolyData->GetBounds(b);
 
   // calculate statistics
   for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int p = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
     float length = 0;
     for (int j=0; j<p-1; j++)
     {
       double p1[3];
       points->GetPoint(j, p1);
       double p2[3];
       points->GetPoint(j+1, p2);
 
       float dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2]));
       length += dist;
     }
     m_FiberLengths.push_back(length);
     m_MeanFiberLength += length;
     if (i==0)
     {
       m_MinFiberLength = length;
       m_MaxFiberLength = length;
     }
     else
     {
       if (length<m_MinFiberLength)
         m_MinFiberLength = length;
       if (length>m_MaxFiberLength)
         m_MaxFiberLength = length;
     }
   }
   m_MeanFiberLength /= m_NumFibers;
 
   std::vector< float > sortedLengths = m_FiberLengths;
   std::sort(sortedLengths.begin(), sortedLengths.end());
   for (int i=0; i<m_NumFibers; i++)
     m_LengthStDev += (m_MeanFiberLength-sortedLengths.at(i))*(m_MeanFiberLength-sortedLengths.at(i));
   if (m_NumFibers>1)
     m_LengthStDev /= (m_NumFibers-1);
   else
     m_LengthStDev = 0;
   m_LengthStDev = std::sqrt(m_LengthStDev);
   m_MedianFiberLength = sortedLengths.at(m_NumFibers/2);
 
   mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New();
   geometry->SetFloatBounds(b);
   this->SetGeometry(geometry);
 
   m_UpdateTime3D.Modified();
   m_UpdateTime2D.Modified();
 }
 
 float mitk::FiberBundle::GetFiberWeight(unsigned int fiber) const
 {
   return m_FiberWeights->GetValue(fiber);
 }
 
 void mitk::FiberBundle::SetFiberWeights(float newWeight)
 {
   for (int i=0; i<m_FiberWeights->GetNumberOfValues(); i++)
     m_FiberWeights->SetValue(i, newWeight);
 }
 
 void mitk::FiberBundle::SetFiberWeights(vtkSmartPointer<vtkFloatArray> weights)
 {
   if (m_NumFibers!=weights->GetNumberOfValues())
   {
     MITK_INFO << "Weights array not equal to number of fibers! " << weights->GetNumberOfValues() << " vs " << m_NumFibers;
     return;
   }
 
   for (int i=0; i<weights->GetNumberOfValues(); i++)
     m_FiberWeights->SetValue(i, weights->GetValue(i));
 
   m_FiberWeights->SetName("FIBER_WEIGHTS");
 }
 
 void mitk::FiberBundle::SetFiberWeight(unsigned int fiber, float weight)
 {
   m_FiberWeights->SetValue(fiber, weight);
 }
 
 void mitk::FiberBundle::SetFiberColors(vtkSmartPointer<vtkUnsignedCharArray> fiberColors)
 {
   for(long i=0; i<m_FiberPolyData->GetNumberOfPoints(); ++i)
   {
     unsigned char source[4] = {0,0,0,0};
     fiberColors->GetTypedTuple(i, source);
 
     unsigned char target[4] = {0,0,0,0};
     target[0] = source[0];
     target[1] = source[1];
     target[2] = source[2];
     target[3] = source[3];
     m_FiberColors->InsertTypedTuple(i, target);
   }
   m_UpdateTime3D.Modified();
   m_UpdateTime2D.Modified();
 }
 
 itk::Matrix< double, 3, 3 > mitk::FiberBundle::TransformMatrix(itk::Matrix< double, 3, 3 > m, double rx, double ry, double rz)
 {
   rx = rx*M_PI/180;
   ry = ry*M_PI/180;
   rz = rz*M_PI/180;
 
   itk::Matrix< double, 3, 3 > rotX; rotX.SetIdentity();
   rotX[1][1] = cos(rx);
   rotX[2][2] = rotX[1][1];
   rotX[1][2] = -sin(rx);
   rotX[2][1] = -rotX[1][2];
 
   itk::Matrix< double, 3, 3 > rotY; rotY.SetIdentity();
   rotY[0][0] = cos(ry);
   rotY[2][2] = rotY[0][0];
   rotY[0][2] = sin(ry);
   rotY[2][0] = -rotY[0][2];
 
   itk::Matrix< double, 3, 3 > rotZ; rotZ.SetIdentity();
   rotZ[0][0] = cos(rz);
   rotZ[1][1] = rotZ[0][0];
   rotZ[0][1] = -sin(rz);
   rotZ[1][0] = -rotZ[0][1];
 
   itk::Matrix< double, 3, 3 > rot = rotZ*rotY*rotX;
 
   m = rot*m;
 
   return m;
 }
 
 itk::Point<float, 3> mitk::FiberBundle::TransformPoint(vnl_vector_fixed< double, 3 > point, double rx, double ry, double rz, double tx, double ty, double tz)
 {
   rx = rx*M_PI/180;
   ry = ry*M_PI/180;
   rz = rz*M_PI/180;
 
   vnl_matrix_fixed< double, 3, 3 > rotX; rotX.set_identity();
   rotX[1][1] = cos(rx);
   rotX[2][2] = rotX[1][1];
   rotX[1][2] = -sin(rx);
   rotX[2][1] = -rotX[1][2];
 
   vnl_matrix_fixed< double, 3, 3 > rotY; rotY.set_identity();
   rotY[0][0] = cos(ry);
   rotY[2][2] = rotY[0][0];
   rotY[0][2] = sin(ry);
   rotY[2][0] = -rotY[0][2];
 
   vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity();
   rotZ[0][0] = cos(rz);
   rotZ[1][1] = rotZ[0][0];
   rotZ[0][1] = -sin(rz);
   rotZ[1][0] = -rotZ[0][1];
 
   vnl_matrix_fixed< double, 3, 3 > rot = rotZ*rotY*rotX;
 
   mitk::BaseGeometry::Pointer geom = this->GetGeometry();
   mitk::Point3D center = geom->GetCenter();
 
   point[0] -= center[0];
   point[1] -= center[1];
   point[2] -= center[2];
   point = rot*point;
   point[0] += center[0]+tx;
   point[1] += center[1]+ty;
   point[2] += center[2]+tz;
   itk::Point<float, 3> out; out[0] = point[0]; out[1] = point[1]; out[2] = point[2];
   return out;
 }
 
 void mitk::FiberBundle::TransformFibers(double rx, double ry, double rz, double tx, double ty, double tz)
 {
   rx = rx*M_PI/180;
   ry = ry*M_PI/180;
   rz = rz*M_PI/180;
 
   vnl_matrix_fixed< double, 3, 3 > rotX; rotX.set_identity();
   rotX[1][1] = cos(rx);
   rotX[2][2] = rotX[1][1];
   rotX[1][2] = -sin(rx);
   rotX[2][1] = -rotX[1][2];
 
   vnl_matrix_fixed< double, 3, 3 > rotY; rotY.set_identity();
   rotY[0][0] = cos(ry);
   rotY[2][2] = rotY[0][0];
   rotY[0][2] = sin(ry);
   rotY[2][0] = -rotY[0][2];
 
   vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity();
   rotZ[0][0] = cos(rz);
   rotZ[1][1] = rotZ[0][0];
   rotZ[0][1] = -sin(rz);
   rotZ[1][0] = -rotZ[0][1];
 
   vnl_matrix_fixed< double, 3, 3 > rot = rotZ*rotY*rotX;
 
   mitk::BaseGeometry::Pointer geom = this->GetGeometry();
   mitk::Point3D center = geom->GetCenter();
 
   vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
   vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
 
   for (int i=0; i<m_NumFibers; i++)
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     for (int j=0; j<numPoints; j++)
     {
       double* p = points->GetPoint(j);
       vnl_vector_fixed< double, 3 > dir;
       dir[0] = p[0]-center[0];
       dir[1] = p[1]-center[1];
       dir[2] = p[2]-center[2];
       dir = rot*dir;
       dir[0] += center[0]+tx;
       dir[1] += center[1]+ty;
       dir[2] += center[2]+tz;
       vtkIdType id = vtkNewPoints->InsertNextPoint(dir.data_block());
       container->GetPointIds()->InsertNextId(id);
     }
     vtkNewCells->InsertNextCell(container);
   }
 
   m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
   m_FiberPolyData->SetPoints(vtkNewPoints);
   m_FiberPolyData->SetLines(vtkNewCells);
   this->SetFiberPolyData(m_FiberPolyData, true);
 }
 
 void mitk::FiberBundle::RotateAroundAxis(double x, double y, double z)
 {
   x = x*M_PI/180;
   y = y*M_PI/180;
   z = z*M_PI/180;
 
   vnl_matrix_fixed< double, 3, 3 > rotX; rotX.set_identity();
   rotX[1][1] = cos(x);
   rotX[2][2] = rotX[1][1];
   rotX[1][2] = -sin(x);
   rotX[2][1] = -rotX[1][2];
 
   vnl_matrix_fixed< double, 3, 3 > rotY; rotY.set_identity();
   rotY[0][0] = cos(y);
   rotY[2][2] = rotY[0][0];
   rotY[0][2] = sin(y);
   rotY[2][0] = -rotY[0][2];
 
   vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity();
   rotZ[0][0] = cos(z);
   rotZ[1][1] = rotZ[0][0];
   rotZ[0][1] = -sin(z);
   rotZ[1][0] = -rotZ[0][1];
 
   mitk::BaseGeometry::Pointer geom = this->GetGeometry();
   mitk::Point3D center = geom->GetCenter();
 
   vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
   vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
 
   for (int i=0; i<m_NumFibers; i++)
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     for (int j=0; j<numPoints; j++)
     {
       double* p = points->GetPoint(j);
       vnl_vector_fixed< double, 3 > dir;
       dir[0] = p[0]-center[0];
       dir[1] = p[1]-center[1];
       dir[2] = p[2]-center[2];
       dir = rotZ*rotY*rotX*dir;
       dir[0] += center[0];
       dir[1] += center[1];
       dir[2] += center[2];
       vtkIdType id = vtkNewPoints->InsertNextPoint(dir.data_block());
       container->GetPointIds()->InsertNextId(id);
     }
     vtkNewCells->InsertNextCell(container);
   }
 
   m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
   m_FiberPolyData->SetPoints(vtkNewPoints);
   m_FiberPolyData->SetLines(vtkNewCells);
   this->SetFiberPolyData(m_FiberPolyData, true);
 }
 
 void mitk::FiberBundle::ScaleFibers(double x, double y, double z, bool subtractCenter)
 {
   MITK_INFO << "Scaling fibers";
   boost::progress_display disp(m_NumFibers);
 
   mitk::BaseGeometry* geom = this->GetGeometry();
   mitk::Point3D c = geom->GetCenter();
 
   vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
   vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
 
   for (int i=0; i<m_NumFibers; i++)
   {
     ++disp ;
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     for (int j=0; j<numPoints; j++)
     {
       double* p = points->GetPoint(j);
       if (subtractCenter)
       {
         p[0] -= c[0]; p[1] -= c[1]; p[2] -= c[2];
       }
       p[0] *= x;
       p[1] *= y;
       p[2] *= z;
       if (subtractCenter)
       {
         p[0] += c[0]; p[1] += c[1]; p[2] += c[2];
       }
       vtkIdType id = vtkNewPoints->InsertNextPoint(p);
       container->GetPointIds()->InsertNextId(id);
     }
     vtkNewCells->InsertNextCell(container);
   }
 
   m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
   m_FiberPolyData->SetPoints(vtkNewPoints);
   m_FiberPolyData->SetLines(vtkNewCells);
   this->SetFiberPolyData(m_FiberPolyData, true);
 }
 
 void mitk::FiberBundle::TranslateFibers(double x, double y, double z)
 {
   vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
   vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
 
   for (int i=0; i<m_NumFibers; i++)
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     for (int j=0; j<numPoints; j++)
     {
       double* p = points->GetPoint(j);
       p[0] += x;
       p[1] += y;
       p[2] += z;
       vtkIdType id = vtkNewPoints->InsertNextPoint(p);
       container->GetPointIds()->InsertNextId(id);
     }
     vtkNewCells->InsertNextCell(container);
   }
 
   m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
   m_FiberPolyData->SetPoints(vtkNewPoints);
   m_FiberPolyData->SetLines(vtkNewCells);
   this->SetFiberPolyData(m_FiberPolyData, true);
 }
 
 void mitk::FiberBundle::MirrorFibers(unsigned int axis)
 {
   if (axis>2)
     return;
 
   MITK_INFO << "Mirroring fibers";
   boost::progress_display disp(m_NumFibers);
 
   vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
   vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
 
   for (int i=0; i<m_NumFibers; i++)
   {
     ++disp;
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     for (int j=0; j<numPoints; j++)
     {
       double* p = points->GetPoint(j);
       p[axis] = -p[axis];
       vtkIdType id = vtkNewPoints->InsertNextPoint(p);
       container->GetPointIds()->InsertNextId(id);
     }
     vtkNewCells->InsertNextCell(container);
   }
 
   m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
   m_FiberPolyData->SetPoints(vtkNewPoints);
   m_FiberPolyData->SetLines(vtkNewCells);
   this->SetFiberPolyData(m_FiberPolyData, true);
 }
 
 void mitk::FiberBundle::RemoveDir(vnl_vector_fixed<double,3> dir, double threshold)
 {
   dir.normalize();
   vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
   vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
 
   boost::progress_display disp(m_FiberPolyData->GetNumberOfCells());
   for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
   {
     ++disp ;
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     // calculate curvatures
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     bool discard = false;
     for (int j=0; j<numPoints-1; j++)
     {
       double p1[3];
       points->GetPoint(j, p1);
       double p2[3];
       points->GetPoint(j+1, p2);
 
       vnl_vector_fixed< double, 3 > v1;
       v1[0] = p2[0]-p1[0];
       v1[1] = p2[1]-p1[1];
       v1[2] = p2[2]-p1[2];
       if (v1.magnitude()>0.001)
       {
         v1.normalize();
 
         if (fabs(dot_product(v1,dir))>threshold)
         {
           discard = true;
           break;
         }
       }
     }
     if (!discard)
     {
       for (int j=0; j<numPoints; j++)
       {
         double p1[3];
         points->GetPoint(j, p1);
 
         vtkIdType id = vtkNewPoints->InsertNextPoint(p1);
         container->GetPointIds()->InsertNextId(id);
       }
       vtkNewCells->InsertNextCell(container);
     }
   }
 
   m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
   m_FiberPolyData->SetPoints(vtkNewPoints);
   m_FiberPolyData->SetLines(vtkNewCells);
 
   this->SetFiberPolyData(m_FiberPolyData, true);
 
   //    UpdateColorCoding();
   //    UpdateFiberGeometry();
 }
 
 bool mitk::FiberBundle::ApplyCurvatureThreshold(float minRadius, bool deleteFibers)
 {
   if (minRadius<0)
     return true;
 
   vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
   vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
 
   MITK_INFO << "Applying curvature threshold";
   boost::progress_display disp(m_FiberPolyData->GetNumberOfCells());
   for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
   {
     ++disp ;
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     // calculate curvatures
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     for (int j=0; j<numPoints-2; j++)
     {
       double p1[3];
       points->GetPoint(j, p1);
       double p2[3];
       points->GetPoint(j+1, p2);
       double p3[3];
       points->GetPoint(j+2, p3);
 
       vnl_vector_fixed< float, 3 > v1, v2, v3;
 
       v1[0] = p2[0]-p1[0];
       v1[1] = p2[1]-p1[1];
       v1[2] = p2[2]-p1[2];
 
       v2[0] = p3[0]-p2[0];
       v2[1] = p3[1]-p2[1];
       v2[2] = p3[2]-p2[2];
 
       v3[0] = p1[0]-p3[0];
       v3[1] = p1[1]-p3[1];
       v3[2] = p1[2]-p3[2];
 
       float a = v1.magnitude();
       float b = v2.magnitude();
       float c = v3.magnitude();
       float r = a*b*c/std::sqrt((a+b+c)*(a+b-c)*(b+c-a)*(a-b+c)); // radius of triangle via Heron's formula (area of triangle)
 
       vtkIdType id = vtkNewPoints->InsertNextPoint(p1);
       container->GetPointIds()->InsertNextId(id);
 
       if (deleteFibers && r<minRadius)
         break;
 
       if (r<minRadius)
       {
         j += 2;
         vtkNewCells->InsertNextCell(container);
         container = vtkSmartPointer<vtkPolyLine>::New();
       }
       else if (j==numPoints-3)
       {
         id = vtkNewPoints->InsertNextPoint(p2);
         container->GetPointIds()->InsertNextId(id);
         id = vtkNewPoints->InsertNextPoint(p3);
         container->GetPointIds()->InsertNextId(id);
         vtkNewCells->InsertNextCell(container);
       }
     }
   }
 
   if (vtkNewCells->GetNumberOfCells()<=0)
     return false;
 
   m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
   m_FiberPolyData->SetPoints(vtkNewPoints);
   m_FiberPolyData->SetLines(vtkNewCells);
   this->SetFiberPolyData(m_FiberPolyData, true);
   return true;
 }
 
 bool mitk::FiberBundle::RemoveShortFibers(float lengthInMM)
 {
   MITK_INFO << "Removing short fibers";
   if (lengthInMM<=0 || lengthInMM<m_MinFiberLength)
   {
     MITK_INFO << "No fibers shorter than " << lengthInMM << " mm found!";
     return true;
   }
 
   if (lengthInMM>m_MaxFiberLength)    // can't remove all fibers
   {
     MITK_WARN << "Process aborted. No fibers would be left!";
     return false;
   }
 
   vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
   vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
   float min = m_MaxFiberLength;
 
   boost::progress_display disp(m_NumFibers);
   for (int i=0; i<m_NumFibers; i++)
   {
     ++disp;
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     if (m_FiberLengths.at(i)>=lengthInMM)
     {
       vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
       for (int j=0; j<numPoints; j++)
       {
         double* p = points->GetPoint(j);
         vtkIdType id = vtkNewPoints->InsertNextPoint(p);
         container->GetPointIds()->InsertNextId(id);
       }
       vtkNewCells->InsertNextCell(container);
       if (m_FiberLengths.at(i)<min)
         min = m_FiberLengths.at(i);
     }
   }
 
   if (vtkNewCells->GetNumberOfCells()<=0)
     return false;
 
   m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
   m_FiberPolyData->SetPoints(vtkNewPoints);
   m_FiberPolyData->SetLines(vtkNewCells);
   this->SetFiberPolyData(m_FiberPolyData, true);
   return true;
 }
 
 bool mitk::FiberBundle::RemoveLongFibers(float lengthInMM)
 {
   if (lengthInMM<=0 || lengthInMM>m_MaxFiberLength)
     return true;
 
   if (lengthInMM<m_MinFiberLength)    // can't remove all fibers
     return false;
 
   vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
   vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
 
   MITK_INFO << "Removing long fibers";
   boost::progress_display disp(m_NumFibers);
   for (int i=0; i<m_NumFibers; i++)
   {
     ++disp;
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     if (m_FiberLengths.at(i)<=lengthInMM)
     {
       vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
       for (int j=0; j<numPoints; j++)
       {
         double* p = points->GetPoint(j);
         vtkIdType id = vtkNewPoints->InsertNextPoint(p);
         container->GetPointIds()->InsertNextId(id);
       }
       vtkNewCells->InsertNextCell(container);
     }
   }
 
   if (vtkNewCells->GetNumberOfCells()<=0)
     return false;
 
   m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
   m_FiberPolyData->SetPoints(vtkNewPoints);
   m_FiberPolyData->SetLines(vtkNewCells);
   this->SetFiberPolyData(m_FiberPolyData, true);
   return true;
 }
 
 void mitk::FiberBundle::ResampleSpline(float pointDistance, double tension, double continuity, double bias )
 {
   if (pointDistance<=0)
     return;
 
   vtkSmartPointer<vtkPoints> vtkSmoothPoints = vtkSmartPointer<vtkPoints>::New(); //in smoothpoints the interpolated points representing a fiber are stored.
 
   //in vtkcells all polylines are stored, actually all id's of them are stored
   vtkSmartPointer<vtkCellArray> vtkSmoothCells = vtkSmartPointer<vtkCellArray>::New(); //cellcontainer for smoothed lines
   vtkIdType pointHelperCnt = 0;
 
   MITK_INFO << "Smoothing fibers";
   vtkSmartPointer<vtkFloatArray> newFiberWeights = vtkSmartPointer<vtkFloatArray>::New();
   newFiberWeights->SetName("FIBER_WEIGHTS");
   newFiberWeights->SetNumberOfValues(m_NumFibers);
 
   boost::progress_display disp(m_NumFibers);
 #pragma omp parallel for
   for (int i=0; i<m_NumFibers; i++)
   {
     vtkSmartPointer<vtkPoints> newPoints = vtkSmartPointer<vtkPoints>::New();
     float length = 0;
     float weight = 1;
 #pragma omp critical
     {
       length = m_FiberLengths.at(i);
       weight = m_FiberWeights->GetValue(i);
       ++disp;
       vtkCell* cell = m_FiberPolyData->GetCell(i);
       int numPoints = cell->GetNumberOfPoints();
       vtkPoints* points = cell->GetPoints();
       for (int j=0; j<numPoints; j++)
         newPoints->InsertNextPoint(points->GetPoint(j));
     }
 
     int sampling = std::ceil(length/pointDistance);
 
     vtkSmartPointer<vtkKochanekSpline> xSpline = vtkSmartPointer<vtkKochanekSpline>::New();
     vtkSmartPointer<vtkKochanekSpline> ySpline = vtkSmartPointer<vtkKochanekSpline>::New();
     vtkSmartPointer<vtkKochanekSpline> zSpline = vtkSmartPointer<vtkKochanekSpline>::New();
     xSpline->SetDefaultBias(bias); xSpline->SetDefaultTension(tension); xSpline->SetDefaultContinuity(continuity);
     ySpline->SetDefaultBias(bias); ySpline->SetDefaultTension(tension); ySpline->SetDefaultContinuity(continuity);
     zSpline->SetDefaultBias(bias); zSpline->SetDefaultTension(tension); zSpline->SetDefaultContinuity(continuity);
 
     vtkSmartPointer<vtkParametricSpline> spline = vtkSmartPointer<vtkParametricSpline>::New();
     spline->SetXSpline(xSpline);
     spline->SetYSpline(ySpline);
     spline->SetZSpline(zSpline);
     spline->SetPoints(newPoints);
 
     vtkSmartPointer<vtkParametricFunctionSource> functionSource = vtkSmartPointer<vtkParametricFunctionSource>::New();
     functionSource->SetParametricFunction(spline);
     functionSource->SetUResolution(sampling);
     functionSource->SetVResolution(sampling);
     functionSource->SetWResolution(sampling);
     functionSource->Update();
 
     vtkPolyData* outputFunction = functionSource->GetOutput();
     vtkPoints* tmpSmoothPnts = outputFunction->GetPoints(); //smoothPoints of current fiber
 
     vtkSmartPointer<vtkPolyLine> smoothLine = vtkSmartPointer<vtkPolyLine>::New();
     smoothLine->GetPointIds()->SetNumberOfIds(tmpSmoothPnts->GetNumberOfPoints());
 
 #pragma omp critical
     {
       for (int j=0; j<smoothLine->GetNumberOfPoints(); j++)
       {
         smoothLine->GetPointIds()->SetId(j, j+pointHelperCnt);
         vtkSmoothPoints->InsertNextPoint(tmpSmoothPnts->GetPoint(j));
       }
 
       newFiberWeights->SetValue(vtkSmoothCells->GetNumberOfCells(), weight);
       vtkSmoothCells->InsertNextCell(smoothLine);
       pointHelperCnt += tmpSmoothPnts->GetNumberOfPoints();
     }
   }
 
   SetFiberWeights(newFiberWeights);
   m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
   m_FiberPolyData->SetPoints(vtkSmoothPoints);
   m_FiberPolyData->SetLines(vtkSmoothCells);
   this->SetFiberPolyData(m_FiberPolyData, true);
 }
 
 void mitk::FiberBundle::ResampleSpline(float pointDistance)
 {
   ResampleSpline(pointDistance, 0, 0, 0 );
 }
 
 unsigned long mitk::FiberBundle::GetNumberOfPoints() const
 {
   unsigned long points = 0;
   for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     points += cell->GetNumberOfPoints();
   }
   return points;
 }
 
 void mitk::FiberBundle::Compress(float error)
 {
   vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
   vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
 
   MITK_INFO << "Compressing fibers";
   unsigned long numRemovedPoints = 0;
   boost::progress_display disp(m_FiberPolyData->GetNumberOfCells());
   vtkSmartPointer<vtkFloatArray> newFiberWeights = vtkSmartPointer<vtkFloatArray>::New();
   newFiberWeights->SetName("FIBER_WEIGHTS");
   newFiberWeights->SetNumberOfValues(m_NumFibers);
 
 #pragma omp parallel for
   for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
   {
 
     std::vector< vnl_vector_fixed< double, 3 > > vertices;
     float weight = 1;
 
 #pragma omp critical
     {
       ++disp;
       weight = m_FiberWeights->GetValue(i);
       vtkCell* cell = m_FiberPolyData->GetCell(i);
       int numPoints = cell->GetNumberOfPoints();
       vtkPoints* points = cell->GetPoints();
 
       for (int j=0; j<numPoints; j++)
       {
         double cand[3];
         points->GetPoint(j, cand);
         vnl_vector_fixed< double, 3 > candV;
         candV[0]=cand[0]; candV[1]=cand[1]; candV[2]=cand[2];
         vertices.push_back(candV);
       }
     }
 
     // calculate curvatures
     int numPoints = vertices.size();
     std::vector< int > removedPoints; removedPoints.resize(numPoints, 0);
     removedPoints[0]=-1; removedPoints[numPoints-1]=-1;
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     int remCounter = 0;
 
     bool pointFound = true;
     while (pointFound)
     {
       pointFound = false;
       double minError = error;
       int removeIndex = -1;
 
       for (unsigned int j=0; j<vertices.size(); j++)
       {
         if (removedPoints[j]==0)
         {
           vnl_vector_fixed< double, 3 > candV = vertices.at(j);
 
           int validP = -1;
           vnl_vector_fixed< double, 3 > pred;
           for (int k=j-1; k>=0; k--)
             if (removedPoints[k]<=0)
             {
               pred = vertices.at(k);
               validP = k;
               break;
             }
           int validS = -1;
           vnl_vector_fixed< double, 3 > succ;
           for (int k=j+1; k<numPoints; k++)
             if (removedPoints[k]<=0)
             {
               succ = vertices.at(k);
               validS = k;
               break;
             }
 
           if (validP>=0 && validS>=0)
           {
             double a = (candV-pred).magnitude();
             double b = (candV-succ).magnitude();
             double c = (pred-succ).magnitude();
             double s=0.5*(a+b+c);
             double hc=(2.0/c)*sqrt(fabs(s*(s-a)*(s-b)*(s-c)));
 
             if (hc<minError)
             {
               removeIndex = j;
               minError = hc;
               pointFound = true;
             }
           }
         }
       }
 
       if (pointFound)
       {
         removedPoints[removeIndex] = 1;
         remCounter++;
       }
     }
 
     for (int j=0; j<numPoints; j++)
     {
       if (removedPoints[j]<=0)
       {
 #pragma omp critical
         {
           vtkIdType id = vtkNewPoints->InsertNextPoint(vertices.at(j).data_block());
           container->GetPointIds()->InsertNextId(id);
         }
       }
     }
 
 #pragma omp critical
     {
       newFiberWeights->SetValue(vtkNewCells->GetNumberOfCells(), weight);
       numRemovedPoints += remCounter;
       vtkNewCells->InsertNextCell(container);
     }
   }
 
   if (vtkNewCells->GetNumberOfCells()>0)
   {
     MITK_INFO << "Removed points: " << numRemovedPoints;
     SetFiberWeights(newFiberWeights);
     m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
     m_FiberPolyData->SetPoints(vtkNewPoints);
     m_FiberPolyData->SetLines(vtkNewCells);
     this->SetFiberPolyData(m_FiberPolyData, true);
   }
 }
 
 void mitk::FiberBundle::ResampleLinear(double pointDistance)
 {
   vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
   vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
 
   MITK_INFO << "Resampling fibers (linear)";
   boost::progress_display disp(m_FiberPolyData->GetNumberOfCells());
   vtkSmartPointer<vtkFloatArray> newFiberWeights = vtkSmartPointer<vtkFloatArray>::New();
   newFiberWeights->SetName("FIBER_WEIGHTS");
   newFiberWeights->SetNumberOfValues(m_NumFibers);
 
 #pragma omp parallel for
   for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
   {
 
     std::vector< vnl_vector_fixed< double, 3 > > vertices;
     float weight = 1;
 
 #pragma omp critical
     {
       ++disp;
       weight = m_FiberWeights->GetValue(i);
       vtkCell* cell = m_FiberPolyData->GetCell(i);
       int numPoints = cell->GetNumberOfPoints();
       vtkPoints* points = cell->GetPoints();
 
       for (int j=0; j<numPoints; j++)
       {
         double cand[3];
         points->GetPoint(j, cand);
         vnl_vector_fixed< double, 3 > candV;
         candV[0]=cand[0]; candV[1]=cand[1]; candV[2]=cand[2];
         vertices.push_back(candV);
       }
     }
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     vnl_vector_fixed< double, 3 > lastV = vertices.at(0);
 
 #pragma omp critical
     {
       vtkIdType id = vtkNewPoints->InsertNextPoint(lastV.data_block());
       container->GetPointIds()->InsertNextId(id);
     }
     for (unsigned int j=1; j<vertices.size(); j++)
     {
       vnl_vector_fixed< double, 3 > vec = vertices.at(j) - lastV;
       double new_dist = vec.magnitude();
 
       if (new_dist >= pointDistance)
       {
         vnl_vector_fixed< double, 3 > newV = lastV;
         if ( new_dist-pointDistance <= mitk::eps )
         {
           vec.normalize();
           newV += vec * pointDistance;
         }
         else
         {
           // intersection between sphere (radius 'pointDistance', center 'lastV') and line (direction 'd' and point 'p')
           vnl_vector_fixed< double, 3 > p = vertices.at(j-1);
           vnl_vector_fixed< double, 3 > d = vertices.at(j) - p;
 
           double a = d[0]*d[0] + d[1]*d[1] + d[2]*d[2];
           double b = 2 * (d[0] * (p[0] - lastV[0]) + d[1] * (p[1] - lastV[1]) + d[2] * (p[2] - lastV[2]));
           double c = (p[0] - lastV[0])*(p[0] - lastV[0]) + (p[1] - lastV[1])*(p[1] - lastV[1]) + (p[2] - lastV[2])*(p[2] - lastV[2]) - pointDistance*pointDistance;
 
           double v1 =(-b + std::sqrt(b*b-4*a*c))/(2*a);
           double v2 =(-b - std::sqrt(b*b-4*a*c))/(2*a);
 
           if (v1>0)
             newV = p + d * v1;
           else if (v2>0)
             newV = p + d * v2;
           else
             MITK_INFO << "ERROR1 - linear resampling";
 
           j--;
         }
 
 #pragma omp critical
         {
           vtkIdType id = vtkNewPoints->InsertNextPoint(newV.data_block());
           container->GetPointIds()->InsertNextId(id);
         }
         lastV = newV;
       }
       else if (j==vertices.size()-1 && new_dist>0.0001)
       {
 #pragma omp critical
         {
           vtkIdType id = vtkNewPoints->InsertNextPoint(vertices.at(j).data_block());
           container->GetPointIds()->InsertNextId(id);
         }
       }
     }
 
 #pragma omp critical
     {
       newFiberWeights->SetValue(vtkNewCells->GetNumberOfCells(), weight);
       vtkNewCells->InsertNextCell(container);
     }
   }
 
   if (vtkNewCells->GetNumberOfCells()>0)
   {
     SetFiberWeights(newFiberWeights);
     m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
     m_FiberPolyData->SetPoints(vtkNewPoints);
     m_FiberPolyData->SetLines(vtkNewCells);
     this->SetFiberPolyData(m_FiberPolyData, true);
   }
 }
 
 // reapply selected colorcoding in case PolyData structure has changed
 bool mitk::FiberBundle::Equals(mitk::FiberBundle* fib, double eps)
 {
   if (fib==nullptr)
   {
     MITK_INFO << "Reference bundle is nullptr!";
     return false;
   }
 
   if (m_NumFibers!=fib->GetNumFibers())
   {
     MITK_INFO << "Unequal number of fibers!";
     MITK_INFO << m_NumFibers << " vs. " << fib->GetNumFibers();
     return false;
   }
 
   for (int i=0; i<m_NumFibers; i++)
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     vtkCell* cell2 = fib->GetFiberPolyData()->GetCell(i);
     int numPoints2 = cell2->GetNumberOfPoints();
     vtkPoints* points2 = cell2->GetPoints();
 
     if (numPoints2!=numPoints)
     {
       MITK_INFO << "Unequal number of points in fiber " << i << "!";
       MITK_INFO << numPoints2 << " vs. " << numPoints;
       return false;
     }
 
     for (int j=0; j<numPoints; j++)
     {
       double* p1 = points->GetPoint(j);
       double* p2 = points2->GetPoint(j);
       if (fabs(p1[0]-p2[0])>eps || fabs(p1[1]-p2[1])>eps || fabs(p1[2]-p2[2])>eps)
       {
         MITK_INFO << "Unequal points in fiber " << i << " at position " << j << "!";
         MITK_INFO << "p1: " << p1[0] << ", " << p1[1] << ", " << p1[2];
         MITK_INFO << "p2: " << p2[0] << ", " << p2[1] << ", " << p2[2];
         return false;
       }
     }
   }
 
   return true;
 }
 
 void mitk::FiberBundle::PrintSelf(std::ostream &os, itk::Indent indent) const
 {
   os << indent << this->GetNameOfClass() << ":\n";
   os << indent << "Number of fibers: " << this->GetNumFibers() << std::endl;
   os << indent << "Min. fiber length: " << this->GetMinFiberLength() << std::endl;
   os << indent << "Max. fiber length: " << this->GetMaxFiberLength() << std::endl;
   os << indent << "Mean fiber length: " << this->GetMeanFiberLength() << std::endl;
   os << indent << "Median fiber length: " << this->GetMedianFiberLength() << std::endl;
   os << indent << "STDEV fiber length: " << this->GetLengthStDev() << std::endl;
   os << indent << "Number of points: " << this->GetNumberOfPoints() << std::endl;
   Superclass::PrintSelf(os, indent);
 }
 
 /* ESSENTIAL IMPLEMENTATION OF SUPERCLASS METHODS */
 void mitk::FiberBundle::UpdateOutputInformation()
 {
 
 }
 void mitk::FiberBundle::SetRequestedRegionToLargestPossibleRegion()
 {
 
 }
 bool mitk::FiberBundle::RequestedRegionIsOutsideOfTheBufferedRegion()
 {
   return false;
 }
 bool mitk::FiberBundle::VerifyRequestedRegion()
 {
   return true;
 }
 void mitk::FiberBundle::SetRequestedRegion(const itk::DataObject* )
 {
 
 }
diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibilityFit.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/AnchorBasedScoring.cpp
similarity index 78%
rename from Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibilityFit.cpp
rename to Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/AnchorBasedScoring.cpp
index 19f16dc350..16e127144f 100755
--- a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibilityFit.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/AnchorBasedScoring.cpp
@@ -1,470 +1,444 @@
 /*===================================================================
 
 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>
 
 using namespace std;
 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< string > get_file_list(const std::string& path, std::vector< string > extensions={".fib", ".trk"})
 {
   std::vector< 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("");
+  parser.setTitle("Anchor Based Scoring");
   parser.setCategory("Fiber Tracking Evaluation");
   parser.setDescription("");
   parser.setContributor("MIC");
 
   parser.setArgumentPrefix("--", "-");
-  parser.addArgument("", "i1", mitkCommandLineParser::InputFile, "Anchor tracts:", "anchor tracts in one tractogram file", us::Any(), false);
+  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 masks for accuracy evaluation", false);
-
-  parser.addArgument("mask", "", mitkCommandLineParser::InputFile, "", "", us::Any());
-
-  parser.addArgument("greedy_add", "", mitkCommandLineParser::Bool, "", "", false);
-  parser.addArgument("normalize_gain", "", mitkCommandLineParser::Bool, "", "", 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);
 
   map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
   if (parsedArgs.size()==0)
     return EXIT_FAILURE;
 
-  string fib_file = us::any_cast<string>(parsedArgs["i1"]);
+  string anchors_file = us::any_cast<string>(parsedArgs["i1"]);
   string peak_file_name = us::any_cast<string>(parsedArgs["i2"]);
   string candidate_folder = us::any_cast<string>(parsedArgs["i3"]);
   string out_folder = us::any_cast<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 normalize_gain = false;
-  if (parsedArgs.count("normalize_gain"))
-    normalize_gain = us::any_cast<bool>(parsedArgs["normalize_gain"]);
-
   bool filter_outliers = true;
   if (parsedArgs.count("dont_filter_outliers"))
     filter_outliers = !us::any_cast<bool>(parsedArgs["dont_filter_outliers"]);
 
-
   string mask_file = "";
   if (parsedArgs.count("mask"))
     mask_file = us::any_cast<string>(parsedArgs["mask"]);
 
-
   string reference_mask_folder = "";
   if (parsedArgs.count("reference_mask_folder"))
     reference_mask_folder = us::any_cast<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";
 
     streambuf *old = cout.rdbuf(); // <-- save
     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< 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 reference tractogram consisting of all known tracts
-    std::vector< mitk::FiberBundle::Pointer > input_reference;
-    mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(mitk::IOUtil::Load(fib_file)[0].GetPointer());
-    if (fib.IsNull())
-      return EXIT_FAILURE;
-    fib->ResampleLinear(minSpacing/10.0);
-    input_reference.push_back(fib);
-
     // Load all candidate tracts
     std::vector< mitk::FiberBundle::Pointer > input_candidates;
     std::vector< string > candidate_tract_files = get_file_list(candidate_folder);
     for (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
 
-    // Fit known tracts to peak image to obtain underexplained image
-    MITK_INFO << "Fit anchor tracts";
+    double rmse = 0.0;
     int iteration = 0;
-    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();
-    std::string name = ist::GetFilenameWithoutExtension(fib_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<string>(iteration) + "_" + name + ".nrrd");
-    peak_image_writer->Update();
+    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) )
+    {
+      streambuf *old = cout.rdbuf(); // <-- save
+      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<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);
 
         string bundle_name = ist::GetFilenameWithoutExtension(candidate_tract_files.at(c));
-        try
-        {
-          streambuf *old = cout.rdbuf(); // <-- save
-          stringstream ss;
-          std::cout.rdbuf (ss.rdbuf());       // <-- redirect
-          mitk::IOUtil::Save(fib, out_folder + boost::lexical_cast<string>((int)(10000*rms_diff[c])) + "_" + bundle_name + ".fib");
-          std::cout.rdbuf (old);              // <-- restore
-        }
-        catch(...)
-        {
-          MITK_INFO << fib->GetNumFibers();
-          MITK_INFO << "ERROR!!!!!!!!!!!!!!";
-        }
 
         streambuf *old = cout.rdbuf(); // <-- save
         stringstream ss;
         std::cout.rdbuf (ss.rdbuf());       // <-- redirect
+        mitk::IOUtil::Save(fib, out_folder + boost::lexical_cast<string>((int)(10000*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
     {
-
-      double coverage = fitter->GetCoverage();
-      MITK_INFO << "Iteration: " << iteration << " - " << std::fixed << "Coverage: " << setprecision(2) << 100.0*coverage << "%";
-      logfile << "Iteration: " << iteration << " - " << std::fixed << "Coverage: " << setprecision(2) << 100.0*coverage << "%\n\n";
+      MITK_INFO << "RMSE: " << setprecision(5) << rmse;
       //    fitter->SetPeakImage(peak_image);
 
-
-      // Iteratively add/subtract candidate bundles in a greedy manner
+      // Iteratively add candidate bundles in a greedy manner
       while (!input_candidates.empty())
       {
-        double next_coverage = 0;
-        double next_covered_directions = 1.0;
+        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)});
 
           streambuf *old = cout.rdbuf(); // <-- save
           stringstream ss;
           std::cout.rdbuf (ss.rdbuf());       // <-- redirect
           fitter->Update();
           std::cout.rdbuf (old);              // <-- restore
 
-          double candidate_coverage = fitter->GetCoverage();
-          if (normalize_gain)
-            candidate_coverage /= fitter->GetNumCoveredDirections();
+          double candidate_rmse = fitter->GetRMSE();
 
-          if (candidate_coverage>next_coverage)
+          if (candidate_rmse<next_rmse)
           {
-            next_covered_directions = fitter->GetNumCoveredDirections();
-            next_coverage = candidate_coverage;
-
+            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< 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++;
         streambuf *old = cout.rdbuf(); // <-- save
         stringstream ss;
         std::cout.rdbuf (ss.rdbuf());       // <-- redirect
         // Save winning candidate
         mitk::IOUtil::Save(best_candidate, out_folder + boost::lexical_cast<string>(iteration) + "_" + name + ".fib");
         peak_image_writer->SetInput(peak_image);
         peak_image_writer->SetFileName(out_folder + boost::lexical_cast<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
 
-        if (normalize_gain)
-          next_coverage *= next_covered_directions;
-
-        MITK_INFO << "Iteration: " << iteration << " - " << std::fixed << "Coverage: " << setprecision(2) << 100.0*(coverage + (1.0-coverage) * next_coverage) << "% (+" << 100*(1.0-coverage) * next_coverage << "%)";
-        logfile << "Iteration: " << iteration << " - " << std::fixed << "Coverage: " << setprecision(2) << 100.0*(coverage + (1.0-coverage) * next_coverage) << "% (+" << 100*(1.0-coverage) * next_coverage << "%)\n";
-        coverage += (1.0-coverage) * next_coverage;
-
+        logfile << "RMSE: " << setprecision(5) << rmse << " " << name << " " << num_peaks << "\n";
         if (best_overlap_index>=0)
-        {
-          MITK_INFO << "Best overlap: " << best_overlap << " - " << ist::GetFilenameWithoutExtension(reference_mask_filenames.at(best_overlap_index));
-          logfile << "Best overlap: " << best_overlap << " - " << ist::GetFilenameWithoutExtension(reference_mask_filenames.at(best_overlap_index)) << "\n";
-        }
-        logfile << "\n";
+          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 760a65a867..a72860b7ab 100755
--- a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/CMakeLists.txt
+++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/CMakeLists.txt
@@ -1,44 +1,44 @@
 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
-    MergeDuplicateTracts^^MitkFiberTracking
-    TractPlausibility^^MitkFiberTracking
-    TractPlausibilityFit^^MitkFiberTracking
+    MergeOverlappingTracts^^MitkFiberTracking
+    ExtractAnchorTracts^^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/EvaluateLiFE.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/EvaluateLiFE.cpp
index 89d7ef5e9d..80e65bf60a 100755
--- a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/EvaluateLiFE.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/EvaluateLiFE.cpp
@@ -1,217 +1,227 @@
 /*===================================================================
 
 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>
 
 using namespace std;
 typedef itksys::SystemTools ist;
 typedef itk::Image<unsigned char, 3>    ItkUcharImgType;
 typedef std::tuple< ItkUcharImgType::Pointer, std::string > MaskType;
 
 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)
 {
   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();
 
   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 c = -1;
     for (int r : im_indices)
     {
       c++;
       const char *filename = dir->GetFile(r);
 
       MITK_INFO << "Loading " << ist::GetFilenameWithoutExtension(filename);
       streambuf *old = cout.rdbuf(); // <-- save
       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;
 }
 
 /*!
 \brief
 */
 int main(int argc, char* argv[])
 {
   mitkCommandLineParser parser;
 
   parser.setTitle("Evaluate LiFE results");
   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", false);
-  parser.addArgument("overlap", "", mitkCommandLineParser::Float, "", "", 0.7);
-  parser.addArgument("steps", "", mitkCommandLineParser::Int, "", "", 10);
+  parser.addArgument("overlap", "", mitkCommandLineParser::Float, "Overlap threshold:", "Overlap threshold used to identify true positives", 0.8);
+  parser.addArgument("steps", "", mitkCommandLineParser::Int, "Threshold steps:", "number of weight thresholds used to calculate the ROC curve", 100);
+  parser.addArgument("pre_filter_zeros", "", mitkCommandLineParser::Bool, "Remove zero weights:", "remove fibers with zero weights before starting the evaluation");
 
   map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
   if (parsedArgs.size()==0)
     return EXIT_FAILURE;
 
   string fibFile = us::any_cast<string>(parsedArgs["input"]);
   string reference_mask_folder = us::any_cast<string>(parsedArgs["reference_mask_folder"]);
   string out_folder = us::any_cast<string>(parsedArgs["out"]);
 
-  float overlap = 0.7;
+  float overlap = 0.8;
   if (parsedArgs.count("overlap"))
     overlap = us::any_cast<float>(parsedArgs["overlap"]);
 
   int steps = 10;
   if (parsedArgs.count("steps"))
     steps = us::any_cast<int>(parsedArgs["steps"]);
 
+  bool pre_filter_zeros = false;
+  if (parsedArgs.count("pre_filter_zeros"))
+    pre_filter_zeros = us::any_cast<bool>(parsedArgs["pre_filter_zeros"]);
+
   try
   {
     std::vector< MaskType > known_tract_masks = get_file_list(reference_mask_folder);
     if (known_tract_masks.empty())
       return EXIT_FAILURE;
 
     mitk::FiberBundle::Pointer inputTractogram = dynamic_cast<mitk::FiberBundle*>(mitk::IOUtil::Load(fibFile)[0].GetPointer());
 
     // resample fibers
     float minSpacing = 1;
     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);
 
     std::vector< float > weights;
     for (int i=0; i<inputTractogram->GetNumFibers(); ++i)
       weights.push_back(inputTractogram->GetFiberWeight(i));
     std::sort(weights.begin(), weights.end());
 
-    inputTractogram = inputTractogram->FilterByWeights(0.0);
+    if (pre_filter_zeros)
+      inputTractogram = inputTractogram->FilterByWeights(0.0);
     mitk::FiberBundle::Pointer pred_positives = inputTractogram->GetDeepCopy();
     mitk::FiberBundle::Pointer pred_negatives = mitk::FiberBundle::New(nullptr);
 
     ofstream logfile;
     logfile.open (out_folder + "LiFE_ROC.txt");
 
+    float fpr = 1.0;
+    float tpr = 1.0;
     float step = weights.back()/steps;
     float w = 0;
-    while (pred_positives->GetNumFibers()>0)
+    if (!pre_filter_zeros)
+      w -= step;
+    while (pred_positives->GetNumFibers()>0 && fpr>0.001 && tpr>0.001)
     {
       w += step;
 
       streambuf *old = cout.rdbuf(); // <-- save
       stringstream ss;
       std::cout.rdbuf (ss.rdbuf());       // <-- redirect
 
       mitk::FiberBundle::Pointer tp_tracts = mitk::FiberBundle::New(nullptr);
       mitk::FiberBundle::Pointer fn_tracts = mitk::FiberBundle::New(nullptr);
       for ( MaskType mask : known_tract_masks )
       {
         ItkUcharImgType::Pointer mask_image = std::get<0>(mask);
 
         mitk::FiberBundle::Pointer a = pred_positives->ExtractFiberSubset(mask_image, true, false, false, overlap, false);
         tp_tracts = tp_tracts->AddBundle(a);
 
         mitk::FiberBundle::Pointer b = pred_negatives->ExtractFiberSubset(mask_image, true, false, false, overlap, false);
         fn_tracts = fn_tracts->AddBundle(b);
       }
       mitk::FiberBundle::Pointer fp_tracts = pred_positives->SubtractBundle(tp_tracts);
       mitk::FiberBundle::Pointer tn_tracts = pred_negatives->SubtractBundle(fn_tracts);
 
       std::cout.rdbuf (old);              // <-- restore
 
       float positives = tp_tracts->GetNumFibers() + fn_tracts->GetNumFibers();
       float negatives = tn_tracts->GetNumFibers() + fp_tracts->GetNumFibers();
-      float fpr = (float)fp_tracts->GetNumFibers() / negatives;
-      float tpr = (float)tp_tracts->GetNumFibers() / positives;
+      fpr = (float)fp_tracts->GetNumFibers() / negatives;
+      tpr = (float)tp_tracts->GetNumFibers() / positives;
 
       float accuracy = 1.0;
       if (pred_positives->GetNumFibers()>0)
         accuracy = (float)tp_tracts->GetNumFibers()/pred_positives->GetNumFibers();
       logfile << w << " " << fpr << " " << tpr << " " << accuracy << " \n";
       MITK_INFO << "#Fibers: " << pred_positives->GetNumFibers();
       MITK_INFO << "FPR/TPR: " << fpr << "/" << tpr;
       MITK_INFO << "Acc: " << accuracy;
 
       pred_positives = inputTractogram->FilterByWeights(w);
       pred_negatives = inputTractogram->FilterByWeights(w, true);
     }
     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/TractPlausibility.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/ExtractAnchorTracts.cpp
similarity index 70%
rename from Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibility.cpp
rename to Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/ExtractAnchorTracts.cpp
index e6781a202b..8479b8344a 100755
--- a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibility.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/ExtractAnchorTracts.cpp
@@ -1,271 +1,269 @@
 /*===================================================================
 
 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>
 
 using namespace std;
 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 + "/known_tracts/");
+  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, int dropout, const std::string& skipped_path)
+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 = 0;
-    if (dropout>1)
-      skipping_num = (float)num_images/dropout;
-    if (dropout>=num_images)
-      skipping_num = (float)num_images/2.0;
-
+    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 && ist::GetFilenameWithoutExtension(filename)!="CC")
+      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);
       streambuf *old = cout.rdbuf(); // <-- save
       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("Tract Plausibility");
+  parser.setTitle("Extract Anchor 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", false);
-  parser.addArgument("gray_matter_mask", "gm", mitkCommandLineParser::String, "", "");
-  parser.addArgument("dropout", "", mitkCommandLineParser::Int, "", "", false);
-  parser.addArgument("overlap", "", mitkCommandLineParser::Float, "", "", false, 0.7);
-  parser.addArgument("subsample", "", mitkCommandLineParser::Float, "", "", false, 1.0);
+  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);
 
   map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
   if (parsedArgs.size()==0)
     return EXIT_FAILURE;
 
   string fibFile = us::any_cast<string>(parsedArgs["input"]);
   string reference_mask_folder = us::any_cast<string>(parsedArgs["reference_mask_folder"]);
   string out_folder = us::any_cast<string>(parsedArgs["out"]);
 
   string gray_matter_mask = "";
   if (parsedArgs.count("gray_matter_mask"))
     gray_matter_mask = us::any_cast<string>(parsedArgs["gray_matter_mask"]);
 
-  int dropout = 1;
-  if (parsedArgs.count("dropout"))
-    dropout = us::any_cast<int>(parsedArgs["dropout"]);
+  float anchor_fraction = 0.5;
+  if (parsedArgs.count("anchor_fraction"))
+    anchor_fraction = us::any_cast<float>(parsedArgs["anchor_fraction"]);
 
-  float overlap = 0.7;
+  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, dropout, out_folder + "/skipped_masks/");
-    if (known_tract_masks.empty())
-      return EXIT_FAILURE;
-
-
+    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)
     {
       streambuf *old = cout.rdbuf(); // <-- save
       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(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];
+    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_INFO << "Find known tracts via overlap match";
-    boost::progress_display disp(known_tract_masks.size());
-    mitk::FiberBundle::Pointer all_known_tracts = nullptr;
-    for ( MaskType mask : known_tract_masks )
+    mitk::FiberBundle::Pointer all_anchor_tracts = mitk::FiberBundle::New(nullptr);
+    if (!known_tract_masks.empty())
     {
-      ++disp;
-      streambuf *old = cout.rdbuf(); // <-- save
-      stringstream ss;
-      std::cout.rdbuf (ss.rdbuf());       // <-- redirect
+      MITK_INFO << "Find known tracts via overlap match";
+      boost::progress_display disp(known_tract_masks.size());
+      for ( MaskType mask : known_tract_masks )
+      {
+        ++disp;
+        streambuf *old = cout.rdbuf(); // <-- save
+        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 + "/known_tracts/" + mask_name + ".trk");
+        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");
 
-      if (all_known_tracts.IsNull())
-        all_known_tracts = mitk::FiberBundle::New(known_tract->GetFiberPolyData());
-      else
-        all_known_tracts = all_known_tracts->AddBundle(known_tract);
+        all_anchor_tracts = all_anchor_tracts->AddBundle(known_tract);
 
-      std::cout.rdbuf (old);              // <-- restore
+        std::cout.rdbuf (old);              // <-- restore
+      }
     }
-    mitk::IOUtil::Save(all_known_tracts, out_folder + "/known_tracts/all_known_tracts.trk");
-    mitk::IOUtil::Save(TransformToMRtrixSpace(all_known_tracts), out_folder + "/known_tracts/all_known_tracts_mrtrixspace.fib");
 
-    mitk::FiberBundle::Pointer remaining_tracts = inputTractogram->SubtractBundle(all_known_tracts);
-    mitk::IOUtil::Save(remaining_tracts, out_folder + "/candidate_tracts/remaining_tracts.trk");
-    mitk::IOUtil::Save(TransformToMRtrixSpace(remaining_tracts), out_folder + "/candidate_tracts/remaining_tracts_mrtrixspace.fib");
+    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/MergeDuplicateTracts.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/MergeOverlappingTracts.cpp
similarity index 88%
rename from Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/MergeDuplicateTracts.cpp
rename to Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/MergeOverlappingTracts.cpp
index ee7cb77cdf..673202ff03 100755
--- a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/MergeDuplicateTracts.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/MergeOverlappingTracts.cpp
@@ -1,250 +1,252 @@
 /*===================================================================
 
 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>
 #include <itkTractsToFiberEndingsImageFilter.h>
 #include <itkTractDensityImageFilter.h>
 #include <itkImageRegionConstIterator.h>
 #include <itkImageFileWriter.h>
 
 using namespace std;
 typedef itksys::SystemTools ist;
 typedef itk::Image<unsigned char, 3>    ItkUcharImgType;
 typedef itk::Image<unsigned int, 3>    ItkUIntImgType;
 
 std::vector< string > get_file_list(const std::string& path, std::vector< string > extensions={".fib", ".trk"})
 {
   std::vector< 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;
 }
 
 /*!
 \brief
 */
 int main(int argc, char* argv[])
 {
   mitkCommandLineParser parser;
 
-  parser.setTitle("Merge Duplicate Tracts");
+  parser.setTitle("Merge Overlapping Tracts");
   parser.setCategory("Fiber Tracking Evaluation");
   parser.setDescription("");
   parser.setContributor("MIC");
 
   parser.setArgumentPrefix("--", "-");
   parser.addArgument("in", "i", mitkCommandLineParser::InputFile, "Input Folder:", "input folder", us::Any(), false);
   parser.addArgument("out", "o", mitkCommandLineParser::OutputDirectory, "Output Folder:", "output folder", us::Any(), false);
-  parser.addArgument("overlap", "", mitkCommandLineParser::Float, "", "", false, 0.8);
+  parser.addArgument("overlap", "", mitkCommandLineParser::Float, "Overlap threshold:", "Tracts with overlap larger than this threshold are merged", false, 0.8);
 
   map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
   if (parsedArgs.size()==0)
     return EXIT_FAILURE;
 
   string input_folder = us::any_cast<string>(parsedArgs["in"]);
   string out_folder = us::any_cast<string>(parsedArgs["out"]);
 
   float overlap = 0.8;
   if (parsedArgs.count("overlap"))
     overlap = us::any_cast<float>(parsedArgs["overlap"]);
 
   try
   {
     if (ist::PathExists(out_folder))
       ist::RemoveADirectory(out_folder);
     ist::MakeDirectory(out_folder);
 
     std::vector< string > fib_files = get_file_list(input_folder, {".fib", ".trk", ".tck"});
     if (fib_files.empty())
       return EXIT_FAILURE;
 
+    streambuf *old = cout.rdbuf(); // <-- save
+    stringstream ss;
+    std::cout.rdbuf (ss.rdbuf());       // <-- redirect
+
     std::vector< mitk::FiberBundle::Pointer > fibs;
     for (string f : fib_files)
     {
       mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(mitk::IOUtil::Load(f)[0].GetPointer());
       fibs.push_back(fib);
     }
 
     mitk::FiberBundle::Pointer combined = mitk::FiberBundle::New();
     combined = combined->AddBundles(fibs);
 
     itk::TractsToFiberEndingsImageFilter< ItkUcharImgType >::Pointer endings = itk::TractsToFiberEndingsImageFilter< ItkUcharImgType >::New();
     endings->SetFiberBundle(combined);
     endings->SetUpsamplingFactor(0.25);
     endings->Update();
     ItkUcharImgType::Pointer ref_image = endings->GetOutput();
 
-    int original_size = fibs.size();
+    std::cout.rdbuf (old);              // <-- restore
+
     for (int its = 0; its<3; its++)
     {
+      streambuf *old = cout.rdbuf(); // <-- save
+      stringstream ss;
+      std::cout.rdbuf (ss.rdbuf());       // <-- redirect
+
       std::vector< ItkUcharImgType::Pointer > mask_images;
-      MITK_INFO << "Generating tract envelopes";
       for (auto fib : fibs)
       {
         itk::TractDensityImageFilter< ItkUcharImgType >::Pointer masks = itk::TractDensityImageFilter< ItkUcharImgType >::New();
         masks->SetInputImage(ref_image);
         masks->SetBinaryOutput(true);
         masks->SetFiberBundle(fib);
         masks->SetUseImageGeometry(true);
         masks->Update();
         mask_images.push_back(masks->GetOutput());
       }
 
-      MITK_INFO << "Calculating overlaps";
       int r=0;
       vnl_matrix< int > mat; mat.set_size(mask_images.size(), mask_images.size()); mat.fill(0);
       for (auto m1 : mask_images)
       {
-        //      itk::ImageFileWriter<ItkUcharImgType>::Pointer writer = itk::ImageFileWriter<ItkUcharImgType>::New();
-        //      writer->SetInput(m1);
-        //      writer->SetFileName(out_folder + "/mask_" + boost::lexical_cast<string>(r) + ".nii.gz");
-        //      writer->Update();
-
         float max_overlap = overlap;
         int c = 0;
         for (auto m2 : mask_images)
         {
           if (c<=r)
           {
             ++c;
             continue;
           }
 
           itk::ImageRegionConstIterator<ItkUcharImgType> it1(m1, m1->GetLargestPossibleRegion());
           itk::ImageRegionConstIterator<ItkUcharImgType> it2(m2, m2->GetLargestPossibleRegion());
 
           unsigned int c1 = 0;
           unsigned int c2 = 0;
           unsigned int intersect = 0;
           while( !it1.IsAtEnd() )
           {
             if( it1.Get()>0 && it2.Get()>0)
               ++intersect;
             if(it1.Get()>0)
               ++c1;
             if(it2.Get()>0)
               ++c2;
             ++it1;
             ++it2;
           }
 
-          //        MITK_INFO << "******************************";
-          //        MITK_INFO << r << "," << c;
-          //        MITK_INFO << intersect << "/" << c1+c2;
           if ( (float)intersect/c1>max_overlap )
           {
             max_overlap = (float)intersect/c1;
             mat.put(r,c, 1);
           }
           if ( (float)intersect/c2>max_overlap )
           {
             max_overlap = (float)intersect/c2;
             mat.put(r,c, 1);
           }
           ++c;
         }
 
         ++r;
       }
 
-      MITK_INFO << "Merging tracts";
       std::vector< mitk::FiberBundle::Pointer > out_fibs;
       std::vector< bool > used;
       for (unsigned int i=0; i<fibs.size(); i++)
         used.push_back(false);
       for (unsigned int r=0; r<mask_images.size(); r++)
       {
         if (used.at(r))
           continue;
 
-        MITK_INFO << "**********************";
-        MITK_INFO << r;
         mitk::FiberBundle::Pointer fib = fibs.at(r);
         for (unsigned int c=r+1; c<mask_images.size(); c++)
         {
           if (mat.get(r,c)>0)
           {
             fib = fib->AddBundle(fibs.at(c));
             MITK_INFO << c;
             used[c] = true;
           }
         }
 
         out_fibs.push_back(fib);
       }
+      std::cout.rdbuf (old);              // <-- restore
+
+      MITK_INFO << fibs.size() << " --> " << out_fibs.size();
 
       if (fibs.size()==out_fibs.size())
         break;
       fibs = out_fibs;
     }
-    int final_size = fibs.size();
 
-    MITK_INFO << original_size << " --> " << final_size;
     int c = 0;
     for (auto fib : fibs)
     {
+      streambuf *old = cout.rdbuf(); // <-- save
+      stringstream ss;
+      std::cout.rdbuf (ss.rdbuf());       // <-- redirect
       mitk::IOUtil::Save(fib, out_folder + "/bundle_" + boost::lexical_cast<string>(c) + ".trk");
+      std::cout.rdbuf (old);              // <-- restore
+
       ++c;
     }
   }
   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;
 }