diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.cpp
index 94bfc6ffd5..9b0a9ed25e 100644
--- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.cpp
@@ -1,425 +1,427 @@
 #include "itkFitFibersToImageFilter.h"
 
 #include <boost/progress.hpp>
 
 namespace itk{
 
 FitFibersToImageFilter::FitFibersToImageFilter()
   : m_FitIndividualFibers(true)
   , m_GradientTolerance(1e-5)
   , m_Lambda(0.1)
   , m_MaxIterations(20)
   , m_FiberSampling(10)
   , m_Coverage(0)
   , m_Overshoot(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)
 {
   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];
 
   unsigned int num_unknowns = m_Tractograms.size();
   if (m_FitIndividualFibers)
   {
     num_unknowns = 0;
     for (unsigned int bundle=0; bundle<m_Tractograms.size(); bundle++)
       num_unknowns += m_Tractograms.at(bundle)->GetNumFibers();
   }
 
   for (unsigned int bundle=0; bundle<m_Tractograms.size(); bundle++)
   {
     if (m_DeepCopy)
       m_Tractograms[bundle] = m_Tractograms.at(bundle)->GetDeepCopy();
     m_Tractograms.at(bundle)->ResampleLinear(minSpacing/m_FiberSampling);
   }
 
   unsigned int number_of_residuals = num_voxels * sz_peaks;
 
   MITK_INFO << "Num. unknowns: " << num_unknowns;
   MITK_INFO << "Num. residuals: " << number_of_residuals;
   MITK_INFO << "Creating system ...";
 
   vnl_sparse_matrix<double> A;
   vnl_vector<double> b;
   A.set_size(number_of_residuals, num_unknowns);
   b.set_size(number_of_residuals); b.fill(0.0);
 
   double TD = 0;
   double FD = 0;
   unsigned int dir_count = 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);
         if (!m_PeakImage->GetLargestPossibleRegion().IsInside(idx4))
           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)
         {
           dir_count++;
           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 /= (dir_count*fiber_count);
   FD /= dir_count;
   A /= TD;
   b *= 100.0/FD;  // times 100 because we want to avoid too small values for computational reasons
 
-  double init_lambda = 1e5;  // initialization for lambda estimation
+  double init_lambda = fiber_count;  // initialization for lambda estimation
 
   itk::TimeProbe clock;
   clock.Start();
 
   VnlCostFunction cost(num_unknowns);
   cost.SetProblem(A, b, init_lambda);
-  m_Weights.set_size(num_unknowns); m_Weights.fill( TD/100.0 * FD/2.0 );
+  m_Weights.set_size(num_unknowns); // 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(num_unknowns); l.fill(0);
   vnl_vector<long> bound_selection; bound_selection.set_size(num_unknowns); 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(1);
+  minimizer.set_max_function_evals(2);
   minimizer.minimize(m_Weights);
   vnl_vector<double> dx; dx.set_size(num_unknowns); 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(num_unknowns*0.95);
-    vnl_vector<double> u; u.set_size(num_unknowns); u.fill(weights.at(num_unknowns*0.95));
+    MITK_INFO << "Setting upper weight bound to " << weights.at(num_unknowns*0.99);
+    vnl_vector<double> u; u.set_size(num_unknowns); u.fill(weights.at(num_unknowns*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(num_unknowns*0.5);
   m_MinWeight = weights.at(0);
   m_MaxWeight = weights.at(num_unknowns-1);
 
   MITK_INFO << "*************************";
   MITK_INFO << "Weight statistics";
   MITK_INFO << "Mean: " << m_MeanWeight;
   MITK_INFO << "Median: " << m_MedianWeight;
   MITK_INFO << "75% quantile: " << weights.at(num_unknowns*0.75);
   MITK_INFO << "95% quantile: " << weights.at(num_unknowns*0.95);
   MITK_INFO << "99% quantile: " << weights.at(num_unknowns*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();
   MITK_INFO << "Final RMS: " << cost.S->get_rms_error(m_Weights);
 
   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";
 
   // transform back for peak image creation
   A *= FD/100.0;
   b *= FD/100.0;
 
   MITK_INFO << "Weighting fibers";
   if (m_FitIndividualFibers)
   {
     unsigned int fiber_count = 0;
     for (unsigned int bundle=0; bundle<m_Tractograms.size(); bundle++)
     {
       for (int i=0; i<m_Tractograms.at(bundle)->GetNumFibers(); i++)
       {
         m_Tractograms.at(bundle)->SetFiberWeight(i, m_Weights[fiber_count]);
         ++fiber_count;
       }
       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);
     }
   }
 
   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> idx;
     unsigned int linear_index = r;
     idx[0] = linear_index % sz_x; linear_index /= sz_x;
     idx[1] = linear_index % sz_y; linear_index /= sz_y;
     idx[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;
 
       idx[3] = peak_id*3;
       peak_dir[0] = m_PeakImage->GetPixel(idx);
       idx[3] += 1;
       peak_dir[1] = m_PeakImage->GetPixel(idx);
       idx[3] += 1;
       peak_dir[2] = m_PeakImage->GetPixel(idx);
 
       peak_dir.normalize();
       peak_dir *= fitted_b[r];
 
       idx[3] = peak_id*3;
       m_FittedImage->SetPixel(idx, peak_dir[0]);
 
       idx[3] += 1;
       m_FittedImage->SetPixel(idx, peak_dir[1]);
 
       idx[3] += 1;
       m_FittedImage->SetPixel(idx, peak_dir[2]);
     }
   }
 
   FD = 0;
   m_Coverage = 0;
   m_Overshoot = 0;
 
   itk::Index<4> idx;
   for (idx[0]=0; idx[0]<sz_x; ++idx[0])
     for (idx[1]=0; idx[1]<sz_y; ++idx[1])
       for (idx[2]=0; idx[2]<sz_z; ++idx[2])
       {
         vnl_vector_fixed<float,3> peak_dir;
         vnl_vector_fixed<float,3> fitted_dir;
         vnl_vector_fixed<float,3> overshoot_dir;
         for (idx[3]=0; idx[3]<(itk::IndexValueType)m_PeakImage->GetLargestPossibleRegion().GetSize(3); ++idx[3])
         {
           peak_dir[idx[3]%3] = m_PeakImage->GetPixel(idx);
           fitted_dir[idx[3]%3] = m_FittedImage->GetPixel(idx);
           m_ResidualImage->SetPixel(idx, m_PeakImage->GetPixel(idx) - m_FittedImage->GetPixel(idx));
 
           if (idx[3]%3==2)
           {
             FD += peak_dir.magnitude();
 
             itk::Index<4> tidx= idx;
             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(1) << 100.0*m_Coverage << "%";
   MITK_INFO << std::fixed << "Overshoot: " << setprecision(1) << 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.8;
 
   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/cmdapps/FiberProcessing/FitFibersToImage.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FitFibersToImage.cpp
index a8a0539732..e58629356e 100755
--- a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FitFibersToImage.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FitFibersToImage.cpp
@@ -1,203 +1,203 @@
 /*===================================================================
 
 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>
 
 using namespace std;
 typedef itksys::SystemTools ist;
 typedef itk::Point<float, 4> PointType4;
 typedef itk::Image< float, 4 >  PeakImgType;
 
 std::vector< string > get_file_list(const std::string& path)
 {
   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);
       if (ext==".fib" || ext==".trk")
         file_list.push_back(path + '/' + filename);
     }
   }
   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("Fit Fibers To Image");
   parser.setCategory("Fiber Tracking and Processing Methods");
   parser.setDescription("Assigns a weight to each fiber in order to optimally explain the input peak image");
   parser.setContributor("MIC");
 
   parser.setArgumentPrefix("--", "-");
   parser.addArgument("", "i1", mitkCommandLineParser::StringList, "Input tractograms:", "input tractograms (.fib, vtk ascii file format)", us::Any(), false);
   parser.addArgument("", "i2", mitkCommandLineParser::InputFile, "Input peaks:", "input peak image", us::Any(), false);
   parser.addArgument("", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output root", us::Any(), false);
 
   parser.addArgument("max_iter", "", mitkCommandLineParser::Int, "Max. iterations:", "maximum number of optimizer iterations", 20);
   parser.addArgument("bundle_based", "", mitkCommandLineParser::Bool, "Bundle based fit:", "fit one weight per input tractogram/bundle, not for each fiber", false);
   parser.addArgument("min_g", "", mitkCommandLineParser::Float, "Min. g:", "lower termination threshold for gradient magnitude", 1e-5);
   parser.addArgument("lambda", "", mitkCommandLineParser::Float, "Lambda:", "modifier for regularization", 0.1);
   parser.addArgument("save_res", "", mitkCommandLineParser::Bool, "Residuals:", "save residual images", false);
-  parser.addArgument("filter_outliers", "", mitkCommandLineParser::Bool, "Filter outliers:", "perform second optimization run with an upper weight bound based on the first weight estimation (95% quantile)", true);
+  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;
 
   mitkCommandLineParser::StringContainerType fib_files = us::any_cast<mitkCommandLineParser::StringContainerType>(parsedArgs["i1"]);
   string peak_file_name = us::any_cast<string>(parsedArgs["i2"]);
   string outRoot = us::any_cast<string>(parsedArgs["o"]);
 
   bool single_fib = true;
   if (parsedArgs.count("bundle_based"))
     single_fib = !us::any_cast<bool>(parsedArgs["bundle_based"]);
 
   bool save_residuals = false;
   if (parsedArgs.count("save_res"))
     save_residuals = us::any_cast<bool>(parsedArgs["save_res"]);
 
   int max_iter = 20;
   if (parsedArgs.count("max_iter"))
     max_iter = us::any_cast<int>(parsedArgs["max_iter"]);
 
   float g_tol = 1e-5;
   if (parsedArgs.count("min_g"))
     g_tol = us::any_cast<float>(parsedArgs["min_g"]);
 
   float lambda = 0.1;
   if (parsedArgs.count("lambda"))
     lambda = us::any_cast<float>(parsedArgs["lambda"]);
 
   bool filter_outliers = true;
-  if (parsedArgs.count("filter_outliers"))
-    filter_outliers = us::any_cast<bool>(parsedArgs["filter_outliers"]);
+  if (parsedArgs.count("dont_filter_outliers"))
+    filter_outliers = !us::any_cast<bool>(parsedArgs["dont_filter_outliers"]);
 
   try
   {
     std::vector< mitk::FiberBundle::Pointer > input_tracts;
 
     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());
 
     typedef mitk::ImageToItk< PeakImgType > CasterType;
     CasterType::Pointer caster = CasterType::New();
     caster->SetInput(inputImage);
     caster->Update();
     PeakImgType::Pointer peak_image = caster->GetOutput();
 
     std::vector< std::string > fib_names;
     for (auto item : fib_files)
     {
       if ( ist::FileIsDirectory(item) )
       {
         for ( auto fibFile : get_file_list(item) )
         {
           mitk::FiberBundle::Pointer inputTractogram = dynamic_cast<mitk::FiberBundle*>(mitk::IOUtil::Load(fibFile)[0].GetPointer());
           if (inputTractogram.IsNull())
             continue;
           input_tracts.push_back(inputTractogram);
           fib_names.push_back(fibFile);
         }
       }
       else
       {
         mitk::FiberBundle::Pointer inputTractogram = dynamic_cast<mitk::FiberBundle*>(mitk::IOUtil::Load(item)[0].GetPointer());
         if (inputTractogram.IsNull())
           continue;
         input_tracts.push_back(inputTractogram);
         fib_names.push_back(item);
       }
     }
 
     itk::FitFibersToImageFilter::Pointer fitter = itk::FitFibersToImageFilter::New();
     fitter->SetPeakImage(peak_image);
     fitter->SetTractograms(input_tracts);
     fitter->SetFitIndividualFibers(single_fib);
     fitter->SetMaxIterations(max_iter);
     fitter->SetGradientTolerance(g_tol);
     fitter->SetLambda(lambda);
     fitter->SetFilterOutliers(filter_outliers);
     fitter->Update();
 
     if (save_residuals)
     {
       itk::ImageFileWriter< PeakImgType >::Pointer writer = itk::ImageFileWriter< PeakImgType >::New();
       writer->SetInput(fitter->GetFittedImage());
       writer->SetFileName(outRoot + "fitted_image.nrrd");
       writer->Update();
 
       writer->SetInput(fitter->GetResidualImage());
       writer->SetFileName(outRoot + "residual_image.nrrd");
       writer->Update();
 
       writer->SetInput(fitter->GetOverexplainedImage());
       writer->SetFileName(outRoot + "overexplained_image.nrrd");
       writer->Update();
 
       writer->SetInput(fitter->GetUnderexplainedImage());
       writer->SetFileName(outRoot + "underexplained_image.nrrd");
       writer->Update();
     }
 
     std::vector< mitk::FiberBundle::Pointer > output_tracts = fitter->GetTractograms();
     for (unsigned int bundle=0; bundle<output_tracts.size(); bundle++)
     {
       std::string name = fib_names.at(bundle);
       name = ist::GetFilenameWithoutExtension(name);
       mitk::IOUtil::Save(output_tracts.at(bundle), outRoot + name + "_fitted.fib");
     }
   }
   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/TractPlausibility.cpp
index 2b74fd5a94..899d63b403 100755
--- a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibility.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibility.cpp
@@ -1,187 +1,201 @@
 /*===================================================================
 
 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>
 
 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 + "/plausible_tracts/");
   ist::MakeDirectory(path + "/implausible_tracts/");
 }
 
 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::vector< MaskType > mask_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);
+
+      if (std::rand()%3==0 && ist::GetFilenameWithoutExtension(filename)!="CC")
+      {
+        MITK_INFO << "Dismissing " << ist::GetFilenameWithoutExtension(filename);
+        continue;
+      }
+
       std::string ext = ist::GetFilenameExtension(filename);
       if (ext==".nii" || ext==".nii.gz" || ext==".nrrd")
       {
+        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::GetFilenameWithoutExtension(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.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 root", 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, "Reference Mask Folder:", "reference masks of known bundles", false);
+  parser.addArgument("gray_matter_mask", "gm", mitkCommandLineParser::String, "", "", false);
 
   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 gray_matter_mask = us::any_cast<string>(parsedArgs["gray_matter_mask"]);
   string reference_mask_folder = us::any_cast<string>(parsedArgs["reference_mask_folder"]);
   string outRoot = us::any_cast<string>(parsedArgs["out"]);
 
   try
   {
     CreateFolderStructure(outRoot);
 
     std::vector< MaskType > known_tract_masks = get_file_list(reference_mask_folder);
     if (known_tract_masks.empty())
       return EXIT_FAILURE;
 
     ItkUcharImgType::Pointer gm_image = LoadItkMaskImage(gray_matter_mask);
 
     mitk::FiberBundle::Pointer inputTractogram = dynamic_cast<mitk::FiberBundle*>(mitk::IOUtil::Load(fibFile)[0].GetPointer());
 
     // filter gray matter fibers
     mitk::FiberBundle::Pointer not_gm_fibers = inputTractogram->ExtractFiberSubset(gm_image, false, true, true);
     mitk::IOUtil::Save(not_gm_fibers, outRoot + "/implausible_tracts/no_gm_endings.trk");
     inputTractogram = inputTractogram->ExtractFiberSubset(gm_image, false, false, true);
 
     // 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);
 
     // find known tracts via overlap match
     mitk::FiberBundle::Pointer all_known_tracts = nullptr;
     for ( MaskType mask : known_tract_masks )
     {
       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, 0.9, false);
       mitk::IOUtil::Save(known_tract, outRoot + "/known_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);
     }
     mitk::IOUtil::Save(all_known_tracts, outRoot + "/known_tracts/all_known_tracts.trk");
-    mitk::IOUtil::Save(TransformToMRtrixSpace(all_known_tracts), outRoot + "/known_tracts/all_known_tracts.fib");
+    mitk::IOUtil::Save(TransformToMRtrixSpace(all_known_tracts), outRoot + "/known_tracts/all_known_tracts_mrtrixspace.fib");
 
     mitk::FiberBundle::Pointer remaining_tracts = inputTractogram->SubtractBundle(all_known_tracts);
     mitk::IOUtil::Save(remaining_tracts, outRoot + "/plausible_tracts/remaining_tracts.trk");
-    mitk::IOUtil::Save(TransformToMRtrixSpace(remaining_tracts), outRoot + "/plausible_tracts/remaining_tracts.fib");
+    mitk::IOUtil::Save(TransformToMRtrixSpace(remaining_tracts), outRoot + "/plausible_tracts/remaining_tracts_mrtrixspace.fib");
     MITK_INFO << "step_size: " << minSpacing/5;
   }
   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/TractPlausibilityFit.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibilityFit.cpp
index b1eecb53d6..bceac56af1 100755
--- a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibilityFit.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibilityFit.cpp
@@ -1,271 +1,277 @@
 /*===================================================================
 
 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>
 
 using namespace std;
 typedef itksys::SystemTools ist;
 typedef itk::Point<float, 4> PointType4;
 typedef itk::Image< float, 4 >  PeakImgType;
 
 std::vector< string > get_file_list(const std::string& path)
 {
   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);
       if (ext==".fib" || ext==".trk")
         file_list.push_back(path + '/' + filename);
     }
   }
   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("");
   parser.setCategory("Fiber Tracking Evaluation");
   parser.setDescription("");
   parser.setContributor("MIC");
 
   parser.setArgumentPrefix("--", "-");
   parser.addArgument("", "i1", mitkCommandLineParser::InputFile, "Input tractogram:", "input tractogram (.fib, vtk ascii file format)", us::Any(), false);
   parser.addArgument("", "i2", mitkCommandLineParser::InputFile, "Input peaks:", "input peak image", us::Any(), false);
   parser.addArgument("", "i3", mitkCommandLineParser::InputFile, "", "", us::Any(), false);
   parser.addArgument("min_gain", "", mitkCommandLineParser::Float, "Min. gain:", "process stops if remaining bundles don't contribute enough", 0.05);
 
   parser.addArgument("", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output root", us::Any(), false);
 
   parser.addArgument("max_iter", "", mitkCommandLineParser::Int, "Max. iterations:", "maximum number of optimizer iterations", 20);
   parser.addArgument("bundle_based", "", mitkCommandLineParser::Bool, "Bundle based fit:", "fit one weight per input tractogram/bundle, not for each fiber", false);
   parser.addArgument("min_g", "", mitkCommandLineParser::Float, "Min. gradient:", "lower termination threshold for gradient magnitude", 1e-5);
   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 (95% quantile)", true);
+  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 peak_file_name = us::any_cast<string>(parsedArgs["i2"]);
   string candidate_folder = us::any_cast<string>(parsedArgs["i3"]);
   string outRoot = us::any_cast<string>(parsedArgs["o"]);
 
   bool single_fib = true;
   if (parsedArgs.count("bundle_based"))
     single_fib = !us::any_cast<bool>(parsedArgs["bundle_based"]);
 
   int max_iter = 20;
   if (parsedArgs.count("max_iter"))
     max_iter = us::any_cast<int>(parsedArgs["max_iter"]);
 
   float g_tol = 1e-5;
   if (parsedArgs.count("min_g"))
     g_tol = us::any_cast<float>(parsedArgs["min_g"]);
 
   float min_gain = 0.05;
   if (parsedArgs.count("min_gain"))
     min_gain = us::any_cast<float>(parsedArgs["min_gain"]);
 
   float lambda = 0.1;
   if (parsedArgs.count("lambda"))
     lambda = us::any_cast<float>(parsedArgs["lambda"]);
 
   bool filter_outliers = true;
-  if (parsedArgs.count("filter_outliers"))
-    filter_outliers = us::any_cast<bool>(parsedArgs["filter_outliers"]);
+  if (parsedArgs.count("dont_filter_outliers"))
+    filter_outliers = !us::any_cast<bool>(parsedArgs["dont_filter_outliers"]);
 
   try
   {
+    itk::TimeProbe clock;
+    clock.Start();
 
     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());
 
     typedef mitk::ImageToItk< PeakImgType > CasterType;
     CasterType::Pointer caster = CasterType::New();
     caster->SetInput(inputImage);
     caster->Update();
     PeakImgType::Pointer peak_image = caster->GetOutput();
 
     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;
     input_reference.push_back(fib);
 
     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;
       input_candidates.push_back(fib);
     }
 
     int iteration = 0;
     std::string name = ist::GetFilenameWithoutExtension(fib_file);
 
     itk::FitFibersToImageFilter::Pointer fitter = itk::FitFibersToImageFilter::New();
     fitter->SetTractograms(input_reference);
     fitter->SetFitIndividualFibers(single_fib);
     fitter->SetMaxIterations(max_iter);
     fitter->SetGradientTolerance(g_tol);
     fitter->SetLambda(lambda);
     fitter->SetFilterOutliers(filter_outliers);
     fitter->SetPeakImage(peak_image);
-    fitter->SetVerbose(false);
+    fitter->SetVerbose(true);
     fitter->SetDeepCopy(false);
     fitter->Update();
 
-
-    fitter->GetTractograms().at(0)->SetFiberWeights(fitter->GetCoverage());
-    fitter->GetTractograms().at(0)->ColorFibersByFiberWeights(false, false);
+//    fitter->GetTractograms().at(0)->SetFiberWeights(fitter->GetCoverage());
+//    fitter->GetTractograms().at(0)->ColorFibersByFiberWeights(false, false);
 
     mitk::IOUtil::Save(fitter->GetTractograms().at(0), outRoot + "0_" + name + ".fib");
     peak_image = fitter->GetUnderexplainedImage();
 
     itk::ImageFileWriter< PeakImgType >::Pointer writer = itk::ImageFileWriter< PeakImgType >::New();
     writer->SetInput(peak_image);
     writer->SetFileName(outRoot + boost::lexical_cast<string>(iteration) + "_" + name + ".nrrd");
     writer->Update();
 
     double coverage = fitter->GetCoverage();
     MITK_INFO << "Iteration: " << iteration;
     MITK_INFO << std::fixed << "Coverage: " << setprecision(1) << 100.0*coverage << "%";
 //    fitter->SetPeakImage(peak_image);
 
     while (!input_candidates.empty())
     {
       streambuf *old = cout.rdbuf(); // <-- save
       stringstream ss;
       std::cout.rdbuf (ss.rdbuf());       // <-- redirect
 
       double next_coverage = 0;
       mitk::FiberBundle::Pointer best_candidate = nullptr;
       for (auto fib : input_candidates)
       {
 
         // WHY NECESSARY AGAIN??
         itk::FitFibersToImageFilter::Pointer fitter = itk::FitFibersToImageFilter::New();
         fitter->SetFitIndividualFibers(single_fib);
         fitter->SetMaxIterations(max_iter);
         fitter->SetGradientTolerance(g_tol);
         fitter->SetLambda(lambda);
         fitter->SetFilterOutliers(filter_outliers);
         fitter->SetVerbose(false);
         fitter->SetPeakImage(peak_image);
         fitter->SetDeepCopy(false);
         // ******************************
         fitter->SetTractograms({fib});
         fitter->Update();
 
         double candidate_coverage = fitter->GetCoverage();
 
         if (candidate_coverage>next_coverage)
         {
           next_coverage = candidate_coverage;
           if ((1.0-coverage) * next_coverage >= min_gain)
           {
             best_candidate = fitter->GetTractograms().at(0);
             peak_image = fitter->GetUnderexplainedImage();
           }
         }
       }
 
       if (best_candidate.IsNull())
       {
         std::cout.rdbuf (old);              // <-- restore
         break;
       }
 
 //      fitter->SetPeakImage(peak_image);
-      best_candidate->SetFiberWeights((1.0-coverage) * next_coverage);
-      best_candidate->ColorFibersByFiberWeights(false, false);
+//      best_candidate->SetFiberWeights((1.0-coverage) * next_coverage);
+//      best_candidate->ColorFibersByFiberWeights(false, false);
 
       coverage += (1.0-coverage) * next_coverage;
 
       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++;
       mitk::IOUtil::Save(best_candidate, outRoot + boost::lexical_cast<string>(iteration) + "_" + name + ".fib");
       writer->SetInput(peak_image);
       writer->SetFileName(outRoot + boost::lexical_cast<string>(iteration) + "_" + name + ".nrrd");
       writer->Update();
       std::cout.rdbuf (old);              // <-- restore
 
       MITK_INFO << "Iteration: " << iteration;
       MITK_INFO << std::fixed << "Coverage: " << setprecision(1) << 100.0*coverage << "% (+" << 100*(1.0-coverage) * next_coverage << "%)";
     }
-    MITK_INFO << "DONE";
+
+    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";
   }
   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/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationView.cpp
index 4ff98294ba..885c90dc8f 100644
--- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationView.cpp
+++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationView.cpp
@@ -1,450 +1,453 @@
 /*===================================================================
 
 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.
 
 ===================================================================*/
 
 
 // Blueberry
 #include <berryISelectionService.h>
 #include <berryIWorkbenchWindow.h>
 
 // Qmitk
 #include "QmitkFiberQuantificationView.h"
 
 // Qt
 #include <QMessageBox>
 
 // MITK
 #include <mitkNodePredicateProperty.h>
 #include <mitkImageCast.h>
 #include <mitkPointSet.h>
 #include <mitkPlanarCircle.h>
 #include <mitkPlanarPolygon.h>
 #include <mitkPlanarRectangle.h>
 #include <mitkPlanarFigureInteractor.h>
 #include <mitkImageAccessByItk.h>
 #include <mitkDataNodeObject.h>
 #include <mitkTensorImage.h>
 #include <mitkPeakImage.h>
 
 // ITK
 #include <itkResampleImageFilter.h>
 #include <itkGaussianInterpolateImageFunction.h>
 #include <itkImageRegionIteratorWithIndex.h>
 #include <itkTractsToFiberEndingsImageFilter.h>
 #include <itkTractDensityImageFilter.h>
 #include <itkImageRegion.h>
 #include <itkTractsToRgbaImageFilter.h>
 #include <itkTractsToVectorImageFilter.h>
 
 #include <math.h>
 #include <boost/lexical_cast.hpp>
 
 const std::string QmitkFiberQuantificationView::VIEW_ID = "org.mitk.views.fiberquantification";
 const std::string id_DataManager = "org.mitk.views.datamanager";
 using namespace mitk;
 
 QmitkFiberQuantificationView::QmitkFiberQuantificationView()
   : QmitkAbstractView()
   , m_Controls( 0 )
   , m_UpsamplingFactor(5)
 {
 
 }
 
 // Destructor
 QmitkFiberQuantificationView::~QmitkFiberQuantificationView()
 {
 
 }
 
 void QmitkFiberQuantificationView::CreateQtPartControl( QWidget *parent )
 {
   // build up qt view, unless already done
   if ( !m_Controls )
   {
     // create GUI widgets from the Qt Designer's .ui file
     m_Controls = new Ui::QmitkFiberQuantificationViewControls;
     m_Controls->setupUi( parent );
 
     connect( m_Controls->m_ProcessFiberBundleButton, SIGNAL(clicked()), this, SLOT(ProcessSelectedBundles()) );
     connect( m_Controls->m_ExtractFiberPeaks, SIGNAL(clicked()), this, SLOT(CalculateFiberDirections()) );
   }
 }
 
 void QmitkFiberQuantificationView::SetFocus()
 {
   m_Controls->m_ProcessFiberBundleButton->setFocus();
 }
 
 void QmitkFiberQuantificationView::CalculateFiberDirections()
 {
   typedef itk::Image<unsigned char, 3>                                            ItkUcharImgType;
 
   // load fiber bundle
   mitk::FiberBundle::Pointer inputTractogram = dynamic_cast<mitk::FiberBundle*>(m_SelectedFB.back()->GetData());
 
   itk::TractsToVectorImageFilter<float>::Pointer fOdfFilter = itk::TractsToVectorImageFilter<float>::New();
   if (m_SelectedImage.IsNotNull())
   {
     ItkUcharImgType::Pointer itkMaskImage = ItkUcharImgType::New();
     mitk::CastToItkImage<ItkUcharImgType>(m_SelectedImage, itkMaskImage);
     fOdfFilter->SetMaskImage(itkMaskImage);
   }
 
   // extract directions from fiber bundle
   fOdfFilter->SetFiberBundle(inputTractogram);
   fOdfFilter->SetAngularThreshold(cos(m_Controls->m_AngularThreshold->value()*M_PI/180));
   switch (m_Controls->m_FiberDirNormBox->currentIndex())
   {
   case 0:
     fOdfFilter->SetNormalizationMethod(itk::TractsToVectorImageFilter<float>::NormalizationMethods::GLOBAL_MAX);
     break;
   case 1:
     fOdfFilter->SetNormalizationMethod(itk::TractsToVectorImageFilter<float>::NormalizationMethods::SINGLE_VEC_NORM);
     break;
   case 2:
     fOdfFilter->SetNormalizationMethod(itk::TractsToVectorImageFilter<float>::NormalizationMethods::MAX_VEC_NORM);
     break;
   }
   fOdfFilter->SetUseWorkingCopy(true);
   fOdfFilter->SetSizeThreshold(m_Controls->m_PeakThreshold->value());
   fOdfFilter->SetMaxNumDirections(m_Controls->m_MaxNumDirections->value());
   fOdfFilter->Update();
 
   QString name = m_SelectedFB.back()->GetName().c_str();
 
   if (m_Controls->m_NumDirectionsBox->isChecked())
   {
     mitk::Image::Pointer mitkImage = mitk::Image::New();
     mitkImage->InitializeByItk( fOdfFilter->GetNumDirectionsImage().GetPointer() );
     mitkImage->SetVolume( fOdfFilter->GetNumDirectionsImage()->GetBufferPointer() );
 
     mitk::DataNode::Pointer node = mitk::DataNode::New();
     node->SetData(mitkImage);
     node->SetName((name+"_NUM_DIRECTIONS").toStdString().c_str());
     GetDataStorage()->Add(node, m_SelectedFB.back());
   }
 
   Image::Pointer mitkImage = dynamic_cast<Image*>(PeakImage::New().GetPointer());
   mitk::CastToMitkImage(fOdfFilter->GetDirectionImage(), mitkImage);
   mitkImage->SetVolume(fOdfFilter->GetDirectionImage()->GetBufferPointer());
 
   mitk::DataNode::Pointer node = mitk::DataNode::New();
   node->SetData(mitkImage);
   node->SetName( (name+"_DIRECTIONS").toStdString().c_str());
   GetDataStorage()->Add(node, m_SelectedFB.back());
 }
 
 void QmitkFiberQuantificationView::UpdateGui()
 {
   m_Controls->m_ProcessFiberBundleButton->setEnabled(!m_SelectedFB.empty());
   m_Controls->m_ExtractFiberPeaks->setEnabled(!m_SelectedFB.empty());
 }
 
 void QmitkFiberQuantificationView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList<mitk::DataNode::Pointer>& nodes)
 {
   //reset existing Vectors containing FiberBundles and PlanarFigures from a previous selection
   m_SelectedFB.clear();
   m_SelectedSurfaces.clear();
   m_SelectedImage = nullptr;
 
   for (mitk::DataNode::Pointer node: nodes)
   {
     if ( dynamic_cast<mitk::FiberBundle*>(node->GetData()) )
     {
       m_SelectedFB.push_back(node);
     }
     else if (dynamic_cast<mitk::Image*>(node->GetData()))
       m_SelectedImage = dynamic_cast<mitk::Image*>(node->GetData());
     else if (dynamic_cast<mitk::Surface*>(node->GetData()))
     {
       m_SelectedSurfaces.push_back(dynamic_cast<mitk::Surface*>(node->GetData()));
     }
   }
   UpdateGui();
   GenerateStats();
 }
 
 void QmitkFiberQuantificationView::GenerateStats()
 {
   if ( m_SelectedFB.empty() )
     return;
 
   QString stats("");
 
   for( unsigned int i=0; i<m_SelectedFB.size(); i++ )
   {
     mitk::DataNode::Pointer node = m_SelectedFB[i];
     if (node.IsNotNull() && dynamic_cast<mitk::FiberBundle*>(node->GetData()))
     {
       if (i>0)
         stats += "\n-----------------------------\n";
       stats += QString(node->GetName().c_str()) + "\n";
       mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
       stats += "Number of fibers: "+ QString::number(fib->GetNumFibers()) + "\n";
       stats += "Number of points: "+ QString::number(fib->GetNumberOfPoints()) + "\n";
       stats += "Min. length:         "+ QString::number(fib->GetMinFiberLength(),'f',1) + " mm\n";
       stats += "Max. length:         "+ QString::number(fib->GetMaxFiberLength(),'f',1) + " mm\n";
       stats += "Mean length:         "+ QString::number(fib->GetMeanFiberLength(),'f',1) + " mm\n";
       stats += "Median length:       "+ QString::number(fib->GetMedianFiberLength(),'f',1) + " mm\n";
       stats += "Standard deviation:  "+ QString::number(fib->GetLengthStDev(),'f',1) + " mm\n";
 
       vtkSmartPointer<vtkFloatArray> weights = fib->GetFiberWeights();
+
       if (weights!=nullptr)
       {
-        float weight = -1;
-        int c = 0;
+        std::vector< float > weights2;
         for (int i=0; i<weights->GetSize(); i++)
-          if (!mitk::Equal(weights->GetValue(i),weight,0.0001))
-          {
-            weight = weights->GetValue(i);
-                c++;
-            if (c>1)
-              break;
-          }
-        if (c>1)
-          stats += "Detected fiber weights. Fibers are not weighted uniformly.\n";
-        else
-          stats += "Fibers are weighted equally.\n";
+          weights2.push_back(weights->GetValue(i));
+
+        std::sort(weights2.begin(), weights2.end());
+
+        stats += "\nFiber weight statistics\n";
+        stats += "Min: " + QString::number(weights2.front()) + "\n";
+        stats += "1% quantile: " + QString::number(weights2.at(weights2.size()*0.01)) + "\n";
+        stats += "5% quantile: " + QString::number(weights2.at(weights2.size()*0.05)) + "\n";
+        stats += "25% quantile: " + QString::number(weights2.at(weights2.size()*0.25)) + "\n";
+        stats += "Median: " + QString::number(weights2.at(weights2.size()*0.5)) + "\n";
+        stats += "75% quantile: " + QString::number(weights2.at(weights2.size()*0.75)) + "\n";
+        stats += "95% quantile: " + QString::number(weights2.at(weights2.size()*0.95)) + "\n";
+        stats += "99% quantile: " + QString::number(weights2.at(weights2.size()*0.99)) + "\n";
+        stats += "Max: " + QString::number(weights2.back()) + "\n";
       }
       else
         stats += "No fiber weight array found.\n";
     }
   }
   this->m_Controls->m_StatsTextEdit->setText(stats);
 }
 
 void QmitkFiberQuantificationView::ProcessSelectedBundles()
 {
   if ( m_SelectedFB.empty() ){
     QMessageBox::information( nullptr, "Warning", "No fibe bundle selected!");
     MITK_WARN("QmitkFiberQuantificationView") << "no fibe bundle selected";
     return;
   }
 
   int generationMethod = m_Controls->m_GenerationBox->currentIndex();
 
   for( unsigned int i=0; i<m_SelectedFB.size(); i++ )
   {
     mitk::DataNode::Pointer node = m_SelectedFB[i];
     if (node.IsNotNull() && dynamic_cast<mitk::FiberBundle*>(node->GetData()))
     {
       mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
       QString name(node->GetName().c_str());
       DataNode::Pointer newNode = nullptr;
       switch(generationMethod){
       case 0:
         newNode = GenerateTractDensityImage(fib, false, true);
         name += "_TDI";
         break;
       case 1:
         newNode = GenerateTractDensityImage(fib, false, false);
         name += "_TDI";
         break;
       case 2:
         newNode = GenerateTractDensityImage(fib, true, false);
         name += "_envelope";
         break;
       case 3:
         newNode = GenerateColorHeatmap(fib);
         break;
       case 4:
         newNode = GenerateFiberEndingsImage(fib);
         name += "_fiber_endings";
         break;
       case 5:
         newNode = GenerateFiberEndingsPointSet(fib);
         name += "_fiber_endings";
         break;
       }
       if (newNode.IsNotNull())
       {
         newNode->SetName(name.toStdString());
         GetDataStorage()->Add(newNode);
       }
     }
   }
 }
 
 // generate pointset displaying the fiber endings
 mitk::DataNode::Pointer QmitkFiberQuantificationView::GenerateFiberEndingsPointSet(mitk::FiberBundle::Pointer fib)
 {
   mitk::PointSet::Pointer pointSet = mitk::PointSet::New();
   vtkSmartPointer<vtkPolyData> fiberPolyData = fib->GetFiberPolyData();
 
   int count = 0;
   int numFibers = fib->GetNumFibers();
   for( int i=0; i<numFibers; i++ )
   {
     vtkCell* cell = fiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     if (numPoints>0)
     {
       double* point = points->GetPoint(0);
       itk::Point<float,3> itkPoint;
       itkPoint[0] = point[0];
       itkPoint[1] = point[1];
       itkPoint[2] = point[2];
       pointSet->InsertPoint(count, itkPoint);
       count++;
     }
     if (numPoints>2)
     {
       double* point = points->GetPoint(numPoints-1);
       itk::Point<float,3> itkPoint;
       itkPoint[0] = point[0];
       itkPoint[1] = point[1];
       itkPoint[2] = point[2];
       pointSet->InsertPoint(count, itkPoint);
       count++;
     }
   }
 
   mitk::DataNode::Pointer node = mitk::DataNode::New();
   node->SetData( pointSet );
   return node;
 }
 
 // generate image displaying the fiber endings
 mitk::DataNode::Pointer QmitkFiberQuantificationView::GenerateFiberEndingsImage(mitk::FiberBundle::Pointer fib)
 {
   typedef unsigned int OutPixType;
   typedef itk::Image<OutPixType,3> OutImageType;
 
   typedef itk::TractsToFiberEndingsImageFilter< OutImageType > ImageGeneratorType;
   ImageGeneratorType::Pointer generator = ImageGeneratorType::New();
   generator->SetFiberBundle(fib);
   generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value());
   if (m_SelectedImage.IsNotNull())
   {
     OutImageType::Pointer itkImage = OutImageType::New();
     CastToItkImage(m_SelectedImage, itkImage);
     generator->SetInputImage(itkImage);
     generator->SetUseImageGeometry(true);
   }
   generator->Update();
 
   // get output image
   OutImageType::Pointer outImg = generator->GetOutput();
   mitk::Image::Pointer img = mitk::Image::New();
   img->InitializeByItk(outImg.GetPointer());
   img->SetVolume(outImg->GetBufferPointer());
 
   // init data node
   mitk::DataNode::Pointer node = mitk::DataNode::New();
   node->SetData(img);
   return node;
 }
 
 // generate rgba heatmap from fiber bundle
 mitk::DataNode::Pointer QmitkFiberQuantificationView::GenerateColorHeatmap(mitk::FiberBundle::Pointer fib)
 {
   typedef itk::RGBAPixel<unsigned char> OutPixType;
   typedef itk::Image<OutPixType, 3> OutImageType;
   typedef itk::TractsToRgbaImageFilter< OutImageType > ImageGeneratorType;
   ImageGeneratorType::Pointer generator = ImageGeneratorType::New();
   generator->SetFiberBundle(fib);
   generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value());
   if (m_SelectedImage.IsNotNull())
   {
     itk::Image<unsigned char, 3>::Pointer itkImage = itk::Image<unsigned char, 3>::New();
     CastToItkImage(m_SelectedImage, itkImage);
     generator->SetInputImage(itkImage);
     generator->SetUseImageGeometry(true);
   }
   generator->Update();
 
   // get output image
   typedef itk::Image<OutPixType,3> OutType;
   OutType::Pointer outImg = generator->GetOutput();
   mitk::Image::Pointer img = mitk::Image::New();
   img->InitializeByItk(outImg.GetPointer());
   img->SetVolume(outImg->GetBufferPointer());
 
   // init data node
   mitk::DataNode::Pointer node = mitk::DataNode::New();
   node->SetData(img);
   return node;
 }
 
 // generate tract density image from fiber bundle
 mitk::DataNode::Pointer QmitkFiberQuantificationView::GenerateTractDensityImage(mitk::FiberBundle::Pointer fib, bool binary, bool absolute)
 {
   mitk::DataNode::Pointer node = mitk::DataNode::New();
   if (binary)
   {
     typedef unsigned char OutPixType;
     typedef itk::Image<OutPixType, 3> OutImageType;
 
     itk::TractDensityImageFilter< OutImageType >::Pointer generator = itk::TractDensityImageFilter< OutImageType >::New();
     generator->SetFiberBundle(fib);
     generator->SetBinaryOutput(binary);
     generator->SetOutputAbsoluteValues(absolute);
     generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value());
     if (m_SelectedImage.IsNotNull())
     {
       OutImageType::Pointer itkImage = OutImageType::New();
       CastToItkImage(m_SelectedImage, itkImage);
       generator->SetInputImage(itkImage);
       generator->SetUseImageGeometry(true);
 
     }
     generator->Update();
 
     // get output image
     typedef itk::Image<OutPixType,3> OutType;
     OutType::Pointer outImg = generator->GetOutput();
     mitk::Image::Pointer img = mitk::Image::New();
     img->InitializeByItk(outImg.GetPointer());
     img->SetVolume(outImg->GetBufferPointer());
 
     // init data node
     node->SetData(img);
   }
   else
   {
     typedef float OutPixType;
     typedef itk::Image<OutPixType, 3> OutImageType;
 
     itk::TractDensityImageFilter< OutImageType >::Pointer generator = itk::TractDensityImageFilter< OutImageType >::New();
     generator->SetFiberBundle(fib);
     generator->SetBinaryOutput(binary);
     generator->SetOutputAbsoluteValues(absolute);
     generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value());
     if (m_SelectedImage.IsNotNull())
     {
       OutImageType::Pointer itkImage = OutImageType::New();
       CastToItkImage(m_SelectedImage, itkImage);
       generator->SetInputImage(itkImage);
       generator->SetUseImageGeometry(true);
 
     }
     //generator->SetDoFiberResampling(false);
     generator->Update();
 
     // get output image
     typedef itk::Image<OutPixType,3> OutType;
     OutType::Pointer outImg = generator->GetOutput();
     mitk::Image::Pointer img = mitk::Image::New();
     img->InitializeByItk(outImg.GetPointer());
     img->SetVolume(outImg->GetBufferPointer());
 
     // init data node
     node->SetData(img);
   }
   return node;
 }