diff --git a/Modules/DiffusionCore/IODataStructures/mitkFiberBundle.cpp b/Modules/DiffusionCore/IODataStructures/mitkFiberBundle.cpp
index 3c27ede..2da5848 100644
--- a/Modules/DiffusionCore/IODataStructures/mitkFiberBundle.cpp
+++ b/Modules/DiffusionCore/IODataStructures/mitkFiberBundle.cpp
@@ -1,2756 +1,2762 @@
 /*===================================================================
 
 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 "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 <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";
 
 mitk::FiberBundle::FiberBundle( vtkPolyData* fiberPolyData )
   : m_NumFibers(0)
 {
   m_TrackVisHeader.hdr_size = 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);
   newFib->SetTrackVisHeader(this->GetTrackVisHeader());
   return newFib;
 }
 
 vtkSmartPointer<vtkPolyData> mitk::FiberBundle::GeneratePolyDataByIds(std::vector<unsigned int> 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>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();
 
   auto num_weights = this->GetNumFibers();
   for (auto fib : fibs)
     num_weights += fib->GetNumFibers();
   weights->SetNumberOfValues(num_weights);
 
   unsigned int counter = 0;
   for (unsigned int i=0; i<m_FiberPolyData->GetNumberOfCells(); ++i)
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     auto numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     for (unsigned 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 (unsigned int i=0; i<fib->GetFiberPolyData()->GetNumberOfCells(); i++)
     {
       vtkCell* cell = fib->GetFiberPolyData()->GetCell(i);
       auto numPoints = cell->GetNumberOfPoints();
       vtkPoints* points = cell->GetPoints();
 
       vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
       for (unsigned int j=0; j<numPoints; j++)
       {
         double p[3];
         points->GetPoint(j, p);
 
         vtkIdType id = vNewPoints->InsertNextPoint(p);
         container->GetPointIds()->InsertNextId(id);
       }
       weights->InsertValue(counter, fib->GetFiberWeight(i));
       vNewLines->InsertNextCell(container);
       counter++;
     }
   }
 
   // initialize PolyData
   vNewPolyData->SetPoints(vNewPoints);
   vNewPolyData->SetLines(vNewLines);
 
   // initialize fiber bundle
   mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(vNewPolyData);
   newFib->SetFiberWeights(weights);
   return newFib;
 }
 
 // merge two fiber bundles
 mitk::FiberBundle::Pointer mitk::FiberBundle::AddBundle(mitk::FiberBundle* fib)
 {
   if (fib==nullptr)
     return this->GetDeepCopy();
 
   MITK_INFO << "Adding fibers";
 
   vtkSmartPointer<vtkPolyData> vNewPolyData = vtkSmartPointer<vtkPolyData>::New();
   vtkSmartPointer<vtkCellArray> vNewLines = vtkSmartPointer<vtkCellArray>::New();
   vtkSmartPointer<vtkPoints> vNewPoints = vtkSmartPointer<vtkPoints>::New();
 
   // add current fiber bundle
   vtkSmartPointer<vtkFloatArray> weights = vtkSmartPointer<vtkFloatArray>::New();
   weights->SetNumberOfValues(this->GetNumFibers()+fib->GetNumFibers());
 
   unsigned int counter = 0;
   for (unsigned int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     auto numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     for (unsigned 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 (unsigned int i=0; i<fib->GetFiberPolyData()->GetNumberOfCells(); i++)
   {
     vtkCell* cell = fib->GetFiberPolyData()->GetCell(i);
     auto numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     for (unsigned int j=0; j<numPoints; j++)
     {
       double p[3];
       points->GetPoint(j, p);
 
       vtkIdType id = vNewPoints->InsertNextPoint(p);
       container->GetPointIds()->InsertNextId(id);
     }
     weights->InsertValue(counter, fib->GetFiberWeight(i));
     vNewLines->InsertNextCell(container);
     counter++;
   }
 
   // initialize PolyData
   vNewPolyData->SetPoints(vNewPoints);
   vNewPolyData->SetLines(vNewLines);
 
   // initialize fiber bundle
   mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(vNewPolyData);
   newFib->SetFiberWeights(weights);
   return newFib;
 }
 
 // Only retain fibers with a weight larger than the specified threshold
 mitk::FiberBundle::Pointer mitk::FiberBundle::FilterByWeights(float weight_thr, bool invert)
 {
   vtkSmartPointer<vtkPolyData> vNewPolyData = vtkSmartPointer<vtkPolyData>::New();
   vtkSmartPointer<vtkCellArray> vNewLines = vtkSmartPointer<vtkCellArray>::New();
   vtkSmartPointer<vtkPoints> vNewPoints = vtkSmartPointer<vtkPoints>::New();
 
   std::vector<float> weights;
 
   for (unsigned int i=0; i<this->GetNumFibers(); i++)
   {
     if ( (invert && this->GetFiberWeight(i)>weight_thr) || (!invert && this->GetFiberWeight(i)<=weight_thr))
       continue;
 
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     auto 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));
+  newFib->SetTrackVisHeader(this->GetTrackVisHeader());
   return newFib;
 }
 
 // Only retain a subsample of the fibers
 mitk::FiberBundle::Pointer mitk::FiberBundle::SubsampleFibers(float factor, bool random_seed)
 {
   vtkSmartPointer<vtkPolyData> vNewPolyData = vtkSmartPointer<vtkPolyData>::New();
   vtkSmartPointer<vtkCellArray> vNewLines = vtkSmartPointer<vtkCellArray>::New();
   vtkSmartPointer<vtkPoints> vNewPoints = vtkSmartPointer<vtkPoints>::New();
 
   unsigned int new_num_fibs = static_cast<unsigned int>(std::round(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< unsigned int > ids;
   for (unsigned int i=0; i<this->GetNumFibers(); i++)
     ids.push_back(i);
   if (random_seed)
     std::srand(static_cast<unsigned int>(std::time(nullptr)));
   else
     std::srand(0);
   std::random_shuffle(ids.begin(), ids.end());
 
   unsigned int counter = 0;
   for (unsigned int i=0; i<new_num_fibs; i++)
   {
     vtkCell* cell = m_FiberPolyData->GetCell(ids.at(i));
     auto 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);
+  newFib->SetTrackVisHeader(this->GetTrackVisHeader());
   return newFib;
 }
 
 // subtract two fiber bundles
 mitk::FiberBundle::Pointer mitk::FiberBundle::SubtractBundle(mitk::FiberBundle* fib)
 {
   if (fib==nullptr)
     return this->GetDeepCopy();
 
   MITK_INFO << "Subtracting fibers";
   vtkSmartPointer<vtkPolyData> vNewPolyData = vtkSmartPointer<vtkPolyData>::New();
   vtkSmartPointer<vtkCellArray> vNewLines = vtkSmartPointer<vtkCellArray>::New();
   vtkSmartPointer<vtkPoints> vNewPoints = vtkSmartPointer<vtkPoints>::New();
 
   std::vector< std::vector< itk::Point<float, 3> > > points1;
   for(unsigned int i=0; i<m_NumFibers; i++ )
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     auto numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     if (points==nullptr || numPoints<=0)
       continue;
 
     itk::Point<float, 3> start = mitk::imv::GetItkPoint(points->GetPoint(0));
     itk::Point<float, 3> end = mitk::imv::GetItkPoint(points->GetPoint(numPoints-1));
 
     points1.push_back( {start, end} );
   }
 
   std::vector< std::vector< itk::Point<float, 3> > > points2;
   for(unsigned int i=0; i<fib->GetNumFibers(); i++ )
   {
     vtkCell* cell = fib->GetFiberPolyData()->GetCell(i);
     auto numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     if (points==nullptr || numPoints<=0)
       continue;
 
     itk::Point<float, 3> start =mitk::imv::GetItkPoint(points->GetPoint(0));
     itk::Point<float, 3> end =mitk::imv::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<static_cast<int>(points1.size()); i++)
   {
     bool match = false;
     for (unsigned int j=0; j<points2.size(); j++)
     {
       auto v1 = points1.at(static_cast<unsigned int>(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<0.000001f)
       {
         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<0.000001f)
       {
         match = true;
         break;
       }
     }
 
 #pragma omp critical
     if (!match)
       ids.push_back(i);
   }
 
   for( int i : ids )
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     auto numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     if (points==nullptr || numPoints<=0)
       continue;
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     for( int j=0; j<numPoints; j++)
     {
       vtkIdType id = vNewPoints->InsertNextPoint(points->GetPoint(j));
       container->GetPointIds()->InsertNextId(id);
     }
     vNewLines->InsertNextCell(container);
   }
   if(vNewLines->GetNumberOfCells()==0)
     return mitk::FiberBundle::New();
   // initialize PolyData
   vNewPolyData->SetPoints(vNewPoints);
   vNewPolyData->SetLines(vNewLines);
 
   // initialize fiber bundle
   return mitk::FiberBundle::New(vNewPolyData);
 }
 
 /*
  * 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 = static_cast<unsigned int>(m_FiberPolyData->GetNumberOfLines());
 
   if (updateGeometry)
     UpdateFiberGeometry();
   GenerateFiberIds();
   ColorFibersByOrientation();
 }
 
 /*
  * return vtkPolyData
  */
 vtkSmartPointer<vtkPolyData> mitk::FiberBundle::GetFiberPolyData() const
 {
   return m_FiberPolyData;
 }
 
 void mitk::FiberBundle::ColorFibersByLength(bool opacity, bool normalize)
 {
   if (m_MaxFiberLength<=0)
     return;
 
   auto numOfPoints = this->GetNumberOfPoints();
 
   //colors and alpha value for each single point, RGBA = 4 components
   unsigned char rgba[4] = {0,0,0,0};
   m_FiberColors = vtkSmartPointer<vtkUnsignedCharArray>::New();
   m_FiberColors->Allocate(numOfPoints * 4);
   m_FiberColors->SetNumberOfComponents(4);
   m_FiberColors->SetName("FIBER_COLORS");
 
   auto numOfFibers = m_FiberPolyData->GetNumberOfLines();
   if (numOfFibers < 1)
     return;
 
   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 int count = 0;
   for (unsigned int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     auto numPoints = cell->GetNumberOfPoints();
 
     float l = m_FiberLengths.at(i)/m_MaxFiberLength;
     if (!normalize)
     {
       l = m_FiberLengths.at(i)/255.0f;
       if (l > 1.0f)
         l = 1.0;
     }
     for (int j=0; j<numPoints; j++)
     {
       double color[3];
       lookupTable->GetColor(1.0 - static_cast<double>(l), color);
 
       rgba[0] = static_cast<unsigned char>(255.0 * color[0]);
       rgba[1] = static_cast<unsigned char>(255.0 * color[1]);
       rgba[2] = static_cast<unsigned char>(255.0 * color[2]);
       if (opacity)
         rgba[3] = static_cast<unsigned char>(255.0f * l);
       else
         rgba[3] = static_cast<unsigned char>(255.0);
       m_FiberColors->InsertTypedTuple(cell->GetPointId(j), rgba);
       count++;
     }
   }
   m_UpdateTime3D.Modified();
   m_UpdateTime2D.Modified();
 }
 
 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 = m_FiberPolyData->GetPoints();
   vtkIdType 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};
   m_FiberColors = vtkSmartPointer<vtkUnsignedCharArray>::New();
   m_FiberColors->Allocate(numOfPoints * 4);
   m_FiberColors->SetNumberOfComponents(4);
   m_FiberColors->SetName("FIBER_COLORS");
 
   auto 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] = static_cast<unsigned char>(255.0 * std::fabs(diff[0]));
           rgba[1] = static_cast<unsigned char>(255.0 * std::fabs(diff[1]));
           rgba[2] = static_cast<unsigned char>(255.0 * std::fabs(diff[2]));
           rgba[3] = static_cast<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] = static_cast<unsigned char>(255.0 * std::fabs(diff1[0]));
           rgba[1] = static_cast<unsigned char>(255.0 * std::fabs(diff1[1]));
           rgba[2] = static_cast<unsigned char>(255.0 * std::fabs(diff1[2]));
           rgba[3] = static_cast<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] = static_cast<unsigned char>(255.0 * std::fabs(diff2[0]));
           rgba[1] = static_cast<unsigned char>(255.0 * std::fabs(diff2[1]));
           rgba[2] = static_cast<unsigned char>(255.0 * std::fabs(diff2[2]));
           rgba[3] = static_cast<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};
   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);
 
   std::vector< double > values;
   double min = 1;
   double max = 0;
   MITK_INFO << "Coloring fibers by curvature";
   boost::progress_display disp(static_cast<unsigned long>(m_FiberPolyData->GetNumberOfCells()));
   for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
   {
     ++disp;
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     auto 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< double, 3 > > vectors;
       vnl_vector_fixed< double, 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< double, 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);
         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< double, 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);
         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/itk::Math::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);
     auto 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] = static_cast<unsigned char>(255.0 * color[0]);
       rgba[1] = static_cast<unsigned char>(255.0 * color[1]);
       rgba[2] = static_cast<unsigned char>(255.0 * color[2]);
       rgba[3] = static_cast<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, static_cast<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 = 999999;
   double max = -999999;
   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];
     auto pixelValue = static_cast<double>(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];
     auto pixelValue = static_cast<double>(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] = static_cast<unsigned char>(255.0 * color[0]);
     rgba[1] = static_cast<unsigned char>(255.0 * color[1]);
     rgba[2] = static_cast<unsigned char>(255.0 * color[2]);
     if (opacity)
       rgba[3] = static_cast<unsigned char>(255.0 * pixelValue);
     else
       rgba[3] = static_cast<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 (unsigned int i=0; i<m_NumFibers; i++)
   {
     float weight = this->GetFiberWeight(i);
     if (weight>max)
       max = weight;
     if (weight<min)
       min = weight;
   }
   if (fabs(max-min)<0.00001f)
   {
     max = 1;
     min = 0;
   }
 
   for (unsigned int i=0; i<m_NumFibers; i++)
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     auto numPoints = cell->GetNumberOfPoints();
     auto 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(static_cast<double>(1-v), color);
 
       rgba[0] = static_cast<unsigned char>(255.0 * color[0]);
       rgba[1] = static_cast<unsigned char>(255.0 * color[1]);
       rgba[2] = static_cast<unsigned char>(255.0 * color[2]);
       if (opacity)
         rgba[3] = static_cast<unsigned char>(255.0f * v);
       else
         rgba[3] = static_cast<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] = static_cast<unsigned char>(r);
     rgba[1] = static_cast<unsigned char>(g);
     rgba[2] = static_cast<unsigned char>(b);
     rgba[3] = static_cast<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::GetNumEpFractionInMask(ItkUcharImgType* mask, bool different_label)
 {
   vtkSmartPointer<vtkPolyData> PolyData = m_FiberPolyData;
 
   MITK_INFO << "Calculating EP-Fraction";
 
   boost::progress_display disp(m_NumFibers);
   unsigned int in_mask = 0;
 
   for (unsigned int i=0; i<m_NumFibers; i++)
   {
     ++disp;
     vtkCell* cell = PolyData->GetCell(i);
     auto numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     itk::Point<float, 3> startVertex =mitk::imv::GetItkPoint(points->GetPoint(0));
     itk::Index<3> startIndex;
     mask->TransformPhysicalPointToIndex(startVertex, startIndex);
 
     itk::Point<float, 3> endVertex =mitk::imv::GetItkPoint(points->GetPoint(numPoints-1));
     itk::Index<3> endIndex;
     mask->TransformPhysicalPointToIndex(endVertex, endIndex);
 
     if (mask->GetLargestPossibleRegion().IsInside(startIndex) && mask->GetLargestPossibleRegion().IsInside(endIndex))
     {
       float v1 = mask->GetPixel(startIndex);
       if (v1 < 0.5f)
         continue;
       float v2 = mask->GetPixel(startIndex);
       if (v2 < 0.5f)
         continue;
 
       if (!different_label)
         ++in_mask;
       else if (fabs(v1-v2)>0.00001f)
         ++in_mask;
     }
   }
   return float(in_mask)/m_NumFibers;
 }
 
 std::tuple<float, float> mitk::FiberBundle::GetDirectionalOverlap(ItkUcharImgType* mask, mitk::PeakImage::ItkPeakImageType* peak_image)
 {
   vtkSmartPointer<vtkPolyData> PolyData = m_FiberPolyData;
 
   MITK_INFO << "Calculating overlap";
   auto spacing = mask->GetSpacing();
   boost::progress_display disp(m_NumFibers);
   double length_sum = 0;
   double in_mask_length = 0;
   double aligned_length = 0;
   for (unsigned int i=0; i<m_NumFibers; i++)
   {
     ++disp;
     vtkCell* cell = PolyData->GetCell(i);
     auto numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     for (int j=0; j<numPoints-1; j++)
     {
       itk::Point<float, 3> startVertex =mitk::imv::GetItkPoint(points->GetPoint(j));
       itk::Index<3> startIndex;
       itk::ContinuousIndex<float, 3> startIndexCont;
       mask->TransformPhysicalPointToIndex(startVertex, startIndex);
       mask->TransformPhysicalPointToContinuousIndex(startVertex, startIndexCont);
 
       itk::Point<float, 3> endVertex =mitk::imv::GetItkPoint(points->GetPoint(j + 1));
       itk::Index<3> endIndex;
       itk::ContinuousIndex<float, 3> endIndexCont;
       mask->TransformPhysicalPointToIndex(endVertex, endIndex);
       mask->TransformPhysicalPointToContinuousIndex(endVertex, endIndexCont);
       
       vnl_vector_fixed< float, 3 > fdir;
       fdir[0] = endVertex[0] - startVertex[0];
       fdir[1] = endVertex[1] - startVertex[1];
       fdir[2] = endVertex[2] - startVertex[2];
       fdir.normalize();
 
       std::vector< std::pair< itk::Index<3>, double > > segments = mitk::imv::IntersectImage(spacing, startIndex, endIndex, startIndexCont, endIndexCont);
       for (std::pair< itk::Index<3>, double > segment : segments)
       {
         if ( mask->GetLargestPossibleRegion().IsInside(segment.first) && mask->GetPixel(segment.first) > 0 )
         {
           in_mask_length += segment.second;
           
           mitk::PeakImage::ItkPeakImageType::IndexType idx4; 
           idx4[0] = segment.first[0]; 
           idx4[1] = segment.first[1]; 
           idx4[2] = segment.first[2];
           
           vnl_vector_fixed< float, 3 > peak;
           idx4[3] = 0;
           peak[0] = peak_image->GetPixel(idx4);
           idx4[3] = 1;
           peak[1] = peak_image->GetPixel(idx4);
           idx4[3] = 2;
           peak[2] = peak_image->GetPixel(idx4);
           if (std::isnan(peak[0]) || std::isnan(peak[1]) || std::isnan(peak[2]) || peak.magnitude()<0.0001f)
             continue;
           peak.normalize();
           
           double f = 1.0 - std::acos(std::fabs(static_cast<double>(dot_product(fdir, peak)))) * 2.0/itk::Math::pi;
           aligned_length += segment.second * f;
         }
         length_sum += segment.second;
       }
     }
   }
 
   if (length_sum<=0.0001)
   {
     MITK_INFO << "Fiber length sum is zero!";
     return std::make_tuple(0,0);
   }
   return std::make_tuple(aligned_length/length_sum, in_mask_length/length_sum);
 }
 
 float mitk::FiberBundle::GetOverlap(ItkUcharImgType* mask)
 {
   vtkSmartPointer<vtkPolyData> PolyData = m_FiberPolyData;
 
   MITK_INFO << "Calculating overlap";
   auto spacing = mask->GetSpacing();
   boost::progress_display disp(m_NumFibers);
   double length_sum = 0;
   double in_mask_length = 0;
   for (unsigned int i=0; i<m_NumFibers; i++)
   {
     ++disp;
     vtkCell* cell = PolyData->GetCell(i);
     auto numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     for (int j=0; j<numPoints-1; j++)
     {
       itk::Point<float, 3> startVertex =mitk::imv::GetItkPoint(points->GetPoint(j));
       itk::Index<3> startIndex;
       itk::ContinuousIndex<float, 3> startIndexCont;
       mask->TransformPhysicalPointToIndex(startVertex, startIndex);
       mask->TransformPhysicalPointToContinuousIndex(startVertex, startIndexCont);
 
       itk::Point<float, 3> endVertex =mitk::imv::GetItkPoint(points->GetPoint(j + 1));
       itk::Index<3> endIndex;
       itk::ContinuousIndex<float, 3> endIndexCont;
       mask->TransformPhysicalPointToIndex(endVertex, endIndex);
       mask->TransformPhysicalPointToContinuousIndex(endVertex, endIndexCont);
 
       std::vector< std::pair< itk::Index<3>, double > > segments = mitk::imv::IntersectImage(spacing, startIndex, endIndex, startIndexCont, endIndexCont);
       for (std::pair< itk::Index<3>, double > segment : segments)
       {
         if ( mask->GetLargestPossibleRegion().IsInside(segment.first) && mask->GetPixel(segment.first) > 0 )
           in_mask_length += segment.second;
         length_sum += segment.second;
       }
     }
   }
 
   if (length_sum<=0.000001)
   {
     MITK_INFO << "Fiber length sum is zero!";
     return 0;
   }
   return static_cast<float>(in_mask_length/length_sum);
 }
 
 mitk::FiberBundle::Pointer mitk::FiberBundle::RemoveFibersOutside(ItkUcharImgType* mask, bool invert)
 {
   vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
   vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
 
   std::vector< float > fib_weights;
 
   MITK_INFO << "Cutting fibers";
   boost::progress_display disp(m_NumFibers);
   for (unsigned int i=0; i<m_NumFibers; i++)
   {
     ++disp;
 
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     auto numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     int newNumPoints = 0;
     if (numPoints>1)
     {
       for (int j=0; j<numPoints; j++)
       {
         itk::Point<float, 3> itkP =mitk::imv::GetItkPoint(points->GetPoint(j));
         itk::Index<3> idx;
         mask->TransformPhysicalPointToIndex(itkP, idx);
 
         bool inside = false;
         if ( mask->GetLargestPossibleRegion().IsInside(idx) && mask->GetPixel(idx)!=0 )
           inside = true;
 
         if (inside && !invert)
         {
           vtkIdType id = vtkNewPoints->InsertNextPoint(itkP.GetDataPointer());
           container->GetPointIds()->InsertNextId(id);
           newNumPoints++;
         }
         else if ( !inside && invert )
         {
           vtkIdType id = vtkNewPoints->InsertNextPoint(itkP.GetDataPointer());
           container->GetPointIds()->InsertNextId(id);
           newNumPoints++;
         }
         else if (newNumPoints>1)
         {
           fib_weights.push_back(this->GetFiberWeight(i));
           vtkNewCells->InsertNextCell(container);
           newNumPoints = 0;
           container = vtkSmartPointer<vtkPolyLine>::New();
         }
         else
         {
           newNumPoints = 0;
           container = vtkSmartPointer<vtkPolyLine>::New();
         }
       }
 
       if (newNumPoints>1)
       {
         fib_weights.push_back(this->GetFiberWeight(i));
         vtkNewCells->InsertNextCell(container);
       }
     }
 
   }
 
   vtkSmartPointer<vtkFloatArray> newFiberWeights = vtkSmartPointer<vtkFloatArray>::New();
   newFiberWeights->SetName("FIBER_WEIGHTS");
   newFiberWeights->SetNumberOfValues(static_cast<vtkIdType>(fib_weights.size()));
 
   if (vtkNewCells->GetNumberOfCells()<=0)
     return nullptr;
 
   for (unsigned int i=0; i<newFiberWeights->GetNumberOfValues(); i++)
     newFiberWeights->SetValue(i, fib_weights.at(i));
 
 //  vtkSmartPointer<vtkUnsignedCharArray> newFiberColors = vtkSmartPointer<vtkUnsignedCharArray>::New();
 //  newFiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * 4);
 //  newFiberColors->SetNumberOfComponents(4);
 //  newFiberColors->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);
 //  }
 
   vtkSmartPointer<vtkPolyData> newPolyData = vtkSmartPointer<vtkPolyData>::New();
   newPolyData->SetPoints(vtkNewPoints);
   newPolyData->SetLines(vtkNewCells);
   mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(newPolyData);
   newFib->SetFiberWeights(newFiberWeights);
 //  newFib->Compress(0.1);
+  newFib->SetTrackVisHeader(this->GetTrackVisHeader());
   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<unsigned int> 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);
