diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFiberCurvatureFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFiberCurvatureFilter.h
index a6b0c5cfa6..3101de522d 100644
--- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFiberCurvatureFilter.h
+++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFiberCurvatureFilter.h
@@ -1,83 +1,83 @@
 /*===================================================================
 
 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.
 
 ===================================================================*/
 #ifndef itkFiberCurvatureFilter_h
 #define itkFiberCurvatureFilter_h
 
 // MITK
 #include <mitkPlanarEllipse.h>
 #include <mitkFiberBundle.h>
 #include <mitkFiberfoxParameters.h>
 
 // ITK
 #include <itkProcessObject.h>
 
 // VTK
 #include <vtkSmartPointer.h>
 #include <vtkPolyData.h>
 #include <vtkCellArray.h>
 #include <vtkPoints.h>
 #include <vtkPolyLine.h>
 
 using namespace std;
 
 namespace itk{
 
 /**
-* \brief Generates artificial fibers distributed in and interpolated between the input planar figures.   */
+* \brief    */
 
 class FiberCurvatureFilter : public ProcessObject
 {
 public:
 
     typedef FiberCurvatureFilter Self;
     typedef ProcessObject                                       Superclass;
     typedef SmartPointer< Self >                                Pointer;
     typedef SmartPointer< const Self >                          ConstPointer;
 
     itkFactorylessNewMacro(Self)
     itkCloneMacro(Self)
     itkTypeMacro( FiberCurvatureFilter, ProcessObject )
 
     virtual void Update() override{
         this->GenerateData();
     }
 
     itkSetMacro( Distance, double )
     itkSetMacro( AngularDeviation, double )
     itkSetMacro( RemoveFibers, bool )
     itkSetMacro( InputFiberBundle, FiberBundle::Pointer )
     itkGetMacro( OutputFiberBundle, FiberBundle::Pointer )
 
 protected:
 
     void GenerateData() override;
 
     FiberCurvatureFilter();
     virtual ~FiberCurvatureFilter();
 
     FiberBundle::Pointer m_InputFiberBundle;
     FiberBundle::Pointer m_OutputFiberBundle;
     double  m_AngularDeviation;
     double  m_Distance;
     bool    m_RemoveFibers;
 };
 }
 
 #ifndef ITK_MANUAL_INSTANTIATION
 #include "itkFiberCurvatureFilter.cpp"
 #endif
 
 #endif
diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp
index ef74a08581..d2e8e5b73a 100755
--- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp
@@ -1,2584 +1,2584 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 
 
 #define _USE_MATH_DEFINES
 #include "mitkFiberBundle.h"
 
 #include <mitkPlanarCircle.h>
 #include <mitkPlanarPolygon.h>
 #include <mitkPlanarFigureComposite.h>
 #include "mitkImagePixelReadAccessor.h"
 #include <mitkPixelTypeMultiplex.h>
 
 #include <vtkPointData.h>
 #include <vtkDataArray.h>
 #include <vtkUnsignedCharArray.h>
 #include <vtkPolyLine.h>
 #include <vtkCellArray.h>
 #include <vtkCellData.h>
 #include <vtkIdFilter.h>
 #include <vtkClipPolyData.h>
 #include <vtkPlane.h>
 #include <vtkDoubleArray.h>
 #include <vtkKochanekSpline.h>
 #include <vtkParametricFunctionSource.h>
 #include <vtkParametricSpline.h>
 #include <vtkPolygon.h>
 #include <vtkCleanPolyData.h>
 #include <cmath>
 #include <boost/progress.hpp>
 #include <vtkTransformPolyDataFilter.h>
 #include <mitkTransferFunction.h>
 #include <vtkLookupTable.h>
 #include <mitkLookupTable.h>
 #include <vtkCardinalSpline.h>
 #include <vtkAppendPolyData.h>
 
 const char* mitk::FiberBundle::FIBER_ID_ARRAY = "Fiber_IDs";
 
 using namespace std;
 
 mitk::FiberBundle::FiberBundle( vtkPolyData* fiberPolyData )
   : m_NumFibers(0)
 {
   m_FiberWeights = vtkSmartPointer<vtkFloatArray>::New();
   m_FiberWeights->SetName("FIBER_WEIGHTS");
 
   m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
   if (fiberPolyData != nullptr)
     m_FiberPolyData = fiberPolyData;
   else
   {
     this->m_FiberPolyData->SetPoints(vtkSmartPointer<vtkPoints>::New());
     this->m_FiberPolyData->SetLines(vtkSmartPointer<vtkCellArray>::New());
   }
 
   this->UpdateFiberGeometry();
   this->GenerateFiberIds();
   this->ColorFibersByOrientation();
 }
 
 mitk::FiberBundle::~FiberBundle()
 {
 
 }
 
 mitk::FiberBundle::Pointer mitk::FiberBundle::GetDeepCopy()
 {
   mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(m_FiberPolyData);
   newFib->SetFiberColors(this->m_FiberColors);
   newFib->SetFiberWeights(this->m_FiberWeights);
   return newFib;
 }
 
 vtkSmartPointer<vtkPolyData> mitk::FiberBundle::GeneratePolyDataByIds(std::vector<long> fiberIds, vtkSmartPointer<vtkFloatArray> weights)
 {
   vtkSmartPointer<vtkPolyData> newFiberPolyData = vtkSmartPointer<vtkPolyData>::New();
   vtkSmartPointer<vtkCellArray> newLineSet = vtkSmartPointer<vtkCellArray>::New();
   vtkSmartPointer<vtkPoints> newPointSet = vtkSmartPointer<vtkPoints>::New();
   weights->SetNumberOfValues(fiberIds.size());
 
   int counter = 0;
   auto finIt = fiberIds.begin();
   while ( finIt != fiberIds.end() )
   {
     if (*finIt < 0 || *finIt>GetNumFibers()){
       MITK_INFO << "FiberID can not be negative or >NumFibers!!! check id Extraction!" << *finIt;
       break;
     }
 
     vtkSmartPointer<vtkCell> fiber = m_FiberIdDataSet->GetCell(*finIt);//->DeepCopy(fiber);
     vtkSmartPointer<vtkPoints> fibPoints = fiber->GetPoints();
     vtkSmartPointer<vtkPolyLine> newFiber = vtkSmartPointer<vtkPolyLine>::New();
     newFiber->GetPointIds()->SetNumberOfIds( fibPoints->GetNumberOfPoints() );
 
     for(int i=0; i<fibPoints->GetNumberOfPoints(); i++)
     {
       newFiber->GetPointIds()->SetId(i, newPointSet->GetNumberOfPoints());
       newPointSet->InsertNextPoint(fibPoints->GetPoint(i)[0], fibPoints->GetPoint(i)[1], fibPoints->GetPoint(i)[2]);
     }
 
     weights->InsertValue(counter, this->GetFiberWeight(*finIt));
     newLineSet->InsertNextCell(newFiber);
     ++finIt;
     ++counter;
   }
 
   newFiberPolyData->SetPoints(newPointSet);
   newFiberPolyData->SetLines(newLineSet);
   return newFiberPolyData;
 }
 
 // merge two fiber bundles
 mitk::FiberBundle::Pointer mitk::FiberBundle::AddBundles(std::vector< mitk::FiberBundle::Pointer > fibs)
 {
   vtkSmartPointer<vtkPolyData> vNewPolyData = vtkSmartPointer<vtkPolyData>::New();
   vtkSmartPointer<vtkCellArray> vNewLines = vtkSmartPointer<vtkCellArray>::New();
   vtkSmartPointer<vtkPoints> vNewPoints = vtkSmartPointer<vtkPoints>::New();
 
   // add current fiber bundle
   vtkSmartPointer<vtkFloatArray> weights = vtkSmartPointer<vtkFloatArray>::New();
 
   int num_weights = this->GetNumFibers();
   for (auto fib : fibs)
     num_weights += fib->GetNumFibers();
   weights->SetNumberOfValues(num_weights);
 
   unsigned int counter = 0;
   for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     for (int j=0; j<numPoints; j++)
     {
       double p[3];
       points->GetPoint(j, p);
 
       vtkIdType id = vNewPoints->InsertNextPoint(p);
       container->GetPointIds()->InsertNextId(id);
     }
     weights->InsertValue(counter, this->GetFiberWeight(i));
     vNewLines->InsertNextCell(container);
     counter++;
   }
 
   for (auto fib : fibs)
   {
     // add new fiber bundle
     for (int i=0; i<fib->GetFiberPolyData()->GetNumberOfCells(); i++)
     {
       vtkCell* cell = fib->GetFiberPolyData()->GetCell(i);
       int numPoints = cell->GetNumberOfPoints();
       vtkPoints* points = cell->GetPoints();
 
       vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
       for (int j=0; j<numPoints; j++)
       {
         double p[3];
         points->GetPoint(j, p);
 
         vtkIdType id = vNewPoints->InsertNextPoint(p);
         container->GetPointIds()->InsertNextId(id);
       }
       weights->InsertValue(counter, fib->GetFiberWeight(i));
       vNewLines->InsertNextCell(container);
       counter++;
     }
   }
 
   // initialize PolyData
   vNewPolyData->SetPoints(vNewPoints);
   vNewPolyData->SetLines(vNewLines);
 
   // initialize fiber bundle
   mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(vNewPolyData);
   newFib->SetFiberWeights(weights);
   return newFib;
 }
 
 // merge two fiber bundles
 mitk::FiberBundle::Pointer mitk::FiberBundle::AddBundle(mitk::FiberBundle* fib)
 {
   if (fib==nullptr)
-  {
-    MITK_WARN << "trying to call AddBundle with nullptr argument";
-    return nullptr;
-  }
+    return this->GetDeepCopy();
+
   MITK_INFO << "Adding fibers";
 
   vtkSmartPointer<vtkPolyData> vNewPolyData = vtkSmartPointer<vtkPolyData>::New();
   vtkSmartPointer<vtkCellArray> vNewLines = vtkSmartPointer<vtkCellArray>::New();
   vtkSmartPointer<vtkPoints> vNewPoints = vtkSmartPointer<vtkPoints>::New();
 
   // add current fiber bundle
   vtkSmartPointer<vtkFloatArray> weights = vtkSmartPointer<vtkFloatArray>::New();
   weights->SetNumberOfValues(this->GetNumFibers()+fib->GetNumFibers());
 
   unsigned int counter = 0;
   for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     for (int j=0; j<numPoints; j++)
     {
       double p[3];
       points->GetPoint(j, p);
 
       vtkIdType id = vNewPoints->InsertNextPoint(p);
       container->GetPointIds()->InsertNextId(id);
     }
     weights->InsertValue(counter, this->GetFiberWeight(i));
     vNewLines->InsertNextCell(container);
     counter++;
   }
 
   // add new fiber bundle
   for (int i=0; i<fib->GetFiberPolyData()->GetNumberOfCells(); i++)
   {
     vtkCell* cell = fib->GetFiberPolyData()->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     for (int j=0; j<numPoints; j++)
     {
       double p[3];
       points->GetPoint(j, p);
 
       vtkIdType id = vNewPoints->InsertNextPoint(p);
       container->GetPointIds()->InsertNextId(id);
     }
     weights->InsertValue(counter, fib->GetFiberWeight(i));
     vNewLines->InsertNextCell(container);
     counter++;
   }
 
   // initialize PolyData
   vNewPolyData->SetPoints(vNewPoints);
   vNewPolyData->SetLines(vNewLines);
 
   // initialize fiber bundle
   mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(vNewPolyData);
   newFib->SetFiberWeights(weights);
   return newFib;
 }
 
 // Only retain fibers with a weight larger than the specified threshold
-mitk::FiberBundle::Pointer mitk::FiberBundle::FilterByWeights(float weight_thr)
+mitk::FiberBundle::Pointer mitk::FiberBundle::FilterByWeights(float weight_thr, bool invert)
 {
   vtkSmartPointer<vtkPolyData> vNewPolyData = vtkSmartPointer<vtkPolyData>::New();
   vtkSmartPointer<vtkCellArray> vNewLines = vtkSmartPointer<vtkCellArray>::New();
   vtkSmartPointer<vtkPoints> vNewPoints = vtkSmartPointer<vtkPoints>::New();
 
   std::vector<float> weights;
 
   for (int i=0; i<this->GetNumFibers(); i++)
   {
-    if (this->GetFiberWeight(i)<=weight_thr)
+    if ( (invert && this->GetFiberWeight(i)>weight_thr) || (!invert && this->GetFiberWeight(i)<=weight_thr))
       continue;
 
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     for (int j=0; j<numPoints; j++)
     {
       double p[3];
       points->GetPoint(j, p);
 
       vtkIdType id = vNewPoints->InsertNextPoint(p);
       container->GetPointIds()->InsertNextId(id);
     }
     vNewLines->InsertNextCell(container);
     weights.push_back(this->GetFiberWeight(i));
   }
 
   // initialize PolyData
   vNewPolyData->SetPoints(vNewPoints);
   vNewPolyData->SetLines(vNewLines);
 
   // initialize fiber bundle
   mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(vNewPolyData);
   for (unsigned int i=0; i<weights.size(); ++i)
     newFib->SetFiberWeight(i, weights.at(i));
   return newFib;
 }
 
 // Only retain a subsample of the fibers
 mitk::FiberBundle::Pointer mitk::FiberBundle::SubsampleFibers(float factor)
 {
   vtkSmartPointer<vtkPolyData> vNewPolyData = vtkSmartPointer<vtkPolyData>::New();
   vtkSmartPointer<vtkCellArray> vNewLines = vtkSmartPointer<vtkCellArray>::New();
   vtkSmartPointer<vtkPoints> vNewPoints = vtkSmartPointer<vtkPoints>::New();
 
   int new_num_fibs = this->GetNumFibers()*factor;
   MITK_INFO << "Subsampling fibers with factor " << factor << "(" << new_num_fibs << "/" << this->GetNumFibers() << ")";
 
   // add current fiber bundle
   vtkSmartPointer<vtkFloatArray> weights = vtkSmartPointer<vtkFloatArray>::New();
   weights->SetNumberOfValues(new_num_fibs);
 
   std::vector< int > ids;
   for (int i=0; i<this->GetNumFibers(); i++)
     ids.push_back(i);
   std::random_shuffle(ids.begin(), ids.end());
 
   unsigned int counter = 0;
   for (int i=0; i<new_num_fibs; i++)
   {
     vtkCell* cell = m_FiberPolyData->GetCell(ids.at(i));
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     for (int j=0; j<numPoints; j++)
     {
       double p[3];
       points->GetPoint(j, p);
 
       vtkIdType id = vNewPoints->InsertNextPoint(p);
       container->GetPointIds()->InsertNextId(id);
     }
     weights->InsertValue(counter, this->GetFiberWeight(ids.at(i)));
     vNewLines->InsertNextCell(container);
     counter++;
   }
 
   // initialize PolyData
   vNewPolyData->SetPoints(vNewPoints);
   vNewPolyData->SetLines(vNewLines);
 
   // initialize fiber bundle
   mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(vNewPolyData);
   newFib->SetFiberWeights(weights);
   return newFib;
 }
 
 // subtract two fiber bundles
 mitk::FiberBundle::Pointer mitk::FiberBundle::SubtractBundle(mitk::FiberBundle* fib)
 {
-  MITK_INFO << "Subtracting fibers";
+  if (fib==nullptr)
+    return this->GetDeepCopy();
 
+  MITK_INFO << "Subtracting fibers";
   vtkSmartPointer<vtkPolyData> vNewPolyData = vtkSmartPointer<vtkPolyData>::New();
   vtkSmartPointer<vtkCellArray> vNewLines = vtkSmartPointer<vtkCellArray>::New();
   vtkSmartPointer<vtkPoints> vNewPoints = vtkSmartPointer<vtkPoints>::New();
 
   std::vector< std::vector< itk::Point<float, 3> > > points1;
   for( int i=0; i<m_NumFibers; i++ )
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     if (points==nullptr || numPoints<=0)
       continue;
 
     itk::Point<float, 3> start = GetItkPoint(points->GetPoint(0));
     itk::Point<float, 3> end = GetItkPoint(points->GetPoint(numPoints-1));
 
     points1.push_back( {start, end} );
   }
 
   std::vector< std::vector< itk::Point<float, 3> > > points2;
   for( int i=0; i<fib->GetNumFibers(); i++ )
   {
     vtkCell* cell = fib->GetFiberPolyData()->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     if (points==nullptr || numPoints<=0)
       continue;
 
     itk::Point<float, 3> start = GetItkPoint(points->GetPoint(0));
     itk::Point<float, 3> end = GetItkPoint(points->GetPoint(numPoints-1));
 
     points2.push_back( {start, end} );
   }
 
 //  int progress = 0;
   std::vector< int > ids;
 #pragma omp parallel for
   for (int i=0; i<(int)points1.size(); i++)
   {
 //#pragma omp critical
 //    {
 //      progress++;
 //      std::cout << (int)(100*(float)progress/points1.size()) << "%" << '\r';
 //      cout.flush();
 //    }
 
     bool match = false;
     for (unsigned int j=0; j<points2.size(); j++)
     {
       auto v1 = points1.at(i);
       auto v2 = points2.at(j);
 
       float dist=0;
       for (unsigned int c=0; c<v1.size(); c++)
       {
         float d = v1[c][0]-v2[c][0];
         dist += d*d;
 
         d = v1[c][1]-v2[c][1];
         dist += d*d;
 
         d = v1[c][2]-v2[c][2];
         dist += d*d;
       }
       dist /= v1.size();
 
       if (dist<mitk::eps)
       {
         match = true;
         break;
       }
 
       dist=0;
       for (unsigned int c=0; c<v1.size(); c++)
       {
         float d = v1[v1.size()-1-c][0]-v2[c][0];
         dist += d*d;
 
         d = v1[v1.size()-1-c][1]-v2[c][1];
         dist += d*d;
 
         d = v1[v1.size()-1-c][2]-v2[c][2];
         dist += d*d;
       }
       dist /= v1.size();
 
       if (dist<mitk::eps)
       {
         match = true;
         break;
       }
     }
 
 #pragma omp critical
     if (!match)
       ids.push_back(i);
   }
 
   for( int i : ids )
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     if (points==nullptr || numPoints<=0)
       continue;
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     for( int j=0; j<numPoints; j++)
     {
       vtkIdType id = vNewPoints->InsertNextPoint(points->GetPoint(j));
       container->GetPointIds()->InsertNextId(id);
     }
     vNewLines->InsertNextCell(container);
   }
   if(vNewLines->GetNumberOfCells()==0)
