diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.cpp
index a313be04f9..2d9645a07a 100644
--- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.cpp
@@ -1,971 +1,971 @@
 #include "itkFitFibersToImageFilter.h"
 
 #include <boost/progress.hpp>
 
 namespace itk{
 
 FitFibersToImageFilter::FitFibersToImageFilter()
   : m_PeakImage(nullptr)
   , m_DiffImage(nullptr)
   , m_ScalarImage(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(false)
   , 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)
   , m_SignalModel(nullptr)
   , sz_x(0)
   , sz_y(0)
   , sz_z(0)
   , m_MeanTractDensity(0)
   , m_MeanSignal(0)
   , fiber_count(0)
   , m_Regularization(VnlCostFunction::REGU::VOXEL_VARIANCE)
 {
   this->SetNumberOfRequiredOutputs(3);
 }
 
 FitFibersToImageFilter::~FitFibersToImageFilter()
 {
 
 }
 
 void FitFibersToImageFilter::CreateDiffSystem()
 {
   sz_x = m_DiffImage->GetLargestPossibleRegion().GetSize(0);
   sz_y = m_DiffImage->GetLargestPossibleRegion().GetSize(1);
   sz_z = m_DiffImage->GetLargestPossibleRegion().GetSize(2);
   dim_four_size = m_DiffImage->GetVectorLength();
   int num_voxels = sz_x*sz_y*sz_z;
 
   float minSpacing = 1;
   if(m_DiffImage->GetSpacing()[0]<m_DiffImage->GetSpacing()[1] && m_DiffImage->GetSpacing()[0]<m_DiffImage->GetSpacing()[2])
     minSpacing = m_DiffImage->GetSpacing()[0];
   else if (m_DiffImage->GetSpacing()[1] < m_DiffImage->GetSpacing()[2])
     minSpacing = m_DiffImage->GetSpacing()[1];
   else
     minSpacing = m_DiffImage->GetSpacing()[2];
 
   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 * dim_four_size;
 
   MITK_INFO << "Num. unknowns: " << m_NumUnknowns;
   MITK_INFO << "Num. residuals: " << m_NumResiduals;
   MITK_INFO << "Creating system ...";
 
   A.set_size(m_NumResiduals, m_NumUnknowns);
   b.set_size(m_NumResiduals); b.fill(0.0);
 
   m_MeanTractDensity = 0;
   m_MeanSignal = 0;
   m_NumCoveredDirections = 0;
   fiber_count = 0;
   vnl_vector<int> voxel_indicator; voxel_indicator.set_size(sz_x*sz_y*sz_z); voxel_indicator.fill(0);
 
   m_GroupSizes.clear();
   for (unsigned int bundle=0; bundle<m_Tractograms.size(); bundle++)
   {
     vtkSmartPointer<vtkPolyData> polydata = m_Tractograms.at(bundle)->GetFiberPolyData();
     m_GroupSizes.push_back(m_Tractograms.at(bundle)->GetNumFibers());
     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);
         PointType3 p;
         p[0]=p1[0];
         p[1]=p1[1];
         p[2]=p1[2];
 
         itk::Index<3> idx3;
         m_DiffImage->TransformPhysicalPointToIndex(p, idx3);
         if (!m_DiffImage->GetLargestPossibleRegion().IsInside(idx3) || (m_MaskImage.IsNotNull() && m_MaskImage->GetPixel(idx3)==0))
           continue;
 
         double* p2 = points->GetPoint(j+1);
         mitk::DiffusionSignalModel<>::GradientType 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();
 
         int x = idx3[0];
         int y = idx3[1];
         int z = idx3[2];
 
         mitk::DiffusionSignalModel<>::PixelType simulated_pixel = m_SignalModel->SimulateMeasurement(fiber_dir);
         VectorImgType::PixelType measured_pixel = m_DiffImage->GetPixel(idx3);
 
         double simulated_mean = 0;
         double measured_mean = 0;
         int num_nonzero_g = 0;
         for (int g=0; g<dim_four_size; ++g)
         {
           if( m_SignalModel->GetGradientDirection(g).GetNorm()<mitk::eps )
             continue;
           simulated_mean += simulated_pixel[g];
           measured_mean += (double)measured_pixel[g];
           ++num_nonzero_g;
         }
         simulated_mean /= num_nonzero_g;
         measured_mean /= num_nonzero_g;
         simulated_pixel -= simulated_mean;
 
         if (voxel_indicator[x + sz_x*y + sz_x*sz_y*z]==0)
           m_MeanSignal += measured_mean;
         m_MeanTractDensity += simulated_mean;
         voxel_indicator[x + sz_x*y + sz_x*sz_y*z] = 1;
 
         for (int g=0; g<dim_four_size; ++g)
         {
           unsigned int linear_index = x + sz_x*y + sz_x*sz_y*z + sz_x*sz_y*sz_z*g;
 
           if (m_FitIndividualFibers)
           {
             b[linear_index] = (double)measured_pixel[g] - measured_mean;
             A.put(linear_index, fiber_count, A.get(linear_index, fiber_count) + simulated_pixel[g]);
           }
           else
           {
             b[linear_index] = (double)measured_pixel[g] - measured_mean;
             A.put(linear_index, bundle, A.get(linear_index, bundle) + simulated_pixel[g]);
           }
         }
 
       }
 
       ++fiber_count;
     }
   }
 
   m_NumCoveredDirections = voxel_indicator.sum();
 
   m_MeanTractDensity /= (m_NumCoveredDirections*fiber_count);
   m_MeanSignal /= m_NumCoveredDirections;
   A /= m_MeanTractDensity;
   b *= 100.0/m_MeanSignal;  // times 100 because we want to avoid too small values for computational reasons
 
   // NEW FIT
   //  m_MeanTractDensity /= m_NumCoveredDirections;
   //  m_MeanSignal /= m_NumCoveredDirections;
   //  b /= m_MeanSignal;
   //  b *= m_MeanTractDensity;
 }
 
 void FitFibersToImageFilter::CreatePeakSystem()
 {
   sz_x = m_PeakImage->GetLargestPossibleRegion().GetSize(0);
   sz_y = m_PeakImage->GetLargestPossibleRegion().GetSize(1);
   sz_z = m_PeakImage->GetLargestPossibleRegion().GetSize(2);
   dim_four_size = 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];
 
   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 * dim_four_size;
 
   MITK_INFO << "Num. unknowns: " << m_NumUnknowns;
   MITK_INFO << "Num. residuals: " << m_NumResiduals;
   MITK_INFO << "Creating system ...";
 
   A.set_size(m_NumResiduals, m_NumUnknowns);
   b.set_size(m_NumResiduals); b.fill(0.0);
 
   m_MeanTractDensity = 0;
   m_MeanSignal = 0;
   m_NumCoveredDirections = 0;
   fiber_count = 0;
 
   m_GroupSizes.clear();
   for (unsigned int bundle=0; bundle<m_Tractograms.size(); bundle++)
   {
     vtkSmartPointer<vtkPolyData> polydata = m_Tractograms.at(bundle)->GetFiberPolyData();
 
     m_GroupSizes.push_back(m_Tractograms.at(bundle)->GetNumFibers());
     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 = dim_four_size-1;
 
         double peak_mag = 0;
         GetClosestPeak(idx4, m_PeakImage, fiber_dir, peak_id, w, peak_mag);
 
         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<dim_four_size-1)
         {
           m_NumCoveredDirections++;
           m_MeanSignal += peak_mag;
         }
         m_MeanTractDensity += 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;
     }
   }
   m_MeanTractDensity /= (m_NumCoveredDirections*fiber_count);
   m_MeanSignal /= m_NumCoveredDirections;
   A /= m_MeanTractDensity;
   b *= 100.0/m_MeanSignal;  // times 100 because we want to avoid too small values for computational reasons
 
   // NEW FIT
   //  m_MeanTractDensity /= m_NumCoveredDirections;
   //  m_MeanSignal /= m_NumCoveredDirections;
   //  b /= m_MeanSignal;
   //  b *= m_MeanTractDensity;
 }
 
 void FitFibersToImageFilter::CreateScalarSystem()
 {
   sz_x = m_ScalarImage->GetLargestPossibleRegion().GetSize(0);
   sz_y = m_ScalarImage->GetLargestPossibleRegion().GetSize(1);
   sz_z = m_ScalarImage->GetLargestPossibleRegion().GetSize(2);
   int num_voxels = sz_x*sz_y*sz_z;
 
   float minSpacing = 1;
   if(m_ScalarImage->GetSpacing()[0]<m_ScalarImage->GetSpacing()[1] && m_ScalarImage->GetSpacing()[0]<m_ScalarImage->GetSpacing()[2])
     minSpacing = m_ScalarImage->GetSpacing()[0];
   else if (m_ScalarImage->GetSpacing()[1] < m_ScalarImage->GetSpacing()[2])
     minSpacing = m_ScalarImage->GetSpacing()[1];
   else
     minSpacing = m_ScalarImage->GetSpacing()[2];
 
   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;
 
   MITK_INFO << "Num. unknowns: " << m_NumUnknowns;
   MITK_INFO << "Num. residuals: " << m_NumResiduals;
   MITK_INFO << "Creating system ...";
 
   A.set_size(m_NumResiduals, m_NumUnknowns);
   b.set_size(m_NumResiduals); b.fill(0.0);
 
   m_MeanTractDensity = 0;
   m_MeanSignal = 0;
   int numCoveredVoxels = 0;
   fiber_count = 0;
 
   m_GroupSizes.clear();
   for (unsigned int bundle=0; bundle<m_Tractograms.size(); bundle++)
   {
     vtkSmartPointer<vtkPolyData> polydata = m_Tractograms.at(bundle)->GetFiberPolyData();
 
     m_GroupSizes.push_back(m_Tractograms.at(bundle)->GetNumFibers());
     for (int i=0; i<m_Tractograms.at(bundle)->GetNumFibers(); ++i)
     {
       vtkCell* cell = polydata->GetCell(i);
       int numPoints = cell->GetNumberOfPoints();
       vtkPoints* points = cell->GetPoints();
 
       for (int j=0; j<numPoints; ++j)
       {
         double* p1 = points->GetPoint(j);
         PointType3 p;
         p[0]=p1[0];
         p[1]=p1[1];
         p[2]=p1[2];
 
         itk::Index<3> idx3;
         m_ScalarImage->TransformPhysicalPointToIndex(p, idx3);
         if (!m_ScalarImage->GetLargestPossibleRegion().IsInside(idx3) || (m_MaskImage.IsNotNull() && m_MaskImage->GetPixel(idx3)==0))
           continue;
 
         float image_value = m_ScalarImage->GetPixel(idx3);
         int x = idx3[0];
         int y = idx3[1];
         int z = idx3[2];
 
         unsigned int linear_index = x + sz_x*y + sz_x*sz_y*z;
 
         if (b[linear_index] == 0)
         {
           numCoveredVoxels++;
           m_MeanSignal += image_value;
         }
         m_MeanTractDensity += 1;
 
         if (m_FitIndividualFibers)
         {
           b[linear_index] = image_value;
           A.put(linear_index, fiber_count, A.get(linear_index, fiber_count) + 1);
         }
         else
         {
           b[linear_index] = image_value;
           A.put(linear_index, bundle, A.get(linear_index, bundle) + 1);
         }
       }
 
       ++fiber_count;
     }
   }
   m_MeanTractDensity /= (numCoveredVoxels*fiber_count);
   m_MeanSignal /= numCoveredVoxels;
   A /= m_MeanTractDensity;
   b *= 100.0/m_MeanSignal;  // times 100 because we want to avoid too small values for computational reasons
 
   // NEW FIT
   //  m_MeanTractDensity /= numCoveredVoxels;
   //  m_MeanSignal /= numCoveredVoxels;
   //  b /= m_MeanSignal;
   //  b *= m_MeanTractDensity;
 }
 
 void FitFibersToImageFilter::GenerateData()
 {
   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();
   }
   else
     m_FilterOutliers = false;
   if (m_NumUnknowns<1)
   {
     MITK_INFO << "No fibers in tractogram.";
     return;
   }
 
   fiber_count = 0;
   sz_x = 0;
   sz_y = 0;
   sz_z = 0;
   m_MeanTractDensity = 0;
   m_MeanSignal = 0;
   if (m_PeakImage.IsNotNull())
     CreatePeakSystem();
   else if (m_DiffImage.IsNotNull())
     CreateDiffSystem();
   else if (m_ScalarImage.IsNotNull())
     CreateScalarSystem();
   else
     mitkThrow() << "No input image set!";
 
   double init_lambda = fiber_count;  // initialization for lambda estimation
 
   itk::TimeProbe clock;
   clock.Start();
 
   cost = VnlCostFunction(m_NumUnknowns);
   cost.SetProblem(A, b, init_lambda, m_Regularization);
   cost.SetGroupSizes(m_GroupSizes);
   m_Weights.set_size(m_NumUnknowns);
   m_Weights.fill( 1.0/m_NumUnknowns );
   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);
 
   if (m_Regularization==VnlCostFunction::REGU::MSM)
     MITK_INFO << "Regularization type: MSM";
   else if (m_Regularization==VnlCostFunction::REGU::VARIANCE)
     MITK_INFO << "Regularization type: VARIANCE";
   else if (m_Regularization==VnlCostFunction::REGU::LASSO)
     MITK_INFO << "Regularization type: LASSO";
   else if (m_Regularization==VnlCostFunction::REGU::VOXEL_VARIANCE)
     MITK_INFO << "Regularization type: VOXEL_VARIANCE";
   else if (m_Regularization==VnlCostFunction::REGU::GROUP_LASSO)
     MITK_INFO << "Regularization type: GROUP_LASSO";
   else if (m_Regularization==VnlCostFunction::REGU::GROUP_VARIANCE)
     MITK_INFO << "Regularization type: GROUP_VARIANCE";
   else if (m_Regularization==VnlCostFunction::REGU::NONE)
     MITK_INFO << "Regularization type: NONE";
 
   if (m_Regularization!=VnlCostFunction::REGU::NONE)  // REMOVE FOR NEW FIT AND SET cost.m_Lambda = m_Lambda
   {
     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.calc_regularization_gradient(m_Weights, dx);
 
     if (m_Weights.magnitude()==0)
     {
       MITK_INFO << "Regularization estimation failed. Using default value.";
       cost.m_Lambda = fiber_count*m_Lambda;
     }
     else
     {
       double r = dx.magnitude()/m_Weights.magnitude();  // wtf???
       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*m_Lambda;
       }
     }
   }
   else
     cost.m_Lambda = 0;
   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 << "Sum: " << m_Weights.sum();
   MITK_INFO << "Mean: " << m_MeanWeight;
   MITK_INFO << "1% quantile: " << weights.at(m_NumUnknowns*0.01);
   MITK_INFO << "5% quantile: " << weights.at(m_NumUnknowns*0.05);
   MITK_INFO << "25% quantile: " << weights.at(m_NumUnknowns*0.25);
   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();
   m_RMSE = cost.S->get_rms_error(m_Weights);