+  fib->SetTrackVisHeader(this->GetTrackVisHeader());
   return fib;
 }
 
 std::vector<unsigned int> mitk::FiberBundle::ExtractFiberIdSubset(DataNode *roi, DataStorage* storage)
 {
   std::vector<unsigned int> 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<unsigned int>::iterator it;
       for (unsigned int i=1; i<children->Size(); ++i)
       {
         std::vector<unsigned int> inRoi = this->ExtractFiberIdSubset(children->ElementAt(i), storage);
 
         std::vector<unsigned int> rest(std::min(result.size(),inRoi.size()));
         it = std::set_intersection(result.begin(), result.end(), inRoi.begin(), inRoi.end(), rest.begin() );
         rest.resize( static_cast<unsigned int>(it - rest.begin()) );
         result = rest;
       }
       break;
     }
     case 1: // OR
     {
       MITK_INFO << "OR";
       result = ExtractFiberIdSubset(children->ElementAt(0), storage);
       std::vector<unsigned int>::iterator it;
       for (unsigned int i=1; i<children->Size(); ++i)
       {
         it = result.end();
         std::vector<unsigned int> 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( static_cast<unsigned int>(it - result.begin()) );
       break;
     }
     case 2: // NOT
     {
       MITK_INFO << "NOT";
       for(unsigned int i=0; i<this->GetNumFibers(); i++)
         result.push_back(i);
 
       std::vector<unsigned int>::iterator it;
       for (unsigned int i=0; i<children->Size(); ++i)
       {
         std::vector<unsigned int> inRoi = ExtractFiberIdSubset(children->ElementAt(i), storage);
 
         std::vector<unsigned int> rest(result.size()-inRoi.size());
         it = std::set_difference(result.begin(), result.end(), inRoi.begin(), inRoi.end(), rest.begin() );
         rest.resize( static_cast<unsigned int>(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 (unsigned int i=0; i<m_NumFibers; i++)
       {
         ++disp ;
         vtkCell* cell = m_FiberPolyData->GetCell(i);
         auto 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 (unsigned int i=0; i<m_NumFibers; i++)
       {
         ++disp ;
         vtkCell* cell = m_FiberPolyData->GetCell(i);
         auto 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 = static_cast<unsigned int>(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);
     auto 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);
 
       double 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 += static_cast<float>(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 (unsigned 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);
 
+  GetTrackVisHeader();
+
   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*itk::Math::pi/180;
   ry = ry*itk::Math::pi/180;
   rz = rz*itk::Math::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;
 }
 
 void mitk::FiberBundle::TransformFibers(itk::ScalableAffineTransform< mitk::ScalarType >::Pointer transform)
 {
   vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
   vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
 
   for (unsigned int i=0; i<m_NumFibers; i++)
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     auto numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     for (int j=0; j<numPoints; j++)
     {
       itk::Point<float, 3> p =mitk::imv::GetItkPoint(points->GetPoint(j));
       p = transform->TransformPoint(p);
       vtkIdType id = vtkNewPoints->InsertNextPoint(p.GetDataPointer());
       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::TransformFibers(double rx, double ry, double rz, double tx, double ty, double tz)
 {
   vnl_matrix_fixed< double, 3, 3 > rot = mitk::imv::GetRotationMatrixVnl(rx, ry, rz);
 
   mitk::BaseGeometry::Pointer geom = this->GetGeometry();
   mitk::Point3D center = geom->GetCenter();
 
   vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
   vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
 
   for (unsigned int i=0; i<m_NumFibers; i++)
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     auto 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*itk::Math::pi/180;
   y = y*itk::Math::pi/180;
   z = z*itk::Math::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 (unsigned int i=0; i<m_NumFibers; i++)
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     auto 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 (unsigned int i=0; i<m_NumFibers; i++)
   {
     ++disp ;
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     auto 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 (unsigned int i=0; i<m_NumFibers; i++)
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     auto 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 (unsigned int i=0; i<m_NumFibers; i++)
   {
     ++disp;
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     auto 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(static_cast<unsigned long>(m_FiberPolyData->GetNumberOfCells()));
   for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
   {
     ++disp ;
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     auto 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(static_cast<unsigned long>(m_FiberPolyData->GetNumberOfCells()));
   for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
   {
     ++disp ;
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     auto 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] = static_cast<float>(p2[0]-p1[0]);
       v1[1] = static_cast<float>(p2[1]-p1[1]);
       v1[2] = static_cast<float>(p2[2]-p1[2]);
 
       v2[0] = static_cast<float>(p3[0]-p2[0]);
       v2[1] = static_cast<float>(p3[1]-p2[1]);
       v2[2] = static_cast<float>(p3[2]-p2[2]);
 
       v3[0] = static_cast<float>(p1[0]-p3[0]);
       v3[1] = static_cast<float>(p1[1]-p3[1]);
       v3[2] = static_cast<float>(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 (unsigned int i=0; i<m_NumFibers; i++)
   {
     ++disp;
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     auto 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 (unsigned int i=0; i<m_NumFibers; i++)
   {
     ++disp;
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     auto 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
 
   MITK_INFO << "Smoothing fibers";
   vtkSmartPointer<vtkFloatArray> newFiberWeights = vtkSmartPointer<vtkFloatArray>::New();
   newFiberWeights->SetName("FIBER_WEIGHTS");
   newFiberWeights->SetNumberOfValues(m_NumFibers);
 
   std::vector< vtkSmartPointer<vtkPolyLine> > resampled_streamlines;
   resampled_streamlines.resize(m_NumFibers);
 
   boost::progress_display disp(m_NumFibers);
 #pragma omp parallel for
   for (int i=0; i<static_cast<int>(m_NumFibers); i++)
   {
     vtkSmartPointer<vtkPoints> newPoints = vtkSmartPointer<vtkPoints>::New();
     float length = 0;
 #pragma omp critical
     {
       length = m_FiberLengths.at(static_cast<unsigned int>(i));
       ++disp;
       vtkCell* cell = m_FiberPolyData->GetCell(i);
       auto numPoints = cell->GetNumberOfPoints();
       vtkPoints* points = cell->GetPoints();
       for (int j=0; j<numPoints; j++)
         newPoints->InsertNextPoint(points->GetPoint(j));
     }
 
     int sampling = static_cast<int>(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();
 
 #pragma omp critical
     {
       for (int j=0; j<tmpSmoothPnts->GetNumberOfPoints(); j++)
       {
         vtkIdType id = vtkSmoothPoints->InsertNextPoint(tmpSmoothPnts->GetPoint(j));
         smoothLine->GetPointIds()->InsertNextId(id);
       }
 
       resampled_streamlines[static_cast<unsigned long>(i)] = smoothLine;
     }
   }
 
   for (auto container : resampled_streamlines)
   {
     vtkSmoothCells->InsertNextCell(container);
   }
 
   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 int mitk::FiberBundle::GetNumberOfPoints() const
 {
   unsigned int 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 int numRemovedPoints = 0;
   boost::progress_display disp(static_cast<unsigned long>(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<static_cast<int>(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);
       auto 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
     auto numPoints = vertices.size();
     std::vector< int > removedPoints; removedPoints.resize(numPoints, 0);
     removedPoints[0]=-1; removedPoints[numPoints-1]=-1;
 
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     unsigned int remCounter = 0;
 
     bool pointFound = true;
     while (pointFound)
     {
       pointFound = false;
       double minError = static_cast<double>(error);
       unsigned int removeIndex = 0;
 
       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=static_cast<int>(j)-1; k>=0; k--)
             if (removedPoints[static_cast<unsigned int>(k)]<=0)
             {
               pred = vertices.at(static_cast<unsigned int>(k));
               validP = k;
               break;
             }
           int validS = -1;
           vnl_vector_fixed< double, 3 > succ;
           for (unsigned int k=j+1; k<numPoints; k++)
             if (removedPoints[k]<=0)
             {
               succ = vertices.at(k);
               validS = static_cast<int>(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 (unsigned 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::ResampleToNumPoints(unsigned int targetPoints)
 {
   if (targetPoints<2)
     mitkThrow() << "Minimum two points required for resampling!";
 
   MITK_INFO << "Resampling fibers (number of points " << targetPoints << ")";
 
   bool unequal_fibs = true;
   while (unequal_fibs)
   {
     vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
     vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
 
     vtkSmartPointer<vtkFloatArray> newFiberWeights = vtkSmartPointer<vtkFloatArray>::New();
     newFiberWeights->SetName("FIBER_WEIGHTS");
     newFiberWeights->SetNumberOfValues(m_NumFibers);
 
     unequal_fibs = false;
     for (unsigned int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
     {
 
       std::vector< vnl_vector_fixed< double, 3 > > vertices;
       float weight = 1;
       double seg_len = 0;
 
       {
         weight = m_FiberWeights->GetValue(i);
         vtkCell* cell = m_FiberPolyData->GetCell(i);
         auto numPoints = cell->GetNumberOfPoints();
         if (numPoints!=targetPoints)
           seg_len = static_cast<double>(this->GetFiberLength(i)/(targetPoints-1));
         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);
 
       {
         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 >= seg_len && seg_len>0)
         {
           vnl_vector_fixed< double, 3 > newV = lastV;
           if ( new_dist-seg_len <= mitk::eps )
           {
             vec.normalize();
             newV += vec * seg_len;
           }
           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]) - seg_len*seg_len;
 
             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) || seg_len<=0.0000001)
         {
 //#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 (container->GetNumberOfPoints()!=targetPoints)
           unequal_fibs = true;
       }
     }
 
     if (vtkNewCells->GetNumberOfCells()>0)
     {
       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(static_cast<unsigned long>(m_FiberPolyData->GetNumberOfCells()));
   vtkSmartPointer<vtkFloatArray> newFiberWeights = vtkSmartPointer<vtkFloatArray>::New();
   newFiberWeights->SetName("FIBER_WEIGHTS");
   newFiberWeights->SetNumberOfValues(m_NumFibers);
 
   std::vector< vtkSmartPointer<vtkPolyLine> > resampled_streamlines;
   resampled_streamlines.resize(static_cast<unsigned long>(m_FiberPolyData->GetNumberOfCells()));
 
 #pragma omp parallel for
   for (int i=0; i<static_cast<int>(m_FiberPolyData->GetNumberOfCells()); i++)
   {
 
     std::vector< vnl_vector_fixed< double, 3 > > vertices;
 
 #pragma omp critical
     {
       ++disp;
       vtkCell* cell = m_FiberPolyData->GetCell(i);
       auto 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
     {
       resampled_streamlines[static_cast<unsigned int>(i)] = container;
     }
   }
 
   for (auto container : resampled_streamlines)
   {
     vtkNewCells->InsertNextCell(container);
   }
 
   if (vtkNewCells->GetNumberOfCells()>0)
   {
     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 (unsigned int i=0; i<m_NumFibers; i++)
   {
     vtkCell* cell = m_FiberPolyData->GetCell(i);
     auto numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     vtkCell* cell2 = fib->GetFiberPolyData()->GetCell(i);
     auto 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 << 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;
 
   os << indent << "Extent x: " << this->GetGeometry()->GetExtentInMM(0) << "mm" << std::endl;
   os << indent << "Extent y: " << this->GetGeometry()->GetExtentInMM(1) << "mm" << std::endl;
   os << indent << "Extent z: " << this->GetGeometry()->GetExtentInMM(2) << "mm" << std::endl;
   os << indent << "Diagonal: " << this->GetGeometry()->GetDiagonalLength()  << "mm" << std::endl;
 
   os << "\nReference geometry:" << std::endl;
   os << indent << "Size: [" << std::defaultfloat << m_TrackVisHeader.dim[0] << " " << m_TrackVisHeader.dim[1] << " " << m_TrackVisHeader.dim[2] << "]" << std::endl;
   os << indent << "Voxel size: [" << m_TrackVisHeader.voxel_size[0] << " " << m_TrackVisHeader.voxel_size[1] << " " << m_TrackVisHeader.voxel_size[2] << "]" << std::endl;
   os << indent << "Origin: [" << m_TrackVisHeader.origin[0] << " " << m_TrackVisHeader.origin[1] << " " << m_TrackVisHeader.origin[2] << "]" << std::endl;
   os << indent << "Matrix: " << std::scientific << std::endl;
   os << indent << "[[" << m_TrackVisHeader.vox_to_ras[0][0] << ", " << m_TrackVisHeader.vox_to_ras[0][1] << ", " << m_TrackVisHeader.vox_to_ras[0][2] << ", " << m_TrackVisHeader.vox_to_ras[0][3] << "]" << std::endl;
   os << indent << " [" << m_TrackVisHeader.vox_to_ras[1][0] << ", " << m_TrackVisHeader.vox_to_ras[1][1] << ", " << m_TrackVisHeader.vox_to_ras[1][2] << ", " << m_TrackVisHeader.vox_to_ras[1][3] << "]" << std::endl;
   os << indent << " [" << m_TrackVisHeader.vox_to_ras[2][0] << ", " << m_TrackVisHeader.vox_to_ras[2][1] << ", " << m_TrackVisHeader.vox_to_ras[2][2] << ", " << m_TrackVisHeader.vox_to_ras[2][3] << "]" << std::endl;
   os << indent << " [" << m_TrackVisHeader.vox_to_ras[3][0] << ", " << m_TrackVisHeader.vox_to_ras[3][1] << ", " << m_TrackVisHeader.vox_to_ras[3][2] << ", " << m_TrackVisHeader.vox_to_ras[3][3] << "]]" << std::defaultfloat << std::endl;
 
   if (m_FiberWeights!=nullptr)
   {
     std::vector< float > weights;
     for (int i=0; i<m_FiberWeights->GetSize(); i++)
       weights.push_back(m_FiberWeights->GetValue(i));
 
     std::sort(weights.begin(), weights.end());
 
     os << "\nFiber weight statistics" << std::endl;
     os << indent << "Min: " << weights.front() << std::endl;
     os << indent << "1% quantile: " << weights.at(static_cast<unsigned long>(weights.size()*0.01)) << std::endl;
     os << indent << "5% quantile: " << weights.at(static_cast<unsigned long>(weights.size()*0.05)) << std::endl;
     os << indent << "25% quantile: " << weights.at(static_cast<unsigned long>(weights.size()*0.25)) << std::endl;
     os << indent << "Median: " << weights.at(static_cast<unsigned long>(weights.size()*0.5)) << std::endl;
     os << indent << "75% quantile: " << weights.at(static_cast<unsigned long>(weights.size()*0.75)) << std::endl;
     os << indent << "95% quantile: " << weights.at(static_cast<unsigned long>(weights.size()*0.95)) << std::endl;
     os << indent << "99% quantile: " << weights.at(static_cast<unsigned long>(weights.size()*0.99)) << std::endl;
     os << indent << "Max: " << weights.back() << std::endl;
   }
   else
     os << indent << "\n\nNo fiber weight array found." << std::endl;
 
   Superclass::PrintSelf(os, indent);
 }
 
 mitk::FiberBundle::TrackVis_header mitk::FiberBundle::GetTrackVisHeader()
 {
   if (m_TrackVisHeader.hdr_size==0)
   {
     mitk::Geometry3D::Pointer geom = dynamic_cast<mitk::Geometry3D*>(this->GetGeometry());
     SetTrackVisHeader(geom);
   }
   return m_TrackVisHeader;
 }
 
 void mitk::FiberBundle::SetTrackVisHeader(const mitk::FiberBundle::TrackVis_header &TrackVisHeader)
 {
   m_TrackVisHeader = TrackVisHeader;
 }
 
 void mitk::FiberBundle::SetTrackVisHeader(mitk::BaseGeometry* geometry)
 {
   vtkSmartPointer< vtkMatrix4x4 > matrix = vtkSmartPointer< vtkMatrix4x4 >::New();
   matrix->Identity();
 
   if (geometry==nullptr)
     return;
 
   for(int i=0; i<3 ;i++)
   {
     m_TrackVisHeader.dim[i]            = geometry->GetExtent(i);
     m_TrackVisHeader.voxel_size[i]     = geometry->GetSpacing()[i];
     m_TrackVisHeader.origin[i]         = geometry->GetOrigin()[i];
     matrix = geometry->GetVtkMatrix();
   }
 
   for (int i=0; i<4; ++i)
     for (int j=0; j<4; ++j)
       m_TrackVisHeader.vox_to_ras[i][j] = matrix->GetElement(i, j);
 
   m_TrackVisHeader.n_scalars = 0;
   m_TrackVisHeader.n_properties = 0;
   sprintf(m_TrackVisHeader.voxel_order,"LPS");
   m_TrackVisHeader.image_orientation_patient[0] = 1.0;
   m_TrackVisHeader.image_orientation_patient[1] = 0.0;
   m_TrackVisHeader.image_orientation_patient[2] = 0.0;
   m_TrackVisHeader.image_orientation_patient[3] = 0.0;
   m_TrackVisHeader.image_orientation_patient[4] = 1.0;
   m_TrackVisHeader.image_orientation_patient[5] = 0.0;
   m_TrackVisHeader.pad1[0] = 0;
   m_TrackVisHeader.pad1[1] = 0;
   m_TrackVisHeader.pad2[0] = 0;
   m_TrackVisHeader.pad2[1] = 0;
   m_TrackVisHeader.invert_x = 0;
   m_TrackVisHeader.invert_y = 0;
   m_TrackVisHeader.invert_z = 0;
   m_TrackVisHeader.swap_xy = 0;
   m_TrackVisHeader.swap_yz = 0;
   m_TrackVisHeader.swap_zx = 0;
   m_TrackVisHeader.n_count = 0;
   m_TrackVisHeader.version = 2;
   m_TrackVisHeader.hdr_size = 1000;
   std::string id = "TRACK";
   strcpy(m_TrackVisHeader.id_string, id.c_str());
 }
 
 /* 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/FiberTracking/Algorithms/itkFiberCurvatureFilter.cpp b/Modules/FiberTracking/Algorithms/itkFiberCurvatureFilter.cpp
index 6ed4255..d4e2777 100644
--- a/Modules/FiberTracking/Algorithms/itkFiberCurvatureFilter.cpp
+++ b/Modules/FiberTracking/Algorithms/itkFiberCurvatureFilter.cpp
@@ -1,166 +1,167 @@
 /*===================================================================
 
 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 "itkFiberCurvatureFilter.h"
 
 #include <vtkDoubleArray.h>
 #include <vtkPointData.h>
 #include <boost/progress.hpp>
 
 namespace itk{
 
 FiberCurvatureFilter::FiberCurvatureFilter()
     : m_AngularDeviation(30)
     , m_Distance(10.0)
     , m_RemoveFibers(false)
     , m_UseMedian(false)
 {
 
 }
 
 FiberCurvatureFilter::~FiberCurvatureFilter()
 {
 
 }
 
 void FiberCurvatureFilter::GenerateData()
 {
     vtkSmartPointer<vtkPolyData> inputPoly = m_InputFiberBundle->GetFiberPolyData();
     vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
     vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
 
     MITK_INFO << "Applying curvature threshold";
     boost::progress_display disp(inputPoly->GetNumberOfCells());
 #pragma omp parallel for
     for (int i=0; i<inputPoly->GetNumberOfCells(); i++)
     {
         std::vector< vnl_vector_fixed< double, 3 > > vertices;
 
 #pragma omp critical
         {
             ++disp;
             vtkCell* cell = inputPoly->GetCell(i);
             int numPoints = cell->GetNumberOfPoints();
             vtkPoints* points = cell->GetPoints();
 
             for (int j=0; j<numPoints; j++)
             {
                 double p[3];
                 points->GetPoint(j, p);
                 vnl_vector_fixed< double, 3 > p_vec;
                 p_vec[0]=p[0]; p_vec[1]=p[1]; p_vec[2]=p[2];
                 vertices.push_back(p_vec);
             }
         }
 
         // calculate curvatures
         int numPoints = vertices.size();
         vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
         for (int j=0; j<numPoints; j++)
         {
             double dist = 0;
             int c = j;
             std::vector< vnl_vector_fixed< double, 3 > > vectors;
             vnl_vector_fixed< double, 3 > meanV; meanV.fill(0.0);
             while(dist<m_Distance/2 && c>1)
             {
               vnl_vector_fixed< double, 3 > p1 = vertices.at(c-1);
               vnl_vector_fixed< double, 3 > p2 = vertices.at(c);
 
                 vnl_vector_fixed< double, 3 > v = p2-p1;
                 dist += v.magnitude();
                 v.normalize();
                 vectors.push_back(v);
                 if (m_UseMedian && c==j)
                     meanV += v;
                 else if (!m_UseMedian)
                   meanV += v;
                 c--;
             }
             c = j;
             dist = 0;
             while(dist<m_Distance/2 && c<numPoints-1)
             {
               vnl_vector_fixed< double, 3 > p1 = vertices.at(c);
               vnl_vector_fixed< double, 3 > p2 = vertices.at(c+1);
 
                 vnl_vector_fixed< double, 3 > v = p2-p1;
                 dist += v.magnitude();
                 v.normalize();
                 vectors.push_back(v);
                 if (m_UseMedian && c==j)
                     meanV += v;
                 else if (!m_UseMedian)
                   meanV += v;
                 c++;
             }
             meanV.normalize();
 
             double dev = 0;
             for (auto vec : vectors)
             {
                 double angle = dot_product(meanV, vec);
                 if (angle>1.0)
                     angle = 1.0;
                 if (angle<-1.0)
                     angle = -1.0;
                 dev += acos(angle)*180/itk::Math::pi;
             }
             if (vectors.size()>0)
                 dev /= vectors.size();
 
             if (dev<m_AngularDeviation)
             {
 #pragma omp critical
               {
                 vtkIdType id = vtkNewPoints->InsertNextPoint(vertices.at(j).data_block());
                 container->GetPointIds()->InsertNextId(id);
               }
             }
             else
             {
                 if (m_RemoveFibers)
                 {
                     container = vtkSmartPointer<vtkPolyLine>::New();
                     break;
                 }
 
                 if (container->GetNumberOfPoints()>0)
                 {
 #pragma omp critical
                     vtkNewCells->InsertNextCell(container);
                 }
                 container = vtkSmartPointer<vtkPolyLine>::New();
             }
         }
 
 #pragma omp critical
         {
           if (container->GetNumberOfPoints()>0)
               vtkNewCells->InsertNextCell(container);
         }
     }
 
     vtkSmartPointer<vtkPolyData> outputPoly = vtkSmartPointer<vtkPolyData>::New();
     outputPoly->SetPoints(vtkNewPoints);
     outputPoly->SetLines(vtkNewCells);
     m_OutputFiberBundle = mitk::FiberBundle::New(outputPoly);
+    m_OutputFiberBundle->SetTrackVisHeader(m_InputFiberBundle->GetTrackVisHeader());
 }
 
 }
 
 
 
 
diff --git a/Modules/FiberTracking/Algorithms/itkFiberExtractionFilter.cpp b/Modules/FiberTracking/Algorithms/itkFiberExtractionFilter.cpp
index 50ee12d..f89782f 100644
--- a/Modules/FiberTracking/Algorithms/itkFiberExtractionFilter.cpp
+++ b/Modules/FiberTracking/Algorithms/itkFiberExtractionFilter.cpp
@@ -1,500 +1,505 @@
 /*===================================================================
 
 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 __itkFiberExtractionFilter_cpp
 #define __itkFiberExtractionFilter_cpp
 
 #include "itkFiberExtractionFilter.h"
 
 #define _USE_MATH_DEFINES
 #include <cmath>
 #include <boost/progress.hpp>
 #include <mitkDiffusionFunctionCollection.h>
 #include <boost/lexical_cast.hpp>
 
 namespace itk{
 
 
 template< class PixelType >
 FiberExtractionFilter< PixelType >::FiberExtractionFilter()
   : m_DontResampleFibers(false)
   , m_Mode(MODE::OVERLAP)
   , m_InputType(INPUT::SCALAR_MAP)
   , m_BothEnds(true)
   , m_OverlapFraction(0.8)
   , m_NoNegatives(false)
   , m_NoPositives(false)
   , m_Interpolate(false)
   , m_Threshold(0.5)
   , m_Labels()
   , m_SkipSelfConnections(false)
   , m_OnlySelfConnections(false)
   , m_SplitByRoi(false)
   , m_SplitLabels(false)
   , m_MinFibersPerTract(0)
   , m_PairedStartEndLabels(false)
 {
   m_Interpolator = itk::LinearInterpolateImageFunction< itk::Image< PixelType, 3 >, float >::New();
 }
 
 template< class PixelType >
 FiberExtractionFilter< PixelType >::~FiberExtractionFilter()
 {
 
 }
 
 template< class PixelType >
 void FiberExtractionFilter< PixelType >::SetRoiImages(const std::vector<ItkInputImgType*> &rois)
 {
   for (auto roi : rois)
   {
     if (roi==nullptr)
     {
       MITK_INFO << "ROI image is NULL!";
       return;
     }
   }
   m_RoiImages = rois;
 }
 
 template< class PixelType >
 void FiberExtractionFilter< PixelType >::SetRoiImageNames(const std::vector< std::string > roi_names)
 {
   m_RoiImageNames = roi_names;
 }
 
 template< class PixelType >
 mitk::FiberBundle::Pointer FiberExtractionFilter< PixelType >::CreateFib(std::vector< unsigned int >& ids)
 {
   vtkSmartPointer<vtkFloatArray> weights = vtkSmartPointer<vtkFloatArray>::New();
   vtkSmartPointer<vtkPolyData> pTmp = m_InputFiberBundle->GeneratePolyDataByIds(ids, weights);
   mitk::FiberBundle::Pointer fib = mitk::FiberBundle::New(pTmp);
   fib->SetFiberWeights(weights);
   return fib;
 }
 
 template< class PixelType >
 bool FiberExtractionFilter< PixelType >::IsPositive(const itk::Point<float, 3>& itkP)
 {
   if( m_InputType == INPUT::SCALAR_MAP )
     return mitk::imv::IsInsideMask<PixelType>(itkP, m_Interpolate, m_Interpolator, m_Threshold);
   else if( m_InputType == INPUT::LABEL_MAP )
   {
     auto val = mitk::imv::GetImageValue<PixelType>(itkP, false, m_Interpolator);
     for (auto l : m_Labels)
       if (l==val)
         return true;
   }
   else
     mitkThrow() << "No valid input type selected!";
   return false;
 }
 
 template< class PixelType >
 std::vector< std::string > FiberExtractionFilter< PixelType >::GetPositiveLabels() const
 {
   return m_PositiveLabels;
 }
 
 template< class PixelType >
 void FiberExtractionFilter< PixelType >::ExtractOverlap(mitk::FiberBundle::Pointer fib)
 {
   MITK_INFO << "Extracting fibers (min. overlap " << m_OverlapFraction << ")";
   vtkSmartPointer<vtkPolyData> polydata = fib->GetFiberPolyData();
 
   std::vector< std::vector< unsigned int > > positive_ids;  // one ID vector per ROI
   positive_ids.resize(m_RoiImages.size());
 
   std::vector< unsigned int > negative_ids; // fibers not overlapping with ANY mask
 
   boost::progress_display disp(m_InputFiberBundle->GetNumFibers());
   for (unsigned int i=0; i<m_InputFiberBundle->GetNumFibers(); i++)
   {
     ++disp;
     vtkCell* cell = polydata->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     float best_ol = 0;
     int best_ol_idx = -1;
     for (unsigned int m=0; m<m_RoiImages.size(); ++m)
     {
       auto roi = m_RoiImages.at(m);
       m_Interpolator->SetInputImage(roi);
       PixelType inside = 0;
       PixelType outside = 0;
       for (int j=0; j<numPoints-1; j++)
       {
         itk::Point<float, 3> startVertex = mitk::imv::GetItkPoint(points->GetPoint(j));
         itk::Index<3> startIndex;
         itk::ContinuousIndex<float, 3> startIndexCont;
         roi->TransformPhysicalPointToIndex(startVertex, startIndex);
         roi->TransformPhysicalPointToContinuousIndex(startVertex, startIndexCont);
 
         itk::Point<float, 3> endVertex = mitk::imv::GetItkPoint(points->GetPoint(j + 1));
         itk::Index<3> endIndex;
         itk::ContinuousIndex<float, 3> endIndexCont;
         roi->TransformPhysicalPointToIndex(endVertex, endIndex);
         roi->TransformPhysicalPointToContinuousIndex(endVertex, endIndexCont);
 
         std::vector< std::pair< itk::Index<3>, double > > segments = mitk::imv::IntersectImage(roi->GetSpacing(), startIndex, endIndex, startIndexCont, endIndexCont);
         for (std::pair< itk::Index<3>, double > segment : segments)
         {
           if (!roi->GetLargestPossibleRegion().IsInside(segment.first))
             continue;
           if (roi->GetPixel(segment.first)>=m_Threshold)
             inside += segment.second;
           else
             outside += segment.second;
         }
       }
 
       float overlap = (float)inside/(inside+outside);
       if (overlap > best_ol && overlap >= m_OverlapFraction)
       {
         best_ol_idx = m;
         best_ol = overlap;
       }
     }
     if (best_ol_idx<0)
       negative_ids.push_back(i);
     else
       positive_ids.at(best_ol_idx).push_back(i);
   }
 
   if (!m_NoNegatives)
     m_Negatives.push_back(CreateFib(negative_ids));
   if (!m_NoPositives)
   {
     for (unsigned int i=0; i<positive_ids.size(); ++i)
     {
       std::string name;
       if (i<m_RoiImageNames.size())
         name = m_RoiImageNames.at(i);
       else
         name = "Roi" + boost::lexical_cast<std::string>(i);
 
       if (positive_ids.at(i).size()>=m_MinFibersPerTract)
       {
         m_Positives.push_back(CreateFib(positive_ids.at(i)));
         m_PositiveLabels.push_back(name);
       }
     }
     if (!m_SplitByRoi)
     {
       mitk::FiberBundle::Pointer output = mitk::FiberBundle::New(nullptr);
       output = output->AddBundles(m_Positives);
       m_Positives.clear();
       m_Positives.push_back(output);
       m_PositiveLabels.clear();
       m_PositiveLabels.push_back("");
     }
   }
 }
 
 template< class PixelType >
 void FiberExtractionFilter< PixelType >::ExtractEndpoints(mitk::FiberBundle::Pointer fib)
 {
   MITK_INFO << "Extracting fibers (endpoints in mask)";
   vtkSmartPointer<vtkPolyData> polydata = fib->GetFiberPolyData();
 
   std::vector< std::vector< unsigned int > > positive_ids;  // one ID vector per ROI
   positive_ids.resize(m_RoiImages.size());
 
   std::vector< unsigned int > negative_ids; // fibers not overlapping with ANY mask
 
   boost::progress_display disp(m_InputFiberBundle->GetNumFibers());
   for (unsigned int i=0; i<m_InputFiberBundle->GetNumFibers(); i++)
   {
     ++disp;
     vtkCell* cell = polydata->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     bool positive = false;
     if (numPoints>1)
       for (unsigned int m=0; m<m_RoiImages.size(); ++m)
       {
         auto roi = m_RoiImages.at(m);
         m_Interpolator->SetInputImage(roi);
 
         int inside = 0;
 
         // check first fiber point
         {
           double* p = points->GetPoint(0);
           itk::Point<float, 3> itkP = mitk::imv::GetItkPoint(p);
 
           if ( IsPositive(itkP) )
             inside++;
         }
 
         // check second fiber point
         {
           double* p = points->GetPoint(numPoints-1);
           itk::Point<float, 3> itkP = mitk::imv::GetItkPoint(p);
 
           if ( IsPositive(itkP) )
             inside++;
         }
 
         if (inside==2 || (inside==1 && !m_BothEnds))
         {
           positive = true;
           positive_ids[m].push_back(i);
         }
       }
     if (!positive)
       negative_ids.push_back(i);
   }
 
   if (!m_NoNegatives)
     m_Negatives.push_back(CreateFib(negative_ids));
   if (!m_NoPositives)
   {
     for (unsigned int i=0; i<positive_ids.size(); ++i)
     {
       std::string name;
       if (i<m_RoiImageNames.size())
         name = m_RoiImageNames.at(i);
       else
         name = "Roi" + boost::lexical_cast<std::string>(i);
 
       if (positive_ids.at(i).size()>=m_MinFibersPerTract)
       {
         m_Positives.push_back(CreateFib(positive_ids.at(i)));
         m_PositiveLabels.push_back(name);
       }
     }
     if (!m_SplitByRoi)
     {
       mitk::FiberBundle::Pointer output = mitk::FiberBundle::New(nullptr);
       output = output->AddBundles(m_Positives);
       m_Positives.clear();
       m_Positives.push_back(output);
       m_PositiveLabels.clear();
       m_PositiveLabels.push_back("");
     }
   }
 }
 
 template< class PixelType >
 void FiberExtractionFilter< PixelType >::ExtractLabels(mitk::FiberBundle::Pointer fib)
 {
   MITK_INFO << "Extracting fibers by labels";
   vtkSmartPointer<vtkPolyData> polydata = fib->GetFiberPolyData();
 
   std::map< std::string, std::vector< unsigned int > > positive_ids;
 
   std::vector< unsigned int > negative_ids; // fibers not overlapping with ANY label
 
   boost::progress_display disp(m_InputFiberBundle->GetNumFibers());
   for (unsigned int i=0; i<m_InputFiberBundle->GetNumFibers(); i++)
   {
     ++disp;
     vtkCell* cell = polydata->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     bool positive = false;
     if (numPoints>1)
       for (unsigned int m=0; m<m_RoiImages.size(); ++m)
       {
         auto roi = m_RoiImages.at(m);
         m_Interpolator->SetInputImage(roi);
 
         int inside = 0;
 
         double* p1 = points->GetPoint(0);
         itk::Point<float, 3> itkP1 = mitk::imv::GetItkPoint(p1);
         short label1 = mitk::imv::GetImageValue<PixelType>(itkP1, false, m_Interpolator);
 
         double* p2 = points->GetPoint(numPoints-1);
         itk::Point<float, 3> itkP2 = mitk::imv::GetItkPoint(p2);
         short label2 = mitk::imv::GetImageValue<PixelType>(itkP2, false, m_Interpolator);
 
         if (!m_Labels.empty())  // extract fibers from all pairwise label combinations
         {
           for (auto l : m_Labels)
           {
             if (l==label1)
               inside++;
             if (l==label2)
               inside++;
             if (inside==2)
               break;
           }
         }
         else if (!m_StartLabels.empty() || !m_EndLabels.empty()) // extract fibers between start and end labels
         {
           if (!m_StartLabels.empty() && !m_EndLabels.empty())
             m_BothEnds = true;  // if we have start and end labels it does not make sense to not use both endpoints
           if (m_PairedStartEndLabels)
           {
             if (m_StartLabels.size()!=m_EndLabels.size())
               mitkThrow() << "Start and end label lists must have same size if paired labels are used";
             for (unsigned int ii=0; ii<m_StartLabels.size(); ++ii)
               if ( (m_StartLabels.at(ii)==label1 && m_EndLabels.at(ii)==label2) || (m_StartLabels.at(ii)==label2 && m_EndLabels.at(ii)==label1) )
               {
                 inside = 2;
                 break;
               }
           }
           else
           {
             for (unsigned int ii=0; ii<m_StartLabels.size(); ++ii)
             {
               if ( m_StartLabels.at(ii)==label1 || m_StartLabels.at(ii)==label2 )
               {
                 ++inside;
                 break;
               }
             }
             for (unsigned int jj=0; jj<m_EndLabels.size(); ++jj)
             {
               if ( m_EndLabels.at(jj)==label1 || m_EndLabels.at(jj)==label2 )
               {
                 ++inside;
                 break;
               }
             }
           }
         }
         else  // use all labels
         {
           if (label1!=0)
             ++inside;
           if (label2!=0)
             ++inside;
         }
 
         std::string key = "";
 
         if (m_SplitLabels)
         {
           if (label1<label2)
             key = boost::lexical_cast<std::string>(label1) + "-" + boost::lexical_cast<std::string>(label2);
           else
             key = boost::lexical_cast<std::string>(label2) + "-" + boost::lexical_cast<std::string>(label1);
         }
 
         if (m_SplitByRoi)
         {
           if (m<m_RoiImageNames.size())
             key = m_RoiImageNames.at(m) + "_" + key;
           else
             key = "Roi" + boost::lexical_cast<std::string>(m) + "_" + key;
         }
 
         if (m_BothEnds)
         {
           if ( (inside==2 && (!m_SkipSelfConnections || label1!=label2)) || (inside==2 && m_OnlySelfConnections && label1==label2) )
           {
             positive = true;
             if ( positive_ids.count(key)==0 )
               positive_ids.insert( std::pair< std::string, std::vector< unsigned int > >( key, {i}) );
             else
               positive_ids[key].push_back(i);
           }
         }
         else
         {
           if ( (inside>=1 && (!m_SkipSelfConnections || label1!=label2)) || (inside==2 && m_OnlySelfConnections && label1==label2) )
           {
             positive = true;
             if ( positive_ids.count(key)==0 )
               positive_ids.insert( std::pair< std::string, std::vector< unsigned int > >( key, {i}) );
             else
               positive_ids[key].push_back(i);
           }
         }
       }
     if (!positive)
       negative_ids.push_back(i);
   }
 
   if (!m_NoNegatives)
     m_Negatives.push_back(CreateFib(negative_ids));
   if (!m_NoPositives)
   {
     for (auto label : positive_ids)
     {
       if (label.second.size()>=m_MinFibersPerTract)
       {
         m_Positives.push_back(CreateFib(label.second));
         m_PositiveLabels.push_back(label.first);
       }
     }
   }
 }
 
 template< class PixelType >
 void FiberExtractionFilter< PixelType >::SetLabels(const std::vector<unsigned short> &Labels)
 {
   m_Labels = Labels;
 }
 
 template< class PixelType >
 void FiberExtractionFilter< PixelType >::SetStartLabels(const std::vector<unsigned short> &Labels)
 {
   m_StartLabels = Labels;
 }
 
 template< class PixelType >
 void FiberExtractionFilter< PixelType >::SetEndLabels(const std::vector<unsigned short> &Labels)
 {
   m_EndLabels = Labels;
 }
 
 template< class PixelType >
 std::vector<mitk::FiberBundle::Pointer> FiberExtractionFilter< PixelType >::GetNegatives() const
 {
   return m_Negatives;
 }
 
 template< class PixelType >
 std::vector<mitk::FiberBundle::Pointer> FiberExtractionFilter< PixelType >::GetPositives() const
 {
   return m_Positives;
 }
 
 template< class PixelType >
 void FiberExtractionFilter< PixelType >::GenerateData()
 {
   mitk::FiberBundle::Pointer fib = m_InputFiberBundle;
   if (fib->GetNumFibers()<=0)
   {
     MITK_INFO << "No fibers in tractogram!";
     return;
   }
 
   if (m_Mode == MODE::OVERLAP)
     ExtractOverlap(fib);
   else if (m_Mode == MODE::ENDPOINTS)
   {
     if (m_InputType==INPUT::LABEL_MAP)
       ExtractLabels(fib);
     else
       ExtractEndpoints(fib);
   }
+
+  for (auto f : m_Negatives)
+    f->SetTrackVisHeader(fib->GetTrackVisHeader());
+  for (auto f : m_Positives)
+    f->SetTrackVisHeader(fib->GetTrackVisHeader());
 }
 
 }
 
 #endif // __itkFiberExtractionFilter_cpp
 
 
 
diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp
index 9d41f0c..c0fac5d 100644
--- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp
+++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp
@@ -1,1730 +1,1728 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 
 
 // Blueberry
 #include <berryISelectionService.h>
 #include <berryIWorkbenchPart.h>
 #include <berryIWorkbenchWindow.h>
 
 // Qmitk
 #include "QmitkFiberProcessingView.h"
 
 #include <QMessageBox>
 
 #include <mitkNodePredicateProperty.h>
 #include <mitkImageCast.h>
 #include <mitkPointSet.h>
 #include <mitkPlanarCircle.h>
 #include <mitkPlanarPolygon.h>
 #include <mitkPlanarRectangle.h>
 #include <mitkPlanarFigureInteractor.h>
 #include <mitkImageAccessByItk.h>
 #include <mitkDataNodeObject.h>
 #include <mitkTensorImage.h>
 #include "usModuleRegistry.h"
 #include <itkFiberCurvatureFilter.h>
 
 #include "mitkNodePredicateDataType.h"
 #include <mitkNodePredicateProperty.h>
 #include <mitkNodePredicateAnd.h>
 #include <mitkNodePredicateNot.h>
 #include <mitkNodePredicateOr.h>
 
 #include <itkResampleImageFilter.h>
 #include <itkGaussianInterpolateImageFunction.h>
 #include <itkImageRegionIteratorWithIndex.h>
 #include <itkTractsToFiberEndingsImageFilter.h>
 #include <itkTractDensityImageFilter.h>
 #include <itkImageRegion.h>
 #include <itkTractsToRgbaImageFilter.h>
 #include <itkFiberExtractionFilter.h>
 
 #include <boost/algorithm/string.hpp>
 #include <boost/lexical_cast.hpp>
 
 const std::string QmitkFiberProcessingView::VIEW_ID = "org.mitk.views.fiberprocessing";
 const std::string id_DataManager = "org.mitk.views.datamanager";
 using namespace mitk;
 
 QmitkFiberProcessingView::QmitkFiberProcessingView()
   : QmitkAbstractView()
   , m_Controls( 0 )
   , m_CircleCounter(0)
   , m_PolygonCounter(0)
   , m_UpsamplingFactor(1)
 {
 
 }
 
 // Destructor
 QmitkFiberProcessingView::~QmitkFiberProcessingView()
 {
   RemoveObservers();
 }
 
 void QmitkFiberProcessingView::CreateQtPartControl( QWidget *parent )
 {
   // build up qt view, unless already done
   if ( !m_Controls )
   {
     // create GUI widgets from the Qt Designer's .ui file
     m_Controls = new Ui::QmitkFiberProcessingViewControls;
     m_Controls->setupUi( parent );
 
     connect( m_Controls->m_CircleButton, SIGNAL( clicked() ), this, SLOT( OnDrawCircle() ) );
     connect( m_Controls->m_PolygonButton, SIGNAL( clicked() ), this, SLOT( OnDrawPolygon() ) );
     connect(m_Controls->PFCompoANDButton, SIGNAL(clicked()), this, SLOT(GenerateAndComposite()) );
     connect(m_Controls->PFCompoORButton, SIGNAL(clicked()), this, SLOT(GenerateOrComposite()) );
     connect(m_Controls->PFCompoNOTButton, SIGNAL(clicked()), this, SLOT(GenerateNotComposite()) );
     connect(m_Controls->m_GenerateRoiImage, SIGNAL(clicked()), this, SLOT(GenerateRoiImage()) );
 
     connect(m_Controls->m_JoinBundles, SIGNAL(clicked()), this, SLOT(JoinBundles()) );
     connect(m_Controls->m_SubstractBundles, SIGNAL(clicked()), this, SLOT(SubstractBundles()) );
     connect(m_Controls->m_CopyBundle, SIGNAL(clicked()), this, SLOT(CopyBundles()) );
 
     connect(m_Controls->m_ExtractFibersButton, SIGNAL(clicked()), this, SLOT(Extract()));
     connect(m_Controls->m_RemoveButton, SIGNAL(clicked()), this, SLOT(Remove()));
     connect(m_Controls->m_ModifyButton, SIGNAL(clicked()), this, SLOT(Modify()));
 
     connect(m_Controls->m_ExtractionMethodBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui()));
     connect(m_Controls->m_RemovalMethodBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui()));
     connect(m_Controls->m_ModificationMethodBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui()));
     connect(m_Controls->m_ExtractionBoxMask, SIGNAL(currentIndexChanged(int)), this, SLOT(OnMaskExtractionChanged()));
 
     m_Controls->m_ColorMapBox->SetDataStorage(this->GetDataStorage());
     mitk::TNodePredicateDataType<mitk::Image>::Pointer isMitkImage = mitk::TNodePredicateDataType<mitk::Image>::New();
     mitk::NodePredicateDataType::Pointer isDwi = mitk::NodePredicateDataType::New("DiffusionImage");
     mitk::NodePredicateDataType::Pointer isDti = mitk::NodePredicateDataType::New("TensorImage");
     mitk::NodePredicateDataType::Pointer isOdf = mitk::NodePredicateDataType::New("OdfImage");
     mitk::NodePredicateOr::Pointer isDiffusionImage = mitk::NodePredicateOr::New(isDwi, isDti);
     isDiffusionImage = mitk::NodePredicateOr::New(isDiffusionImage, isOdf);
     mitk::NodePredicateNot::Pointer noDiffusionImage = mitk::NodePredicateNot::New(isDiffusionImage);
     mitk::NodePredicateAnd::Pointer finalPredicate = mitk::NodePredicateAnd::New(isMitkImage, noDiffusionImage);
     m_Controls->m_ColorMapBox->SetPredicate(finalPredicate);
   }
 
   UpdateGui();
   OnMaskExtractionChanged();
 }
 
 void QmitkFiberProcessingView::OnMaskExtractionChanged()
 {
   m_Controls->m_FiberExtractionFractionLabel->setVisible(false);
   m_Controls->m_FiberExtractionFractionBox->setVisible(false);
 
   m_Controls->m_FiberExtractionThresholdLabel->setVisible(false);
   m_Controls->m_FiberExtractionThresholdBox->setVisible(false);
 
   m_Controls->m_InterpolateRoiBox->setVisible(false);
   m_Controls->m_BothEnds->setVisible(false);
 
   m_Controls->m_LabelsBox->setVisible(false);
   m_Controls->m_LabelsLabel->setVisible(false);
 
   if (m_Controls->m_ExtractionBoxMask->currentIndex() == 2 || m_Controls->m_ExtractionBoxMask->currentIndex() == 3)
   {
     m_Controls->m_FiberExtractionFractionLabel->setVisible(true);
     m_Controls->m_FiberExtractionFractionBox->setVisible(true);
 
     m_Controls->m_FiberExtractionThresholdLabel->setVisible(true);
     m_Controls->m_FiberExtractionThresholdBox->setVisible(true);
 
     m_Controls->m_InterpolateRoiBox->setVisible(true);
   }
   else if (m_Controls->m_ExtractionBoxMask->currentIndex() == 0 || m_Controls->m_ExtractionBoxMask->currentIndex() == 1)
   {
     if (m_Controls->m_ExtractionBoxMask->currentIndex() != 3)
       m_Controls->m_BothEnds->setVisible(true);
     m_Controls->m_InterpolateRoiBox->setVisible(true);
 
     m_Controls->m_FiberExtractionThresholdLabel->setVisible(true);
     m_Controls->m_FiberExtractionThresholdBox->setVisible(true);
   }
   else if (m_Controls->m_ExtractionBoxMask->currentIndex() == 4)
   {
     m_Controls->m_BothEnds->setVisible(true);
 
     m_Controls->m_LabelsBox->setVisible(true);
     m_Controls->m_LabelsLabel->setVisible(true);
   }
 }
 
 void QmitkFiberProcessingView::SetFocus()
 {
   m_Controls->toolBoxx->setFocus();
 }
 
 void QmitkFiberProcessingView::Modify()
 {
   switch (m_Controls->m_ModificationMethodBox->currentIndex())
   {
   case 0:
   {
     ResampleSelectedBundlesSpline();
     break;
   }
   case 1:
   {
     ResampleSelectedBundlesLinear();
     break;
   }
   case 2:
   {
     CompressSelectedBundles();
     break;
   }
   case 3:
   {
     DoImageColorCoding();
     break;
   }
   case 4:
   {
     MirrorFibers();
     break;
   }
   case 5:
   {
     WeightFibers();
     break;
   }
   case 6:
   {
     DoCurvatureColorCoding();
     break;
   }
   case 7:
   {
     DoWeightColorCoding();
     break;
   }
   case 8:
   {
     DoLengthColorCoding();
     break;
   }
   }
 }
 
 void QmitkFiberProcessingView::WeightFibers()
 {
   float weight = this->m_Controls->m_BundleWeightBox->value();
   for (auto node : m_SelectedFB)
   {
     mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
     fib->SetFiberWeights(weight);
   }
 
 }
 
 void QmitkFiberProcessingView::Remove()
 {
   switch (m_Controls->m_RemovalMethodBox->currentIndex())
   {
   case 0:
   {
     RemoveDir();
     break;
   }
   case 1:
   {
     PruneBundle();
     break;
   }
   case 2:
   {
     ApplyCurvatureThreshold();
     break;
   }
   case 3:
   {
     RemoveWithMask(false);
     break;
   }
   case 4:
   {
     RemoveWithMask(true);
     break;
   }
   case 5:
   {
     ApplyWeightThreshold();
     break;
   }
   case 6:
   {
     ApplyDensityThreshold();
     break;
   }
   }
 }
 
 void QmitkFiberProcessingView::Extract()
 {
   switch (m_Controls->m_ExtractionMethodBox->currentIndex())
   {
   case 0:
   {
     ExtractWithPlanarFigure();
     break;
   }
   case 1:
   {
     switch (m_Controls->m_ExtractionBoxMask->currentIndex())
     {
     {
     case 0:
         ExtractWithMask(true, false, false);
         break;
     }
     {
     case 1:
         ExtractWithMask(true, true, false);
         break;
     }
     {
     case 2:
         ExtractWithMask(false, false, false);
         break;
     }
     {
     case 3:
         ExtractWithMask(false, true, false);
         break;
     }
     {
     case 4:
         ExtractWithMask(true, false, true);
         break;
     }
     }
     break;
   }
   }
 }
 
 void QmitkFiberProcessingView::PruneBundle()
 {
   int minLength = this->m_Controls->m_PruneFibersMinBox->value();
   int maxLength = this->m_Controls->m_PruneFibersMaxBox->value();
   for (auto node : m_SelectedFB)
   {
     mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
     if (!fib->RemoveShortFibers(minLength))
       QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers.");
     else if (!fib->RemoveLongFibers(maxLength))
       QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers.");
   }
   RenderingManager::GetInstance()->RequestUpdateAll();
 }
 
 void QmitkFiberProcessingView::ApplyWeightThreshold()
 {
   float thr = this->m_Controls->m_WeightThresholdBox->value();
   std::vector< DataNode::Pointer > nodes = m_SelectedFB;
   for (auto node : nodes)
   {
     mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
 
     mitk::FiberBundle::Pointer newFib = fib->FilterByWeights(thr);
     if (newFib->GetNumFibers()>0)
     {
       newFib->ColorFibersByFiberWeights(false, true);
       node->SetData(newFib);
     }
     else
       QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers.");
   }
   RenderingManager::GetInstance()->RequestUpdateAll();
 }
 
 void QmitkFiberProcessingView::ApplyDensityThreshold()
 {
   float thr = this->m_Controls->m_DensityThresholdBox->value();
   float ol = this->m_Controls->m_DensityOverlapBox->value();
   std::vector< DataNode::Pointer > nodes = m_SelectedFB;
   for (auto node : nodes)
   {
     mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
 
     itk::TractDensityImageFilter< ItkFloatImageType >::Pointer generator = itk::TractDensityImageFilter< ItkFloatImageType >::New();
     generator->SetFiberBundle(fib);
     generator->SetBinaryOutput(false);
     generator->SetOutputAbsoluteValues(false);
     generator->Update();
 
     itk::FiberExtractionFilter<float>::Pointer extractor = itk::FiberExtractionFilter<float>::New();
     extractor->SetRoiImages({generator->GetOutput()});
     extractor->SetInputFiberBundle(fib);
     extractor->SetOverlapFraction(ol);
     extractor->SetInterpolate(true);
     extractor->SetThreshold(thr);
     extractor->SetNoNegatives(true);
     extractor->Update();
 
     if (extractor->GetPositives().empty())
     {
       QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers.");
       continue;
     }
 
     mitk::FiberBundle::Pointer newFib = extractor->GetPositives().at(0);
     if (newFib->GetNumFibers()>0)
       node->SetData(newFib);
     else
       QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers.");
   }
   RenderingManager::GetInstance()->RequestUpdateAll();
 }
 
 void QmitkFiberProcessingView::ApplyCurvatureThreshold()
 {
   int angle = this->m_Controls->m_CurvSpinBox->value();
   int dist = this->m_Controls->m_CurvDistanceSpinBox->value();
   std::vector< DataNode::Pointer > nodes = m_SelectedFB;
   for (auto node : nodes)
   {
     mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
 
     itk::FiberCurvatureFilter::Pointer filter = itk::FiberCurvatureFilter::New();
     filter->SetInputFiberBundle(fib);
     filter->SetAngularDeviation(angle);
     filter->SetDistance(dist);
     filter->SetRemoveFibers(m_Controls->m_RemoveCurvedFibersBox->isChecked());
     filter->Update();
     mitk::FiberBundle::Pointer newFib = filter->GetOutputFiberBundle();
     if (newFib->GetNumFibers()>0)
     {
       newFib->ColorFibersByOrientation();
       node->SetData(newFib);
     }
     else
       QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers.");
   }
   RenderingManager::GetInstance()->RequestUpdateAll();
 }
 
 void QmitkFiberProcessingView::RemoveDir()
 {
   for (auto node : m_SelectedFB)
   {
     mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
     vnl_vector_fixed<double,3> dir;
     dir[0] = m_Controls->m_ExtractDirX->value();
     dir[1] = m_Controls->m_ExtractDirY->value();
     dir[2] = m_Controls->m_ExtractDirZ->value();
     fib->RemoveDir(dir,cos((float)m_Controls->m_ExtractAngle->value()*itk::Math::pi/180));
   }
   RenderingManager::GetInstance()->RequestUpdateAll();
 }
 
 void QmitkFiberProcessingView::RemoveWithMask(bool removeInside)
 {
   if (m_RoiImageNode.IsNull())
     return;
 
   mitk::Image::Pointer mitkMask = dynamic_cast<mitk::Image*>(m_RoiImageNode->GetData());
   for (auto node : m_SelectedFB)
   {
     mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
 
     ItkUCharImageType::Pointer mask = ItkUCharImageType::New();
     mitk::CastToItkImage(mitkMask, mask);
     mitk::FiberBundle::Pointer newFib = fib->RemoveFibersOutside(mask, removeInside);
     if (newFib->GetNumFibers()<=0)
     {
       QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers.");
       continue;
     }
     node->SetData(newFib);
   }
   RenderingManager::GetInstance()->RequestUpdateAll();
 }
 
 void QmitkFiberProcessingView::ExtractWithMask(bool onlyEnds, bool invert, bool labelmap)
 {
   if (m_RoiImageNode.IsNull())
     return;
 
   mitk::Image::Pointer mitkMask = dynamic_cast<mitk::Image*>(m_RoiImageNode->GetData());
   for (auto node : m_SelectedFB)
   {
     std::string roi_name = m_RoiImageNode->GetName();
     mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
     ItkFloatImageType::Pointer mask = ItkFloatImageType::New();
     mitk::CastToItkImage(mitkMask, mask);
 
     itk::FiberExtractionFilter<float>::Pointer extractor = itk::FiberExtractionFilter<float>::New();
     extractor->SetInputFiberBundle(fib);
     extractor->SetRoiImages({mask});
     extractor->SetRoiImageNames({roi_name});
     extractor->SetThreshold(m_Controls->m_FiberExtractionThresholdBox->value());
     extractor->SetOverlapFraction(m_Controls->m_FiberExtractionFractionBox->value());
     extractor->SetBothEnds(m_Controls->m_BothEnds->isChecked());
     extractor->SetInterpolate(m_Controls->m_InterpolateRoiBox->isChecked());
     extractor->SetMinFibersPerTract(m_Controls->m_MinExtractedFibersBox->value());
     extractor->SetSplitByRoi(true);
     if (invert)
       extractor->SetNoPositives(true);
     else
       extractor->SetNoNegatives(true);
 
     if (labelmap)
     {
       std::string labels_string = m_Controls->m_LabelsBox->text().toStdString();
       if (labels_string!="ALL")
       {
         std::vector<std::string> strs;
         boost::split(strs,labels_string,boost::is_any_of(" ,;\t"));
 
         std::vector< unsigned short > labels_vector;
         for (auto v : strs)
         {
           try{
             unsigned short l = boost::lexical_cast<unsigned short>(v);
             labels_vector.push_back(l);
           }
           catch(...)
           {
 
           }
         }
         extractor->SetLabels(labels_vector);
       }
 
       extractor->SetInterpolate(false);
       extractor->SetInputType(itk::FiberExtractionFilter<float>::INPUT::LABEL_MAP);
       extractor->SetSplitLabels(true);
       onlyEnds = true;
     }
 
     if (onlyEnds)
       extractor->SetMode(itk::FiberExtractionFilter<float>::MODE::ENDPOINTS);
     extractor->Update();
 
     std::vector< mitk::FiberBundle::Pointer > newFibs;
     if (invert)
       newFibs = extractor->GetNegatives();
     else
       newFibs = extractor->GetPositives();
 
     if (newFibs.empty())
     {
       QMessageBox::information(nullptr, "No output generated:", "No fibers could be extracted.");
       continue;
     }
 
     auto labels = extractor->GetPositiveLabels();
     for (unsigned int i=0; i<newFibs.size(); ++i)
     {
       DataNode::Pointer newNode = DataNode::New();
       auto fib = newFibs.at(i);
 
       newNode->SetData(fib);
       std::string name = roi_name;
       if (i<labels.size())
         name = labels.at(i);
 
       if (invert)
       {
         if (onlyEnds)
           name = "not-ending-in-" + name;
         else
           name = "not-passing-" + name;
       }
       else if (!labelmap)
       {
         if (onlyEnds)
           name = "ending-in-" + name;
         else
           name = "passing-" + name;
       }
 
       newNode->SetName(name);
       GetDataStorage()->Add(newNode, node);
     }
     node->SetVisibility(false);
   }
 }
 
 void QmitkFiberProcessingView::GenerateRoiImage()
 {
   if (m_SelectedPF.empty())
     return;
 
   mitk::BaseGeometry::Pointer geometry;
   if (!m_SelectedFB.empty())
   {
     mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(m_SelectedFB.front()->GetData());
     geometry = fib->GetGeometry();
   }
   else if (m_SelectedImage)
     geometry = m_SelectedImage->GetGeometry();
   else
     return;
 
   itk::Vector<double,3> spacing = geometry->GetSpacing();
   spacing /= m_UpsamplingFactor;
 
   mitk::Point3D newOrigin = geometry->GetOrigin();
   mitk::Geometry3D::BoundsArrayType bounds = geometry->GetBounds();
   newOrigin[0] += bounds.GetElement(0);
   newOrigin[1] += bounds.GetElement(2);
   newOrigin[2] += bounds.GetElement(4);
 
   itk::Matrix<double, 3, 3> direction;
   itk::ImageRegion<3> imageRegion;
   for (int i=0; i<3; i++)
     for (int j=0; j<3; j++)
       direction[j][i] = geometry->GetMatrixColumn(i)[j]/spacing[j];
   imageRegion.SetSize(0, geometry->GetExtent(0)*m_UpsamplingFactor);
   imageRegion.SetSize(1, geometry->GetExtent(1)*m_UpsamplingFactor);
   imageRegion.SetSize(2, geometry->GetExtent(2)*m_UpsamplingFactor);
 
   m_PlanarFigureImage = ItkUCharImageType::New();
   m_PlanarFigureImage->SetSpacing( spacing );   // Set the image spacing
   m_PlanarFigureImage->SetOrigin( newOrigin );     // Set the image origin
   m_PlanarFigureImage->SetDirection( direction );  // Set the image direction
   m_PlanarFigureImage->SetRegions( imageRegion );
   m_PlanarFigureImage->Allocate();
   m_PlanarFigureImage->FillBuffer( 0 );
 
   Image::Pointer tmpImage = Image::New();
   tmpImage->InitializeByItk(m_PlanarFigureImage.GetPointer());
   tmpImage->SetVolume(m_PlanarFigureImage->GetBufferPointer());
 
   std::string name = m_SelectedPF.at(0)->GetName();
   WritePfToImage(m_SelectedPF.at(0), tmpImage);
   for (unsigned int i=1; i<m_SelectedPF.size(); i++)
   {
     name += "+";
     name += m_SelectedPF.at(i)->GetName();
     WritePfToImage(m_SelectedPF.at(i), tmpImage);
   }
 
   DataNode::Pointer node = DataNode::New();
   tmpImage = Image::New();
   tmpImage->InitializeByItk(m_PlanarFigureImage.GetPointer());
   tmpImage->SetVolume(m_PlanarFigureImage->GetBufferPointer());
   node->SetData(tmpImage);
   node->SetName(name);
   this->GetDataStorage()->Add(node);
 }
 
 void QmitkFiberProcessingView::WritePfToImage(mitk::DataNode::Pointer node, mitk::Image* image)
 {
   if (dynamic_cast<mitk::PlanarFigure*>(node->GetData()))
   {
     m_PlanarFigure = dynamic_cast<mitk::PlanarFigure*>(node->GetData());
     AccessFixedDimensionByItk_2(
           image,
           InternalReorientImagePlane, 3,
           m_PlanarFigure->GetGeometry(), -1);
 
     AccessFixedDimensionByItk_2(
           m_InternalImage,
           InternalCalculateMaskFromPlanarFigure,
           3, 2, node->GetName() );
   }
   else if (dynamic_cast<mitk::PlanarFigureComposite*>(node->GetData()))
   {
     DataStorage::SetOfObjects::ConstPointer children = GetDataStorage()->GetDerivations(node);
     for (unsigned int i=0; i<children->Size(); i++)
     {
       WritePfToImage(children->at(i), image);
     }
   }
 }
 
 template < typename TPixel, unsigned int VImageDimension >
 void QmitkFiberProcessingView::InternalReorientImagePlane( const itk::Image< TPixel, VImageDimension > *image, mitk::BaseGeometry* planegeo3D, int additionalIndex )
 {
   typedef itk::Image< TPixel, VImageDimension > ImageType;
   typedef itk::Image< float, VImageDimension > FloatImageType;
 
   typedef itk::ResampleImageFilter<ImageType, FloatImageType, double> ResamplerType;
   typename ResamplerType::Pointer resampler = ResamplerType::New();
 
   mitk::PlaneGeometry* planegeo = dynamic_cast<mitk::PlaneGeometry*>(planegeo3D);
 
   float upsamp = m_UpsamplingFactor;
   float gausssigma = 0.5;
 
   // Spacing
   typename ResamplerType::SpacingType spacing = planegeo->GetSpacing();
   spacing[0] = image->GetSpacing()[0] / upsamp;
   spacing[1] = image->GetSpacing()[1] / upsamp;
   spacing[2] = image->GetSpacing()[2];
   resampler->SetOutputSpacing( spacing );
 
   // Size
   typename ResamplerType::SizeType size;
   size[0] = planegeo->GetExtentInMM(0) / spacing[0];
   size[1] = planegeo->GetExtentInMM(1) / spacing[1];
   size[2] = 1;
   resampler->SetSize( size );
 
   // Origin
   typename mitk::Point3D orig = planegeo->GetOrigin();
   typename mitk::Point3D corrorig;
   planegeo3D->WorldToIndex(orig,corrorig);
   corrorig[0] += 0.5/upsamp;
   corrorig[1] += 0.5/upsamp;
   corrorig[2] += 0;
   planegeo3D->IndexToWorld(corrorig,corrorig);
   resampler->SetOutputOrigin(corrorig );
 
   // Direction
   typename ResamplerType::DirectionType direction;
   typename mitk::AffineTransform3D::MatrixType matrix = planegeo->GetIndexToWorldTransform()->GetMatrix();
   for(unsigned int c=0; c<matrix.ColumnDimensions; c++)
   {
     double sum = 0;
     for(unsigned int r=0; r<matrix.RowDimensions; r++)
       sum += matrix(r,c)*matrix(r,c);
     for(unsigned int r=0; r<matrix.RowDimensions; r++)
       direction(r,c) = matrix(r,c)/sqrt(sum);
   }
   resampler->SetOutputDirection( direction );
 
   // Gaussian interpolation
   if(gausssigma != 0)
   {
     double sigma[3];
     for( unsigned int d = 0; d < 3; d++ )
       sigma[d] = gausssigma * image->GetSpacing()[d];
     double alpha = 2.0;
     typedef itk::GaussianInterpolateImageFunction<ImageType, double> GaussianInterpolatorType;
     typename GaussianInterpolatorType::Pointer interpolator = GaussianInterpolatorType::New();
     interpolator->SetInputImage( image );
     interpolator->SetParameters( sigma, alpha );
     resampler->SetInterpolator( interpolator );
   }
   else
   {
     typedef typename itk::LinearInterpolateImageFunction<ImageType, double> InterpolatorType;
     typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
     interpolator->SetInputImage( image );
     resampler->SetInterpolator( interpolator );
   }
 
   resampler->SetInput( image );
   resampler->SetDefaultPixelValue(0);
   resampler->Update();
 
   if(additionalIndex < 0)
   {
     this->m_InternalImage = mitk::Image::New();
     this->m_InternalImage->InitializeByItk( resampler->GetOutput() );
     this->m_InternalImage->SetVolume( resampler->GetOutput()->GetBufferPointer() );
   }
 }
 
 template < typename TPixel, unsigned int VImageDimension >
 void QmitkFiberProcessingView::InternalCalculateMaskFromPlanarFigure( itk::Image< TPixel, VImageDimension > *image, unsigned int axis, std::string )
 {
   typedef itk::Image< TPixel, VImageDimension > ImageType;
 
   // Generate mask image as new image with same header as input image and
   // initialize with "1".
   ItkUCharImageType::Pointer newMaskImage = ItkUCharImageType::New();
   newMaskImage->SetSpacing( image->GetSpacing() );   // Set the image spacing
   newMaskImage->SetOrigin( image->GetOrigin() );     // Set the image origin
   newMaskImage->SetDirection( image->GetDirection() );  // Set the image direction
   newMaskImage->SetRegions( image->GetLargestPossibleRegion() );
   newMaskImage->Allocate();
   newMaskImage->FillBuffer( 1 );
 
   // Generate VTK polygon from (closed) PlanarFigure polyline
   // (The polyline points are shifted by -0.5 in z-direction to make sure
   // that the extrusion filter, which afterwards elevates all points by +0.5
   // in z-direction, creates a 3D object which is cut by the the plane z=0)
   const PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry();
   const PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 );
   const BaseGeometry *imageGeometry3D = m_InternalImage->GetGeometry( 0 );
 
   vtkPolyData *polyline = vtkPolyData::New();
   polyline->Allocate( 1, 1 );
 
   // Determine x- and y-dimensions depending on principal axis
   int i0, i1;
   switch ( axis )
   {
   case 0:
     i0 = 1;
     i1 = 2;
     break;
   case 1:
     i0 = 0;
     i1 = 2;
     break;
   case 2:
   default:
     i0 = 0;
     i1 = 1;
     break;
   }
 
   // Create VTK polydata object of polyline contour
   vtkPoints *points = vtkPoints::New();
   PlanarFigure::PolyLineType::const_iterator it;
   unsigned int numberOfPoints = 0;
 
   for ( it = planarFigurePolyline.begin(); it != planarFigurePolyline.end(); ++it )
   {
     Point3D point3D;
 
     // Convert 2D point back to the local index coordinates of the selected image
     Point2D point2D = *it;
     planarFigurePlaneGeometry->WorldToIndex(point2D, point2D);
     point2D[0] -= 0.5/m_UpsamplingFactor;
     point2D[1] -= 0.5/m_UpsamplingFactor;
     planarFigurePlaneGeometry->IndexToWorld(point2D, point2D);
     planarFigurePlaneGeometry->Map( point2D, point3D );
 
     // Polygons (partially) outside of the image bounds can not be processed further due to a bug in vtkPolyDataToImageStencil
     if ( !imageGeometry3D->IsInside( point3D ) )
     {
       float bounds[2] = {0,0};
       bounds[0] =
           this->m_InternalImage->GetLargestPossibleRegion().GetSize().GetElement(i0);
       bounds[1] =
           this->m_InternalImage->GetLargestPossibleRegion().GetSize().GetElement(i1);
 
       imageGeometry3D->WorldToIndex( point3D, point3D );
 
       if (point3D[i0]<0)
         point3D[i0] = 0.0;
       else if (point3D[i0]>bounds[0])
         point3D[i0] = bounds[0]-0.001;
 
       if (point3D[i1]<0)
         point3D[i1] = 0.0;
       else if (point3D[i1]>bounds[1])
         point3D[i1] = bounds[1]-0.001;
 
       points->InsertNextPoint( point3D[i0], point3D[i1], -0.5 );
       numberOfPoints++;
     }
     else
     {
       imageGeometry3D->WorldToIndex( point3D, point3D );
       // Add point to polyline array
       points->InsertNextPoint( point3D[i0], point3D[i1], -0.5 );
       numberOfPoints++;
     }
   }
   polyline->SetPoints( points );
   points->Delete();
 
   vtkIdType *ptIds = new vtkIdType[numberOfPoints];
   for ( vtkIdType i = 0; i < numberOfPoints; ++i )
     ptIds[i] = i;
   polyline->InsertNextCell( VTK_POLY_LINE, numberOfPoints, ptIds );
 
   // Extrude the generated contour polygon
   vtkLinearExtrusionFilter *extrudeFilter = vtkLinearExtrusionFilter::New();
   extrudeFilter->SetInputData( polyline );
   extrudeFilter->SetScaleFactor( 1 );
   extrudeFilter->SetExtrusionTypeToNormalExtrusion();
   extrudeFilter->SetVector( 0.0, 0.0, 1.0 );
 
   // Make a stencil from the extruded polygon
   vtkPolyDataToImageStencil *polyDataToImageStencil = vtkPolyDataToImageStencil::New();
   polyDataToImageStencil->SetInputConnection( extrudeFilter->GetOutputPort() );
 
   // Export from ITK to VTK (to use a VTK filter)
   typedef itk::VTKImageImport< ItkUCharImageType > ImageImportType;
   typedef itk::VTKImageExport< ItkUCharImageType > ImageExportType;
 
   typename ImageExportType::Pointer itkExporter = ImageExportType::New();
   itkExporter->SetInput( newMaskImage );
 
   vtkImageImport *vtkImporter = vtkImageImport::New();
   this->ConnectPipelines( itkExporter, vtkImporter );
   vtkImporter->Update();
 
   // Apply the generated image stencil to the input image
   vtkImageStencil *imageStencilFilter = vtkImageStencil::New();
   imageStencilFilter->SetInputConnection( vtkImporter->GetOutputPort() );
   imageStencilFilter->SetStencilConnection(polyDataToImageStencil->GetOutputPort() );
   imageStencilFilter->ReverseStencilOff();
   imageStencilFilter->SetBackgroundValue( 0 );
   imageStencilFilter->Update();
 
   // Export from VTK back to ITK
   vtkImageExport *vtkExporter = vtkImageExport::New();
   vtkExporter->SetInputConnection( imageStencilFilter->GetOutputPort() );
   vtkExporter->Update();
 
   typename ImageImportType::Pointer itkImporter = ImageImportType::New();
   this->ConnectPipelines( vtkExporter, itkImporter );
   itkImporter->Update();
 
   // calculate cropping bounding box
   m_InternalImageMask3D = itkImporter->GetOutput();
   m_InternalImageMask3D->SetDirection(image->GetDirection());
 
   itk::ImageRegionConstIterator<ItkUCharImageType>
       itmask(m_InternalImageMask3D, m_InternalImageMask3D->GetLargestPossibleRegion());
   itk::ImageRegionIterator<ImageType>
       itimage(image, image->GetLargestPossibleRegion());
 
   itmask.GoToBegin();
   itimage.GoToBegin();
 
   typename ImageType::SizeType lowersize = {{itk::NumericTraits<unsigned long>::max(),itk::NumericTraits<unsigned long>::max(),itk::NumericTraits<unsigned long>::max()}};
   typename ImageType::SizeType uppersize = {{0,0,0}};
   while( !itmask.IsAtEnd() )
   {
     if(itmask.Get() == 0)
       itimage.Set(0);
     else
     {
       typename ImageType::IndexType index = itimage.GetIndex();
       typename ImageType::SizeType signedindex;
       signedindex[0] = index[0];
       signedindex[1] = index[1];
       signedindex[2] = index[2];
 
       lowersize[0] = signedindex[0] < lowersize[0] ? signedindex[0] : lowersize[0];
       lowersize[1] = signedindex[1] < lowersize[1] ? signedindex[1] : lowersize[1];
       lowersize[2] = signedindex[2] < lowersize[2] ? signedindex[2] : lowersize[2];
 
       uppersize[0] = signedindex[0] > uppersize[0] ? signedindex[0] : uppersize[0];
       uppersize[1] = signedindex[1] > uppersize[1] ? signedindex[1] : uppersize[1];
       uppersize[2] = signedindex[2] > uppersize[2] ? signedindex[2] : uppersize[2];
     }
 
     ++itmask;
     ++itimage;
   }
 
   typename ImageType::IndexType index;
   index[0] = lowersize[0];
   index[1] = lowersize[1];
   index[2] = lowersize[2];
 
   typename ImageType::SizeType size;
   size[0] = uppersize[0] - lowersize[0] + 1;
   size[1] = uppersize[1] - lowersize[1] + 1;
   size[2] = uppersize[2] - lowersize[2] + 1;
 
   itk::ImageRegion<3> cropRegion = itk::ImageRegion<3>(index, size);
 
   // crop internal mask
   typedef itk::RegionOfInterestImageFilter< ItkUCharImageType, ItkUCharImageType > ROIMaskFilterType;
   typename ROIMaskFilterType::Pointer roi2 = ROIMaskFilterType::New();
   roi2->SetRegionOfInterest(cropRegion);
   roi2->SetInput(m_InternalImageMask3D);
   roi2->Update();
   m_InternalImageMask3D = roi2->GetOutput();
 
   Image::Pointer tmpImage = Image::New();
   tmpImage->InitializeByItk(m_InternalImageMask3D.GetPointer());
   tmpImage->SetVolume(m_InternalImageMask3D->GetBufferPointer());
 
   Image::Pointer tmpImage2 = Image::New();
   tmpImage2->InitializeByItk(m_PlanarFigureImage.GetPointer());
   const BaseGeometry *pfImageGeometry3D = tmpImage2->GetGeometry( 0 );
 
   const BaseGeometry *intImageGeometry3D = tmpImage->GetGeometry( 0 );
 
   typedef itk::ImageRegionIteratorWithIndex<ItkUCharImageType> IteratorType;
   IteratorType imageIterator (m_InternalImageMask3D, m_InternalImageMask3D->GetRequestedRegion());
   imageIterator.GoToBegin();
   while ( !imageIterator.IsAtEnd() )
   {
     unsigned char val = imageIterator.Value();
     if (val>0)
     {
       itk::Index<3> index = imageIterator.GetIndex();
       Point3D point;
       point[0] = index[0];
       point[1] = index[1];
       point[2] = index[2];
 
       intImageGeometry3D->IndexToWorld(point, point);
       pfImageGeometry3D->WorldToIndex(point, point);
 
       point[i0] += 0.5;
       point[i1] += 0.5;
 
       index[0] = point[0];
       index[1] = point[1];
       index[2] = point[2];
 
       if (pfImageGeometry3D->IsIndexInside(index))
         m_PlanarFigureImage->SetPixel(index, 1);
     }
     ++imageIterator;
   }
 
   // Clean up VTK objects
   polyline->Delete();
   extrudeFilter->Delete();
   polyDataToImageStencil->Delete();
   vtkImporter->Delete();
   imageStencilFilter->Delete();
   //vtkExporter->Delete(); // TODO: crashes when outcommented; memory leak??
   delete[] ptIds;
 }
 
 void QmitkFiberProcessingView::UpdateGui()
 {
   m_Controls->m_FibLabel->setText("<font color='red'>mandatory</font>");
   m_Controls->m_PfLabel->setText("<font color='grey'>needed for extraction</font>");
   m_Controls->m_InputData->setTitle("Please Select Input Data");
 
   m_Controls->m_RemoveButton->setEnabled(false);
 
   m_Controls->m_PlanarFigureButtonsFrame->setEnabled(false);
   m_Controls->PFCompoANDButton->setEnabled(false);
   m_Controls->PFCompoORButton->setEnabled(false);
   m_Controls->PFCompoNOTButton->setEnabled(false);
 
   m_Controls->m_GenerateRoiImage->setEnabled(false);
   m_Controls->m_ExtractFibersButton->setEnabled(false);
   m_Controls->m_ModifyButton->setEnabled(false);
 
   m_Controls->m_CopyBundle->setEnabled(false);
   m_Controls->m_JoinBundles->setEnabled(false);
   m_Controls->m_SubstractBundles->setEnabled(false);
 
   // disable alle frames
   m_Controls->m_BundleWeightFrame->setVisible(false);
   m_Controls->m_ExtactionFramePF->setVisible(false);
   m_Controls->m_RemoveDirectionFrame->setVisible(false);
   m_Controls->m_RemoveLengthFrame->setVisible(false);
   m_Controls->m_RemoveCurvatureFrame->setVisible(false);
   m_Controls->m_RemoveByWeightFrame->setVisible(false);
   m_Controls->m_RemoveByDensityFrame->setVisible(false);
   m_Controls->m_SmoothFibersFrame->setVisible(false);
   m_Controls->m_CompressFibersFrame->setVisible(false);
   m_Controls->m_ColorFibersFrame->setVisible(false);
   m_Controls->m_MirrorFibersFrame->setVisible(false);
   m_Controls->m_MaskExtractionFrame->setVisible(false);
   m_Controls->m_ColorMapBox->setVisible(false);
 
   bool pfSelected = !m_SelectedPF.empty();
   bool fibSelected = !m_SelectedFB.empty();
   bool multipleFibsSelected = (m_SelectedFB.size()>1);
   bool maskSelected = m_RoiImageNode.IsNotNull();
   bool imageSelected = m_SelectedImage.IsNotNull();
 
   // toggle visibility of elements according to selected method
   switch ( m_Controls->m_ExtractionMethodBox->currentIndex() )
   {
   case 0:
     m_Controls->m_ExtactionFramePF->setVisible(true);
     break;
   case 1:
     m_Controls->m_MaskExtractionFrame->setVisible(true);
     break;
   }
 
   switch ( m_Controls->m_RemovalMethodBox->currentIndex() )
   {
   case 0:
     m_Controls->m_RemoveDirectionFrame->setVisible(true);
     if ( fibSelected )
       m_Controls->m_RemoveButton->setEnabled(true);
     break;
   case 1:
     m_Controls->m_RemoveLengthFrame->setVisible(true);
     if ( fibSelected )
       m_Controls->m_RemoveButton->setEnabled(true);
     break;
   case 2:
     m_Controls->m_RemoveCurvatureFrame->setVisible(true);
     if ( fibSelected )
       m_Controls->m_RemoveButton->setEnabled(true);
     break;
   case 3:
     break;
   case 4:
     break;
   case 5:
     m_Controls->m_RemoveByWeightFrame->setVisible(true);
     if ( fibSelected )
       m_Controls->m_RemoveButton->setEnabled(true);
     break;
   case 6:
     m_Controls->m_RemoveByDensityFrame->setVisible(true);
     if ( fibSelected )
       m_Controls->m_RemoveButton->setEnabled(true);
     break;
   }
 
   switch ( m_Controls->m_ModificationMethodBox->currentIndex() )
   {
   case 0:
     m_Controls->m_SmoothFibersFrame->setVisible(true);
     break;
   case 1:
     m_Controls->m_SmoothFibersFrame->setVisible(true);
     break;
   case 2:
     m_Controls->m_CompressFibersFrame->setVisible(true);
     break;
   case 3:
     m_Controls->m_ColorFibersFrame->setVisible(true);
     m_Controls->m_ColorMapBox->setVisible(true);
     break;
   case 4:
     m_Controls->m_MirrorFibersFrame->setVisible(true);
     if (m_SelectedSurfaces.size()>0)
       m_Controls->m_ModifyButton->setEnabled(true);
     break;
   case 5:
     m_Controls->m_BundleWeightFrame->setVisible(true);
     break;
   case 6:
     m_Controls->m_ColorFibersFrame->setVisible(true);
     break;
   case 7:
     m_Controls->m_ColorFibersFrame->setVisible(true);
     break;
   case 8:
     m_Controls->m_ColorFibersFrame->setVisible(true);
     break;
   }
 
   // are fiber bundles selected?
   if ( fibSelected )
   {
     m_Controls->m_CopyBundle->setEnabled(true);
     m_Controls->m_ModifyButton->setEnabled(true);
     m_Controls->m_PlanarFigureButtonsFrame->setEnabled(true);
     m_Controls->m_FibLabel->setText(QString(m_SelectedFB.at(0)->GetName().c_str()));
 
     // one bundle and one planar figure needed to extract fibers
     if (pfSelected && m_Controls->m_ExtractionMethodBox->currentIndex()==0)
     {
       m_Controls->m_InputData->setTitle("Input Data");
       m_Controls->m_PfLabel->setText(QString(m_SelectedPF.at(0)->GetName().c_str()));
       m_Controls->m_ExtractFibersButton->setEnabled(true);
     }
 
     // more than two bundles needed to join/subtract
     if (multipleFibsSelected)
     {
       m_Controls->m_FibLabel->setText("multiple bundles selected");
 
       m_Controls->m_JoinBundles->setEnabled(true);
       m_Controls->m_SubstractBundles->setEnabled(true);
     }
 
     if (maskSelected && m_Controls->m_ExtractionMethodBox->currentIndex()==1)
     {
       m_Controls->m_InputData->setTitle("Input Data");
       m_Controls->m_PfLabel->setText(QString(m_RoiImageNode->GetName().c_str()));
       m_Controls->m_ExtractFibersButton->setEnabled(true);
     }
 
     if (maskSelected && (m_Controls->m_RemovalMethodBox->currentIndex()==3 || m_Controls->m_RemovalMethodBox->currentIndex()==4) )
     {
       m_Controls->m_InputData->setTitle("Input Data");
       m_Controls->m_PfLabel->setText(QString(m_RoiImageNode->GetName().c_str()));
       m_Controls->m_RemoveButton->setEnabled(true);
     }
   }
 
   // are planar figures selected?
   if (pfSelected)
   {
     if ( fibSelected || m_SelectedImage.IsNotNull())
       m_Controls->m_GenerateRoiImage->setEnabled(true);
 
     if (m_SelectedPF.size() > 1)
     {
       m_Controls->PFCompoANDButton->setEnabled(true);
       m_Controls->PFCompoORButton->setEnabled(true);
     }
     else
       m_Controls->PFCompoNOTButton->setEnabled(true);
   }
 
   // is image selected
   if (imageSelected || maskSelected)
   {
     m_Controls->m_PlanarFigureButtonsFrame->setEnabled(true);
   }
 }
 
 void QmitkFiberProcessingView::NodeRemoved(const mitk::DataNode* node )
 {
   for (auto fnode: m_SelectedFB)
     if (node == fnode)
     {
       m_SelectedFB.clear();
       break;
     }
 
   berry::IWorkbenchPart::Pointer nullPart;
   QList<mitk::DataNode::Pointer> nodes;
   OnSelectionChanged(nullPart, nodes);
 }
 
 void QmitkFiberProcessingView::NodeAdded(const mitk::DataNode* )
 {
   if (!m_Controls->m_InteractiveBox->isChecked())
   {
     berry::IWorkbenchPart::Pointer nullPart;
     QList<mitk::DataNode::Pointer> nodes;
     OnSelectionChanged(nullPart, nodes);
   }
 }
 
 void QmitkFiberProcessingView::OnEndInteraction()
 {
   if (m_Controls->m_InteractiveBox->isChecked())
     ExtractWithPlanarFigure(true);
 }
 
 void QmitkFiberProcessingView::AddObservers()
 {
   typedef itk::SimpleMemberCommand< QmitkFiberProcessingView > SimpleCommandType;
   for (auto node : m_SelectedPF)
   {
     mitk::PlanarFigure* figure = dynamic_cast<mitk::PlanarFigure*>(node->GetData());
     if (figure!=nullptr)
     {
       figure->RemoveAllObservers();
 
       // add observer for event when interaction with figure starts
       SimpleCommandType::Pointer endInteractionCommand = SimpleCommandType::New();
       endInteractionCommand->SetCallbackFunction( this, &QmitkFiberProcessingView::OnEndInteraction);
       m_EndInteractionObserverTag = figure->AddObserver( mitk::EndInteractionPlanarFigureEvent(), endInteractionCommand );
     }
   }
 }
 
 void QmitkFiberProcessingView::RemoveObservers()
 {
   for (auto node : m_SelectedPF)
   {
     mitk::PlanarFigure* figure = dynamic_cast<mitk::PlanarFigure*>(node->GetData());
     if (figure!=nullptr)
       figure->RemoveAllObservers();
   }
 }
 
 void QmitkFiberProcessingView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList<mitk::DataNode::Pointer>& nodes)
 {
   RemoveObservers();
 
   //reset existing Vectors containing FiberBundles and PlanarFigures from a previous selection
   std::vector<mitk::DataNode::Pointer>  lastSelectedFB = m_SelectedFB;
   m_SelectedFB.clear();
   m_SelectedPF.clear();
   m_SelectedSurfaces.clear();
   m_SelectedImage = nullptr;
   m_RoiImageNode = nullptr;
 
   for (auto node: nodes)
   {
     if ( dynamic_cast<mitk::FiberBundle*>(node->GetData()) )
       m_SelectedFB.push_back(node);
     else if (dynamic_cast<mitk::PlanarPolygon*>(node->GetData()) || dynamic_cast<mitk::PlanarFigureComposite*>(node->GetData()) || dynamic_cast<mitk::PlanarCircle*>(node->GetData()))
       m_SelectedPF.push_back(node);
     else if (dynamic_cast<mitk::Image*>(node->GetData()))
     {
       m_SelectedImage = dynamic_cast<mitk::Image*>(node->GetData());
       if (m_SelectedImage->GetDimension()==3)
         m_RoiImageNode = node;
     }
     else if (dynamic_cast<mitk::Surface*>(node->GetData()))
       m_SelectedSurfaces.push_back(dynamic_cast<mitk::Surface*>(node->GetData()));
   }
 
   // if we perform interactive fiber extraction, we want to avoid auto-selection of the extracted bundle
   if (m_SelectedFB.empty() && m_Controls->m_InteractiveBox->isChecked())
     m_SelectedFB = lastSelectedFB;
 
   // if no fibers or surfaces are selected, select topmost
   if (m_SelectedFB.empty() && m_SelectedSurfaces.empty())
   {
     int maxLayer = 0;
     itk::VectorContainer<unsigned int, mitk::DataNode::Pointer>::ConstPointer nodes = this->GetDataStorage()->GetAll();
     for (unsigned int i=0; i<nodes->Size(); i++)
       if (dynamic_cast<mitk::FiberBundle*>(nodes->at(i)->GetData()))
       {
         mitk::DataStorage::SetOfObjects::ConstPointer sources = GetDataStorage()->GetSources(nodes->at(i));
         if (sources->Size()>0)
           continue;
 
         int layer = 0;
         nodes->at(i)->GetPropertyValue("layer", layer);
         if (layer>=maxLayer)
         {
           maxLayer = layer;
           m_SelectedFB.clear();
           m_SelectedFB.push_back(nodes->at(i));
         }
       }
   }
 
   // if no plar figure is selected, select topmost
   if (m_SelectedPF.empty())
   {
     int maxLayer = 0;
     itk::VectorContainer<unsigned int, mitk::DataNode::Pointer>::ConstPointer nodes = this->GetDataStorage()->GetAll();
     for (unsigned int i=0; i<nodes->Size(); i++)
       if (dynamic_cast<mitk::PlanarPolygon*>(nodes->at(i)->GetData()) || dynamic_cast<mitk::PlanarFigureComposite*>(nodes->at(i)->GetData()) || dynamic_cast<mitk::PlanarCircle*>(nodes->at(i)->GetData()))
       {
         mitk::DataStorage::SetOfObjects::ConstPointer sources = GetDataStorage()->GetSources(nodes->at(i));
         if (sources->Size()>0)
           continue;
 
         int layer = 0;
         nodes->at(i)->GetPropertyValue("layer", layer);
         if (layer>=maxLayer)
         {
           maxLayer = layer;
           m_SelectedPF.clear();
           m_SelectedPF.push_back(nodes->at(i));
         }
       }
   }
 
   AddObservers();
   UpdateGui();
 }
 
 void QmitkFiberProcessingView::OnDrawPolygon()
 {
   mitk::PlanarPolygon::Pointer figure = mitk::PlanarPolygon::New();
   figure->ClosedOn();
   this->AddFigureToDataStorage(figure, QString("Polygon%1").arg(++m_PolygonCounter));
 }
 
 void QmitkFiberProcessingView::OnDrawCircle()
 {
   mitk::PlanarCircle::Pointer figure = mitk::PlanarCircle::New();
   this->AddFigureToDataStorage(figure, QString("Circle%1").arg(++m_CircleCounter));
 }
 
 void QmitkFiberProcessingView::AddFigureToDataStorage(mitk::PlanarFigure* figure, const QString& name, const char *, mitk::BaseProperty* )
 {
   // initialize figure's geometry with empty geometry
   mitk::PlaneGeometry::Pointer emptygeometry = mitk::PlaneGeometry::New();
   figure->SetPlaneGeometry( emptygeometry );
 
   //set desired data to DataNode where Planarfigure is stored
   mitk::DataNode::Pointer newNode = mitk::DataNode::New();
   newNode->SetName(name.toStdString());
   newNode->SetData(figure);
   newNode->SetBoolProperty("planarfigure.3drendering", true);
   newNode->SetBoolProperty("planarfigure.3drendering.fill", true);
 
   mitk::PlanarFigureInteractor::Pointer figureInteractor = dynamic_cast<mitk::PlanarFigureInteractor*>(newNode->GetDataInteractor().GetPointer());
   if(figureInteractor.IsNull())
   {
     figureInteractor = mitk::PlanarFigureInteractor::New();
     us::Module* planarFigureModule = us::ModuleRegistry::GetModule( "MitkPlanarFigure" );
     figureInteractor->LoadStateMachine("PlanarFigureInteraction.xml", planarFigureModule );
     figureInteractor->SetEventConfig( "PlanarFigureConfig.xml", planarFigureModule );
     figureInteractor->SetDataNode(newNode);
   }
 
   // figure drawn on the topmost layer / image
   GetDataStorage()->Add(newNode );
 
   RemoveObservers();
   for(unsigned int i = 0; i < m_SelectedPF.size(); i++)
     m_SelectedPF[i]->SetSelected(false);
 
   newNode->SetSelected(true);
   m_SelectedPF.clear();
   m_SelectedPF.push_back(newNode);
   AddObservers();
   UpdateGui();
 }
 
 void QmitkFiberProcessingView::ExtractWithPlanarFigure(bool interactive)
 {
   if ( m_SelectedFB.empty() || m_SelectedPF.empty() ){
     QMessageBox::information( nullptr, "Warning", "No fibe bundle selected!");
     return;
   }
 
   try
   {
     std::vector<mitk::DataNode::Pointer> fiberBundles = m_SelectedFB;
     mitk::DataNode::Pointer planarFigure = m_SelectedPF.at(0);
     for (unsigned int i=0; i<fiberBundles.size(); i++)
     {
       mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(fiberBundles.at(i)->GetData());
       mitk::FiberBundle::Pointer extFB = fib->ExtractFiberSubset(planarFigure, GetDataStorage());
 
       if (interactive && m_Controls->m_InteractiveBox->isChecked())
       {
         if (m_InteractiveNode.IsNull())
         {
           m_InteractiveNode = mitk::DataNode::New();
           QString name("Interactive");
           m_InteractiveNode->SetName(name.toStdString());
           GetDataStorage()->Add(m_InteractiveNode);
         }
         float op = 5.0/sqrt(fib->GetNumFibers());
         float currentOp = 0;
         fiberBundles.at(i)->GetFloatProperty("opacity", currentOp);
 
         if (currentOp!=op)
         {
           fib->SetFiberColors(255, 255, 255);
           fiberBundles.at(i)->SetFloatProperty("opacity", op);
           fiberBundles.at(i)->SetBoolProperty("Fiber2DfadeEFX", false);
         }
         m_InteractiveNode->SetData(extFB);
       }
       else
       {
         if (extFB->GetNumFibers()<=0)
         {
           QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers.");
           continue;
         }
 
         mitk::DataNode::Pointer node;
         node = mitk::DataNode::New();
         node->SetData(extFB);
         QString name(fiberBundles.at(i)->GetName().c_str());
         name += "*";
         node->SetName(name.toStdString());
         fiberBundles.at(i)->SetVisibility(false);
         GetDataStorage()->Add(node);
       }
     }
   }
   catch(const std::out_of_range& )
   {
     QMessageBox::warning( nullptr, "Fiber extraction failed", "Did you only create the planar figure, using the circle or polygon button, but forgot to actually place it in the image afterwards? \nAfter creating a planar figure, simply left-click at the desired position in the image or on the tractogram to place it.");
   }
 }
 
 void QmitkFiberProcessingView::GenerateAndComposite()
 {
   mitk::PlanarFigureComposite::Pointer PFCAnd = mitk::PlanarFigureComposite::New();
   PFCAnd->setOperationType(mitk::PlanarFigureComposite::AND);
 
   mitk::DataNode::Pointer newPFCNode;
   newPFCNode = mitk::DataNode::New();
   newPFCNode->SetName("AND");
   newPFCNode->SetData(PFCAnd);
 
   AddCompositeToDatastorage(newPFCNode, m_SelectedPF);
   RemoveObservers();
   m_SelectedPF.clear();
   m_SelectedPF.push_back(newPFCNode);
   AddObservers();
   UpdateGui();
 }
 
 void QmitkFiberProcessingView::GenerateOrComposite()
 {
   mitk::PlanarFigureComposite::Pointer PFCOr = mitk::PlanarFigureComposite::New();
   PFCOr->setOperationType(mitk::PlanarFigureComposite::OR);
 
   mitk::DataNode::Pointer newPFCNode;
   newPFCNode = mitk::DataNode::New();
   newPFCNode->SetName("OR");
   newPFCNode->SetData(PFCOr);
 
   RemoveObservers();
   AddCompositeToDatastorage(newPFCNode, m_SelectedPF);
   m_SelectedPF.clear();
   m_SelectedPF.push_back(newPFCNode);
   UpdateGui();
 }
 
 void QmitkFiberProcessingView::GenerateNotComposite()
 {
   mitk::PlanarFigureComposite::Pointer PFCNot = mitk::PlanarFigureComposite::New();
   PFCNot->setOperationType(mitk::PlanarFigureComposite::NOT);
 
   mitk::DataNode::Pointer newPFCNode;
   newPFCNode = mitk::DataNode::New();
   newPFCNode->SetName("NOT");
   newPFCNode->SetData(PFCNot);
 
   RemoveObservers();
   AddCompositeToDatastorage(newPFCNode, m_SelectedPF);
   m_SelectedPF.clear();
   m_SelectedPF.push_back(newPFCNode);
   AddObservers();
   UpdateGui();
 }
 
 void QmitkFiberProcessingView::AddCompositeToDatastorage(mitk::DataNode::Pointer pfc, std::vector<mitk::DataNode::Pointer> children, mitk::DataNode::Pointer parentNode )
 {
   pfc->SetSelected(true);
   if (parentNode.IsNotNull())
     GetDataStorage()->Add(pfc, parentNode);
   else
     GetDataStorage()->Add(pfc);
 
   for (auto child : children)
   {
     if (dynamic_cast<PlanarFigure*>(child->GetData()))
     {
       mitk::DataNode::Pointer newChild;
       newChild = mitk::DataNode::New();
       newChild->SetData(dynamic_cast<PlanarFigure*>(child->GetData()));
       newChild->SetName( child->GetName() );
       newChild->SetBoolProperty("planarfigure.3drendering", true);
       newChild->SetBoolProperty("planarfigure.3drendering.fill", true);
 
       GetDataStorage()->Add(newChild, pfc);
       GetDataStorage()->Remove(child);
     }
     else if (dynamic_cast<PlanarFigureComposite*>(child->GetData()))
     {
       mitk::DataNode::Pointer newChild;
       newChild = mitk::DataNode::New();
       newChild->SetData(dynamic_cast<PlanarFigureComposite*>(child->GetData()));
       newChild->SetName( child->GetName() );
 
       std::vector< mitk::DataNode::Pointer > grandChildVector;
       mitk::DataStorage::SetOfObjects::ConstPointer grandchildren = GetDataStorage()->GetDerivations(child);
       for( mitk::DataStorage::SetOfObjects::const_iterator it = grandchildren->begin(); it != grandchildren->end(); ++it )
         grandChildVector.push_back(*it);
 
       AddCompositeToDatastorage(newChild, grandChildVector, pfc);
       GetDataStorage()->Remove(child);
     }
   }
   UpdateGui();
 }
 
 void QmitkFiberProcessingView::CopyBundles()
 {
   if ( m_SelectedFB.empty() ){
     QMessageBox::information( nullptr, "Warning", "Select at least one fiber bundle!");
     MITK_WARN("QmitkFiberProcessingView") << "Select at least one fiber bundle!";
     return;
   }
 
   for (auto node : m_SelectedFB)
   {
     mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
     mitk::FiberBundle::Pointer newFib = fib->GetDeepCopy();
 
     node->SetVisibility(false);
     QString name("");
     name += QString(m_SelectedFB.at(0)->GetName().c_str());
     name += "_copy";
 
     mitk::DataNode::Pointer fbNode = mitk::DataNode::New();
     fbNode->SetData(newFib);
     fbNode->SetName(name.toStdString());
     fbNode->SetVisibility(true);
     GetDataStorage()->Add(fbNode);
   }
   UpdateGui();
 }
 
 void QmitkFiberProcessingView::JoinBundles()
 {
   if ( m_SelectedFB.size()<2 ){
     QMessageBox::information( nullptr, "Warning", "Select at least two fiber bundles!");
     MITK_WARN("QmitkFiberProcessingView") << "Select at least two fiber bundles!";
     return;
   }
 
   m_SelectedFB.at(0)->SetVisibility(false);
   mitk::FiberBundle::Pointer newBundle = dynamic_cast<mitk::FiberBundle*>(m_SelectedFB.at(0)->GetData());
 
   std::vector< mitk::FiberBundle::Pointer > tractograms;
   for (unsigned int i=1; i<m_SelectedFB.size(); i++)
   {
     m_SelectedFB.at(i)->SetVisibility(false);
     tractograms.push_back(dynamic_cast<mitk::FiberBundle*>(m_SelectedFB.at(i)->GetData()));
   }
   newBundle = newBundle->AddBundles(tractograms);
 
   mitk::DataNode::Pointer fbNode = mitk::DataNode::New();
   fbNode->SetData(newBundle);
   fbNode->SetName("Joined_Tractograms");
   fbNode->SetVisibility(true);
   GetDataStorage()->Add(fbNode);
   UpdateGui();
 }
 
 void QmitkFiberProcessingView::SubstractBundles()
 {
   if ( m_SelectedFB.size()<2 ){
     QMessageBox::information( nullptr, "Warning", "Select at least two fiber bundles!");
     MITK_WARN("QmitkFiberProcessingView") << "Select at least two fiber bundles!";
     return;
   }
 
   mitk::FiberBundle::Pointer newBundle = dynamic_cast<mitk::FiberBundle*>(m_SelectedFB.at(0)->GetData());
   m_SelectedFB.at(0)->SetVisibility(false);
   QString name("");
   name += QString(m_SelectedFB.at(0)->GetName().c_str());
   for (unsigned int i=1; i<m_SelectedFB.size(); i++)
   {
     newBundle = newBundle->SubtractBundle(dynamic_cast<mitk::FiberBundle*>(m_SelectedFB.at(i)->GetData()));
     if (newBundle.IsNull())
       break;
     name += "-"+QString(m_SelectedFB.at(i)->GetName().c_str());
     m_SelectedFB.at(i)->SetVisibility(false);
   }
   if (newBundle.IsNull())
   {
     QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers. Did you select the fiber bundles in the correct order? X-Y is not equal to Y-X!");
     return;
   }
 
   mitk::DataNode::Pointer fbNode = mitk::DataNode::New();
   fbNode->SetData(newBundle);
   fbNode->SetName(name.toStdString());
   fbNode->SetVisibility(true);
   GetDataStorage()->Add(fbNode);
   UpdateGui();
 }
 
 void QmitkFiberProcessingView::ResampleSelectedBundlesSpline()
 {
   double factor = this->m_Controls->m_SmoothFibersBox->value();
   for (auto node : m_SelectedFB)
   {
     mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
     fib->ResampleSpline(factor);
   }
   RenderingManager::GetInstance()->RequestUpdateAll();
 }
 
 void QmitkFiberProcessingView::ResampleSelectedBundlesLinear()
 {
   double factor = this->m_Controls->m_SmoothFibersBox->value();
   for (auto node : m_SelectedFB)
   {
     mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
     fib->ResampleLinear(factor);
   }
   RenderingManager::GetInstance()->RequestUpdateAll();
 }
 
 void QmitkFiberProcessingView::CompressSelectedBundles()
 {
   double factor = this->m_Controls->m_ErrorThresholdBox->value();
   for (auto node : m_SelectedFB)
   {
     mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
     fib->Compress(factor);
     fib->ColorFibersByOrientation();
   }
   RenderingManager::GetInstance()->RequestUpdateAll();
 }
 
 void QmitkFiberProcessingView::DoImageColorCoding()
 {
   if (m_Controls->m_ColorMapBox->GetSelectedNode().IsNull())
   {
     QMessageBox::information(nullptr, "Bundle coloring aborted:", "No image providing the scalar values for coloring the selected bundle available.");
     return;
   }
 
   for (auto node : m_SelectedFB)
   {
     mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
     fib->ColorFibersByScalarMap(dynamic_cast<mitk::Image*>(m_Controls->m_ColorMapBox->GetSelectedNode()->GetData()), m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked());
   }
 
   if (auto renderWindowPart = this->GetRenderWindowPart())
   {
     renderWindowPart->RequestUpdate();
   }
 }
 
 
 void QmitkFiberProcessingView::DoCurvatureColorCoding()
 {
   for (auto node : m_SelectedFB)
   {
     mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
     fib->ColorFibersByCurvature(m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked());
   }
 
   if (auto renderWindowPart = this->GetRenderWindowPart())
   {
     renderWindowPart->RequestUpdate();
   }
 }
 
 void QmitkFiberProcessingView::DoLengthColorCoding()
 {
   for (auto node : m_SelectedFB)
   {
     mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
     fib->ColorFibersByLength(m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked());
   }
 
   if (auto renderWindowPart = this->GetRenderWindowPart())
   {
     renderWindowPart->RequestUpdate();
   }
 }
 
 void QmitkFiberProcessingView::DoWeightColorCoding()
 {
   for (auto node : m_SelectedFB)
   {
     mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
     fib->ColorFibersByFiberWeights(m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked());
   }
 
   if (auto renderWindowPart = this->GetRenderWindowPart())
   {
     renderWindowPart->RequestUpdate();
   }
 }
 
 void QmitkFiberProcessingView::MirrorFibers()
 {
   unsigned int axis = this->m_Controls->m_MirrorFibersBox->currentIndex();
   for (auto node : m_SelectedFB)
   {
     mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
-    if (m_SelectedImage.IsNotNull())
-      fib->SetTrackVisHeader(m_SelectedImage->GetGeometry());
     fib->MirrorFibers(axis);
   }
 
   for (auto surf : m_SelectedSurfaces)
   {
     vtkSmartPointer<vtkPolyData> poly = surf->GetVtkPolyData();
     vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
 
     for (int i=0; i<poly->GetNumberOfPoints(); i++)
     {
       double* point = poly->GetPoint(i);
       point[axis] *= -1;
       vtkNewPoints->InsertNextPoint(point);
     }
     poly->SetPoints(vtkNewPoints);
     surf->CalculateBoundingBox();
   }
 
   if (auto renderWindowPart = this->GetRenderWindowPart())
   {
     renderWindowPart->RequestUpdate();
   }
 }
diff --git a/README.md b/README.md
index 74442a3..129a285 100644
--- a/README.md
+++ b/README.md
@@ -1,246 +1,251 @@
 MITK Diffusion
 ==============
 
 Copyright © German Cancer Research Center (DKFZ), [Division of Medical Image Computing (MIC)](https://www.dkfz.de/en/mic/index.php). Please make sure that your usage of this code is in compliance with the code [license](https://github.com/MIC-DKFZ/MITK-Diffusion/blob/master/LICENSE.txt).
 
 The MITK Diffusion application [[1,2]](#References) offers a selection of image analysis algorithms for the processing of diffusion-weighted MR images. It encompasses the research of the Division of Medical Image Computing at the German Cancer Research Center (DKFZ).
 
 * [Downloads](#Downloads)
 * [Requirements](#Requirements)
 * [Features](#Features)
 * [Related Links](#Related-Links)
 * [Image Gallery](#Image-Gallery)
 * [Building MITK Diffusion from source](#Building-MITK-Diffusion-from-source)
 * [User Manual](http://docs.mitk.org/nightly/org_mitk_gui_qt_diffusionimaging.html)
 * [Report a Bug](https://phabricator.mitk.org/maniphest/task/edit/form/29/)
 * [References](#References)
 * [Contact](#Contact)
 
 ## Downloads
 
 **Nightly Ubuntu and Windows installers: [ftp://ftp.dkfz-heidelberg.de/outgoing/MitkDiffusion](https://bit.ly/2S1QfC8)**
 
 Please also have a look at the [requirements](#Requirements) for running MITK Diffusion with all its features successfully!
 
 The installers come as executable setup wizards that install MITK Diffusion on your system or alternatively as simple .tar.gz or .zip archive where you can execute MITK Diffusion and the command line apps "manually". Should there be no new installer for a while, please [contact](#Contact) us and report the issue.
 
 If you encounter any bugs, please report them in our [bugtracking](https://phabricator.mitk.org/maniphest/task/edit/form/29/) system or use the [MITK-users mailing list](http://mitk.org/wiki/MITK_Mailinglist). We are grateful for any feedback!
 
 
 ### Requirements
 
 * **For Ubuntu users**:
     * Install Python 3.X: `sudo apt install python3 python3-pip`
     * Download Python requirements file: [PythonRequirements.txt](https://github.com/MIC-DKFZ/MITK-Diffusion/tree/master/Modules/DiffusionCore/PythonRequirements.txt)
     * Install Python requirements: `pip3 install -r PythonRequirements.txt`
     * If your are behind a proxy use `pip3 --proxy <proxy> install -r PythonRequirements.txt`
     * **Requirements for all deep-learning based functionalities**:
         * Affected functionalities:
             * Brain extraction
             * TractSeg
         * Pytorch: https://pytorch.org/ (version 0.4.0)
 
 * **For Windows users**:
     * MITK Diffusion requires the Microsoft Visual C++ 2017 Redistributable to be installed on the system. The MITK Diffusion installer automatically installs this redistributable for you if not already present on the system, but it needs administrative privileges to do so. So to install the redistributable, **run the MITK Diffusion installer as administrator**.
-    * **Python modules are disabled under Windows!**
+    * Install python 3 from https://www.python.org/downloads/windows/
+    * **Requirements for all deep-learning based functionalities**:
+        * Affected functionalities:
+            * Brain extraction
+            * TractSeg
+        * Pytorch: https://pytorch.org/ (version 0.4.0)
 
 ## Features
 
 **Support for most established image formats**
 * Images: DICOM, NIFTI, NRRD (peak and SH images compatible with MRtrix)
 * Tractograms: fib/vtk, tck and trk.
 
 **Image preprocessing**
 * Registration
 * Head-motion correction
 * Denoising
 * Skull stripping and brain mask segmentation (Linux only)
 * Resampling, cropping, flipping and merging
 * Header modifications
 * Single volume extraction
 
 **Diffusion gradient/b-value processing**
 * b-value rounding
 * Gradient direction flipping
 * Gradient direction subsampling
 * Averaging of gradient directions/volumes
 * Gradient direction and b-value visualization
 
 **ODF reconstruction and signal modelling**
 * Tensor and Q-ball reconstruction
 * Other reconstructions via Dipy wrapping (CSD, 3D SHORE, SFM) (Linux only)
 * ODF peak calculation
 * MRtrix or camino results can be imported
 
 **Quantification of diffusion-weighted/tensor/ODF images**
 * Intravoxel Incoherent Motion (IVIM) and diffusion kurtosis analysis
 * Calculation of many other derived indices such as ADC, MD, GFA, FA, RA, AD, RD
 * Image statistics
 
 **Segmentation**
 * Automatic white matter bundle segmentation (TractSeg) [[3]](#References) (Linux only)
 * Automatic brain mask segmentation (Linux only)
 * Manual image segmentation and operations on segmentations
 * SOON: automatic brain tissue segmentation
 
 **Fiber tractography**
 * Global tractography [[4]](#References)
 * Streamline tractography
     * Interactive (similar to [[5]](#References)) or seed image based
     * Deterministic or probabilistic
     * Peak, ODF, tensor and raw dMRI based. The latter one in conjunction with machine learning based tractography [[6]](#References)
     * Various possibilities for anatomical constraints.
     * Tractography priors in form of additional peak images, e.g. obtained using TractSeg
 
 **Fiber processing**
 * Tract dissection (parcellation or ROI based)
 * Tract filtering by
     * length
     * curvature
     * direction
     * weight
     * density
 * Tract resampling and compression
 * Tract transformation
     * Mirroring
     * Rotating and translating
     * Registration (apply transform of previously performed image registration)
 * Tract coloring
     * Curvature
     * Length
     * Weight
     * Scalar map (e.g. FA)
 * Other operations
     * Join
     * Subtract
     * Copy
 * Fiber clustering [[7]](#References)
 * Fiber fitting and weighting similar to SIFT2 and LiFE [[8,9]](#References)
 * Principal direction extraction (fibers --> peaks)
 * Tract derived images:
     * Tract density images
     * Tract endpoint images
     * Tract envelopes
 
 **Fiberfox dMRI simulations** [[10]](#References)
 * Multi-compartment signal modeling
 * Simulation of the k-space acquisition including
     * Compartment specific relaxation effects
     * Artifacts such as noise, spikes, ghosts, aliasing, distortions, signal drift, head motion, eddy currents and Gibbs ringing
     * Definition of important acquisition parameters such as bvalues and gradient directions, TE, TR, dwell time, partial Fourier, ...
 * Manual definition of fiber configurations, e.g. for evaluation purposes
 * Automatic generation of random fiber configurations
 
 **Other features**
 * Brain network statistics and visualization (connectomics)
 * Interactive Python console (Linux only)
 * Integrated screenshot maker
 * Command line tools for most functionalities
 
 ## Related Links
 * Great python package for logging your (MITK) command line experiments:
     * https://github.com/MIC-DKFZ/cmdint
     * `pip3 install cmdint`
 * TractSeg reference data of 72 semiautomatically defined bundles in 105 HCP subjects: https://zenodo.org/record/1285152
 * TractSeg python package: https://github.com/MIC-DKFZ/TractSeg
 * Simulated dMRI images and ground truth of random fiber phantoms in various configurations: https://doi.org/10.5281/zenodo.2533250
 * ISMRM 2015 Tractography Challenge Data: https://doi.org/10.5281/zenodo.572345 & https://doi.org/10.5281/zenodo.1007149
 
 
 ## Image Gallery
 
 ![](http://mitk.org/images/8/8f/MitkDiffusion.png) Screenshot of the MITK Diffusion Welcome Screen
 <br/><br/><br/>
 
 ![](http://mitk.org/images/0/09/ScalarMaps.png) Scalar map visualization
 <br/><br/><br/>
 
 ![](http://mitk.org/images/3/3b/Data_Tensors.png) Tensor Visualization
 <br/><br/><br/>
 
 ![](http://mitk.org/images/5/5c/Data_ODF.png) ODF visualization
 <br/><br/><br/>
 
 ![](http://mitk.org/images/7/73/Data_Peaks.png) Peak visualization (uniform white coloring)
 <br/><br/><br/>
 
 ![](http://mitk.org/images/6/68/StreamlineTractography.png) Interactive tractography in MITK Diffusion. The tractogram updates automatically on parameter change and movement of the spherical seed region.
 <br/><br/><br/>
 
 ![](http://mitk.org/images/3/3a/Extraction_1.png) Tract dissection using manually drawn ROIs.
 <br/><br/><br/>
 
 ![](http://mitk.org/images/f/f5/FiberFit.png) Automatic streamline weighting (similar to SIFT2 or LiFE)
 <br/><br/><br/>
 
 ![](http://mitk.org/images/3/33/Fiberfox.png) Illustration of the dMRI phantom simulation process using Fiberfox.
 <br/><br/><br/>
 
 <img src="http://mitk.org/images/e/ec/Eddy.gif" width="250" height="200" /><img src="http://mitk.org/images/a/ad/Spike_and_motion.gif" width="250" height="200" /><img src="http://mitk.org/images/8/8e/Drift.gif" width="250" height="200" />
 <img src="http://mitk.org/images/0/08/645551_axial_motion_ghost_eddy.gif" width="250" height="200" /><img src="http://mitk.org/images/5/55/Ringing.gif" width="250" height="200" /><img src="http://mitk.org/images/3/37/683256_Distortions_Sagittal.gif" width="250" height="200" />\
 Illustration of simulated dMRI images with various artifacts (a bit excessive for illustration purposes): eddy current distortions (1), motion and spike (2), intensity drift (3), motion, eddy and noise (4), ringing (5), B0 inhomogeneity distortions (6), from left to right.
 <br/><br/><br/>
 
 ![](http://mitk.org/images/0/08/RandomFibers_Example.png)
 Automatically generated random fiber configuration for Fiberfox simulations.
 <br/><br/><br/>
 
 
 ## Building MITK Diffusion from source
 * Install [Qt](https://www.qt.io/ Qt) on your system (>= 5.11.1).
 * Clone MITK from [github](https://github.com/MIC-DKFZ/MITK-Diffusion.git) using [Git version control](https://git-scm.com/).
 * Clone MITK Diffusion from [github](https://github.com/MITK/MITK.git).
 * Configure the MITK Superbuild using [CMake](https://cmake.org/) (>= 3.14.5).
     * Choose the MITK source code directory and an empty binary directory.
     * Click "Configure".
     * Set the option MITK_EXTENSION_DIRS to "/path/to/my/mitk-diffusion-repository".
     * Click "Configure".
     * Set the option MITK_BUILD_CONFIGURATION to "DiffusionRelease".
     * Click "Generate".
     * macOS specifics:
         * Use python 3.**6**, since python 3.**7** leads to build errors on macOS.
         * The cmake variables for python 3 might need to be set manually. It is probably enough to specify PYTHON_EXECUTABLE.
         * Openmp needs to be installed manually since it is not included in apple clang anymore: "brew install libomp" should do the trick. It might be necessary to set the corresponding make variables later in the MITK build manually:
             * OpenMP_CXX_FLAGS: -Xpreprocessor -fopenmp -I"/path/to/python3/includes/"
             * OpenMP_C_FLAGS: -Xpreprocessor -fopenmp -I"/path/to/python3/includes/"
             * OpenMPCXX_LIB_NAMES: libomp
             * OpenMPC_LIB_NAMES: libomp
             * OpenMP_libomp_LIBRARY: /path/to/libomp.dylib
 * Build the project
     * Linux/maxOS: Open a console window, navigate to the build folder and type "make -j8" (optionally supply the number threads to be used for a parallel build with -j).
     * Windows (requires visual studio): Open the MITK Superbuild solution file and build all projects.
 * The build may take some time and should yield the binaries in "your_build_folder/MITK-build/bin"
 
 More detailed build instructions can be found in the [documentation](http://docs.mitk.org/nightly/BuildInstructionsPage.html).
 
 
 Continuous integration: http://cdash.mitk.org/index.php?project=MITK&display=project
 
 
 ## References
 
 All publications of the Division of Medical Image Computing can be found [https://www.dkfz.de/en/mic/publications/ here].
 
 [1] Fritzsche, Klaus H., Peter F. Neher, Ignaz Reicht, Thomas van Bruggen, Caspar Goch, Marco Reisert, Marco Nolden, et al. “MITK Diffusion Imaging.” Methods of Information in Medicine 51, no. 5 (2012): 441.
 
 [2] Fritzsche, K., and H.-P. Meinzer. “MITK-DI A New Diffusion Imaging Component for MITK.” In Bildverarbeitung Für Die Medizin, n.d.
 
 [3] Wasserthal, Jakob, Peter Neher, and Klaus H. Maier-Hein. “TractSeg - Fast and Accurate White Matter Tract Segmentation.” NeuroImage 183 (August 4, 2018): 239–53.
 
 [4] Neher, P. F., B. Stieltjes, M. Reisert, I. Reicht, H.P. Meinzer, and K. Maier-Hein. “MITK Global Tractography.” In SPIE Medical Imaging: Image Processing, 2012.
 
 [5] Chamberland, M., K. Whittingstall, D. Fortin, D. Mathieu, und M. Descoteaux. „Real-time multi-peak tractography for instantaneous connectivity display“. Front Neuroinform 8 (2014): 59. doi:10.3389/fninf.2014.00059.
 
 [6] Neher, Peter F., Marc-Alexandre Côté, Jean-Christophe Houde, Maxime Descoteaux, and Klaus H. Maier-Hein. “Fiber Tractography Using Machine Learning.” NeuroImage. Accessed July 17, 2017. doi:10.1016/j.neuroimage.2017.07.028.
 
 [7] Garyfallidis, Eleftherios, Matthew Brett, Marta Morgado Correia, Guy B. Williams, and Ian Nimmo-Smith. “QuickBundles, a Method for Tractography Simplification.” Frontiers in Neuroscience 6 (2012).
 
 [8] Smith, Robert E., Jacques-Donald Tournier, Fernando Calamante, and Alan Connelly. “SIFT2: Enabling Dense Quantitative Assessment of Brain White Matter Connectivity Using Streamlines Tractography.” NeuroImage 119, no. Supplement C (October 1, 2015): 338–51.
 
 [9] Pestilli, Franco, Jason D. Yeatman, Ariel Rokem, Kendrick N. Kay, and Brian A. Wandell. “Evaluation and Statistical Inference for Human Connectomes.” Nature Methods 11, no. 10 (October 2014): 1058–63.
 
 [10] Neher, Peter F., Frederik B. Laun, Bram Stieltjes, and Klaus H. Maier-Hein. “Fiberfox: Facilitating the Creation of Realistic White Matter Software Phantoms.” Magnetic Resonance in Medicine 72, no. 5 (November 2014): 1460–70. doi:10.1002/mrm.25045.
 
 
 ## Contact
 If you have questions about the application or if you would like to give us feedback, don't hesitate to contact us using [our mailing list](http://mitk.org/wiki/MITK_Mailinglist) or, for questions that are of no interest for the community, [directly](https://www.dkfz.de/en/mic/team/people/Peter_Neher.html).