-    return nullptr;
+    return mitk::FiberBundle::New();
   // initialize PolyData
   vNewPolyData->SetPoints(vNewPoints);
   vNewPolyData->SetLines(vNewLines);
 
   // initialize fiber bundle
   return mitk::FiberBundle::New(vNewPolyData);
 }
 
 itk::Point<float, 3> mitk::FiberBundle::GetItkPoint(double point[3])
 {
   itk::Point<float, 3> itkPoint;
   itkPoint[0] = point[0];
   itkPoint[1] = point[1];
   itkPoint[2] = point[2];
   return itkPoint;
 }
 
 /*
  * set PolyData (additional flag to recompute fiber geometry, default = true)
  */
 void mitk::FiberBundle::SetFiberPolyData(vtkSmartPointer<vtkPolyData> fiberPD, bool updateGeometry)
 {
   if (fiberPD == nullptr)
     this->m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
   else
     m_FiberPolyData->DeepCopy(fiberPD);
 
   m_NumFibers = m_FiberPolyData->GetNumberOfLines();
 
   if (updateGeometry)
     UpdateFiberGeometry();
   GenerateFiberIds();
   ColorFibersByOrientation();
 }
 
 /*
  * return vtkPolyData
  */
 vtkSmartPointer<vtkPolyData> mitk::FiberBundle::GetFiberPolyData() const
 {
   return m_FiberPolyData;
 }
 
 void mitk::FiberBundle::ColorFibersByOrientation()
 {
   //===== FOR WRITING A TEST ========================
   //  colorT size == tupelComponents * tupelElements
   //  compare color results
   //  to cover this code 100% also PolyData needed, where colorarray already exists
   //  + one fiber with exactly 1 point
   //  + one fiber with 0 points
   //=================================================
 
   vtkPoints* extrPoints = nullptr;
   extrPoints = m_FiberPolyData->GetPoints();
   int numOfPoints = 0;
   if (extrPoints!=nullptr)
     numOfPoints = extrPoints->GetNumberOfPoints();
 
   //colors and alpha value for each single point, RGBA = 4 components
   unsigned char rgba[4] = {0,0,0,0};
   int componentSize = 4;
   m_FiberColors = vtkSmartPointer<vtkUnsignedCharArray>::New();
   m_FiberColors->Allocate(numOfPoints * componentSize);
   m_FiberColors->SetNumberOfComponents(componentSize);
   m_FiberColors->SetName("FIBER_COLORS");
 
   int numOfFibers = m_FiberPolyData->GetNumberOfLines();
   if (numOfFibers < 1)
     return;
 
   /* extract single fibers of fiberBundle */
   vtkCellArray* fiberList = m_FiberPolyData->GetLines();
   fiberList->InitTraversal();
   for (int fi=0; fi<numOfFibers; ++fi) {
 
     vtkIdType* idList; // contains the point id's of the line
     vtkIdType pointsPerFiber; // number of points for current line
     fiberList->GetNextCell(pointsPerFiber, idList);
 
     /* single fiber checkpoints: is number of points valid */
     if (pointsPerFiber > 1)
     {
       /* operate on points of single fiber */
       for (int i=0; i <pointsPerFiber; ++i)
       {
         /* process all points elastV[0]ept starting and endpoint for calculating color value take current point, previous point and next point */
         if (i<pointsPerFiber-1 && i > 0)
         {
           /* The color value of the current point is influenced by the previous point and next point. */
           vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]);
           vnl_vector_fixed< double, 3 > nextPntvtk(extrPoints->GetPoint(idList[i+1])[0], extrPoints->GetPoint(idList[i+1])[1], extrPoints->GetPoint(idList[i+1])[2]);
           vnl_vector_fixed< double, 3 > prevPntvtk(extrPoints->GetPoint(idList[i-1])[0], extrPoints->GetPoint(idList[i-1])[1], extrPoints->GetPoint(idList[i-1])[2]);
 
           vnl_vector_fixed< double, 3 > diff1;
           diff1 = currentPntvtk - nextPntvtk;
 
           vnl_vector_fixed< double, 3 > diff2;
           diff2 = currentPntvtk - prevPntvtk;
 
           vnl_vector_fixed< double, 3 > diff;
           diff = (diff1 - diff2) / 2.0;
           diff.normalize();
 
           rgba[0] = (unsigned char) (255.0 * std::fabs(diff[0]));
           rgba[1] = (unsigned char) (255.0 * std::fabs(diff[1]));
           rgba[2] = (unsigned char) (255.0 * std::fabs(diff[2]));
           rgba[3] = (unsigned char) (255.0);
         }
         else if (i==0)
         {
           /* First point has no previous point, therefore only diff1 is taken */
 
           vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]);
           vnl_vector_fixed< double, 3 > nextPntvtk(extrPoints->GetPoint(idList[i+1])[0], extrPoints->GetPoint(idList[i+1])[1], extrPoints->GetPoint(idList[i+1])[2]);
 
           vnl_vector_fixed< double, 3 > diff1;
           diff1 = currentPntvtk - nextPntvtk;
           diff1.normalize();
 
           rgba[0] = (unsigned char) (255.0 * std::fabs(diff1[0]));
           rgba[1] = (unsigned char) (255.0 * std::fabs(diff1[1]));
           rgba[2] = (unsigned char) (255.0 * std::fabs(diff1[2]));
           rgba[3] = (unsigned char) (255.0);
         }
         else if (i==pointsPerFiber-1)
         {
           /* Last point has no next point, therefore only diff2 is taken */
           vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]);
           vnl_vector_fixed< double, 3 > prevPntvtk(extrPoints->GetPoint(idList[i-1])[0], extrPoints->GetPoint(idList[i-1])[1], extrPoints->GetPoint(idList[i-1])[2]);
 
           vnl_vector_fixed< double, 3 > diff2;
           diff2 = currentPntvtk - prevPntvtk;
           diff2.normalize();
 
           rgba[0] = (unsigned char) (255.0 * std::fabs(diff2[0]));
           rgba[1] = (unsigned char) (255.0 * std::fabs(diff2[1]));
           rgba[2] = (unsigned char) (255.0 * std::fabs(diff2[2]));
           rgba[3] = (unsigned char) (255.0);
         }
         m_FiberColors->InsertTypedTuple(idList[i], rgba);
       }
     }
     else if (pointsPerFiber == 1)
     {
       /* a single point does not define a fiber (use vertex mechanisms instead */
       continue;
     }
     else
     {
       MITK_DEBUG << "Fiber with 0 points detected... please check your tractography algorithm!" ;
       continue;
     }
   }
   m_UpdateTime3D.Modified();
   m_UpdateTime2D.Modified();
 }
 
 void mitk::FiberBundle::ColorFibersByCurvature(bool, bool normalize)
 {
   double window = 5;
 
   //colors and alpha value for each single point, RGBA = 4 components
   unsigned char rgba[4] = {0,0,0,0};
   int componentSize = 4;
   m_FiberColors = vtkSmartPointer<vtkUnsignedCharArray>::New();
   m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * componentSize);
   m_FiberColors->SetNumberOfComponents(componentSize);
   m_FiberColors->SetName("FIBER_COLORS");
 
   mitk::LookupTable::Pointer mitkLookup = mitk::LookupTable::New();
   vtkSmartPointer<vtkLookupTable> lookupTable = vtkSmartPointer<vtkLookupTable>::New();
   lookupTable->SetTableRange(0.0, 0.8);
   lookupTable->Build();
   mitkLookup->SetVtkLookupTable(lookupTable);
   mitkLookup->SetType(mitk::LookupTable::JET);
 
   vector< double > values;
   double min = 1;
   double max = 0;
   MITK_INFO << "Coloring fibers by curvature";
   boost::progress_display disp(m_FiberPolyData->GetNumberOfCells());
   for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
   {
     ++disp;
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     // calculate curvatures
     for (int j=0; j<numPoints; j++)
     {
       double dist = 0;
       int c = j;
       std::vector< vnl_vector_fixed< float, 3 > > vectors;
       vnl_vector_fixed< float, 3 > meanV; meanV.fill(0.0);
       while(dist<window/2 && c>1)
       {
         double p1[3];
         points->GetPoint(c-1, p1);
         double p2[3];
         points->GetPoint(c, p2);
 
         vnl_vector_fixed< float, 3 > v;
         v[0] = p2[0]-p1[0];
         v[1] = p2[1]-p1[1];
         v[2] = p2[2]-p1[2];
         dist += v.magnitude();
         v.normalize();
         vectors.push_back(v);
         if (c==j)
           meanV += v;
         c--;
       }
       c = j;
       dist = 0;
       while(dist<window/2 && c<numPoints-1)
       {
         double p1[3];
         points->GetPoint(c, p1);
         double p2[3];
         points->GetPoint(c+1, p2);
 
         vnl_vector_fixed< float, 3 > v;
         v[0] = p2[0]-p1[0];
         v[1] = p2[1]-p1[1];
         v[2] = p2[2]-p1[2];
         dist += v.magnitude();
         v.normalize();
         vectors.push_back(v);
         if (c==j)
           meanV += v;
         c++;
       }
       meanV.normalize();
 
       double dev = 0;
       for (unsigned int c=0; c<vectors.size(); c++)
       {
         double angle = dot_product(meanV, vectors.at(c));
         if (angle>1.0)
           angle = 1.0;
         if (angle<-1.0)
           angle = -1.0;
         dev += acos(angle)*180/M_PI;
       }
       if (vectors.size()>0)
         dev /= vectors.size();
 
       dev = 1.0-dev/180.0;
       values.push_back(dev);
       if (dev<min)
         min = dev;
       if (dev>max)
         max = dev;
     }
   }
   unsigned int count = 0;
   for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     for (int j=0; j<numPoints; j++)
     {
       double color[3];
       double dev = values.at(count);
       if (normalize)
         dev = (dev-min)/(max-min);
       else if (dev>1)
         dev = 1;
       lookupTable->GetColor(dev, color);
 
       rgba[0] = (unsigned char) (255.0 * color[0]);
       rgba[1] = (unsigned char) (255.0 * color[1]);
       rgba[2] = (unsigned char) (255.0 * color[2]);
       rgba[3] = (unsigned char) (255.0);
       m_FiberColors->InsertTypedTuple(cell->GetPointId(j), rgba);
       count++;
     }
   }
   m_UpdateTime3D.Modified();
   m_UpdateTime2D.Modified();
 }
 
 void mitk::FiberBundle::SetFiberOpacity(vtkDoubleArray* FAValArray)
 {
   for(long i=0; i<m_FiberColors->GetNumberOfTuples(); i++)
   {
     double faValue = FAValArray->GetValue(i);
     faValue = faValue * 255.0;
     m_FiberColors->SetComponent(i,3, (unsigned char) faValue );
   }
   m_UpdateTime3D.Modified();
   m_UpdateTime2D.Modified();
 }
 
 void mitk::FiberBundle::ResetFiberOpacity()
 {
   for(long i=0; i<m_FiberColors->GetNumberOfTuples(); i++)
     m_FiberColors->SetComponent(i,3, 255.0 );
   m_UpdateTime3D.Modified();
   m_UpdateTime2D.Modified();
 }
 
 void mitk::FiberBundle::ColorFibersByScalarMap(mitk::Image::Pointer FAimage, bool opacity, bool normalize)
 {
   mitkPixelTypeMultiplex3( ColorFibersByScalarMap, FAimage->GetPixelType(), FAimage, opacity, normalize );
   m_UpdateTime3D.Modified();
   m_UpdateTime2D.Modified();
 }
 
 template <typename TPixel>
 void mitk::FiberBundle::ColorFibersByScalarMap(const mitk::PixelType, mitk::Image::Pointer image, bool opacity, bool normalize)
 {
   m_FiberColors = vtkSmartPointer<vtkUnsignedCharArray>::New();
   m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * 4);
   m_FiberColors->SetNumberOfComponents(4);
   m_FiberColors->SetName("FIBER_COLORS");
 
   mitk::ImagePixelReadAccessor<TPixel,3> readimage(image, image->GetVolumeData(0));
 
   unsigned char rgba[4] = {0,0,0,0};
   vtkPoints* pointSet = m_FiberPolyData->GetPoints();
 
   mitk::LookupTable::Pointer mitkLookup = mitk::LookupTable::New();
   vtkSmartPointer<vtkLookupTable> lookupTable = vtkSmartPointer<vtkLookupTable>::New();
   lookupTable->SetTableRange(0.0, 0.8);
   lookupTable->Build();
   mitkLookup->SetVtkLookupTable(lookupTable);
   mitkLookup->SetType(mitk::LookupTable::JET);
 
   double min = 9999999;
   double max = -9999999;
   for(long i=0; i<m_FiberPolyData->GetNumberOfPoints(); ++i)
   {
     Point3D px;
     px[0] = pointSet->GetPoint(i)[0];
     px[1] = pointSet->GetPoint(i)[1];
     px[2] = pointSet->GetPoint(i)[2];
     double pixelValue = readimage.GetPixelByWorldCoordinates(px);
     if (pixelValue>max)
       max = pixelValue;
     if (pixelValue<min)
       min = pixelValue;
   }
 
   for(long i=0; i<m_FiberPolyData->GetNumberOfPoints(); ++i)
   {
     Point3D px;
     px[0] = pointSet->GetPoint(i)[0];
     px[1] = pointSet->GetPoint(i)[1];
     px[2] = pointSet->GetPoint(i)[2];
     double pixelValue = readimage.GetPixelByWorldCoordinates(px);
 
     if (normalize)
       pixelValue = (pixelValue-min)/(max-min);
     else if (pixelValue>1)
       pixelValue = 1;
 
     double color[3];
     lookupTable->GetColor(1-pixelValue, color);
 
     rgba[0] = (unsigned char) (255.0 * color[0]);
     rgba[1] = (unsigned char) (255.0 * color[1]);
     rgba[2] = (unsigned char) (255.0 * color[2]);
     if (opacity)
       rgba[3] = (unsigned char) (255.0 * pixelValue);
     else
       rgba[3] = (unsigned char) (255.0);
     m_FiberColors->InsertTypedTuple(i, rgba);
   }
   m_UpdateTime3D.Modified();
   m_UpdateTime2D.Modified();
 }
 
 
 void mitk::FiberBundle::ColorFibersByFiberWeights(bool opacity, bool normalize)
 {
   m_FiberColors = vtkSmartPointer<vtkUnsignedCharArray>::New();
   m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * 4);
   m_FiberColors->SetNumberOfComponents(4);
   m_FiberColors->SetName("FIBER_COLORS");
 
   mitk::LookupTable::Pointer mitkLookup = mitk::LookupTable::New();
   vtkSmartPointer<vtkLookupTable> lookupTable = vtkSmartPointer<vtkLookupTable>::New();
   lookupTable->SetTableRange(0.0, 0.8);
   lookupTable->Build();
   mitkLookup->SetVtkLookupTable(lookupTable);
   mitkLookup->SetType(mitk::LookupTable::JET);
 
   unsigned char rgba[4] = {0,0,0,0};
   unsigned int counter = 0;
 
   float max = -999999;
   float min = 999999;
   for (int i=0; i<m_NumFibers; i++)
   {
     double weight = this->GetFiberWeight(i);
     if (weight>max)
       max = weight;
     if (weight<min)
       min = weight;
   }
   if (fabs(max-min)<0.00001)
   {
     max = 1;
     min = 0;
   }
 
   for (int i=0; i<m_NumFibers; i++)
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     double weight = this->GetFiberWeight(i);
 
     for (int j=0; j<numPoints; j++)
     {
       float v = weight;
       if (normalize)
         v = (v-min)/(max-min);
       else if (v>1)
         v = 1;
       double color[3];
       lookupTable->GetColor(1-v, color);
 
       rgba[0] = (unsigned char) (255.0 * color[0]);
       rgba[1] = (unsigned char) (255.0 * color[1]);
       rgba[2] = (unsigned char) (255.0 * color[2]);
       if (opacity)
         rgba[3] = (unsigned char) (255.0 * v);
       else
         rgba[3] = (unsigned char) (255.0);
 
       m_FiberColors->InsertTypedTuple(counter, rgba);
       counter++;
     }
   }
 
   m_UpdateTime3D.Modified();
   m_UpdateTime2D.Modified();
 }
 
 void mitk::FiberBundle::SetFiberColors(float r, float g, float b, float alpha)
 {
   m_FiberColors = vtkSmartPointer<vtkUnsignedCharArray>::New();
   m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * 4);
   m_FiberColors->SetNumberOfComponents(4);
   m_FiberColors->SetName("FIBER_COLORS");
 
   unsigned char rgba[4] = {0,0,0,0};
   for(long i=0; i<m_FiberPolyData->GetNumberOfPoints(); ++i)
   {
     rgba[0] = (unsigned char) r;
     rgba[1] = (unsigned char) g;
     rgba[2] = (unsigned char) b;
     rgba[3] = (unsigned char) alpha;
     m_FiberColors->InsertTypedTuple(i, rgba);
   }
   m_UpdateTime3D.Modified();
   m_UpdateTime2D.Modified();
 }
 
 void mitk::FiberBundle::GenerateFiberIds()
 {
   if (m_FiberPolyData == nullptr)
     return;
 
   vtkSmartPointer<vtkIdFilter> idFiberFilter = vtkSmartPointer<vtkIdFilter>::New();
   idFiberFilter->SetInputData(m_FiberPolyData);
   idFiberFilter->CellIdsOn();
   //  idFiberFilter->PointIdsOn(); // point id's are not needed
   idFiberFilter->SetIdsArrayName(FIBER_ID_ARRAY);
   idFiberFilter->FieldDataOn();
   idFiberFilter->Update();
 
   m_FiberIdDataSet = idFiberFilter->GetOutput();
 
 }
 
 float mitk::FiberBundle::GetOverlap(ItkUcharImgType* mask, bool do_resampling)
 {
   vtkSmartPointer<vtkPolyData> PolyData = m_FiberPolyData;
   mitk::FiberBundle::Pointer fibCopy = this;
   if (do_resampling)
   {
     float minSpacing = 1;
     if(mask->GetSpacing()[0]<mask->GetSpacing()[1] && mask->GetSpacing()[0]<mask->GetSpacing()[2])
       minSpacing = mask->GetSpacing()[0];
     else if (mask->GetSpacing()[1] < mask->GetSpacing()[2])
       minSpacing = mask->GetSpacing()[1];
     else
       minSpacing = mask->GetSpacing()[2];
 
     fibCopy = this->GetDeepCopy();
     fibCopy->ResampleLinear(minSpacing/5);
     PolyData = fibCopy->GetFiberPolyData();
   }
 
   MITK_INFO << "Calculating overlap";
   int inside = 0;
   int outside = 0;
   boost::progress_display disp(m_NumFibers);
   for (int i=0; i<m_NumFibers; i++)
   {
     ++disp;
     vtkCell* cell = PolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     for (int j=0; j<numPoints; j++)
     {
       double* p = points->GetPoint(j);
       itk::Point<float, 3> itkP;
       itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2];
       itk::Index<3> idx;
       mask->TransformPhysicalPointToIndex(itkP, idx);
 
       if ( mask->GetLargestPossibleRegion().IsInside(idx) && mask->GetPixel(idx) != 0 )
         inside++;
       else
         outside++;
     }
   }
 
   if (inside+outside==0)
     outside = 1;
   return (float)inside/(inside+outside);
 }
 
 mitk::FiberBundle::Pointer mitk::FiberBundle::ExtractFiberSubset(ItkUcharImgType* mask, bool anyPoint, bool invert, bool bothEnds, float fraction, bool do_resampling)
 {
   vtkSmartPointer<vtkPolyData> PolyData = m_FiberPolyData;
   mitk::FiberBundle::Pointer fibCopy = this;
   if (anyPoint && do_resampling)
   {
     float minSpacing = 1;
     if(mask->GetSpacing()[0]<mask->GetSpacing()[1] && mask->GetSpacing()[0]<mask->GetSpacing()[2])
       minSpacing = mask->GetSpacing()[0];
     else if (mask->GetSpacing()[1] < mask->GetSpacing()[2])
       minSpacing = mask->GetSpacing()[1];
     else
       minSpacing = mask->GetSpacing()[2];
 
     fibCopy = this->GetDeepCopy();
     fibCopy->ResampleLinear(minSpacing/5);
     PolyData = fibCopy->GetFiberPolyData();
   }
   vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
   vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
 
   std::vector< float > new_weights;
 
   MITK_INFO << "Extracting fibers with mask image";
   boost::progress_display disp(m_NumFibers);
   for (int i=0; i<m_NumFibers; i++)
   {
     ++disp;
 
     vtkCell* cell = PolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     vtkCell* cellOriginal = m_FiberPolyData->GetCell(i);
     int numPointsOriginal = cellOriginal->GetNumberOfPoints();
     vtkPoints* pointsOriginal = cellOriginal->GetPoints();
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
 
     if (numPoints>1 && numPointsOriginal)
     {
       if (anyPoint)
       {
         int inside = 0;
         int outside = 0;
         if (!invert)
         {
           for (int j=0; j<numPoints; j++)
           {
             double* p = points->GetPoint(j);
 
             itk::Point<float, 3> itkP;
             itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2];
             itk::Index<3> idx;
             mask->TransformPhysicalPointToIndex(itkP, idx);
 
             if ( mask->GetLargestPossibleRegion().IsInside(idx) && mask->GetPixel(idx) != 0 )
             {
               inside++;
               if (fraction==0)
                 break;
             }
             else
               outside++;
           }
 
           float current_fraction = 0.0;
           if (inside+outside>0)
             current_fraction = (float)inside/(inside+outside);
 
           if (current_fraction>fraction)
           {
             for (int k=0; k<numPoints; k++)
             {
               double* p = points->GetPoint(k);
               vtkIdType id = vtkNewPoints->InsertNextPoint(p);
               container->GetPointIds()->InsertNextId(id);
             }
           }
         }
         else
         {
           bool includeFiber = true;
           for (int j=0; j<numPoints; j++)
           {
             double* p = points->GetPoint(j);
 
             itk::Point<float, 3> itkP;
             itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2];
             itk::Index<3> idx;
             mask->TransformPhysicalPointToIndex(itkP, idx);
 
             if ( mask->GetPixel(idx) != 0 && mask->GetLargestPossibleRegion().IsInside(idx) )
             {
               inside++;
               includeFiber = false;
               break;
             }
             else
               outside++;
           }
           if (includeFiber)
           {
 
             for (int k=0; k<numPoints; k++)
             {
               double* p = points->GetPoint(k);
               vtkIdType id = vtkNewPoints->InsertNextPoint(p);
               container->GetPointIds()->InsertNextId(id);
             }
           }
         }
       }
       else
       {
         double* start = pointsOriginal->GetPoint(0);
         itk::Point<float, 3> itkStart;
         itkStart[0] = start[0]; itkStart[1] = start[1]; itkStart[2] = start[2];
         itk::Index<3> idxStart;
         mask->TransformPhysicalPointToIndex(itkStart, idxStart);
 
         double* end = pointsOriginal->GetPoint(numPointsOriginal-1);
         itk::Point<float, 3> itkEnd;
         itkEnd[0] = end[0]; itkEnd[1] = end[1]; itkEnd[2] = end[2];
         itk::Index<3> idxEnd;
         mask->TransformPhysicalPointToIndex(itkEnd, idxEnd);
 
         if (invert)
         {
           if (bothEnds)
           {
             if ( mask->GetPixel(idxStart) == 0 && mask->GetPixel(idxEnd) == 0 )
             {
               for (int j=0; j<numPointsOriginal; j++)
               {
                 double* p = pointsOriginal->GetPoint(j);
                 vtkIdType id = vtkNewPoints->InsertNextPoint(p);
                 container->GetPointIds()->InsertNextId(id);
               }
             }
           }
           else if ( mask->GetPixel(idxStart) == 0 || mask->GetPixel(idxEnd) == 0 )
           {
             for (int j=0; j<numPointsOriginal; j++)
             {
               double* p = pointsOriginal->GetPoint(j);
               vtkIdType id = vtkNewPoints->InsertNextPoint(p);
               container->GetPointIds()->InsertNextId(id);
             }
           }
         }
         else
         {
           if (bothEnds)
           {
             if ( mask->GetPixel(idxStart) != 0 && mask->GetPixel(idxEnd) != 0 && mask->GetLargestPossibleRegion().IsInside(idxStart) && mask->GetLargestPossibleRegion().IsInside(idxEnd) )
             {
               for (int j=0; j<numPointsOriginal; j++)
               {
                 double* p = pointsOriginal->GetPoint(j);
                 vtkIdType id = vtkNewPoints->InsertNextPoint(p);
                 container->GetPointIds()->InsertNextId(id);
               }
             }
           }
           else if ( (mask->GetPixel(idxStart) != 0 && mask->GetLargestPossibleRegion().IsInside(idxStart)) || (mask->GetPixel(idxEnd) != 0 && mask->GetLargestPossibleRegion().IsInside(idxEnd)) )
           {
             for (int j=0; j<numPointsOriginal; j++)
             {
               double* p = pointsOriginal->GetPoint(j);
               vtkIdType id = vtkNewPoints->InsertNextPoint(p);
               container->GetPointIds()->InsertNextId(id);
             }
           }
         }
       }
     }
 
     if (container->GetNumberOfPoints()>0)
     {
       new_weights.push_back(fibCopy->GetFiberWeight(i));
       vtkNewCells->InsertNextCell(container);
     }
   }
 
   if (vtkNewCells->GetNumberOfCells()<=0)
     return nullptr;
 
   vtkSmartPointer<vtkPolyData> newPolyData = vtkSmartPointer<vtkPolyData>::New();
   newPolyData->SetPoints(vtkNewPoints);
   newPolyData->SetLines(vtkNewCells);
   mitk::FiberBundle::Pointer newfib = mitk::FiberBundle::New(newPolyData);
   for (unsigned int i=0; i<new_weights.size(); i++)
     newfib->SetFiberWeight(i, new_weights.at(i));
   return newfib;
 }
 
 mitk::FiberBundle::Pointer mitk::FiberBundle::RemoveFibersOutside(ItkUcharImgType* mask, bool invert)
 {
   float minSpacing = 1;
   if(mask->GetSpacing()[0]<mask->GetSpacing()[1] && mask->GetSpacing()[0]<mask->GetSpacing()[2])
     minSpacing = mask->GetSpacing()[0];
   else if (mask->GetSpacing()[1] < mask->GetSpacing()[2])
     minSpacing = mask->GetSpacing()[1];
   else
     minSpacing = mask->GetSpacing()[2];
 
   mitk::FiberBundle::Pointer fibCopy = this->GetDeepCopy();
   fibCopy->ResampleLinear(minSpacing/10);
   vtkSmartPointer<vtkPolyData> PolyData =fibCopy->GetFiberPolyData();
 
   vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
   vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
 
   MITK_INFO << "Cutting fibers";
   boost::progress_display disp(m_NumFibers);
   for (int i=0; i<m_NumFibers; i++)
   {
     ++disp;
 
     vtkCell* cell = PolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     if (numPoints>1)
     {
       int newNumPoints = 0;
       for (int j=0; j<numPoints; j++)
       {
         double* p = points->GetPoint(j);
 
         itk::Point<float, 3> itkP;
         itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2];
         itk::Index<3> idx;
         mask->TransformPhysicalPointToIndex(itkP, idx);
 
         if ( mask->GetPixel(idx) != 0 && mask->GetLargestPossibleRegion().IsInside(idx) && !invert )
         {
           vtkIdType id = vtkNewPoints->InsertNextPoint(p);
           container->GetPointIds()->InsertNextId(id);
           newNumPoints++;
         }
         else if ( (mask->GetPixel(idx) == 0 || !mask->GetLargestPossibleRegion().IsInside(idx)) && invert )
         {
           vtkIdType id = vtkNewPoints->InsertNextPoint(p);
           container->GetPointIds()->InsertNextId(id);
           newNumPoints++;
         }
         else if (newNumPoints>0)
         {
           vtkNewCells->InsertNextCell(container);
 
           newNumPoints = 0;
           container = vtkSmartPointer<vtkPolyLine>::New();
         }
       }
 
       if (newNumPoints>0)
         vtkNewCells->InsertNextCell(container);
     }
 
   }
 
   if (vtkNewCells->GetNumberOfCells()<=0)
     return nullptr;
 
   vtkSmartPointer<vtkPolyData> newPolyData = vtkSmartPointer<vtkPolyData>::New();
   newPolyData->SetPoints(vtkNewPoints);
   newPolyData->SetLines(vtkNewCells);
   mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(newPolyData);
   newFib->Compress(0.1);
   return newFib;
 }
 
 mitk::FiberBundle::Pointer mitk::FiberBundle::ExtractFiberSubset(DataNode* roi, DataStorage* storage)
 {
   if (roi==nullptr || !(dynamic_cast<PlanarFigure*>(roi->GetData()) || dynamic_cast<PlanarFigureComposite*>(roi->GetData())) )
     return nullptr;
 
   std::vector<long> tmp = ExtractFiberIdSubset(roi, storage);
 
   if (tmp.size()<=0)
     return mitk::FiberBundle::New();
   vtkSmartPointer<vtkFloatArray> weights = vtkSmartPointer<vtkFloatArray>::New();
   vtkSmartPointer<vtkPolyData> pTmp = GeneratePolyDataByIds(tmp, weights);
   mitk::FiberBundle::Pointer fib = mitk::FiberBundle::New(pTmp);
   fib->SetFiberWeights(weights);
   return fib;
 }
 
 std::vector<long> mitk::FiberBundle::ExtractFiberIdSubset(DataNode *roi, DataStorage* storage)
 {
   std::vector<long> result;
   if (roi==nullptr || roi->GetData()==nullptr)
     return result;
 
   mitk::PlanarFigureComposite::Pointer pfc = dynamic_cast<mitk::PlanarFigureComposite*>(roi->GetData());
   if (!pfc.IsNull()) // handle composite
   {
     DataStorage::SetOfObjects::ConstPointer children = storage->GetDerivations(roi);
     if (children->size()==0)
       return result;
 
     switch (pfc->getOperationType())
     {
     case 0: // AND
     {
       MITK_INFO << "AND";
       result = this->ExtractFiberIdSubset(children->ElementAt(0), storage);
       std::vector<long>::iterator it;
       for (unsigned int i=1; i<children->Size(); ++i)
       {
         std::vector<long> inRoi = this->ExtractFiberIdSubset(children->ElementAt(i), storage);
 
         std::vector<long> rest(std::min(result.size(),inRoi.size()));
         it = std::set_intersection(result.begin(), result.end(), inRoi.begin(), inRoi.end(), rest.begin() );
         rest.resize( it - rest.begin() );
         result = rest;
       }
       break;
     }
     case 1: // OR
     {
       MITK_INFO << "OR";
       result = ExtractFiberIdSubset(children->ElementAt(0), storage);
       std::vector<long>::iterator it;
       for (unsigned int i=1; i<children->Size(); ++i)
       {
         it = result.end();
         std::vector<long> inRoi = ExtractFiberIdSubset(children->ElementAt(i), storage);
         result.insert(it, inRoi.begin(), inRoi.end());
       }
 
       // remove duplicates
       sort(result.begin(), result.end());
       it = unique(result.begin(), result.end());
       result.resize( it - result.begin() );
       break;
     }
     case 2: // NOT
     {
       MITK_INFO << "NOT";
       for(long i=0; i<this->GetNumFibers(); i++)
         result.push_back(i);
 
       std::vector<long>::iterator it;
       for (unsigned int i=0; i<children->Size(); ++i)
       {
         std::vector<long> inRoi = ExtractFiberIdSubset(children->ElementAt(i), storage);
 
         std::vector<long> rest(result.size()-inRoi.size());
         it = std::set_difference(result.begin(), result.end(), inRoi.begin(), inRoi.end(), rest.begin() );
         rest.resize( it - rest.begin() );
         result = rest;
       }
       break;
     }
     }
   }
   else if ( dynamic_cast<mitk::PlanarFigure*>(roi->GetData()) )  // actual extraction
   {
     if ( dynamic_cast<mitk::PlanarPolygon*>(roi->GetData()) )
     {
       mitk::PlanarFigure::Pointer planarPoly = dynamic_cast<mitk::PlanarFigure*>(roi->GetData());
 
       //create vtkPolygon using controlpoints from planarFigure polygon
       vtkSmartPointer<vtkPolygon> polygonVtk = vtkSmartPointer<vtkPolygon>::New();
       for (unsigned int i=0; i<planarPoly->GetNumberOfControlPoints(); ++i)
       {
         itk::Point<double,3> p = planarPoly->GetWorldControlPoint(i);
         vtkIdType id = polygonVtk->GetPoints()->InsertNextPoint(p[0], p[1], p[2] );
         polygonVtk->GetPointIds()->InsertNextId(id);
       }
 
       MITK_INFO << "Extracting with polygon";
       boost::progress_display disp(m_NumFibers);
       for (int i=0; i<m_NumFibers; i++)
       {
         ++disp ;
         vtkCell* cell = m_FiberPolyData->GetCell(i);
         int numPoints = cell->GetNumberOfPoints();
         vtkPoints* points = cell->GetPoints();
 
         for (int j=0; j<numPoints-1; j++)
         {
           // Inputs
           double p1[3] = {0,0,0};
           points->GetPoint(j, p1);
           double p2[3] = {0,0,0};
           points->GetPoint(j+1, p2);
           double tolerance = 0.001;
 
           // Outputs
           double t = 0; // Parametric coordinate of intersection (0 (corresponding to p1) to 1 (corresponding to p2))
           double x[3] = {0,0,0}; // The coordinate of the intersection
           double pcoords[3] = {0,0,0};
           int subId = 0;
 
           int iD = polygonVtk->IntersectWithLine(p1, p2, tolerance, t, x, pcoords, subId);
           if (iD!=0)
           {
             result.push_back(i);
             break;
           }
         }
       }
     }
     else if ( dynamic_cast<mitk::PlanarCircle*>(roi->GetData()) )
     {
       mitk::PlanarFigure::Pointer planarFigure = dynamic_cast<mitk::PlanarFigure*>(roi->GetData());
       Vector3D planeNormal = planarFigure->GetPlaneGeometry()->GetNormal();
       planeNormal.Normalize();
 
       //calculate circle radius
       mitk::Point3D V1w = planarFigure->GetWorldControlPoint(0); //centerPoint
       mitk::Point3D V2w  = planarFigure->GetWorldControlPoint(1); //radiusPoint
 
       double radius = V1w.EuclideanDistanceTo(V2w);
       radius *= radius;
 
       MITK_INFO << "Extracting with circle";
       boost::progress_display disp(m_NumFibers);
       for (int i=0; i<m_NumFibers; i++)
       {
         ++disp ;
         vtkCell* cell = m_FiberPolyData->GetCell(i);
         int numPoints = cell->GetNumberOfPoints();
         vtkPoints* points = cell->GetPoints();
 
         for (int j=0; j<numPoints-1; j++)
         {
           // Inputs
           double p1[3] = {0,0,0};
           points->GetPoint(j, p1);
           double p2[3] = {0,0,0};
           points->GetPoint(j+1, p2);
 
           // Outputs
           double t = 0; // Parametric coordinate of intersection (0 (corresponding to p1) to 1 (corresponding to p2))
           double x[3] = {0,0,0}; // The coordinate of the intersection
 
           int iD = vtkPlane::IntersectWithLine(p1,p2,planeNormal.GetDataPointer(),V1w.GetDataPointer(),t,x);
 
           if (iD!=0)
           {
             double dist = (x[0]-V1w[0])*(x[0]-V1w[0])+(x[1]-V1w[1])*(x[1]-V1w[1])+(x[2]-V1w[2])*(x[2]-V1w[2]);
             if( dist <= radius)
             {
               result.push_back(i);
               break;
             }
           }
         }
       }
     }
     return result;
   }
 
   return result;
 }
 
 void mitk::FiberBundle::UpdateFiberGeometry()
 {
   vtkSmartPointer<vtkCleanPolyData> cleaner = vtkSmartPointer<vtkCleanPolyData>::New();
   cleaner->SetInputData(m_FiberPolyData);
   cleaner->PointMergingOff();
   cleaner->Update();
   m_FiberPolyData = cleaner->GetOutput();
 
   m_FiberLengths.clear();
   m_MeanFiberLength = 0;
   m_MedianFiberLength = 0;
   m_LengthStDev = 0;
   m_NumFibers = m_FiberPolyData->GetNumberOfCells();
 
   if (m_FiberColors==nullptr || m_FiberColors->GetNumberOfTuples()!=m_FiberPolyData->GetNumberOfPoints())
     this->ColorFibersByOrientation();
 
   if (m_FiberWeights->GetNumberOfValues()!=m_NumFibers)
   {
     m_FiberWeights = vtkSmartPointer<vtkFloatArray>::New();
     m_FiberWeights->SetName("FIBER_WEIGHTS");
     m_FiberWeights->SetNumberOfValues(m_NumFibers);
     this->SetFiberWeights(1);
   }
 
   if (m_NumFibers<=0) // no fibers present; apply default geometry
   {
     m_MinFiberLength = 0;
     m_MaxFiberLength = 0;
     mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New();
     geometry->SetImageGeometry(false);
     float b[] = {0, 1, 0, 1, 0, 1};
     geometry->SetFloatBounds(b);
     SetGeometry(geometry);
     return;
   }
   double b[6];
   m_FiberPolyData->GetBounds(b);
 
   // calculate statistics
   for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int p = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
     float length = 0;
     for (int j=0; j<p-1; j++)
     {
       double p1[3];
       points->GetPoint(j, p1);
       double p2[3];
       points->GetPoint(j+1, p2);
 
       float dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2]));
       length += dist;
     }
     m_FiberLengths.push_back(length);
     m_MeanFiberLength += length;
     if (i==0)
     {
       m_MinFiberLength = length;
       m_MaxFiberLength = length;
     }
     else
     {
       if (length<m_MinFiberLength)
         m_MinFiberLength = length;
       if (length>m_MaxFiberLength)
         m_MaxFiberLength = length;
     }
   }
   m_MeanFiberLength /= m_NumFibers;
 
   std::vector< float > sortedLengths = m_FiberLengths;
   std::sort(sortedLengths.begin(), sortedLengths.end());
   for (int i=0; i<m_NumFibers; i++)
     m_LengthStDev += (m_MeanFiberLength-sortedLengths.at(i))*(m_MeanFiberLength-sortedLengths.at(i));
   if (m_NumFibers>1)
     m_LengthStDev /= (m_NumFibers-1);
   else
     m_LengthStDev = 0;
   m_LengthStDev = std::sqrt(m_LengthStDev);
   m_MedianFiberLength = sortedLengths.at(m_NumFibers/2);
 
   mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New();
   geometry->SetFloatBounds(b);
   this->SetGeometry(geometry);
 
   m_UpdateTime3D.Modified();
   m_UpdateTime2D.Modified();
 }
 
 float mitk::FiberBundle::GetFiberWeight(unsigned int fiber) const
 {
   return m_FiberWeights->GetValue(fiber);
 }
 
 void mitk::FiberBundle::SetFiberWeights(float newWeight)
 {
   for (int i=0; i<m_FiberWeights->GetNumberOfValues(); i++)
     m_FiberWeights->SetValue(i, newWeight);
 }
 
 void mitk::FiberBundle::SetFiberWeights(vtkSmartPointer<vtkFloatArray> weights)
 {
   if (m_NumFibers!=weights->GetNumberOfValues())
   {
     MITK_INFO << "Weights array not equal to number of fibers! " << weights->GetNumberOfValues() << " vs " << m_NumFibers;
     return;
   }
 
   for (int i=0; i<weights->GetNumberOfValues(); i++)
     m_FiberWeights->SetValue(i, weights->GetValue(i));
 
   m_FiberWeights->SetName("FIBER_WEIGHTS");
 }
 
 void mitk::FiberBundle::SetFiberWeight(unsigned int fiber, float weight)
 {
   m_FiberWeights->SetValue(fiber, weight);
 }
 
 void mitk::FiberBundle::SetFiberColors(vtkSmartPointer<vtkUnsignedCharArray> fiberColors)
 {
   for(long i=0; i<m_FiberPolyData->GetNumberOfPoints(); ++i)
   {
     unsigned char source[4] = {0,0,0,0};
     fiberColors->GetTypedTuple(i, source);
 
     unsigned char target[4] = {0,0,0,0};
     target[0] = source[0];
     target[1] = source[1];
     target[2] = source[2];
     target[3] = source[3];
     m_FiberColors->InsertTypedTuple(i, target);
   }
   m_UpdateTime3D.Modified();
   m_UpdateTime2D.Modified();
 }
 
 itk::Matrix< double, 3, 3 > mitk::FiberBundle::TransformMatrix(itk::Matrix< double, 3, 3 > m, double rx, double ry, double rz)
 {
   rx = rx*M_PI/180;
   ry = ry*M_PI/180;
   rz = rz*M_PI/180;
 
   itk::Matrix< double, 3, 3 > rotX; rotX.SetIdentity();
   rotX[1][1] = cos(rx);
   rotX[2][2] = rotX[1][1];
   rotX[1][2] = -sin(rx);
   rotX[2][1] = -rotX[1][2];
 
   itk::Matrix< double, 3, 3 > rotY; rotY.SetIdentity();
   rotY[0][0] = cos(ry);
   rotY[2][2] = rotY[0][0];
   rotY[0][2] = sin(ry);
   rotY[2][0] = -rotY[0][2];
 
   itk::Matrix< double, 3, 3 > rotZ; rotZ.SetIdentity();
   rotZ[0][0] = cos(rz);
   rotZ[1][1] = rotZ[0][0];
   rotZ[0][1] = -sin(rz);
   rotZ[1][0] = -rotZ[0][1];
 
   itk::Matrix< double, 3, 3 > rot = rotZ*rotY*rotX;
 
   m = rot*m;
 
   return m;
 }
 
 itk::Point<float, 3> mitk::FiberBundle::TransformPoint(vnl_vector_fixed< double, 3 > point, double rx, double ry, double rz, double tx, double ty, double tz)
 {
   rx = rx*M_PI/180;
   ry = ry*M_PI/180;
   rz = rz*M_PI/180;
 
   vnl_matrix_fixed< double, 3, 3 > rotX; rotX.set_identity();
   rotX[1][1] = cos(rx);
   rotX[2][2] = rotX[1][1];
   rotX[1][2] = -sin(rx);
   rotX[2][1] = -rotX[1][2];
 
   vnl_matrix_fixed< double, 3, 3 > rotY; rotY.set_identity();
   rotY[0][0] = cos(ry);
   rotY[2][2] = rotY[0][0];
   rotY[0][2] = sin(ry);
   rotY[2][0] = -rotY[0][2];
 
   vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity();
   rotZ[0][0] = cos(rz);
   rotZ[1][1] = rotZ[0][0];
   rotZ[0][1] = -sin(rz);
   rotZ[1][0] = -rotZ[0][1];
 
   vnl_matrix_fixed< double, 3, 3 > rot = rotZ*rotY*rotX;
 
   mitk::BaseGeometry::Pointer geom = this->GetGeometry();
   mitk::Point3D center = geom->GetCenter();
 
   point[0] -= center[0];
   point[1] -= center[1];
   point[2] -= center[2];
   point = rot*point;
   point[0] += center[0]+tx;
   point[1] += center[1]+ty;
   point[2] += center[2]+tz;
   itk::Point<float, 3> out; out[0] = point[0]; out[1] = point[1]; out[2] = point[2];
   return out;
 }
 
 void mitk::FiberBundle::TransformFibers(double rx, double ry, double rz, double tx, double ty, double tz)
 {
   rx = rx*M_PI/180;
   ry = ry*M_PI/180;
   rz = rz*M_PI/180;
 
   vnl_matrix_fixed< double, 3, 3 > rotX; rotX.set_identity();
   rotX[1][1] = cos(rx);
   rotX[2][2] = rotX[1][1];
   rotX[1][2] = -sin(rx);
   rotX[2][1] = -rotX[1][2];
 
   vnl_matrix_fixed< double, 3, 3 > rotY; rotY.set_identity();
   rotY[0][0] = cos(ry);
   rotY[2][2] = rotY[0][0];
   rotY[0][2] = sin(ry);
   rotY[2][0] = -rotY[0][2];
 
   vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity();
   rotZ[0][0] = cos(rz);
   rotZ[1][1] = rotZ[0][0];
   rotZ[0][1] = -sin(rz);
   rotZ[1][0] = -rotZ[0][1];
 
   vnl_matrix_fixed< double, 3, 3 > rot = rotZ*rotY*rotX;
 
   mitk::BaseGeometry::Pointer geom = this->GetGeometry();
   mitk::Point3D center = geom->GetCenter();
 
   vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
   vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
 
   for (int i=0; i<m_NumFibers; i++)
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     for (int j=0; j<numPoints; j++)
     {
       double* p = points->GetPoint(j);
       vnl_vector_fixed< double, 3 > dir;
       dir[0] = p[0]-center[0];
       dir[1] = p[1]-center[1];
       dir[2] = p[2]-center[2];
       dir = rot*dir;
       dir[0] += center[0]+tx;
       dir[1] += center[1]+ty;
       dir[2] += center[2]+tz;
       vtkIdType id = vtkNewPoints->InsertNextPoint(dir.data_block());
       container->GetPointIds()->InsertNextId(id);
     }
     vtkNewCells->InsertNextCell(container);
   }
 
   m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
   m_FiberPolyData->SetPoints(vtkNewPoints);
   m_FiberPolyData->SetLines(vtkNewCells);
   this->SetFiberPolyData(m_FiberPolyData, true);
 }
 
 void mitk::FiberBundle::RotateAroundAxis(double x, double y, double z)
 {
   x = x*M_PI/180;
   y = y*M_PI/180;
   z = z*M_PI/180;
 
   vnl_matrix_fixed< double, 3, 3 > rotX; rotX.set_identity();
   rotX[1][1] = cos(x);
   rotX[2][2] = rotX[1][1];
   rotX[1][2] = -sin(x);
   rotX[2][1] = -rotX[1][2];
 
   vnl_matrix_fixed< double, 3, 3 > rotY; rotY.set_identity();
   rotY[0][0] = cos(y);
   rotY[2][2] = rotY[0][0];
   rotY[0][2] = sin(y);
   rotY[2][0] = -rotY[0][2];
 
   vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity();
   rotZ[0][0] = cos(z);
   rotZ[1][1] = rotZ[0][0];
   rotZ[0][1] = -sin(z);
   rotZ[1][0] = -rotZ[0][1];
 
   mitk::BaseGeometry::Pointer geom = this->GetGeometry();
   mitk::Point3D center = geom->GetCenter();
 
   vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
   vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
 
   for (int i=0; i<m_NumFibers; i++)
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     for (int j=0; j<numPoints; j++)
     {
       double* p = points->GetPoint(j);
       vnl_vector_fixed< double, 3 > dir;
       dir[0] = p[0]-center[0];
       dir[1] = p[1]-center[1];
       dir[2] = p[2]-center[2];
       dir = rotZ*rotY*rotX*dir;
       dir[0] += center[0];
       dir[1] += center[1];
       dir[2] += center[2];
       vtkIdType id = vtkNewPoints->InsertNextPoint(dir.data_block());
       container->GetPointIds()->InsertNextId(id);
     }
     vtkNewCells->InsertNextCell(container);
   }
 
   m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
   m_FiberPolyData->SetPoints(vtkNewPoints);
   m_FiberPolyData->SetLines(vtkNewCells);
   this->SetFiberPolyData(m_FiberPolyData, true);
 }
 
 void mitk::FiberBundle::ScaleFibers(double x, double y, double z, bool subtractCenter)
 {
   MITK_INFO << "Scaling fibers";
   boost::progress_display disp(m_NumFibers);
 
   mitk::BaseGeometry* geom = this->GetGeometry();
   mitk::Point3D c = geom->GetCenter();
 
   vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
   vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
 
   for (int i=0; i<m_NumFibers; i++)
   {
     ++disp ;
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     for (int j=0; j<numPoints; j++)
     {
       double* p = points->GetPoint(j);
       if (subtractCenter)
       {
         p[0] -= c[0]; p[1] -= c[1]; p[2] -= c[2];
       }
       p[0] *= x;
       p[1] *= y;
       p[2] *= z;
       if (subtractCenter)
       {
         p[0] += c[0]; p[1] += c[1]; p[2] += c[2];
       }
       vtkIdType id = vtkNewPoints->InsertNextPoint(p);
       container->GetPointIds()->InsertNextId(id);
     }
     vtkNewCells->InsertNextCell(container);
   }
 
   m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
   m_FiberPolyData->SetPoints(vtkNewPoints);
   m_FiberPolyData->SetLines(vtkNewCells);
   this->SetFiberPolyData(m_FiberPolyData, true);
 }
 
 void mitk::FiberBundle::TranslateFibers(double x, double y, double z)
 {
   vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
   vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
 
   for (int i=0; i<m_NumFibers; i++)
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     for (int j=0; j<numPoints; j++)
     {
       double* p = points->GetPoint(j);
       p[0] += x;
       p[1] += y;
       p[2] += z;
       vtkIdType id = vtkNewPoints->InsertNextPoint(p);
       container->GetPointIds()->InsertNextId(id);
     }
     vtkNewCells->InsertNextCell(container);
   }
 
   m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
   m_FiberPolyData->SetPoints(vtkNewPoints);
   m_FiberPolyData->SetLines(vtkNewCells);
   this->SetFiberPolyData(m_FiberPolyData, true);
 }
 
 void mitk::FiberBundle::MirrorFibers(unsigned int axis)
 {
   if (axis>2)
     return;
 
   MITK_INFO << "Mirroring fibers";
   boost::progress_display disp(m_NumFibers);
 
   vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
   vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
 
   for (int i=0; i<m_NumFibers; i++)
   {
     ++disp;
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     for (int j=0; j<numPoints; j++)
     {
       double* p = points->GetPoint(j);
       p[axis] = -p[axis];
       vtkIdType id = vtkNewPoints->InsertNextPoint(p);
       container->GetPointIds()->InsertNextId(id);
     }
     vtkNewCells->InsertNextCell(container);
   }
 
   m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
   m_FiberPolyData->SetPoints(vtkNewPoints);
   m_FiberPolyData->SetLines(vtkNewCells);
   this->SetFiberPolyData(m_FiberPolyData, true);
 }
 
 void mitk::FiberBundle::RemoveDir(vnl_vector_fixed<double,3> dir, double threshold)
 {
   dir.normalize();
   vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
   vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
 
   boost::progress_display disp(m_FiberPolyData->GetNumberOfCells());
   for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
   {
     ++disp ;
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     // calculate curvatures
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     bool discard = false;
     for (int j=0; j<numPoints-1; j++)
     {
       double p1[3];
       points->GetPoint(j, p1);
       double p2[3];
       points->GetPoint(j+1, p2);
 
       vnl_vector_fixed< double, 3 > v1;
       v1[0] = p2[0]-p1[0];
       v1[1] = p2[1]-p1[1];
       v1[2] = p2[2]-p1[2];
       if (v1.magnitude()>0.001)
       {
         v1.normalize();
 
         if (fabs(dot_product(v1,dir))>threshold)
         {
           discard = true;
           break;
         }
       }
     }
     if (!discard)
     {
       for (int j=0; j<numPoints; j++)
       {
         double p1[3];
         points->GetPoint(j, p1);
 
         vtkIdType id = vtkNewPoints->InsertNextPoint(p1);
         container->GetPointIds()->InsertNextId(id);
       }
       vtkNewCells->InsertNextCell(container);
     }
   }
 
   m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
   m_FiberPolyData->SetPoints(vtkNewPoints);
   m_FiberPolyData->SetLines(vtkNewCells);
 
   this->SetFiberPolyData(m_FiberPolyData, true);
 
   //    UpdateColorCoding();
   //    UpdateFiberGeometry();
 }
 
 bool mitk::FiberBundle::ApplyCurvatureThreshold(float minRadius, bool deleteFibers)
 {
   if (minRadius<0)
     return true;
 
   vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
   vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
 
   MITK_INFO << "Applying curvature threshold";
   boost::progress_display disp(m_FiberPolyData->GetNumberOfCells());
   for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
   {
     ++disp ;
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     // calculate curvatures
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     for (int j=0; j<numPoints-2; j++)
     {
       double p1[3];
       points->GetPoint(j, p1);
       double p2[3];
       points->GetPoint(j+1, p2);
       double p3[3];
       points->GetPoint(j+2, p3);
 
       vnl_vector_fixed< float, 3 > v1, v2, v3;
 
       v1[0] = p2[0]-p1[0];
       v1[1] = p2[1]-p1[1];
       v1[2] = p2[2]-p1[2];
 
       v2[0] = p3[0]-p2[0];
       v2[1] = p3[1]-p2[1];
       v2[2] = p3[2]-p2[2];
 
       v3[0] = p1[0]-p3[0];
       v3[1] = p1[1]-p3[1];
       v3[2] = p1[2]-p3[2];
 
       float a = v1.magnitude();
       float b = v2.magnitude();
       float c = v3.magnitude();
       float r = a*b*c/std::sqrt((a+b+c)*(a+b-c)*(b+c-a)*(a-b+c)); // radius of triangle via Heron's formula (area of triangle)
 
       vtkIdType id = vtkNewPoints->InsertNextPoint(p1);
       container->GetPointIds()->InsertNextId(id);
 
       if (deleteFibers && r<minRadius)
         break;
 
       if (r<minRadius)
       {
         j += 2;
         vtkNewCells->InsertNextCell(container);
         container = vtkSmartPointer<vtkPolyLine>::New();
       }
       else if (j==numPoints-3)
       {
         id = vtkNewPoints->InsertNextPoint(p2);
         container->GetPointIds()->InsertNextId(id);
         id = vtkNewPoints->InsertNextPoint(p3);
         container->GetPointIds()->InsertNextId(id);
         vtkNewCells->InsertNextCell(container);
       }
     }
   }
 
   if (vtkNewCells->GetNumberOfCells()<=0)
     return false;
 
   m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
   m_FiberPolyData->SetPoints(vtkNewPoints);
   m_FiberPolyData->SetLines(vtkNewCells);
   this->SetFiberPolyData(m_FiberPolyData, true);
   return true;
 }
 
 bool mitk::FiberBundle::RemoveShortFibers(float lengthInMM)
 {
   MITK_INFO << "Removing short fibers";
   if (lengthInMM<=0 || lengthInMM<m_MinFiberLength)
   {
     MITK_INFO << "No fibers shorter than " << lengthInMM << " mm found!";
     return true;
   }
 
   if (lengthInMM>m_MaxFiberLength)    // can't remove all fibers
   {
     MITK_WARN << "Process aborted. No fibers would be left!";
     return false;
   }
 
   vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
   vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
   float min = m_MaxFiberLength;
 
   boost::progress_display disp(m_NumFibers);
   for (int i=0; i<m_NumFibers; i++)
   {
     ++disp;
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     if (m_FiberLengths.at(i)>=lengthInMM)
     {
       vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
       for (int j=0; j<numPoints; j++)
       {
         double* p = points->GetPoint(j);
         vtkIdType id = vtkNewPoints->InsertNextPoint(p);
         container->GetPointIds()->InsertNextId(id);
       }
       vtkNewCells->InsertNextCell(container);
       if (m_FiberLengths.at(i)<min)
         min = m_FiberLengths.at(i);
     }
   }
 
   if (vtkNewCells->GetNumberOfCells()<=0)
     return false;
 
   m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
   m_FiberPolyData->SetPoints(vtkNewPoints);
   m_FiberPolyData->SetLines(vtkNewCells);
   this->SetFiberPolyData(m_FiberPolyData, true);
   return true;
 }
 
 bool mitk::FiberBundle::RemoveLongFibers(float lengthInMM)
 {
   if (lengthInMM<=0 || lengthInMM>m_MaxFiberLength)
     return true;
 
   if (lengthInMM<m_MinFiberLength)    // can't remove all fibers
     return false;
 
   vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
   vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
 
   MITK_INFO << "Removing long fibers";
   boost::progress_display disp(m_NumFibers);
   for (int i=0; i<m_NumFibers; i++)
   {
     ++disp;
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     if (m_FiberLengths.at(i)<=lengthInMM)
     {
       vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
       for (int j=0; j<numPoints; j++)
       {
         double* p = points->GetPoint(j);
         vtkIdType id = vtkNewPoints->InsertNextPoint(p);
         container->GetPointIds()->InsertNextId(id);
       }
       vtkNewCells->InsertNextCell(container);
     }
   }
 
   if (vtkNewCells->GetNumberOfCells()<=0)
     return false;
 
   m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
   m_FiberPolyData->SetPoints(vtkNewPoints);
   m_FiberPolyData->SetLines(vtkNewCells);
   this->SetFiberPolyData(m_FiberPolyData, true);
   return true;
 }
 
 void mitk::FiberBundle::ResampleSpline(float pointDistance, double tension, double continuity, double bias )
 {
   if (pointDistance<=0)
     return;
 
   vtkSmartPointer<vtkPoints> vtkSmoothPoints = vtkSmartPointer<vtkPoints>::New(); //in smoothpoints the interpolated points representing a fiber are stored.
 
   //in vtkcells all polylines are stored, actually all id's of them are stored
   vtkSmartPointer<vtkCellArray> vtkSmoothCells = vtkSmartPointer<vtkCellArray>::New(); //cellcontainer for smoothed lines
   vtkIdType pointHelperCnt = 0;
 
   MITK_INFO << "Smoothing fibers";
   vtkSmartPointer<vtkFloatArray> newFiberWeights = vtkSmartPointer<vtkFloatArray>::New();
   newFiberWeights->SetName("FIBER_WEIGHTS");
   newFiberWeights->SetNumberOfValues(m_NumFibers);
 
   boost::progress_display disp(m_NumFibers);
 #pragma omp parallel for
   for (int i=0; i<m_NumFibers; i++)
   {
     vtkSmartPointer<vtkPoints> newPoints = vtkSmartPointer<vtkPoints>::New();
     float length = 0;
     float weight = 1;
 #pragma omp critical
     {
       length = m_FiberLengths.at(i);
       weight = m_FiberWeights->GetValue(i);
       ++disp;
       vtkCell* cell = m_FiberPolyData->GetCell(i);
       int numPoints = cell->GetNumberOfPoints();
       vtkPoints* points = cell->GetPoints();
       for (int j=0; j<numPoints; j++)
         newPoints->InsertNextPoint(points->GetPoint(j));
     }
 
     int sampling = std::ceil(length/pointDistance);
 
     vtkSmartPointer<vtkKochanekSpline> xSpline = vtkSmartPointer<vtkKochanekSpline>::New();
     vtkSmartPointer<vtkKochanekSpline> ySpline = vtkSmartPointer<vtkKochanekSpline>::New();
     vtkSmartPointer<vtkKochanekSpline> zSpline = vtkSmartPointer<vtkKochanekSpline>::New();
     xSpline->SetDefaultBias(bias); xSpline->SetDefaultTension(tension); xSpline->SetDefaultContinuity(continuity);
     ySpline->SetDefaultBias(bias); ySpline->SetDefaultTension(tension); ySpline->SetDefaultContinuity(continuity);
     zSpline->SetDefaultBias(bias); zSpline->SetDefaultTension(tension); zSpline->SetDefaultContinuity(continuity);
 
     vtkSmartPointer<vtkParametricSpline> spline = vtkSmartPointer<vtkParametricSpline>::New();
     spline->SetXSpline(xSpline);
     spline->SetYSpline(ySpline);
     spline->SetZSpline(zSpline);
     spline->SetPoints(newPoints);
 
     vtkSmartPointer<vtkParametricFunctionSource> functionSource = vtkSmartPointer<vtkParametricFunctionSource>::New();
     functionSource->SetParametricFunction(spline);
     functionSource->SetUResolution(sampling);
     functionSource->SetVResolution(sampling);
     functionSource->SetWResolution(sampling);
     functionSource->Update();
 
     vtkPolyData* outputFunction = functionSource->GetOutput();
     vtkPoints* tmpSmoothPnts = outputFunction->GetPoints(); //smoothPoints of current fiber
 
     vtkSmartPointer<vtkPolyLine> smoothLine = vtkSmartPointer<vtkPolyLine>::New();
     smoothLine->GetPointIds()->SetNumberOfIds(tmpSmoothPnts->GetNumberOfPoints());
 
 #pragma omp critical
     {
       for (int j=0; j<smoothLine->GetNumberOfPoints(); j++)
       {
         smoothLine->GetPointIds()->SetId(j, j+pointHelperCnt);
         vtkSmoothPoints->InsertNextPoint(tmpSmoothPnts->GetPoint(j));
       }
 
       newFiberWeights->SetValue(vtkSmoothCells->GetNumberOfCells(), weight);
       vtkSmoothCells->InsertNextCell(smoothLine);
       pointHelperCnt += tmpSmoothPnts->GetNumberOfPoints();
     }
   }
 
   SetFiberWeights(newFiberWeights);
   m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
   m_FiberPolyData->SetPoints(vtkSmoothPoints);
   m_FiberPolyData->SetLines(vtkSmoothCells);
   this->SetFiberPolyData(m_FiberPolyData, true);
 }
 
 void mitk::FiberBundle::ResampleSpline(float pointDistance)
 {
   ResampleSpline(pointDistance, 0, 0, 0 );
 }
 
 unsigned long mitk::FiberBundle::GetNumberOfPoints() const
 {
   unsigned long points = 0;
   for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     points += cell->GetNumberOfPoints();
   }
   return points;
 }
 
 void mitk::FiberBundle::Compress(float error)
 {
   vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
   vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
 
   MITK_INFO << "Compressing fibers";
   unsigned long numRemovedPoints = 0;
   boost::progress_display disp(m_FiberPolyData->GetNumberOfCells());
   vtkSmartPointer<vtkFloatArray> newFiberWeights = vtkSmartPointer<vtkFloatArray>::New();
   newFiberWeights->SetName("FIBER_WEIGHTS");
   newFiberWeights->SetNumberOfValues(m_NumFibers);
 
 #pragma omp parallel for
   for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
   {
 
     std::vector< vnl_vector_fixed< double, 3 > > vertices;
     float weight = 1;
 
 #pragma omp critical
     {
       ++disp;
       weight = m_FiberWeights->GetValue(i);
       vtkCell* cell = m_FiberPolyData->GetCell(i);
       int numPoints = cell->GetNumberOfPoints();
       vtkPoints* points = cell->GetPoints();
 
       for (int j=0; j<numPoints; j++)
       {
         double cand[3];
         points->GetPoint(j, cand);
         vnl_vector_fixed< double, 3 > candV;
         candV[0]=cand[0]; candV[1]=cand[1]; candV[2]=cand[2];
         vertices.push_back(candV);
       }
     }
 
     // calculate curvatures
     int numPoints = vertices.size();
     std::vector< int > removedPoints; removedPoints.resize(numPoints, 0);
     removedPoints[0]=-1; removedPoints[numPoints-1]=-1;
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     int remCounter = 0;
 
     bool pointFound = true;
     while (pointFound)
     {
       pointFound = false;
       double minError = error;
       int removeIndex = -1;
 
       for (unsigned int j=0; j<vertices.size(); j++)
       {
         if (removedPoints[j]==0)
         {
           vnl_vector_fixed< double, 3 > candV = vertices.at(j);
 
           int validP = -1;
           vnl_vector_fixed< double, 3 > pred;
           for (int k=j-1; k>=0; k--)
             if (removedPoints[k]<=0)
             {
               pred = vertices.at(k);
               validP = k;
               break;
             }
           int validS = -1;
           vnl_vector_fixed< double, 3 > succ;
           for (int k=j+1; k<numPoints; k++)
             if (removedPoints[k]<=0)
             {
               succ = vertices.at(k);
               validS = k;
               break;
             }
 
           if (validP>=0 && validS>=0)
           {
             double a = (candV-pred).magnitude();
             double b = (candV-succ).magnitude();
             double c = (pred-succ).magnitude();
             double s=0.5*(a+b+c);
             double hc=(2.0/c)*sqrt(fabs(s*(s-a)*(s-b)*(s-c)));
 
             if (hc<minError)
             {
               removeIndex = j;
               minError = hc;
               pointFound = true;
             }
           }
         }
       }
 
       if (pointFound)
       {
         removedPoints[removeIndex] = 1;
         remCounter++;
       }
     }
 
     for (int j=0; j<numPoints; j++)
     {
       if (removedPoints[j]<=0)
       {
 #pragma omp critical
         {
           vtkIdType id = vtkNewPoints->InsertNextPoint(vertices.at(j).data_block());
           container->GetPointIds()->InsertNextId(id);
         }
       }
     }
 
 #pragma omp critical
     {
       newFiberWeights->SetValue(vtkNewCells->GetNumberOfCells(), weight);
       numRemovedPoints += remCounter;
       vtkNewCells->InsertNextCell(container);
     }
   }
 
   if (vtkNewCells->GetNumberOfCells()>0)
   {
     MITK_INFO << "Removed points: " << numRemovedPoints;
     SetFiberWeights(newFiberWeights);
     m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
     m_FiberPolyData->SetPoints(vtkNewPoints);
     m_FiberPolyData->SetLines(vtkNewCells);
     this->SetFiberPolyData(m_FiberPolyData, true);
   }
 }
 
 void mitk::FiberBundle::ResampleLinear(double pointDistance)
 {
   vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
   vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
 
   MITK_INFO << "Resampling fibers (linear)";
   boost::progress_display disp(m_FiberPolyData->GetNumberOfCells());
   vtkSmartPointer<vtkFloatArray> newFiberWeights = vtkSmartPointer<vtkFloatArray>::New();
   newFiberWeights->SetName("FIBER_WEIGHTS");
   newFiberWeights->SetNumberOfValues(m_NumFibers);
 
 #pragma omp parallel for
   for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
   {
 
     std::vector< vnl_vector_fixed< double, 3 > > vertices;
     float weight = 1;
 
 #pragma omp critical
     {
       ++disp;
       weight = m_FiberWeights->GetValue(i);
       vtkCell* cell = m_FiberPolyData->GetCell(i);
       int numPoints = cell->GetNumberOfPoints();
       vtkPoints* points = cell->GetPoints();
 
       for (int j=0; j<numPoints; j++)
       {
         double cand[3];
         points->GetPoint(j, cand);
         vnl_vector_fixed< double, 3 > candV;
         candV[0]=cand[0]; candV[1]=cand[1]; candV[2]=cand[2];
         vertices.push_back(candV);
       }
     }
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     vnl_vector_fixed< double, 3 > lastV = vertices.at(0);
 
 #pragma omp critical
     {
       vtkIdType id = vtkNewPoints->InsertNextPoint(lastV.data_block());
       container->GetPointIds()->InsertNextId(id);
     }
     for (unsigned int j=1; j<vertices.size(); j++)
     {
       vnl_vector_fixed< double, 3 > vec = vertices.at(j) - lastV;
       double new_dist = vec.magnitude();
 
       if (new_dist >= pointDistance)
       {
         vnl_vector_fixed< double, 3 > newV = lastV;
         if ( new_dist-pointDistance <= mitk::eps )
         {
           vec.normalize();
           newV += vec * pointDistance;
         }
         else
         {
           // intersection between sphere (radius 'pointDistance', center 'lastV') and line (direction 'd' and point 'p')
           vnl_vector_fixed< double, 3 > p = vertices.at(j-1);
           vnl_vector_fixed< double, 3 > d = vertices.at(j) - p;
 
           double a = d[0]*d[0] + d[1]*d[1] + d[2]*d[2];
           double b = 2 * (d[0] * (p[0] - lastV[0]) + d[1] * (p[1] - lastV[1]) + d[2] * (p[2] - lastV[2]));
           double c = (p[0] - lastV[0])*(p[0] - lastV[0]) + (p[1] - lastV[1])*(p[1] - lastV[1]) + (p[2] - lastV[2])*(p[2] - lastV[2]) - pointDistance*pointDistance;
 
           double v1 =(-b + std::sqrt(b*b-4*a*c))/(2*a);
           double v2 =(-b - std::sqrt(b*b-4*a*c))/(2*a);
 
           if (v1>0)
             newV = p + d * v1;
           else if (v2>0)
             newV = p + d * v2;
           else
             MITK_INFO << "ERROR1 - linear resampling";
 
           j--;
         }
 
 #pragma omp critical
         {
           vtkIdType id = vtkNewPoints->InsertNextPoint(newV.data_block());
           container->GetPointIds()->InsertNextId(id);
         }
         lastV = newV;
       }
       else if (j==vertices.size()-1 && new_dist>0.0001)
       {
 #pragma omp critical
         {
           vtkIdType id = vtkNewPoints->InsertNextPoint(vertices.at(j).data_block());
           container->GetPointIds()->InsertNextId(id);
         }
       }
     }
 
 #pragma omp critical
     {
       newFiberWeights->SetValue(vtkNewCells->GetNumberOfCells(), weight);
       vtkNewCells->InsertNextCell(container);
     }
   }
 
   if (vtkNewCells->GetNumberOfCells()>0)
   {
     SetFiberWeights(newFiberWeights);
     m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
     m_FiberPolyData->SetPoints(vtkNewPoints);
     m_FiberPolyData->SetLines(vtkNewCells);
     this->SetFiberPolyData(m_FiberPolyData, true);
   }
 }
 
 // reapply selected colorcoding in case PolyData structure has changed
 bool mitk::FiberBundle::Equals(mitk::FiberBundle* fib, double eps)
 {
   if (fib==nullptr)
   {
     MITK_INFO << "Reference bundle is nullptr!";
     return false;
   }
 
   if (m_NumFibers!=fib->GetNumFibers())
   {
     MITK_INFO << "Unequal number of fibers!";
     MITK_INFO << m_NumFibers << " vs. " << fib->GetNumFibers();
     return false;
   }
 
   for (int i=0; i<m_NumFibers; i++)
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     vtkCell* cell2 = fib->GetFiberPolyData()->GetCell(i);
     int numPoints2 = cell2->GetNumberOfPoints();
     vtkPoints* points2 = cell2->GetPoints();
 
     if (numPoints2!=numPoints)
     {
       MITK_INFO << "Unequal number of points in fiber " << i << "!";
       MITK_INFO << numPoints2 << " vs. " << numPoints;
       return false;
     }
 
     for (int j=0; j<numPoints; j++)
     {
       double* p1 = points->GetPoint(j);
       double* p2 = points2->GetPoint(j);
       if (fabs(p1[0]-p2[0])>eps || fabs(p1[1]-p2[1])>eps || fabs(p1[2]-p2[2])>eps)
       {
         MITK_INFO << "Unequal points in fiber " << i << " at position " << j << "!";
         MITK_INFO << "p1: " << p1[0] << ", " << p1[1] << ", " << p1[2];
         MITK_INFO << "p2: " << p2[0] << ", " << p2[1] << ", " << p2[2];
         return false;
       }
     }
   }
 
   return true;
 }
 
 void mitk::FiberBundle::PrintSelf(std::ostream &os, itk::Indent indent) const
 {
   os << indent << this->GetNameOfClass() << ":\n";
   os << indent << "Number of fibers: " << this->GetNumFibers() << std::endl;
   os << indent << "Min. fiber length: " << this->GetMinFiberLength() << std::endl;
   os << indent << "Max. fiber length: " << this->GetMaxFiberLength() << std::endl;
   os << indent << "Mean fiber length: " << this->GetMeanFiberLength() << std::endl;
   os << indent << "Median fiber length: " << this->GetMedianFiberLength() << std::endl;
   os << indent << "STDEV fiber length: " << this->GetLengthStDev() << std::endl;
   os << indent << "Number of points: " << this->GetNumberOfPoints() << std::endl;
   Superclass::PrintSelf(os, indent);
 }
 
 /* ESSENTIAL IMPLEMENTATION OF SUPERCLASS METHODS */
 void mitk::FiberBundle::UpdateOutputInformation()
 {
 
 }
 void mitk::FiberBundle::SetRequestedRegionToLargestPossibleRegion()
 {
 
 }
 bool mitk::FiberBundle::RequestedRegionIsOutsideOfTheBufferedRegion()
 {
   return false;
 }
 bool mitk::FiberBundle::VerifyRequestedRegion()
 {
   return true;
 }
 void mitk::FiberBundle::SetRequestedRegion(const itk::DataObject* )
 {
 
 }
diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.h b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.h
index b57e1ab363..863c50b64d 100644
--- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.h
+++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.h
@@ -1,179 +1,179 @@
 /*===================================================================
 
 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.
 
 ===================================================================*/
 
 
 #ifndef _MITK_FiberBundle_H
 #define _MITK_FiberBundle_H
 
 //includes for MITK datastructure
 #include <mitkBaseData.h>
 #include <MitkFiberTrackingExports.h>
 #include <mitkImage.h>
 #include <mitkDataStorage.h>
 #include <mitkPlanarFigure.h>
 #include <mitkPixelTypeTraits.h>
 #include <mitkPlanarFigureComposite.h>
 
 
 //includes storing fiberdata
 #include <vtkSmartPointer.h>
 #include <vtkPolyData.h>
 #include <vtkPoints.h>
 #include <vtkDataSet.h>
 #include <vtkTransform.h>
 #include <vtkFloatArray.h>
 
 
 namespace mitk {
 
 /**
    * \brief Base Class for Fiber Bundles;   */
 class MITKFIBERTRACKING_EXPORT FiberBundle : public BaseData
 {
 public:
 
     typedef itk::Image<unsigned char, 3> ItkUcharImgType;
 
     // fiber colorcodings
     static const char* FIBER_ID_ARRAY;
 
     virtual void UpdateOutputInformation() override;
     virtual void SetRequestedRegionToLargestPossibleRegion() override;
     virtual bool RequestedRegionIsOutsideOfTheBufferedRegion() override;
     virtual bool VerifyRequestedRegion() override;
     virtual void SetRequestedRegion(const itk::DataObject*) override;
 
     mitkClassMacro( FiberBundle, BaseData )
     itkFactorylessNewMacro(Self)
     itkCloneMacro(Self)
     mitkNewMacro1Param(Self, vtkSmartPointer<vtkPolyData>) // custom constructor
 
     // colorcoding related methods
     void ColorFibersByFiberWeights(bool opacity, bool normalize);
     void ColorFibersByCurvature(bool opacity, bool normalize);
     void ColorFibersByScalarMap(mitk::Image::Pointer, bool opacity, bool normalize);
     template <typename TPixel>
     void ColorFibersByScalarMap(const mitk::PixelType pixelType, mitk::Image::Pointer, bool opacity, bool normalize);
     void ColorFibersByOrientation();
     void SetFiberOpacity(vtkDoubleArray *FAValArray);
     void ResetFiberOpacity();
     void SetFiberColors(vtkSmartPointer<vtkUnsignedCharArray> fiberColors);
     void SetFiberColors(float r, float g, float b, float alpha=255);
     vtkSmartPointer<vtkUnsignedCharArray> GetFiberColors() const { return m_FiberColors; }
 
     // fiber compression
     void Compress(float error = 0.0);
 
     // fiber resampling
     void ResampleSpline(float pointDistance=1);
     void ResampleSpline(float pointDistance, double tension, double continuity, double bias );
     void ResampleLinear(double pointDistance=1);
 
-    mitk::FiberBundle::Pointer FilterByWeights(float weight_thr);
+    mitk::FiberBundle::Pointer FilterByWeights(float weight_thr, bool invert=false);
     bool RemoveShortFibers(float lengthInMM);
     bool RemoveLongFibers(float lengthInMM);
     bool ApplyCurvatureThreshold(float minRadius, bool deleteFibers);
     void MirrorFibers(unsigned int axis);
     void RotateAroundAxis(double x, double y, double z);
     void TranslateFibers(double x, double y, double z);
     void ScaleFibers(double x, double y, double z, bool subtractCenter=true);
     void TransformFibers(double rx, double ry, double rz, double tx, double ty, double tz);
     void RemoveDir(vnl_vector_fixed<double,3> dir, double threshold);
     itk::Point<float, 3> TransformPoint(vnl_vector_fixed< double, 3 > point, double rx, double ry, double rz, double tx, double ty, double tz);
     itk::Matrix< double, 3, 3 > TransformMatrix(itk::Matrix< double, 3, 3 > m, double rx, double ry, double rz);
 
     // add/subtract fibers
     FiberBundle::Pointer AddBundle(FiberBundle* fib);
     mitk::FiberBundle::Pointer AddBundles(std::vector< mitk::FiberBundle::Pointer > fibs);
     FiberBundle::Pointer SubtractBundle(FiberBundle* fib);
 
     // fiber subset extraction
     FiberBundle::Pointer           ExtractFiberSubset(DataNode *roi, DataStorage* storage);
     std::vector<long>              ExtractFiberIdSubset(DataNode* roi, DataStorage* storage);
     FiberBundle::Pointer           ExtractFiberSubset(ItkUcharImgType* mask, bool anyPoint, bool invert=false, bool bothEnds=true, float fraction=0.0, bool do_resampling=true);
     FiberBundle::Pointer           RemoveFibersOutside(ItkUcharImgType* mask, bool invert=false);
     float                          GetOverlap(ItkUcharImgType* mask, bool do_resampling);
     mitk::FiberBundle::Pointer     SubsampleFibers(float factor);
 
     // get/set data
     vtkSmartPointer<vtkFloatArray> GetFiberWeights() const { return m_FiberWeights; }
     float GetFiberWeight(unsigned int fiber) const;
     void SetFiberWeights(float newWeight);
     void SetFiberWeight(unsigned int fiber, float weight);
     void SetFiberWeights(vtkSmartPointer<vtkFloatArray> weights);
     void SetFiberPolyData(vtkSmartPointer<vtkPolyData>, bool updateGeometry = true);
     vtkSmartPointer<vtkPolyData> GetFiberPolyData() const;
     itkGetConstMacro( NumFibers, int)
     //itkGetMacro( FiberSampling, int)
     itkGetConstMacro( MinFiberLength, float )
     itkGetConstMacro( MaxFiberLength, float )
     itkGetConstMacro( MeanFiberLength, float )
     itkGetConstMacro( MedianFiberLength, float )
     itkGetConstMacro( LengthStDev, float )
     itkGetConstMacro( UpdateTime2D, itk::TimeStamp )
     itkGetConstMacro( UpdateTime3D, itk::TimeStamp )
     void RequestUpdate2D(){ m_UpdateTime2D.Modified(); }
     void RequestUpdate3D(){ m_UpdateTime3D.Modified(); }
     void RequestUpdate(){ m_UpdateTime2D.Modified(); m_UpdateTime3D.Modified(); }
 
     unsigned long GetNumberOfPoints() const;
 
     // copy fiber bundle
     mitk::FiberBundle::Pointer GetDeepCopy();
 
     // compare fiber bundles
     bool Equals(FiberBundle* fib, double eps=0.01);
 
     itkSetMacro( ReferenceGeometry, mitk::BaseGeometry::Pointer )
     itkGetConstMacro( ReferenceGeometry, mitk::BaseGeometry::Pointer )
 
 protected:
 
     FiberBundle( vtkPolyData* fiberPolyData = nullptr );
     virtual ~FiberBundle();
 
     vtkSmartPointer<vtkPolyData>    GeneratePolyDataByIds(std::vector<long> fiberIds, vtkSmartPointer<vtkFloatArray> weights);
     void                            GenerateFiberIds();
     itk::Point<float, 3>            GetItkPoint(double point[3]);
     void                            UpdateFiberGeometry();
     virtual void                    PrintSelf(std::ostream &os, itk::Indent indent) const override;
 
 private:
 
     // actual fiber container
     vtkSmartPointer<vtkPolyData>  m_FiberPolyData;
 
     // contains fiber ids
     vtkSmartPointer<vtkDataSet>   m_FiberIdDataSet;
 
     int   m_NumFibers;
 
     vtkSmartPointer<vtkUnsignedCharArray> m_FiberColors;
     vtkSmartPointer<vtkFloatArray> m_FiberWeights;
     std::vector< float > m_FiberLengths;
     float   m_MinFiberLength;
     float   m_MaxFiberLength;
     float   m_MeanFiberLength;
     float   m_MedianFiberLength;
     float   m_LengthStDev;
     itk::TimeStamp m_UpdateTime2D;
     itk::TimeStamp m_UpdateTime3D;
     mitk::BaseGeometry::Pointer m_ReferenceGeometry;
 };
 
 } // namespace mitk
 
 #endif /*  _MITK_FiberBundle_H */
diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FitFibersToImage.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FitFibersToImage.cpp
index e58629356e..421f49ffea 100755
--- a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FitFibersToImage.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FitFibersToImage.cpp
@@ -1,203 +1,217 @@
 /*===================================================================
 
 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("save_res", "", mitkCommandLineParser::Bool, "Save Residuals:", "save residual images", false);
+  parser.addArgument("save_weights", "", mitkCommandLineParser::Bool, "Save Weights:", "save fiber weights in a separate text file", false);
   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"]);
 
+  bool save_weights = false;
+  if (parsedArgs.count("save_weights"))
+    save_weights = us::any_cast<bool>(parsedArgs["save_weights"]);
+
   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("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");
+
+      if (save_weights)
+      {
+        ofstream logfile;
+        logfile.open (outRoot + name + "_weights.txt");
+        for (int f=0; f<output_tracts.at(bundle)->GetNumFibers(); ++f)
+          logfile << output_tracts.at(bundle)->GetFiberWeight(f) << "\n";
+        logfile.close();
+      }
     }
   }
   catch (itk::ExceptionObject e)
   {
     std::cout << e;
     return EXIT_FAILURE;
   }
   catch (std::exception e)
   {
     std::cout << e.what();
     return EXIT_FAILURE;
   }
   catch (...)
   {
     std::cout << "ERROR!?!";
     return EXIT_FAILURE;
   }
   return EXIT_SUCCESS;
 }
diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/CMakeLists.txt b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/CMakeLists.txt
index 0f94e117bf..760a65a867 100755
--- a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/CMakeLists.txt
+++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/CMakeLists.txt
@@ -1,43 +1,44 @@
 option(BUILD_DiffusionTractographyEvaluationCmdApps "Build commandline tools for diffusion fiber tractography evaluation" OFF)
 
 if(BUILD_DiffusionTractographyEvaluationCmdApps OR MITK_BUILD_ALL_APPS)
 
   # needed include directories
   include_directories(
     ${CMAKE_CURRENT_SOURCE_DIR}
     ${CMAKE_CURRENT_BINARY_DIR}
     )
 
     # list of diffusion cmdapps
     # if an app requires additional dependencies
     # they are added after a "^^" and separated by "_"
     set( diffusionTractographyEvaluationcmdapps
     PeaksAngularError^^MitkFiberTracking
     TractometerMetrics^^MitkFiberTracking
     MergeDuplicateTracts^^MitkFiberTracking
     TractPlausibility^^MitkFiberTracking
     TractPlausibilityFit^^MitkFiberTracking
+    EvaluateLiFE^^MitkFiberTracking
     # LocalDirectionalFiberPlausibility^^MitkFiberTracking # HAS TO USE NEW PEAK IMAGE FORMAT
     )
 
     foreach(diffusionTractographyEvaluationcmdapp ${diffusionTractographyEvaluationcmdapps})
       # extract cmd app name and dependencies
       string(REPLACE "^^" "\\;" cmdapp_info ${diffusionTractographyEvaluationcmdapp})
       set(cmdapp_info_list ${cmdapp_info})
       list(GET cmdapp_info_list 0 appname)
       list(GET cmdapp_info_list 1 raw_dependencies)
       string(REPLACE "_" "\\;" dependencies "${raw_dependencies}")
       set(dependencies_list ${dependencies})
 
       mitkFunctionCreateCommandLineApp(
         NAME ${appname}
         DEPENDS MitkCore MitkDiffusionCore ${dependencies_list}
         PACKAGE_DEPENDS ITK
       )
     endforeach()
 
   if(EXECUTABLE_IS_ENABLED)
     MITK_INSTALL_TARGETS(EXECUTABLES ${EXECUTABLE_TARGET})
   endif()
 
   endif()
diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/EvaluateLiFE.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/EvaluateLiFE.cpp
new file mode 100755
index 0000000000..89d7ef5e9d
--- /dev/null
+++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/EvaluateLiFE.cpp
@@ -0,0 +1,217 @@
+/*===================================================================
+
+The Medical Imaging Interaction Toolkit (MITK)
+
+Copyright (c) German Cancer Research Center,
+Division of Medical and Biological Informatics.
+All rights reserved.
+
+This software is distributed WITHOUT ANY WARRANTY; without
+even the implied warranty of MERCHANTABILITY or FITNESS FOR
+A PARTICULAR PURPOSE.
+
+See LICENSE.txt or http://www.mitk.org for details.
+
+===================================================================*/
+
+#include <mitkBaseData.h>
+#include <mitkImageCast.h>
+#include <mitkImageToItk.h>
+#include <metaCommand.h>
+#include <mitkCommandLineParser.h>
+#include <usAny.h>
+#include <mitkIOUtil.h>
+#include <boost/lexical_cast.hpp>
+#include <itksys/SystemTools.hxx>
+#include <itkDirectory.h>
+#include <mitkFiberBundle.h>
+#include <vtkTransformPolyDataFilter.h>
+#include <fstream>
+#include <chrono>
+#include <boost/progress.hpp>
+
+using namespace std;
+typedef itksys::SystemTools ist;
+typedef itk::Image<unsigned char, 3>    ItkUcharImgType;
+typedef std::tuple< ItkUcharImgType::Pointer, std::string > MaskType;
+
+ItkUcharImgType::Pointer LoadItkMaskImage(const std::string& filename)
+{
+  mitk::Image::Pointer img = dynamic_cast<mitk::Image*>(mitk::IOUtil::Load(filename)[0].GetPointer());
+  ItkUcharImgType::Pointer itkMask = ItkUcharImgType::New();
+  mitk::CastToItkImage(img, itkMask);
+  return itkMask;
+}
+
+std::vector< MaskType > get_file_list(const std::string& path)
+{
+  std::chrono::milliseconds ms = std::chrono::duration_cast< std::chrono::milliseconds >(std::chrono::system_clock::now().time_since_epoch());
+  std::srand(ms.count());
+
+  std::vector< MaskType > mask_list;
+
+  itk::Directory::Pointer dir = itk::Directory::New();
+
+  if (dir->Load(path.c_str()))
+  {
+    int n = dir->GetNumberOfFiles();
+    int num_images = 0;
+
+    std::vector< int > im_indices;
+    for (int r = 0; r < n; r++)
+    {
+      const char *filename = dir->GetFile(r);
+      std::string ext = ist::GetFilenameExtension(filename);
+      if (ext==".nii" || ext==".nii.gz" || ext==".nrrd")
+      {
+        ++num_images;
+        im_indices.push_back(r);
+      }
+    }
+
+    int c = -1;
+    for (int r : im_indices)
+    {
+      c++;
+      const char *filename = dir->GetFile(r);
+
+      MITK_INFO << "Loading " << ist::GetFilenameWithoutExtension(filename);
+      streambuf *old = cout.rdbuf(); // <-- save
+      stringstream ss;
+      std::cout.rdbuf (ss.rdbuf());       // <-- redirect
+      MaskType m(LoadItkMaskImage(path + '/' + filename), ist::GetFilenameName(filename));
+      mask_list.push_back(m);
+      std::cout.rdbuf (old);              // <-- restore
+    }
+  }
+
+  return mask_list;
+}
+
+/*!
+\brief
+*/
+int main(int argc, char* argv[])
+{
+  mitkCommandLineParser parser;
+
+  parser.setTitle("Evaluate LiFE results");
+  parser.setCategory("Fiber Tracking Evaluation");
+  parser.setDescription("");
+  parser.setContributor("MIC");
+
+  parser.setArgumentPrefix("--", "-");
+  parser.addArgument("input", "i", mitkCommandLineParser::InputFile, "Input:", "input tractogram (.fib, vtk ascii file format)", us::Any(), false);
+  parser.addArgument("out", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output folder", us::Any(), false);
+  parser.addArgument("reference_mask_folder", "m", mitkCommandLineParser::String, "Reference Mask Folder:", "reference masks of known bundles", false);
+  parser.addArgument("overlap", "", mitkCommandLineParser::Float, "", "", 0.7);
+  parser.addArgument("steps", "", mitkCommandLineParser::Int, "", "", 10);
+
+  map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
+  if (parsedArgs.size()==0)
+    return EXIT_FAILURE;
+
+  string fibFile = us::any_cast<string>(parsedArgs["input"]);
+  string reference_mask_folder = us::any_cast<string>(parsedArgs["reference_mask_folder"]);
+  string out_folder = us::any_cast<string>(parsedArgs["out"]);
+
+  float overlap = 0.7;
+  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"]);
+
+  try
+  {
+    std::vector< MaskType > known_tract_masks = get_file_list(reference_mask_folder);
+    if (known_tract_masks.empty())
+      return EXIT_FAILURE;
+
+    mitk::FiberBundle::Pointer inputTractogram = dynamic_cast<mitk::FiberBundle*>(mitk::IOUtil::Load(fibFile)[0].GetPointer());
+
+    // resample fibers
+    float minSpacing = 1;
+    if(std::get<0>(known_tract_masks.at(0))->GetSpacing()[0]<std::get<0>(known_tract_masks.at(0))->GetSpacing()[1] && std::get<0>(known_tract_masks.at(0))->GetSpacing()[0]<std::get<0>(known_tract_masks.at(0))->GetSpacing()[2])
+      minSpacing = std::get<0>(known_tract_masks.at(0))->GetSpacing()[0];
+    else if (std::get<0>(known_tract_masks.at(0))->GetSpacing()[1] < std::get<0>(known_tract_masks.at(0))->GetSpacing()[2])
+      minSpacing = std::get<0>(known_tract_masks.at(0))->GetSpacing()[1];
+    else
+      minSpacing = std::get<0>(known_tract_masks.at(0))->GetSpacing()[2];
+    inputTractogram->ResampleLinear(minSpacing/5);
+
+    std::vector< float > weights;
+    for (int i=0; i<inputTractogram->GetNumFibers(); ++i)
+      weights.push_back(inputTractogram->GetFiberWeight(i));
+    std::sort(weights.begin(), weights.end());
+
+    inputTractogram = inputTractogram->FilterByWeights(0.0);
+    mitk::FiberBundle::Pointer pred_positives = inputTractogram->GetDeepCopy();
+    mitk::FiberBundle::Pointer pred_negatives = mitk::FiberBundle::New(nullptr);
+
+    ofstream logfile;
+    logfile.open (out_folder + "LiFE_ROC.txt");
+
+    float step = weights.back()/steps;
+    float w = 0;
+    while (pred_positives->GetNumFibers()>0)
+    {
+      w += step;
+
+      streambuf *old = cout.rdbuf(); // <-- save
+      stringstream ss;
+      std::cout.rdbuf (ss.rdbuf());       // <-- redirect
+
+      mitk::FiberBundle::Pointer tp_tracts = mitk::FiberBundle::New(nullptr);
+      mitk::FiberBundle::Pointer fn_tracts = mitk::FiberBundle::New(nullptr);
+      for ( MaskType mask : known_tract_masks )
+      {
+        ItkUcharImgType::Pointer mask_image = std::get<0>(mask);
+
+        mitk::FiberBundle::Pointer a = pred_positives->ExtractFiberSubset(mask_image, true, false, false, overlap, false);
+        tp_tracts = tp_tracts->AddBundle(a);
+
+        mitk::FiberBundle::Pointer b = pred_negatives->ExtractFiberSubset(mask_image, true, false, false, overlap, false);
+        fn_tracts = fn_tracts->AddBundle(b);
+      }
+      mitk::FiberBundle::Pointer fp_tracts = pred_positives->SubtractBundle(tp_tracts);
+      mitk::FiberBundle::Pointer tn_tracts = pred_negatives->SubtractBundle(fn_tracts);
+
+      std::cout.rdbuf (old);              // <-- restore
+
+      float positives = tp_tracts->GetNumFibers() + fn_tracts->GetNumFibers();
+      float negatives = tn_tracts->GetNumFibers() + fp_tracts->GetNumFibers();
+      float fpr = (float)fp_tracts->GetNumFibers() / negatives;
+      float tpr = (float)tp_tracts->GetNumFibers() / positives;
+
+      float accuracy = 1.0;
+      if (pred_positives->GetNumFibers()>0)
+        accuracy = (float)tp_tracts->GetNumFibers()/pred_positives->GetNumFibers();
+      logfile << w << " " << fpr << " " << tpr << " " << accuracy << " \n";
+      MITK_INFO << "#Fibers: " << pred_positives->GetNumFibers();
+      MITK_INFO << "FPR/TPR: " << fpr << "/" << tpr;
+      MITK_INFO << "Acc: " << accuracy;
+
+      pred_positives = inputTractogram->FilterByWeights(w);
+      pred_negatives = inputTractogram->FilterByWeights(w, true);
+    }
+    logfile.close();
+  }
+  catch (itk::ExceptionObject e)
+  {
+    std::cout << e;
+    return EXIT_FAILURE;
+  }
+  catch (std::exception e)
+  {
+    std::cout << e.what();
+    return EXIT_FAILURE;
+  }
+  catch (...)
+  {
+    std::cout << "ERROR!?!";
+    return EXIT_FAILURE;
+  }
+  return EXIT_SUCCESS;
+}
diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibilityFit.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibilityFit.cpp
index 73c4dbefbf..19f16dc350 100755
--- a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibilityFit.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibilityFit.cpp
@@ -1,561 +1,470 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 
 #include <mitkBaseData.h>
 #include <mitkImageCast.h>
 #include <mitkImageToItk.h>
 #include <metaCommand.h>
 #include <mitkCommandLineParser.h>
 #include <usAny.h>
 #include <mitkIOUtil.h>
 #include <itksys/SystemTools.hxx>
 #include <itkDirectory.h>
 #include <mitkFiberBundle.h>
 #include <mitkPreferenceListReaderOptionsFunctor.h>
 #include <itkImageFileWriter.h>
 #include <mitkPeakImage.h>
 #include <itkFitFibersToImageFilter.h>
 #include <boost/lexical_cast.hpp>
 #include <itkTractDensityImageFilter.h>
 
 using namespace std;
 typedef itksys::SystemTools ist;
 typedef itk::Point<float, 4> PointType4;
 typedef itk::Image< float, 4 >  PeakImgType;
 typedef itk::Image< unsigned char, 3 >  ItkUcharImageType;
 
 std::vector< string > get_file_list(const std::string& path, std::vector< string > extensions={".fib", ".trk"})
 {
   std::vector< string > file_list;
   itk::Directory::Pointer dir = itk::Directory::New();
 
   if (dir->Load(path.c_str()))
   {
     int n = dir->GetNumberOfFiles();
     for (int r = 0; r < n; r++)
     {
       const char *filename = dir->GetFile(r);
       std::string ext = ist::GetFilenameExtension(filename);
       for (auto e : extensions)
       {
         if (ext==e)
         {
           file_list.push_back(path + '/' + filename);
           break;
         }
       }
     }
   }
   return file_list;
 }
 
 std::vector< mitk::FiberBundle::Pointer > CombineTractograms(std::vector< mitk::FiberBundle::Pointer > reference, std::vector< mitk::FiberBundle::Pointer > candidates, int skip=-1)
 {
   std::vector< mitk::FiberBundle::Pointer > fib;
   for (auto f : reference)
     fib.push_back(f);
 
   int c = 0;
   for (auto f : candidates)
   {
     if (c!=skip)
       fib.push_back(f);
     ++c;
   }
 
   return fib;
 }
 
 /*!
 \brief Fits the tractogram to the input peak image by assigning a weight to each fiber (similar to https://doi.org/10.1016/j.neuroimage.2015.06.092).
 */
 int main(int argc, char* argv[])
 {
   mitkCommandLineParser parser;
 
   parser.setTitle("");
   parser.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("", "i1", mitkCommandLineParser::InputFile, "Anchor tracts:", "anchor tracts in one tractogram file", us::Any(), false);
   parser.addArgument("", "i2", mitkCommandLineParser::InputFile, "Input peaks:", "input peak image", us::Any(), false);
-  parser.addArgument("", "i3", mitkCommandLineParser::InputFile, "", "", us::Any(), false);
+  parser.addArgument("", "i3", mitkCommandLineParser::InputFile, "Candidate folder:", "folder containing all candidate tracts (separate files)", us::Any(), false);
   parser.addArgument("", "o", mitkCommandLineParser::OutputDirectory, "Output folder:", "output folder", us::Any(), false);
+  parser.addArgument("reference_mask_folder", "", mitkCommandLineParser::String, "Reference Mask Folder:", "reference masks for accuracy evaluation", false);
 
   parser.addArgument("mask", "", mitkCommandLineParser::InputFile, "", "", us::Any());
-  parser.addArgument("min_gain", "", mitkCommandLineParser::Float, "Min. gain:", "process stops if remaining bundles don't contribute enough", 0.05);
-
-  parser.addArgument("reference_mask_folder", "", mitkCommandLineParser::String, "Reference Mask Folder:", "reference masks for accuracy evaluation", false);
 
-  parser.addArgument("subtract", "", mitkCommandLineParser::Bool, "", "", false);
+  parser.addArgument("greedy_add", "", mitkCommandLineParser::Bool, "", "", false);
   parser.addArgument("normalize_gain", "", mitkCommandLineParser::Bool, "", "", 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("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 out_folder = 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 greedy_add = true;
-  if (parsedArgs.count("subtract"))
-    greedy_add = !us::any_cast<bool>(parsedArgs["subtract"]);
-
-  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"]);
+  bool greedy_add = false;
+  if (parsedArgs.count("greedy_add"))
+    greedy_add = us::any_cast<bool>(parsedArgs["greedy_add"]);
 
   float lambda = 0.1;
   if (parsedArgs.count("lambda"))
     lambda = us::any_cast<float>(parsedArgs["lambda"]);
 
   bool normalize_gain = false;
   if (parsedArgs.count("normalize_gain"))
     normalize_gain = us::any_cast<bool>(parsedArgs["normalize_gain"]);
 
   bool filter_outliers = true;
   if (parsedArgs.count("dont_filter_outliers"))
     filter_outliers = !us::any_cast<bool>(parsedArgs["dont_filter_outliers"]);
 
 
   string mask_file = "";
   if (parsedArgs.count("mask"))
     mask_file = us::any_cast<string>(parsedArgs["mask"]);
 
 
   string reference_mask_folder = "";
   if (parsedArgs.count("reference_mask_folder"))
     reference_mask_folder = us::any_cast<string>(parsedArgs["reference_mask_folder"]);
 
   try
   {
     itk::TimeProbe clock;
     clock.Start();
 
     MITK_INFO << "Creating directory structure";
     if (ist::PathExists(out_folder))
       ist::RemoveADirectory(out_folder);
     ist::MakeDirectory(out_folder);
 
     MITK_INFO << "Loading data";
 
     streambuf *old = cout.rdbuf(); // <-- save
     stringstream ss;
     std::cout.rdbuf (ss.rdbuf());       // <-- redirect
 
     ofstream logfile;
     logfile.open (out_folder + "log.txt");
 
     itk::ImageFileWriter< PeakImgType >::Pointer peak_image_writer = itk::ImageFileWriter< PeakImgType >::New();
     mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Peak Image", "Fiberbundles"}, {});
     mitk::Image::Pointer inputImage = dynamic_cast<mitk::PeakImage*>(mitk::IOUtil::Load(peak_file_name, &functor)[0].GetPointer());
 
     float minSpacing = 1;
     if(inputImage->GetGeometry()->GetSpacing()[0]<inputImage->GetGeometry()->GetSpacing()[1] && inputImage->GetGeometry()->GetSpacing()[0]<inputImage->GetGeometry()->GetSpacing()[2])
       minSpacing = inputImage->GetGeometry()->GetSpacing()[0];
     else if (inputImage->GetGeometry()->GetSpacing()[1] < inputImage->GetGeometry()->GetSpacing()[2])
       minSpacing = inputImage->GetGeometry()->GetSpacing()[1];
     else
       minSpacing = inputImage->GetGeometry()->GetSpacing()[2];
 
     // Load mask file. Fit is only performed inside the mask
     itk::FitFibersToImageFilter::UcharImgType::Pointer mask = nullptr;
     if (mask_file.compare("")!=0)
     {
       mitk::Image::Pointer mitk_mask = dynamic_cast<mitk::Image*>(mitk::IOUtil::Load(mask_file)[0].GetPointer());
       mitk::CastToItkImage(mitk_mask, mask);
     }
 
     // Load masks covering the true positives for evaluation purposes
     std::vector< string > reference_mask_filenames;
     std::vector< itk::FitFibersToImageFilter::UcharImgType::Pointer > reference_masks;
     if (reference_mask_folder.compare("")!=0)
     {
       reference_mask_filenames = get_file_list(reference_mask_folder, {".nii", ".nii.gz", ".nrrd"});
       for (auto filename : reference_mask_filenames)
       {
         itk::FitFibersToImageFilter::UcharImgType::Pointer ref_mask = nullptr;
         mitk::Image::Pointer ref_mitk_mask = dynamic_cast<mitk::Image*>(mitk::IOUtil::Load(filename)[0].GetPointer());
         mitk::CastToItkImage(ref_mitk_mask, ref_mask);
         reference_masks.push_back(ref_mask);
       }
     }
 
     // Load peak image
     typedef mitk::ImageToItk< PeakImgType > CasterType;
     CasterType::Pointer caster = CasterType::New();
     caster->SetInput(inputImage);
     caster->Update();
     PeakImgType::Pointer peak_image = caster->GetOutput();
 
     // Load reference tractogram consisting of all known tracts
     std::vector< mitk::FiberBundle::Pointer > input_reference;
     mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(mitk::IOUtil::Load(fib_file)[0].GetPointer());
     if (fib.IsNull())
       return EXIT_FAILURE;
     fib->ResampleLinear(minSpacing/10.0);
     input_reference.push_back(fib);
 
     // Load all candidate tracts
     std::vector< mitk::FiberBundle::Pointer > input_candidates;
     std::vector< string > candidate_tract_files = get_file_list(candidate_folder);
     for (string f : candidate_tract_files)
     {
       mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(mitk::IOUtil::Load(f)[0].GetPointer());
       if (fib.IsNull())
         continue;
       if (fib->GetNumFibers()<=0)
         return 0;
       fib->ResampleLinear(minSpacing/10.0);
       input_candidates.push_back(fib);
     }
     std::cout.rdbuf (old);              // <-- restore
 
     // Fit known tracts to peak image to obtain underexplained image
     MITK_INFO << "Fit anchor tracts";
     int iteration = 0;
     itk::FitFibersToImageFilter::Pointer fitter = itk::FitFibersToImageFilter::New();
-    if (greedy_add)
-      fitter->SetTractograms(input_reference);
-    else
-      fitter->SetTractograms(CombineTractograms(input_reference, input_candidates));
-    fitter->SetFitIndividualFibers(single_fib);
-    fitter->SetMaxIterations(max_iter);
-    fitter->SetGradientTolerance(g_tol);
+    fitter->SetTractograms(input_reference);
     fitter->SetLambda(lambda);
     fitter->SetFilterOutliers(filter_outliers);
     fitter->SetPeakImage(peak_image);
     fitter->SetVerbose(true);
     fitter->SetResampleFibers(false);
     fitter->SetMaskImage(mask);
     fitter->Update();
     std::string name = ist::GetFilenameWithoutExtension(fib_file);
     mitk::FiberBundle::Pointer anchor_tracts = fitter->GetTractograms().at(0);
     anchor_tracts->SetFiberColors(255,255,255);
     mitk::IOUtil::Save(anchor_tracts, out_folder + "0_" + name + ".fib");
-    if (greedy_add)
-    {
-      peak_image = fitter->GetUnderexplainedImage();
-      peak_image_writer->SetInput(peak_image);
-      peak_image_writer->SetFileName(out_folder + boost::lexical_cast<string>(iteration) + "_" + name + ".nrrd");
-      peak_image_writer->Update();
-    }
 
+    peak_image = fitter->GetUnderexplainedImage();
+    peak_image_writer->SetInput(peak_image);
+    peak_image_writer->SetFileName(out_folder + boost::lexical_cast<string>(iteration) + "_" + name + ".nrrd");
+    peak_image_writer->Update();
 
+    if (!greedy_add)
     {
       MITK_INFO << "Fit candidate tracts";
       itk::FitFibersToImageFilter::Pointer fitter = itk::FitFibersToImageFilter::New();
-      fitter->SetFitIndividualFibers(single_fib);
-      fitter->SetMaxIterations(max_iter);
-      fitter->SetGradientTolerance(g_tol);
       fitter->SetLambda(lambda);
       fitter->SetFilterOutliers(filter_outliers);
       fitter->SetVerbose(true);
       fitter->SetPeakImage(peak_image);
       fitter->SetResampleFibers(false);
       fitter->SetMaskImage(mask);
       fitter->SetTractograms(input_candidates);
       fitter->Update();
-      MITK_INFO << "DONE";
       vnl_vector<double> rms_diff = fitter->GetRmsDiffPerBundle();
 
       vnl_vector<double> log_rms_diff = rms_diff-rms_diff.min_value() + 1;
       log_rms_diff = log_rms_diff.apply(std::log);
       log_rms_diff /= log_rms_diff.max_value();
       int c = 0;
       for (auto fib : input_candidates)
       {
         fib->SetFiberWeights( log_rms_diff[c] );
         fib->ColorFibersByFiberWeights(false, true);
 
         string bundle_name = ist::GetFilenameWithoutExtension(candidate_tract_files.at(c));
         try
         {
           streambuf *old = cout.rdbuf(); // <-- save
           stringstream ss;
           std::cout.rdbuf (ss.rdbuf());       // <-- redirect
           mitk::IOUtil::Save(fib, out_folder + boost::lexical_cast<string>((int)(10000*rms_diff[c])) + "_" + bundle_name + ".fib");
           std::cout.rdbuf (old);              // <-- restore
         }
         catch(...)
         {
           MITK_INFO << fib->GetNumFibers();
           MITK_INFO << "ERROR!!!!!!!!!!!!!!";
         }
 
         streambuf *old = cout.rdbuf(); // <-- save
         stringstream ss;
         std::cout.rdbuf (ss.rdbuf());       // <-- redirect
         float best_overlap = 0;
         int best_overlap_index = -1;
         int m_idx = 0;
         for (auto ref_mask : reference_masks)
         {
           float overlap = fib->GetOverlap(ref_mask, false);
           if (overlap>best_overlap)
           {
             best_overlap = overlap;
             best_overlap_index = m_idx;
           }
           ++m_idx;
         }
 
         unsigned int num_voxels = 0;
         {
           itk::TractDensityImageFilter< ItkUcharImageType >::Pointer masks_filter = itk::TractDensityImageFilter< ItkUcharImageType >::New();
           masks_filter->SetInputImage(mask);
           masks_filter->SetBinaryOutput(true);
           masks_filter->SetFiberBundle(fib);
           masks_filter->SetUseImageGeometry(true);
           masks_filter->Update();
           num_voxels = masks_filter->GetNumCoveredVoxels();
         }
 
         std::cout.rdbuf (old);              // <-- restore
 
         logfile << "RMS_DIFF: " << setprecision(5) << rms_diff[c] << " " << bundle_name << " " << num_voxels << "\n";
         if (best_overlap_index>=0)
           logfile << "Best_overlap: " << setprecision(5) << best_overlap << " " << ist::GetFilenameWithoutExtension(reference_mask_filenames.at(best_overlap_index)) << "\n";
         else
           logfile << "No_overlap\n";
 
         ++c;
       }
 
       mitk::FiberBundle::Pointer out_fib = mitk::FiberBundle::New();
       out_fib = out_fib->AddBundles(input_candidates);
       out_fib->ColorFibersByFiberWeights(false, true);
       mitk::IOUtil::Save(out_fib, out_folder + "AllCandidates.fib");
-
-      logfile.close();
-      return EXIT_SUCCESS;
     }
+    else
+    {
 
-    double coverage = fitter->GetCoverage();
-    MITK_INFO << "Iteration: " << iteration << " - " << std::fixed << "Coverage: " << setprecision(2) << 100.0*coverage << "%";
-    logfile << "Iteration: " << iteration << " - " << std::fixed << "Coverage: " << setprecision(2) << 100.0*coverage << "%\n\n";
-    //    fitter->SetPeakImage(peak_image);
+      double coverage = fitter->GetCoverage();
+      MITK_INFO << "Iteration: " << iteration << " - " << std::fixed << "Coverage: " << setprecision(2) << 100.0*coverage << "%";
+      logfile << "Iteration: " << iteration << " - " << std::fixed << "Coverage: " << setprecision(2) << 100.0*coverage << "%\n\n";
+      //    fitter->SetPeakImage(peak_image);
 
 
-    // Iteratively add/subtract candidate bundles in a greedy manner
-    while (!input_candidates.empty())
-    {
-      double next_coverage = 0;
-      double next_covered_directions = 1.0;
-      mitk::FiberBundle::Pointer best_candidate = nullptr;
-      PeakImgType::Pointer best_candidate_peak_image = nullptr;
-      if (greedy_add)
+      // Iteratively add/subtract candidate bundles in a greedy manner
+      while (!input_candidates.empty())
       {
+        double next_coverage = 0;
+        double next_covered_directions = 1.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->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->SetResampleFibers(false);
           fitter->SetMaskImage(mask);
           // ******************************
           fitter->SetTractograms({input_candidates.at(i)});
 
           streambuf *old = cout.rdbuf(); // <-- save
           stringstream ss;
           std::cout.rdbuf (ss.rdbuf());       // <-- redirect
           fitter->Update();
           std::cout.rdbuf (old);              // <-- restore
 
           double candidate_coverage = fitter->GetCoverage();
           if (normalize_gain)
             candidate_coverage /= fitter->GetNumCoveredDirections();
 
           if (candidate_coverage>next_coverage)
           {
             next_covered_directions = fitter->GetNumCoveredDirections();
             next_coverage = candidate_coverage;
-            if ((1.0-coverage) * next_coverage >= min_gain)
-            {
-              best_candidate = fitter->GetTractograms().at(0);
-              best_candidate_peak_image = fitter->GetUnderexplainedImage();
-            }
+
+            best_candidate = fitter->GetTractograms().at(0);
+            best_candidate_peak_image = fitter->GetUnderexplainedImage();
           }
         }
-      }
-      else
-      {
-        next_coverage = 1.0;
-        for (unsigned int i=0; i<input_candidates.size(); ++i)
-        {
-          // 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->SetResampleFibers(false);
-          fitter->SetMaskImage(mask);
-          // ******************************
-          fitter->SetTractograms(CombineTractograms(input_reference, input_candidates, i));
 
-          streambuf *old = cout.rdbuf(); // <-- save
-          stringstream ss;
-          std::cout.rdbuf (ss.rdbuf());       // <-- redirect
-          fitter->Update();
-          std::cout.rdbuf (old);              // <-- restore
+        if (best_candidate.IsNull())
+          break;
 
-          double candidate_coverage = fitter->GetCoverage();
-          if (normalize_gain)
-            candidate_coverage /= fitter->GetNumCoveredDirections();
+        //      fitter->SetPeakImage(peak_image);
+        peak_image = best_candidate_peak_image;
 
-          if (candidate_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)
           {
-            //            MITK_INFO << "Coverage: " << candidate_coverage << " (" << ist::GetFilenameWithoutExtension(candidate_tract_files.at(i)) << ")";
-            next_covered_directions = fitter->GetNumCoveredDirections();
-            next_coverage = candidate_coverage;
-            //            if ((1.0-coverage) * next_coverage >= min_gain)
-            //            {
-            best_candidate = input_candidates.at(i);
-            //            }
+            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;
 
-      if (best_candidate.IsNull())
-        break;
-
-      //      fitter->SetPeakImage(peak_image);
-      peak_image = best_candidate_peak_image;
-
-      int i=0;
-      std::vector< mitk::FiberBundle::Pointer > remaining_candidates;
-      std::vector< string > remaining_candidate_files;
-      for (auto fib : input_candidates)
-      {
-        if (fib!=best_candidate)
-        {
-          remaining_candidates.push_back(fib);
-          remaining_candidate_files.push_back(candidate_tract_files.at(i));
-        }
-        else
-          name = ist::GetFilenameWithoutExtension(candidate_tract_files.at(i));
-        ++i;
-      }
-      input_candidates = remaining_candidates;
-      candidate_tract_files = remaining_candidate_files;
-
-      iteration++;
-      streambuf *old = cout.rdbuf(); // <-- save
-      stringstream ss;
-      std::cout.rdbuf (ss.rdbuf());       // <-- redirect
-      // Save winning candidate
-      mitk::IOUtil::Save(best_candidate, out_folder + boost::lexical_cast<string>(iteration) + "_" + name + ".fib");
-      if (greedy_add)
-      {
+        iteration++;
+        streambuf *old = cout.rdbuf(); // <-- save
+        stringstream ss;
+        std::cout.rdbuf (ss.rdbuf());       // <-- redirect
+        // Save winning candidate
+        mitk::IOUtil::Save(best_candidate, out_folder + boost::lexical_cast<string>(iteration) + "_" + name + ".fib");
         peak_image_writer->SetInput(peak_image);
         peak_image_writer->SetFileName(out_folder + boost::lexical_cast<string>(iteration) + "_" + name + ".nrrd");
         peak_image_writer->Update();
-      }
 
-      // Calculate best overlap with reference masks for evaluation purposes
-      float best_overlap = 0;
-      int best_overlap_index = -1;
-      i = 0;
-      for (auto ref_mask : reference_masks)
-      {
-        float overlap = best_candidate->GetOverlap(ref_mask, false);
-        if (overlap>best_overlap)
+        // 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)
         {
-          best_overlap = overlap;
-          best_overlap_index = i;
+          float overlap = best_candidate->GetOverlap(ref_mask, false);
+          if (overlap>best_overlap)
+          {
+            best_overlap = overlap;
+            best_overlap_index = i;
+          }
+          ++i;
         }
-        ++i;
-      }
-      std::cout.rdbuf (old);              // <-- restore
+        std::cout.rdbuf (old);              // <-- restore
+
+        if (normalize_gain)
+          next_coverage *= next_covered_directions;
 
-      if (normalize_gain)
-        next_coverage *= next_covered_directions;
-      if (greedy_add)
-      {
         MITK_INFO << "Iteration: " << iteration << " - " << std::fixed << "Coverage: " << setprecision(2) << 100.0*(coverage + (1.0-coverage) * next_coverage) << "% (+" << 100*(1.0-coverage) * next_coverage << "%)";
         logfile << "Iteration: " << iteration << " - " << std::fixed << "Coverage: " << setprecision(2) << 100.0*(coverage + (1.0-coverage) * next_coverage) << "% (+" << 100*(1.0-coverage) * next_coverage << "%)\n";
         coverage += (1.0-coverage) * next_coverage;
-      }
-      else
-      {
-        MITK_INFO << "Iteration: " << iteration << " - " << std::fixed << "Coverage: " << setprecision(2) << 100.0*coverage << "% (-" << 100*(coverage-next_coverage) << "%)";
-        logfile << "Iteration: " << iteration << " - " << std::fixed << "Coverage: " << setprecision(2) << 100.0*coverage << "% (-" << 100*(coverage-next_coverage) << "%)\n";
-        coverage = next_coverage;
-      }
 
-      if (best_overlap_index>=0)
-      {
-        MITK_INFO << "Best overlap: " << best_overlap << " - " << ist::GetFilenameWithoutExtension(reference_mask_filenames.at(best_overlap_index));
-        logfile << "Best overlap: " << best_overlap << " - " << ist::GetFilenameWithoutExtension(reference_mask_filenames.at(best_overlap_index)) << "\n";
+        if (best_overlap_index>=0)
+        {
+          MITK_INFO << "Best overlap: " << best_overlap << " - " << ist::GetFilenameWithoutExtension(reference_mask_filenames.at(best_overlap_index));
+          logfile << "Best overlap: " << best_overlap << " - " << ist::GetFilenameWithoutExtension(reference_mask_filenames.at(best_overlap_index)) << "\n";
+        }
+        logfile << "\n";
       }
-      logfile << "\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;
 }