-  MITK_INFO << "Final RMS: " << m_RMSE;
+  MITK_INFO << "Final RMSE: " << 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_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++)
     {
       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;
 
         ++fiber_count;
       }
       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)
     {
       vnl_vector< double > temp_weights;
       temp_weights.set_size(m_Weights.size());
       temp_weights.copy_in(m_Weights.data_block());
       temp_weights[i] = 0;
       double d_rms = cost.S->get_rms_error(temp_weights) - m_RMSE;
       m_RmsDiffPerBundle[i] = d_rms;
 
       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
   A *= m_MeanSignal/100.0;
   b *= m_MeanSignal/100.0;
 
   MITK_INFO << "Generating output images ...";
   if (m_PeakImage.IsNotNull())
     GenerateOutputPeakImages();
   else if (m_DiffImage.IsNotNull())
     GenerateOutputDiffImages();
   else if (m_ScalarImage.IsNotNull())
     GenerateOutputScalarImages();
 
   m_Coverage = m_Coverage/m_MeanSignal;
   m_Overshoot = m_Overshoot/m_MeanSignal;
 
   MITK_INFO << std::fixed << "Coverage: " << setprecision(2) << 100.0*m_Coverage << "%";
   MITK_INFO << std::fixed << "Overshoot: " << setprecision(2) << 100.0*m_Overshoot << "%";
 }
 
 void FitFibersToImageFilter::GenerateOutputDiffImages()
 {
   VectorImgType::PixelType pix; pix.SetSize(m_DiffImage->GetVectorLength()); pix.Fill(0);
   itk::ImageDuplicator< VectorImgType >::Pointer duplicator = itk::ImageDuplicator< VectorImgType >::New();
   duplicator->SetInputImage(m_DiffImage);
   duplicator->Update();
   m_UnderexplainedImageDiff = duplicator->GetOutput();
   m_UnderexplainedImageDiff->FillBuffer(pix);
 
   duplicator->SetInputImage(m_UnderexplainedImageDiff);
   duplicator->Update();
   m_OverexplainedImageDiff = duplicator->GetOutput();
   m_OverexplainedImageDiff->FillBuffer(pix);
 
   duplicator->SetInputImage(m_OverexplainedImageDiff);
   duplicator->Update();
   m_ResidualImageDiff = duplicator->GetOutput();
   m_ResidualImageDiff->FillBuffer(pix);
 
   duplicator->SetInputImage(m_ResidualImageDiff);
   duplicator->Update();
   m_FittedImageDiff = duplicator->GetOutput();
   m_FittedImageDiff->FillBuffer(pix);
 
   vnl_vector<double> fitted_b; fitted_b.set_size(b.size());
   cost.S->multiply(m_Weights, fitted_b);
 
   itk::ImageRegionIterator<VectorImgType> it1 = itk::ImageRegionIterator<VectorImgType>(m_DiffImage, m_DiffImage->GetLargestPossibleRegion());
   itk::ImageRegionIterator<VectorImgType> it2 = itk::ImageRegionIterator<VectorImgType>(m_FittedImageDiff, m_FittedImageDiff->GetLargestPossibleRegion());
   itk::ImageRegionIterator<VectorImgType> it3 = itk::ImageRegionIterator<VectorImgType>(m_ResidualImageDiff, m_ResidualImageDiff->GetLargestPossibleRegion());
   itk::ImageRegionIterator<VectorImgType> it4 = itk::ImageRegionIterator<VectorImgType>(m_UnderexplainedImageDiff, m_UnderexplainedImageDiff->GetLargestPossibleRegion());
   itk::ImageRegionIterator<VectorImgType> it5 = itk::ImageRegionIterator<VectorImgType>(m_OverexplainedImageDiff, m_OverexplainedImageDiff->GetLargestPossibleRegion());
 
   m_MeanSignal = 0;
   m_Coverage = 0;
   m_Overshoot = 0;
 
   while( !it2.IsAtEnd() )
   {
     itk::Index<3> idx3 = it2.GetIndex();
     VectorImgType::PixelType original_pix =it1.Get();
     VectorImgType::PixelType fitted_pix =it2.Get();
     VectorImgType::PixelType residual_pix =it3.Get();
     VectorImgType::PixelType underexplained_pix =it4.Get();
     VectorImgType::PixelType overexplained_pix =it5.Get();
 
     int num_nonzero_g = 0;
     double original_mean = 0;
     for (int g=0; g<dim_four_size; ++g)
     {
       if( m_SignalModel->GetGradientDirection(g).GetNorm()>=mitk::eps )
       {
         original_mean += original_pix[g];
         ++num_nonzero_g;
       }
     }
     original_mean /= num_nonzero_g;
 
     for (int g=0; g<dim_four_size; ++g)
     {
       unsigned int linear_index = idx3[0] + sz_x*idx3[1] + sz_x*sz_y*idx3[2] + sz_x*sz_y*sz_z*g;
 
       fitted_pix[g] = fitted_b[linear_index] + original_mean;
       residual_pix[g] = original_pix[g] - fitted_b[linear_index] - original_mean;
 
       if (residual_pix[g]<0)
       {
         overexplained_pix[g] = residual_pix[g];
         m_Coverage += b[linear_index] + original_mean;
         m_Overshoot -= residual_pix[g];
       }
       else if (residual_pix[g]>=0)
       {
         underexplained_pix[g] = residual_pix[g];
         m_Coverage += fitted_b[linear_index] + original_mean;
       }
       m_MeanSignal += b[linear_index] + original_mean;
     }
 
     it2.Set(fitted_pix);
     it3.Set(residual_pix);
     it4.Set(underexplained_pix);
     it5.Set(overexplained_pix);
 
     ++it1;
     ++it2;
     ++it3;
     ++it4;
     ++it5;
   }
 }
 
 void FitFibersToImageFilter::GenerateOutputScalarImages()
 {
   itk::ImageDuplicator< DoubleImgType >::Pointer duplicator = itk::ImageDuplicator< DoubleImgType >::New();
   duplicator->SetInputImage(m_ScalarImage);
   duplicator->Update();
   m_UnderexplainedImageScalar = duplicator->GetOutput();
   m_UnderexplainedImageScalar->FillBuffer(0);
 
   duplicator->SetInputImage(m_UnderexplainedImageScalar);
   duplicator->Update();
   m_OverexplainedImageScalar = duplicator->GetOutput();
   m_OverexplainedImageScalar->FillBuffer(0);
 
   duplicator->SetInputImage(m_OverexplainedImageScalar);
   duplicator->Update();
   m_ResidualImageScalar = duplicator->GetOutput();
   m_ResidualImageScalar->FillBuffer(0);
 
   duplicator->SetInputImage(m_ResidualImageScalar);
   duplicator->Update();
   m_FittedImageScalar = duplicator->GetOutput();
   m_FittedImageScalar->FillBuffer(0);
 
   vnl_vector<double> fitted_b; fitted_b.set_size(b.size());
   cost.S->multiply(m_Weights, fitted_b);
 
   itk::ImageRegionIterator<DoubleImgType> it1 = itk::ImageRegionIterator<DoubleImgType>(m_ScalarImage, m_ScalarImage->GetLargestPossibleRegion());
   itk::ImageRegionIterator<DoubleImgType> it2 = itk::ImageRegionIterator<DoubleImgType>(m_FittedImageScalar, m_FittedImageScalar->GetLargestPossibleRegion());
   itk::ImageRegionIterator<DoubleImgType> it3 = itk::ImageRegionIterator<DoubleImgType>(m_ResidualImageScalar, m_ResidualImageScalar->GetLargestPossibleRegion());
   itk::ImageRegionIterator<DoubleImgType> it4 = itk::ImageRegionIterator<DoubleImgType>(m_UnderexplainedImageScalar, m_UnderexplainedImageScalar->GetLargestPossibleRegion());
   itk::ImageRegionIterator<DoubleImgType> it5 = itk::ImageRegionIterator<DoubleImgType>(m_OverexplainedImageScalar, m_OverexplainedImageScalar->GetLargestPossibleRegion());
 
   m_MeanSignal = 0;
   m_Coverage = 0;
   m_Overshoot = 0;
 
   while( !it2.IsAtEnd() )
   {
     itk::Index<3> idx3 = it2.GetIndex();
     DoubleImgType::PixelType original_pix =it1.Get();
     DoubleImgType::PixelType fitted_pix =it2.Get();
     DoubleImgType::PixelType residual_pix =it3.Get();
     DoubleImgType::PixelType underexplained_pix =it4.Get();
     DoubleImgType::PixelType overexplained_pix =it5.Get();
 
     unsigned int linear_index = idx3[0] + sz_x*idx3[1] + sz_x*sz_y*idx3[2];
 
     fitted_pix = fitted_b[linear_index];
     residual_pix = original_pix - fitted_pix;
 
     if (residual_pix<0)
     {
       overexplained_pix = residual_pix;
       m_Coverage += b[linear_index];
       m_Overshoot -= residual_pix;
     }
     else if (residual_pix>=0)
     {
       underexplained_pix = residual_pix;
       m_Coverage += fitted_b[linear_index];
     }
     m_MeanSignal += b[linear_index];
 
     it2.Set(fitted_pix);
     it3.Set(residual_pix);
     it4.Set(underexplained_pix);
     it5.Set(overexplained_pix);
 
     ++it1;
     ++it2;
     ++it3;
     ++it4;
     ++it5;
   }
 }
 
 VnlCostFunction::REGU FitFibersToImageFilter::GetRegularization() const
 {
   return m_Regularization;
 }
 
 void FitFibersToImageFilter::SetRegularization(const VnlCostFunction::REGU &Regularization)
 {
   m_Regularization = Regularization;
 }
 
 void FitFibersToImageFilter::GenerateOutputPeakImages()
 {
   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 % dim_four_size;
 
     if (peak_id<dim_four_size-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]);
     }
   }
 
   m_MeanSignal = 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)
           {
             m_MeanSignal += 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]);
             }
           }
         }
       }
 }
 
 void FitFibersToImageFilter::GetClosestPeak(itk::Index<4> idx, PeakImgType::Pointer peak_image , vnl_vector_fixed<float,3> fiber_dir, int& id, double& w, double& peak_mag )
 {
   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;
 
     float a = fabs(dot_product(dir, fiber_dir))/mag;
     if (a>angle)
     {
       angle = a;
       w = angle;
       peak_mag = mag;
       id = i;
     }
   }
 }
 
 std::vector<mitk::FiberBundle::Pointer> FitFibersToImageFilter::GetTractograms() const
 {
   return m_Tractograms;
 }
 
 void FitFibersToImageFilter::SetTractograms(const std::vector<mitk::FiberBundle::Pointer> &tractograms)
 {
   m_Tractograms = tractograms;
 }
 
 void FitFibersToImageFilter::SetSignalModel(mitk::DiffusionSignalModel<> *SignalModel)
 {
   m_SignalModel = SignalModel;
 }
 
 }
 
 
 
diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/AnchorBasedScoring.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/AnchorBasedScoring.cpp
index 35bb3ffafe..357c197a03 100755
--- a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/AnchorBasedScoring.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/AnchorBasedScoring.cpp
@@ -1,473 +1,573 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 
 #include <mitkBaseData.h>
 #include <mitkImageCast.h>
 #include <mitkImageToItk.h>
 #include <metaCommand.h>
 #include <mitkCommandLineParser.h>
 #include <usAny.h>
 #include <mitkIOUtil.h>
 #include <itksys/SystemTools.hxx>
 #include <itkDirectory.h>
 #include <mitkFiberBundle.h>
 #include <mitkPreferenceListReaderOptionsFunctor.h>
 #include <itkImageFileWriter.h>
 #include <mitkPeakImage.h>
 #include <itkFitFibersToImageFilter.h>
 #include <boost/lexical_cast.hpp>
 #include <itkTractDensityImageFilter.h>
 
 typedef itksys::SystemTools ist;
 typedef itk::Point<float, 4> PointType4;
 typedef itk::Image< float, 4 >  PeakImgType;
 typedef itk::Image< unsigned char, 3 >  ItkUcharImageType;
 
 std::vector< 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;
 }
 
 std::vector< std::string > get_file_list(const std::string& path, std::vector< std::string > extensions={".fib", ".trk"})
 {
   std::vector< std::string > file_list;
   itk::Directory::Pointer dir = itk::Directory::New();
 
   if (dir->Load(path.c_str()))
   {
     int n = dir->GetNumberOfFiles();
     for (int r = 0; r < n; r++)
     {
       const char *filename = dir->GetFile(r);
       std::string ext = ist::GetFilenameExtension(filename);
       for (auto e : extensions)
       {
         if (ext==e)
         {
           file_list.push_back(path + '/' + filename);
           break;
         }
       }
     }
   }
   return file_list;
 }
 
 /*!
 \brief Fits the tractogram to the input peak image by assigning a weight to each fiber (similar to https://doi.org/10.1016/j.neuroimage.2015.06.092).
 */
 int main(int argc, char* argv[])
 {
   mitkCommandLineParser parser;
 
   parser.setTitle("Anchor Based Scoring");
   parser.setCategory("Fiber Tracking Evaluation");
   parser.setDescription("");
   parser.setContributor("MIC");
 
   parser.setArgumentPrefix("--", "-");
   parser.addArgument("", "a", mitkCommandLineParser::InputFile, "Anchor tractogram:", "anchor tracts in one tractogram file", us::Any(), false);
   parser.addArgument("", "p", mitkCommandLineParser::InputFile, "Input peaks:", "input peak image", us::Any(), false);
   parser.addArgument("", "c", mitkCommandLineParser::InputDirectory, "Candidates folder:", "folder containing candidate tracts", us::Any(), false);
   parser.addArgument("", "o", mitkCommandLineParser::OutputDirectory, "Output folder:", "output folder", us::Any(), false);
 
   parser.addArgument("anchor_masks", "", mitkCommandLineParser::StringList, "Reference Masks:", "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("filter_outliers", "", mitkCommandLineParser::Bool, "Filter outliers:", "perform second optimization run with an upper weight bound based on the first weight estimation (99% quantile)", false);
   parser.addArgument("regu", "", mitkCommandLineParser::String, "Regularization:", "MSM, Variance, VoxelVariance, Lasso, GroupLasso, GroupVariance, NONE (default)");
+  parser.addArgument("use_num_streamlines", "", mitkCommandLineParser::Bool, "Use number of streamlines as score:", "Don't fit candidates, simply use number of streamlines per candidate as score", false);
+  parser.addArgument("use_weights", "", mitkCommandLineParser::Bool, "Use input weights as score:", "Don't fit candidates, simply use first input streamline weight per candidate as score", false);
+  parser.addArgument("filter_zero_weights", "", mitkCommandLineParser::Bool, "Filter zero-weights", "Remove streamlines with weight 0 from candidates", false);
 
   std::map<std::string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
   if (parsedArgs.size()==0)
     return EXIT_FAILURE;
 
   std::string anchors_file = us::any_cast<std::string>(parsedArgs["a"]);
   std::string peak_file_name = us::any_cast<std::string>(parsedArgs["p"]);
   std::string candidate_tract_folder = us::any_cast<std::string>(parsedArgs["c"]);
   std::string out_folder = us::any_cast<std::string>(parsedArgs["o"]);
 
   bool greedy_add = false;
   if (parsedArgs.count("greedy_add"))
     greedy_add = us::any_cast<bool>(parsedArgs["greedy_add"]);
 
   float lambda = 0.1;
   if (parsedArgs.count("lambda"))
     lambda = us::any_cast<float>(parsedArgs["lambda"]);
 
   bool filter_outliers = false;
   if (parsedArgs.count("filter_outliers"))
     filter_outliers = us::any_cast<bool>(parsedArgs["filter_outliers"]);
 
+  bool filter_zero_weights = false;
+  if (parsedArgs.count("filter_zero_weights"))
+    filter_zero_weights = us::any_cast<bool>(parsedArgs["filter_zero_weights"]);
+
   std::string mask_file = "";
   if (parsedArgs.count("mask"))
     mask_file = us::any_cast<std::string>(parsedArgs["mask"]);
 
-  mitkCommandLineParser::StringContainerType anchor_mask_files;
+  mitkCommandLineParser::StringContainerType anchor_mask_files_folders;
   if (parsedArgs.count("anchor_masks"))
-    anchor_mask_files = us::any_cast<mitkCommandLineParser::StringContainerType>(parsedArgs["anchor_masks"]);
+    anchor_mask_files_folders = us::any_cast<mitkCommandLineParser::StringContainerType>(parsedArgs["anchor_masks"]);
 
   std::string regu = "NONE";
   if (parsedArgs.count("regu"))
     regu = us::any_cast<std::string>(parsedArgs["regu"]);
 
+  bool use_weights = false;
+  if (parsedArgs.count("use_weights"))
+    use_weights = us::any_cast<bool>(parsedArgs["use_weights"]);
+
+  bool use_num_streamlines = false;
+  if (parsedArgs.count("use_num_streamlines"))
+    use_num_streamlines = us::any_cast<bool>(parsedArgs["use_num_streamlines"]);
+
   try
   {
     itk::TimeProbe clock;
     clock.Start();
 
     if (!ist::PathExists(out_folder))
     {
       MITK_INFO << "Creating output directory";
       ist::MakeDirectory(out_folder);
     }
 
     MITK_INFO << "Loading data";
     std::streambuf *old = cout.rdbuf(); // <-- save
     std::stringstream ss;
     std::cout.rdbuf (ss.rdbuf());       // <-- redirect
 
     ofstream logfile;
     logfile.open (out_folder + "log.txt");
 
     itk::ImageFileWriter< PeakImgType >::Pointer peak_image_writer = itk::ImageFileWriter< PeakImgType >::New();
     mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Peak Image", "Fiberbundles"}, {});
-    auto inputImage = mitk::IOUtil::Load<mitk::PeakImage>(peak_file_name, &functor);
+    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 = mitk::IOUtil::Load<mitk::Image>(mask_file);
       mitk::CastToItkImage(mitk_mask, mask);
     }
 
     // Load masks covering the true positives for evaluation purposes
     std::vector< itk::FitFibersToImageFilter::UcharImgType::Pointer > reference_masks;
-    for (auto filename : anchor_mask_files)
+    std::vector< std::string > anchor_mask_files;
+    for (auto filename : anchor_mask_files_folders)
     {
-      itk::FitFibersToImageFilter::UcharImgType::Pointer ref_mask = nullptr;
-      mitk::Image::Pointer ref_mitk_mask = mitk::IOUtil::Load<mitk::Image>(filename);
-      mitk::CastToItkImage(ref_mitk_mask, ref_mask);
-      reference_masks.push_back(ref_mask);
+      if (itksys::SystemTools::PathExists(filename))
+      {
+        auto list = get_file_list(filename, {".nrrd",".nii.gz",".nii"});
+        for (auto f : list)
+        {
+          itk::FitFibersToImageFilter::UcharImgType::Pointer ref_mask = nullptr;
+          mitk::Image::Pointer ref_mitk_mask = mitk::IOUtil::Load<mitk::Image>(f);
+          mitk::CastToItkImage(ref_mitk_mask, ref_mask);
+          reference_masks.push_back(ref_mask);
+          anchor_mask_files.push_back(f);
+        }
+      }
+      else if (itksys::SystemTools::FileExists(filename))
+      {
+        anchor_mask_files.push_back(filename);
+        itk::FitFibersToImageFilter::UcharImgType::Pointer ref_mask = nullptr;
+        mitk::Image::Pointer ref_mitk_mask = mitk::IOUtil::Load<mitk::Image>(filename);
+        mitk::CastToItkImage(ref_mitk_mask, ref_mask);
+        reference_masks.push_back(ref_mask);
+      }
     }
 
     // Load peak image
     typedef mitk::ImageToItk< PeakImgType > CasterType;
     CasterType::Pointer caster = CasterType::New();
     caster->SetInput(inputImage);
     caster->Update();
     PeakImgType::Pointer peak_image = caster->GetOutput();
 
     // Load all candidate tracts
     std::vector< std::string > candidate_tract_files = get_file_list(candidate_tract_folder);
     std::vector< mitk::FiberBundle::Pointer > input_candidates;
     for (std::string f : candidate_tract_files)
     {
       mitk::FiberBundle::Pointer fib = mitk::IOUtil::Load<mitk::FiberBundle>(f);
       if (fib.IsNull())
         continue;
       if (fib->GetNumFibers()<=0)
         continue;
       fib->ResampleLinear(minSpacing/10.0);
       input_candidates.push_back(fib);
     }
     std::cout.rdbuf (old);              // <-- restore
     MITK_INFO << "Loaded " << candidate_tract_files.size() << " candidate tracts.";
 
     double rmse = 0.0;
     int iteration = 0;
     std::string name = "NOANCHOR";
     // Load reference tractogram consisting of all known tracts
     std::vector< mitk::FiberBundle::Pointer > input_reference;
     mitk::FiberBundle::Pointer anchor_tractogram = mitk::IOUtil::Load<mitk::FiberBundle>(anchors_file);
     if ( !(anchor_tractogram.IsNull() || anchor_tractogram->GetNumFibers()==0) )
     {
       std::streambuf *old = cout.rdbuf(); // <-- save
       std::stringstream ss;
       std::cout.rdbuf (ss.rdbuf());       // <-- redirect
       anchor_tractogram->ResampleLinear(minSpacing/10.0);
       std::cout.rdbuf (old);              // <-- restore
       input_reference.push_back(anchor_tractogram);
 
       // Fit known tracts to peak image to obtain underexplained image
       MITK_INFO << "Fit anchor tracts";
       itk::FitFibersToImageFilter::Pointer fitter = itk::FitFibersToImageFilter::New();
       fitter->SetTractograms(input_reference);
       fitter->SetLambda(lambda);
       fitter->SetFilterOutliers(filter_outliers);
       fitter->SetPeakImage(peak_image);
       fitter->SetVerbose(true);
       fitter->SetResampleFibers(false);
       fitter->SetMaskImage(mask);
       fitter->SetRegularization(VnlCostFunction::REGU::NONE);
       fitter->Update();
       rmse = fitter->GetRMSE();
       vnl_vector<double> rms_diff = fitter->GetRmsDiffPerBundle();
       logfile << "RMS_DIFF: " << setprecision(5) << rms_diff[0] << " " << name << " RMSE: " << rmse << "\n";
 
       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 + boost::lexical_cast<std::string>((int)(100000*rms_diff[0])) + "_" + name + ".fib");
 
       peak_image = fitter->GetUnderexplainedImage();
       peak_image_writer->SetInput(peak_image);
       peak_image_writer->SetFileName(out_folder + "Residual_" + name + ".nii.gz");
       peak_image_writer->Update();
     }
 
-    if (!greedy_add)
+    if (use_weights || use_num_streamlines)
+    {
+      MITK_INFO << "Using tract weights as scores";
+      int c = 0;
+      for (auto fib : input_candidates)
+      {
+        int mod = 1;
+        double score = 0;
+        if (use_weights)
+        {
+          score = fib->GetFiberWeight(0);
+          mod = 100000;
+        }
+        else if (use_num_streamlines)
+          score = fib->GetNumFibers();
+        fib->ColorFibersByOrientation();
+
+        std::string bundle_name = ist::GetFilenameWithoutExtension(candidate_tract_files.at(c));
+
+        std::streambuf *old = cout.rdbuf(); // <-- save
+        std::stringstream ss;
+        std::cout.rdbuf (ss.rdbuf());       // <-- redirect
+        mitk::IOUtil::Save(fib, out_folder + boost::lexical_cast<std::string>((int)(mod*score)) + "_" + 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();
+        }
+
+        double weight_sum = 0;
+        for (int i=0; i<fib->GetNumFibers(); i++)
+          weight_sum += fib->GetFiberWeight(i);
+
+        std::cout.rdbuf (old);              // <-- restore
+
+        logfile << "RMS_DIFF: " << setprecision(5) << score << " " << bundle_name << " " << num_voxels << " " << fib->GetNumFibers() << " " << weight_sum << "\n";
+        if (best_overlap_index>=0)
+          logfile << "Best_overlap: " << setprecision(5) << best_overlap << " " << ist::GetFilenameWithoutExtension(anchor_mask_files.at(best_overlap_index)) << "\n";
+        else
+          logfile << "No_overlap\n";
+        ++c;
+      }
+
+    }
+    else 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->SetFitIndividualFibers(true);
 
       if (regu=="MSM")
         fitter->SetRegularization(VnlCostFunction::REGU::MSM);
       else if (regu=="Variance")
         fitter->SetRegularization(VnlCostFunction::REGU::VARIANCE);
       else if (regu=="Lasso")
         fitter->SetRegularization(VnlCostFunction::REGU::LASSO);
       else if (regu=="VoxelVariance")
         fitter->SetRegularization(VnlCostFunction::REGU::VOXEL_VARIANCE);
       else if (regu=="GroupLasso")
         fitter->SetRegularization(VnlCostFunction::REGU::GROUP_LASSO);
       else if (regu=="GroupVariance")
         fitter->SetRegularization(VnlCostFunction::REGU::GROUP_VARIANCE);
       else if (regu=="NONE")
         fitter->SetRegularization(VnlCostFunction::REGU::NONE);
 
       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();
+//      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->ColorFibersByOrientation();
+//        fib->SetFiberWeights( log_rms_diff[c] );
+//        fib->ColorFibersByOrientation();
 
         std::string bundle_name = ist::GetFilenameWithoutExtension(candidate_tract_files.at(c));
 
         std::streambuf *old = cout.rdbuf(); // <-- save
         std::stringstream ss;
         std::cout.rdbuf (ss.rdbuf());       // <-- redirect
+        if (filter_zero_weights)
+          fib = fib->FilterByWeights(0);
         mitk::IOUtil::Save(fib, out_folder + boost::lexical_cast<std::string>((int)(100000*rms_diff[c])) + "_" + bundle_name + ".fib");
 
         float best_overlap = 0;
         int best_overlap_index = -1;
         int m_idx = 0;
         for (auto ref_mask : reference_masks)
         {
           float overlap = fib->GetOverlap(ref_mask, false);
           if (overlap>best_overlap)
           {
             best_overlap = overlap;
             best_overlap_index = m_idx;
           }
           ++m_idx;
         }
 
         unsigned int num_voxels = 0;
         {
           itk::TractDensityImageFilter< ItkUcharImageType >::Pointer masks_filter = itk::TractDensityImageFilter< ItkUcharImageType >::New();
           masks_filter->SetInputImage(mask);
           masks_filter->SetBinaryOutput(true);
           masks_filter->SetFiberBundle(fib);
           masks_filter->SetUseImageGeometry(true);
           masks_filter->Update();
           num_voxels = masks_filter->GetNumCoveredVoxels();
         }
 
         double weight_sum = 0;
         for (int i=0; i<fib->GetNumFibers(); i++)
           weight_sum += fib->GetFiberWeight(i);
 
         std::cout.rdbuf (old);              // <-- restore
 
         logfile << "RMS_DIFF: " << setprecision(5) << rms_diff[c] << " " << bundle_name << " " << num_voxels << " " << fib->GetNumFibers() << " " << weight_sum << "\n";
         if (best_overlap_index>=0)
           logfile << "Best_overlap: " << setprecision(5) << best_overlap << " " << ist::GetFilenameWithoutExtension(anchor_mask_files.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");
 
       peak_image = fitter->GetUnderexplainedImage();
       peak_image_writer->SetInput(peak_image);
       peak_image_writer->SetFileName(out_folder + "Residual_AllCandidates.nii.gz");
       peak_image_writer->Update();
     }
     else
     {
       MITK_INFO << "RMSE: " << setprecision(5) << rmse;
       //    fitter->SetPeakImage(peak_image);
 
       // Iteratively add candidate bundles in a greedy manner
       while (!input_candidates.empty())
       {
         double next_rmse = rmse;
         double num_peaks = 0;
         mitk::FiberBundle::Pointer best_candidate = nullptr;
         PeakImgType::Pointer best_candidate_peak_image = nullptr;
 
         for (int i=0; i<(int)input_candidates.size(); ++i)
         {
           // WHY NECESSARY AGAIN??
           itk::FitFibersToImageFilter::Pointer fitter = itk::FitFibersToImageFilter::New();
           fitter->SetLambda(lambda);
           fitter->SetFilterOutliers(filter_outliers);
           fitter->SetVerbose(false);
           fitter->SetPeakImage(peak_image);
           fitter->SetResampleFibers(false);
           fitter->SetMaskImage(mask);
           // ******************************
           fitter->SetTractograms({input_candidates.at(i)});
 
           std::streambuf *old = cout.rdbuf(); // <-- save
           std::stringstream ss;
           std::cout.rdbuf (ss.rdbuf());       // <-- redirect
           fitter->Update();
           std::cout.rdbuf (old);              // <-- restore
 
           double candidate_rmse = fitter->GetRMSE();
 
           if (candidate_rmse<next_rmse)
           {
             next_rmse = candidate_rmse;
             num_peaks = fitter->GetNumCoveredDirections();
             best_candidate = fitter->GetTractograms().at(0);
             best_candidate_peak_image = fitter->GetUnderexplainedImage();
           }
         }
 
         if (best_candidate.IsNull())
           break;
 
         //      fitter->SetPeakImage(peak_image);
         peak_image = best_candidate_peak_image;
 
         int i=0;
         std::vector< mitk::FiberBundle::Pointer > remaining_candidates;
         std::vector< std::string > remaining_candidate_files;
         for (auto fib : input_candidates)
         {
           if (fib!=best_candidate)
           {
             remaining_candidates.push_back(fib);
             remaining_candidate_files.push_back(candidate_tract_files.at(i));
           }
           else
             name = ist::GetFilenameWithoutExtension(candidate_tract_files.at(i));
           ++i;
         }
         input_candidates = remaining_candidates;
         candidate_tract_files = remaining_candidate_files;
 
         iteration++;
         std::streambuf *old = cout.rdbuf(); // <-- save
         std::stringstream ss;
         std::cout.rdbuf (ss.rdbuf());       // <-- redirect
         // Save winning candidate
+        if (filter_zero_weights)
+          best_candidate = best_candidate->FilterByWeights(0);
         mitk::IOUtil::Save(best_candidate, out_folder + boost::lexical_cast<std::string>(iteration) + "_" + name + ".fib");
         peak_image_writer->SetInput(peak_image);
         peak_image_writer->SetFileName(out_folder + boost::lexical_cast<std::string>(iteration) + "_" + name + ".nrrd");
         peak_image_writer->Update();
 
         // Calculate best overlap with reference masks for evaluation purposes
         float best_overlap = 0;
         int best_overlap_index = -1;
         i = 0;
         for (auto ref_mask : reference_masks)
         {
           float overlap = best_candidate->GetOverlap(ref_mask, false);
           if (overlap>best_overlap)
           {
             best_overlap = overlap;
             best_overlap_index = i;
           }
           ++i;
         }
         std::cout.rdbuf (old);              // <-- restore
 
         logfile << "RMSE: " << setprecision(5) << rmse << " " << name << " " << num_peaks << "\n";
         if (best_overlap_index>=0)
           logfile << "Best_overlap: " << setprecision(5) << best_overlap << " " << ist::GetFilenameWithoutExtension(anchor_mask_files.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/EvaluateLiFE.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/EvaluateLiFE.cpp
index 07ecbbd06a..761fd8daa3 100755
--- a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/EvaluateLiFE.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/EvaluateLiFE.cpp
@@ -1,248 +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 <itkFiberExtractionFilter.h>
 
 typedef itksys::SystemTools ist;
 typedef itk::Image<unsigned char, 3>    ItkFloatImgType;
 typedef std::tuple< ItkFloatImgType::Pointer, std::string > MaskType;
 
 ItkFloatImgType::Pointer LoadItkImage(const std::string& filename)
 {
   mitk::Image::Pointer img = mitk::IOUtil::Load<mitk::Image>(filename);
   ItkFloatImgType::Pointer itkMask = ItkFloatImgType::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);
       std::streambuf *old = cout.rdbuf(); // <-- save
       std::stringstream ss;
       std::cout.rdbuf (ss.rdbuf());       // <-- redirect
       MaskType m(LoadItkImage(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("out", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output text file", 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, "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");
 
   std::map<std::string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
   if (parsedArgs.size()==0)
     return EXIT_FAILURE;
 
   std::string fibFile = us::any_cast<std::string>(parsedArgs["input"]);
   std::string reference_mask_folder = us::any_cast<std::string>(parsedArgs["reference_mask_folder"]);
-  std::string out_folder = us::any_cast<std::string>(parsedArgs["out"]);
+  std::string out_file = us::any_cast<std::string>(parsedArgs["out"]);
 
   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 = mitk::IOUtil::Load<mitk::FiberBundle>(fibFile);
 
     // 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());
 
     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");
+    logfile.open(out_file);
 
     float fpr = 1.0;
     float tpr = 1.0;
     float step = weights.back()/steps;
     float w = 0;
     if (!pre_filter_zeros)
       w -= step;
     while (pred_positives->GetNumFibers()>0 && fpr>0.001 && tpr>0.001)
     {
       w += step;
 
       std::streambuf *old = cout.rdbuf(); // <-- save
       std::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 )
       {
         ItkFloatImgType::Pointer mask_image = std::get<0>(mask);
 
-        mitk::FiberBundle::Pointer a;
+        mitk::FiberBundle::Pointer a=nullptr;
         {
           itk::FiberExtractionFilter<unsigned char>::Pointer extractor = itk::FiberExtractionFilter<unsigned char>::New();
           extractor->SetInputFiberBundle(pred_positives);
           extractor->SetRoiImages({mask_image});
           extractor->SetOverlapFraction(overlap);
           extractor->SetDontResampleFibers(true);
           extractor->SetMode(itk::FiberExtractionFilter<unsigned char>::MODE::OVERLAP);
           extractor->Update();
-          a = extractor->GetPositives().at(0);
+          if (!extractor->GetPositives().empty())
+            a = extractor->GetPositives().at(0);
         }
-        tp_tracts = tp_tracts->AddBundle(a);
+        if (a.IsNotNull())
+          tp_tracts = tp_tracts->AddBundle(a);
 
-        mitk::FiberBundle::Pointer b;
+        mitk::FiberBundle::Pointer b=nullptr;
         {
           itk::FiberExtractionFilter<unsigned char>::Pointer extractor = itk::FiberExtractionFilter<unsigned char>::New();
           extractor->SetInputFiberBundle(pred_negatives);
           extractor->SetRoiImages({mask_image});
           extractor->SetOverlapFraction(overlap);
           extractor->SetDontResampleFibers(true);
           extractor->SetMode(itk::FiberExtractionFilter<unsigned char>::MODE::OVERLAP);
           extractor->Update();
-          b = extractor->GetPositives().at(0);
+          if (!extractor->GetPositives().empty())
+            b = extractor->GetPositives().at(0);
         }
-        fn_tracts = fn_tracts->AddBundle(b);
+        if (b.IsNotNull())
+          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();
       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;
 }