diff --git a/Modules/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp b/Modules/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp
index 130f623..86abf58 100644
--- a/Modules/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp
+++ b/Modules/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp
@@ -1,992 +1,997 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center.
 
 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 <ctime>
 #include <cstdio>
 #include <cstdlib>
 
 #include <omp.h>
 #include "itkStreamlineTrackingFilter.h"
 #include <itkImageRegionConstIterator.h>
 #include <itkImageRegionConstIteratorWithIndex.h>
 #include <itkImageRegionIterator.h>
 #include "itkPointShell.h"
 #include <itkRescaleIntensityImageFilter.h>
 #include <mitkLexicalCast.h>
 #include <TrackingHandlers/mitkTrackingDataHandler.h>
 #include <TrackingHandlers/mitkTrackingHandlerOdf.h>
 #include <TrackingHandlers/mitkTrackingHandlerPeaks.h>
 #include <TrackingHandlers/mitkTrackingHandlerTensor.h>
 #include <TrackingHandlers/mitkTrackingHandlerRandomForest.h>
 #include <mitkDiffusionFunctionCollection.h>
+#include <random>
 
 namespace itk {
 
 
 StreamlineTrackingFilter
 ::StreamlineTrackingFilter()
   : m_PauseTracking(false)
   , m_AbortTracking(false)
   , m_BuildFibersFinished(false)
   , m_BuildFibersReady(0)
   , m_FiberPolyData(nullptr)
   , m_Points(nullptr)
   , m_Cells(nullptr)
   , m_StoppingRegions(nullptr)
   , m_TargetRegions(nullptr)
   , m_SeedImage(nullptr)
   , m_MaskImage(nullptr)
   , m_ExclusionRegions(nullptr)
   , m_OutputProbabilityMap(nullptr)
   , m_Verbose(true)
   , m_DemoMode(false)
   , m_CurrentTracts(0)
   , m_Progress(0)
   , m_StopTracking(false)
   , m_TrackingPriorHandler(nullptr)
 {
   this->SetNumberOfRequiredInputs(0);
 }
 
 std::string StreamlineTrackingFilter::GetStatusText()
 {
   std::string status = "Seedpoints processed: " + boost::lexical_cast<std::string>(m_Progress) + "/" + boost::lexical_cast<std::string>(m_SeedPoints.size());
   if (m_SeedPoints.size()>0)
     status += " (" + boost::lexical_cast<std::string>(100*m_Progress/m_SeedPoints.size()) + "%)";
   if (m_Parameters->m_MaxNumFibers>0)
     status += "\nFibers accepted: " + boost::lexical_cast<std::string>(m_CurrentTracts) + "/" + boost::lexical_cast<std::string>(m_Parameters->m_MaxNumFibers);
   else
     status += "\nFibers accepted: " + boost::lexical_cast<std::string>(m_CurrentTracts);
 
   return status;
 }
 
 void StreamlineTrackingFilter::BeforeTracking()
 {
   m_StopTracking = false;
   m_TrackingHandler->SetParameters(m_Parameters);
   m_TrackingHandler->InitForTracking();
   m_FiberPolyData = PolyDataType::New();
   m_Points = vtkSmartPointer< vtkPoints >::New();
   m_Cells = vtkSmartPointer< vtkCellArray >::New();
 
   if (m_TrackingPriorHandler!=nullptr)
   {
     m_TrackingPriorHandler->InitForTracking();
   }
 
   auto imageSpacing = m_TrackingHandler->GetSpacing();
   if (m_Parameters->m_OutputProbMap)
   {
     m_OutputProbabilityMap = ItkDoubleImgType::New();
     m_OutputProbabilityMap->SetSpacing(imageSpacing);
     m_OutputProbabilityMap->SetOrigin(m_TrackingHandler->GetOrigin());
     m_OutputProbabilityMap->SetDirection(m_TrackingHandler->GetDirection());
     m_OutputProbabilityMap->SetRegions(m_TrackingHandler->GetLargestPossibleRegion());
     m_OutputProbabilityMap->Allocate();
     m_OutputProbabilityMap->FillBuffer(0);
   }
 
   m_MaskInterpolator = itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::New();
   m_StopInterpolator = itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::New();
   m_SeedInterpolator = itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::New();
   m_TargetInterpolator = itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::New();
   m_ExclusionInterpolator = itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::New();
 
   if (m_StoppingRegions.IsNull())
   {
     m_StoppingRegions = ItkFloatImgType::New();
     m_StoppingRegions->SetSpacing( imageSpacing );
     m_StoppingRegions->SetOrigin( m_TrackingHandler->GetOrigin() );
     m_StoppingRegions->SetDirection( m_TrackingHandler->GetDirection() );
     m_StoppingRegions->SetRegions( m_TrackingHandler->GetLargestPossibleRegion() );
     m_StoppingRegions->Allocate();
     m_StoppingRegions->FillBuffer(0);
   }
   else
     std::cout << "StreamlineTracking - Using stopping region image" << std::endl;
   m_StopInterpolator->SetInputImage(m_StoppingRegions);
 
   if (m_ExclusionRegions.IsNotNull())
   {
     std::cout << "StreamlineTracking - Using exclusion region image" << std::endl;
     m_ExclusionInterpolator->SetInputImage(m_ExclusionRegions);
   }
 
   if (m_TargetRegions.IsNull())
   {
     m_TargetImageSet = false;
     m_TargetRegions = ItkFloatImgType::New();
     m_TargetRegions->SetSpacing( imageSpacing );
     m_TargetRegions->SetOrigin( m_TrackingHandler->GetOrigin() );
     m_TargetRegions->SetDirection( m_TrackingHandler->GetDirection() );
     m_TargetRegions->SetRegions( m_TrackingHandler->GetLargestPossibleRegion() );
     m_TargetRegions->Allocate();
     m_TargetRegions->FillBuffer(1);
   }
   else
   {
     m_TargetImageSet = true;
     m_TargetInterpolator->SetInputImage(m_TargetRegions);
     std::cout << "StreamlineTracking - Using target region image" << std::endl;
   }
 
   if (m_SeedImage.IsNull())
   {
     m_SeedImageSet = false;
     m_SeedImage = ItkFloatImgType::New();
     m_SeedImage->SetSpacing( imageSpacing );
     m_SeedImage->SetOrigin( m_TrackingHandler->GetOrigin() );
     m_SeedImage->SetDirection( m_TrackingHandler->GetDirection() );
     m_SeedImage->SetRegions( m_TrackingHandler->GetLargestPossibleRegion() );
     m_SeedImage->Allocate();
     m_SeedImage->FillBuffer(1);
   }
   else
   {
     m_SeedImageSet = true;
     std::cout << "StreamlineTracking - Using seed image" << std::endl;
   }
   m_SeedInterpolator->SetInputImage(m_SeedImage);
 
   if (m_MaskImage.IsNull())
   {
     // initialize mask image
     m_MaskImage = ItkFloatImgType::New();
     m_MaskImage->SetSpacing( imageSpacing );
     m_MaskImage->SetOrigin( m_TrackingHandler->GetOrigin() );
     m_MaskImage->SetDirection( m_TrackingHandler->GetDirection() );
     m_MaskImage->SetRegions( m_TrackingHandler->GetLargestPossibleRegion() );
     m_MaskImage->Allocate();
     m_MaskImage->FillBuffer(1);
   }
   else
     std::cout << "StreamlineTracking - Using mask image" << std::endl;
   m_MaskInterpolator->SetInputImage(m_MaskImage);
 
   // Autosettings for endpoint constraints
   if (m_Parameters->m_EpConstraints==EndpointConstraints::NONE && m_TargetImageSet && m_SeedImageSet)
   {
     MITK_INFO << "No endpoint constraint chosen but seed and target image set --> setting constraint to EPS_IN_SEED_AND_TARGET";
     m_Parameters->m_EpConstraints = EndpointConstraints::EPS_IN_SEED_AND_TARGET;
   }
   else if (m_Parameters->m_EpConstraints==EndpointConstraints::NONE && m_TargetImageSet)
   {
     MITK_INFO << "No endpoint constraint chosen but target image set --> setting constraint to EPS_IN_TARGET";
     m_Parameters->m_EpConstraints = EndpointConstraints::EPS_IN_TARGET;
   }
 
   // Check if endpoint constraints are valid
   FiberType test_fib; itk::Point<float> p; p.Fill(0);
   test_fib.push_back(p); test_fib.push_back(p);
   IsValidFiber(&test_fib);
 
   if (m_SeedPoints.empty())
     GetSeedPointsFromSeedImage();
 
   m_BuildFibersReady = 0;
   m_BuildFibersFinished = false;
   m_Tractogram.clear();
   m_SamplingPointset = mitk::PointSet::New();
   m_AlternativePointset = mitk::PointSet::New();
   m_StopVotePointset = mitk::PointSet::New();
   m_StartTime = std::chrono::system_clock::now();
 
   if (m_DemoMode)
     omp_set_num_threads(1);
 
   if (m_Parameters->m_Mode==mitk::TrackingDataHandler::MODE::DETERMINISTIC)
     std::cout << "StreamlineTracking - Mode: deterministic" << std::endl;
   else if(m_Parameters->m_Mode==mitk::TrackingDataHandler::MODE::PROBABILISTIC)
   {
     std::cout << "StreamlineTracking - Mode: probabilistic" << std::endl;
     std::cout << "StreamlineTracking - Trials per seed: " << m_Parameters->m_TrialsPerSeed << std::endl;
   }
   else
     std::cout << "StreamlineTracking - Mode: ???" << std::endl;
 
   if (m_Parameters->m_EpConstraints==EndpointConstraints::NONE)
     std::cout << "StreamlineTracking - Endpoint constraint: NONE" << std::endl;
   else if (m_Parameters->m_EpConstraints==EndpointConstraints::EPS_IN_TARGET)
     std::cout << "StreamlineTracking - Endpoint constraint: EPS_IN_TARGET" << std::endl;
   else if (m_Parameters->m_EpConstraints==EndpointConstraints::EPS_IN_TARGET_LABELDIFF)
     std::cout << "StreamlineTracking - Endpoint constraint: EPS_IN_TARGET_LABELDIFF" << std::endl;
   else if (m_Parameters->m_EpConstraints==EndpointConstraints::EPS_IN_SEED_AND_TARGET)
     std::cout << "StreamlineTracking - Endpoint constraint: EPS_IN_SEED_AND_TARGET" << std::endl;
   else if (m_Parameters->m_EpConstraints==EndpointConstraints::MIN_ONE_EP_IN_TARGET)
     std::cout << "StreamlineTracking - Endpoint constraint: MIN_ONE_EP_IN_TARGET" << std::endl;
   else if (m_Parameters->m_EpConstraints==EndpointConstraints::ONE_EP_IN_TARGET)
     std::cout << "StreamlineTracking - Endpoint constraint: ONE_EP_IN_TARGET" << std::endl;
   else if (m_Parameters->m_EpConstraints==EndpointConstraints::NO_EP_IN_TARGET)
     std::cout << "StreamlineTracking - Endpoint constraint: NO_EP_IN_TARGET" << std::endl;
 
   std::cout << "StreamlineTracking - Angular threshold: " << m_Parameters->GetAngularThresholdDeg() << "°" << std::endl;
   std::cout << "StreamlineTracking - Stepsize: " << m_Parameters->GetStepSizeMm() << "mm (" << m_Parameters->GetStepSizeMm()/m_Parameters->GetMinVoxelSizeMm() << "*vox)" << std::endl;
   std::cout << "StreamlineTracking - Seeds per voxel: " << m_Parameters->m_SeedsPerVoxel << std::endl;
   std::cout << "StreamlineTracking - Max. tract length: " << m_Parameters->m_MaxTractLengthMm << "mm" << std::endl;
   std::cout << "StreamlineTracking - Min. tract length: " << m_Parameters->m_MinTractLengthMm << "mm" << std::endl;
   std::cout << "StreamlineTracking - Max. num. tracts: " << m_Parameters->m_MaxNumFibers << std::endl;
   std::cout << "StreamlineTracking - Loop check: " << m_Parameters->GetLoopCheckDeg() << "°" << std::endl;
 
   std::cout << "StreamlineTracking - Num. neighborhood samples: " << m_Parameters->m_NumSamples << std::endl;
   std::cout << "StreamlineTracking - Max. sampling distance: " << m_Parameters->GetSamplingDistanceMm() << "mm (" << m_Parameters->GetSamplingDistanceMm()/m_Parameters->GetMinVoxelSizeMm() << "*vox)" << std::endl;
   std::cout << "StreamlineTracking - Deflection modifier: " << m_Parameters->m_DeflectionMod << std::endl;
 
   std::cout << "StreamlineTracking - Use stop votes: " << m_Parameters->m_StopVotes << std::endl;
   std::cout << "StreamlineTracking - Only frontal samples: " << m_Parameters->m_OnlyForwardSamples << std::endl;
 
   if (m_TrackingPriorHandler!=nullptr)
     std::cout << "StreamlineTracking - Using directional prior for tractography (w=" << m_Parameters->m_Weight << ")" << std::endl;
 
   if (m_DemoMode)
   {
     std::cout << "StreamlineTracking - Running in demo mode";
     std::cout << "StreamlineTracking - Starting streamline tracking using 1 thread" << std::endl;
   }
   else
     std::cout << "StreamlineTracking - Starting streamline tracking using " << omp_get_max_threads() << " threads" << std::endl;
 }
 
 void StreamlineTrackingFilter::CalculateNewPosition(itk::Point<float, 3>& pos, vnl_vector_fixed<float, 3>& dir)
 {
   pos[0] += dir[0]*m_Parameters->GetStepSizeMm();
   pos[1] += dir[1]*m_Parameters->GetStepSizeMm();
   pos[2] += dir[2]*m_Parameters->GetStepSizeMm();
 }
 
 std::vector< vnl_vector_fixed<float,3> > StreamlineTrackingFilter::CreateDirections(unsigned int NPoints)
 {
   std::vector< vnl_vector_fixed<float,3> > pointshell;
 
   if (NPoints<2)
     return pointshell;
 
   std::vector< double > theta; theta.resize(NPoints);
 
   std::vector< double > phi; phi.resize(NPoints);
 
   auto C = sqrt(4*itk::Math::pi);
 
   phi[0] = 0.0;
   phi[NPoints-1] = 0.0;
 
   for(unsigned int i=0; i<NPoints; i++)
   {
     theta[i] = acos(-1.0+2.0*i/(NPoints-1.0)) - itk::Math::pi / 2.0;
     if( i>0 && i<NPoints-1)
     {
       phi[i] = (phi[i-1] + C / sqrt(NPoints*(1-(-1.0+2.0*i/(NPoints-1.0))*(-1.0+2.0*i/(NPoints-1.0)))));
       // % (2*DIST_POINTSHELL_PI);
     }
   }
 
 
   for(unsigned int i=0; i<NPoints; i++)
   {
     vnl_vector_fixed<float,3> d;
     d[0] = static_cast<float>(cos(theta[i]) * cos(phi[i]));
     d[1] = static_cast<float>(cos(theta[i]) * sin(phi[i]));
     d[2] = static_cast<float>(sin(theta[i]));
     pointshell.push_back(d);
   }
 
   return pointshell;
 }
 
 
 vnl_vector_fixed<float,3> StreamlineTrackingFilter::GetNewDirection(const itk::Point<float, 3> &pos, std::deque<vnl_vector_fixed<float, 3> >& olddirs, itk::Index<3> &oldIndex)
 {
   if (m_DemoMode)
   {
     m_SamplingPointset->Clear();
     m_AlternativePointset->Clear();
     m_StopVotePointset->Clear();
   }
   vnl_vector_fixed<float,3> direction; direction.fill(0);
 
   if (mitk::imv::IsInsideMask<float>(pos, m_Parameters->m_InterpolateRoiImages, m_MaskInterpolator) && !mitk::imv::IsInsideMask<float>(pos, m_Parameters->m_InterpolateRoiImages, m_StopInterpolator))
     direction = m_TrackingHandler->ProposeDirection(pos, olddirs, oldIndex); // get direction proposal at current streamline position
   else
     return direction;
 
   int stop_votes = 0;
   int possible_stop_votes = 0;
   if (!olddirs.empty())
   {
     vnl_vector_fixed<float,3> olddir = olddirs.back();
     std::vector< vnl_vector_fixed<float,3> > probeVecs = CreateDirections(m_Parameters->m_NumSamples);
     itk::Point<double, 3> sample_pos;
     unsigned int alternatives = 1;
     for (unsigned int i=0; i<probeVecs.size(); i++)
     {
       vnl_vector_fixed<float,3> d;
       bool is_stop_voter = false;
       if (!m_Parameters->m_FixRandomSeed && m_Parameters->m_RandomSampling)
       {
         d[0] = static_cast<float>(m_TrackingHandler->GetRandDouble(-0.5, 0.5));
         d[1] = static_cast<float>(m_TrackingHandler->GetRandDouble(-0.5, 0.5));
         d[2] = static_cast<float>(m_TrackingHandler->GetRandDouble(-0.5, 0.5));
         d.normalize();
         d *= static_cast<float>(m_TrackingHandler->GetRandDouble(0, static_cast<double>(m_Parameters->GetSamplingDistanceMm())));
       }
       else
       {
         d = probeVecs.at(i);
         float dot = dot_product(d, olddir);
         if (m_Parameters->m_StopVotes && dot>0.7f)
         {
           is_stop_voter = true;
           possible_stop_votes++;
         }
         else if (m_Parameters->m_OnlyForwardSamples && dot<0)
           continue;
         d *= m_Parameters->GetSamplingDistanceMm();
       }
 
       sample_pos[0] = pos[0] + d[0];
       sample_pos[1] = pos[1] + d[1];
       sample_pos[2] = pos[2] + d[2];
 
       vnl_vector_fixed<float,3> tempDir; tempDir.fill(0.0);
       if (mitk::imv::IsInsideMask<float>(sample_pos, m_Parameters->m_InterpolateRoiImages, m_MaskInterpolator))
         tempDir = m_TrackingHandler->ProposeDirection(sample_pos, olddirs, oldIndex); // sample neighborhood
       if (tempDir.magnitude()>static_cast<float>(mitk::eps))
       {
         direction += tempDir;
 
         if(m_DemoMode)
           m_SamplingPointset->InsertPoint(i, sample_pos);
       }
       else if (m_Parameters->m_AvoidStop && olddir.magnitude()>0.5f) // out of white matter
       {
         if (is_stop_voter)
           stop_votes++;
         if (m_DemoMode)
           m_StopVotePointset->InsertPoint(i, sample_pos);
 
         float dot = dot_product(d, olddir);
         if (dot >= 0.0f) // in front of plane defined by pos and olddir
           d = -d + 2*dot*olddir; // reflect
         else
           d = -d; // invert
 
         // look a bit further into the other direction
         sample_pos[0] = pos[0] + d[0];
         sample_pos[1] = pos[1] + d[1];
         sample_pos[2] = pos[2] + d[2];
         alternatives++;
         vnl_vector_fixed<float,3> tempDir; tempDir.fill(0.0);
         if (mitk::imv::IsInsideMask<float>(sample_pos, m_Parameters->m_InterpolateRoiImages, m_MaskInterpolator))
           tempDir = m_TrackingHandler->ProposeDirection(sample_pos, olddirs, oldIndex); // sample neighborhood
 
         if (tempDir.magnitude()>static_cast<float>(mitk::eps))  // are we back in the white matter?
         {
           direction += d * m_Parameters->m_DeflectionMod;         // go into the direction of the white matter
           direction += tempDir;  // go into the direction of the white matter direction at this location
 
           if(m_DemoMode)
             m_AlternativePointset->InsertPoint(alternatives, sample_pos);
         }
         else
         {
           if (m_DemoMode)
             m_StopVotePointset->InsertPoint(i, sample_pos);
         }
       }
       else
       {
         if (m_DemoMode)
           m_StopVotePointset->InsertPoint(i, sample_pos);
 
         if (is_stop_voter)
           stop_votes++;
       }
     }
   }
 
   bool valid = false;
   if (direction.magnitude()>0.001f && (possible_stop_votes==0 || static_cast<float>(stop_votes)/possible_stop_votes<0.5f) )
   {
     direction.normalize();
     valid = true;
   }
   else
     direction.fill(0);
 
   if (m_TrackingPriorHandler!=nullptr && (m_Parameters->m_NewDirectionsFromPrior || valid))
   {
     vnl_vector_fixed<float,3> prior = m_TrackingPriorHandler->ProposeDirection(pos, olddirs, oldIndex);
     if (prior.magnitude()>0.001f)
     {
       prior.normalize();
       if (dot_product(prior,direction)<0)
         prior *= -1;
       direction = (1.0f-m_Parameters->m_Weight) * direction + m_Parameters->m_Weight * prior;
       direction.normalize();
     }
     else if (m_Parameters->m_RestrictToPrior)
       direction.fill(0.0);
   }
 
   return direction;
 }
 
 
 float StreamlineTrackingFilter::FollowStreamline(itk::Point<float, 3> pos, vnl_vector_fixed<float,3> dir, FiberType* fib, DirectionContainer* container, float tractLength, bool front, bool &exclude)
 {
   vnl_vector_fixed<float,3> zero_dir; zero_dir.fill(0.0);
   std::deque< vnl_vector_fixed<float,3> > last_dirs;
   for (unsigned int i=0; i<m_Parameters->m_NumPreviousDirections-1; i++)
     last_dirs.push_back(zero_dir);
 
   for (int step=0; step< 5000; step++)
   {
     itk::Index<3> oldIndex;
     m_TrackingHandler->WorldToIndex(pos, oldIndex);
 
     // get new position
     CalculateNewPosition(pos, dir);
 
     if (m_ExclusionRegions.IsNotNull() && mitk::imv::IsInsideMask<float>(pos, m_Parameters->m_InterpolateRoiImages, m_ExclusionInterpolator))
     {
       exclude = true;
       return tractLength;
     }
 
     if (m_AbortTracking)
       return tractLength;
 
     // if yes, add new point to streamline
     dir.normalize();
     if (front)
     {
       fib->push_front(pos);
       container->push_front(dir);
     }
     else
     {
       fib->push_back(pos);
       container->push_back(dir);
     }
     tractLength +=  m_Parameters->GetStepSizeMm();
 
     if (m_Parameters->GetLoopCheckDeg()>=0 && CheckCurvature(container, front)>m_Parameters->GetLoopCheckDeg())
       return tractLength;
 
     if (tractLength>m_Parameters->m_MaxTractLengthMm)
       return tractLength;
 
     if (m_DemoMode && !m_Parameters->m_OutputProbMap) // CHECK: warum sind die samplingpunkte der streamline in der visualisierung immer einen schritt voras?
     {
 #pragma omp critical
       {
         m_BuildFibersReady++;
         m_Tractogram.push_back(*fib);
         BuildFibers(true);
         m_Stop = true;
 
         while (m_Stop){
         }
       }
     }
 
     last_dirs.push_back(dir);
     if (last_dirs.size()>m_Parameters->m_NumPreviousDirections)
       last_dirs.pop_front();
     dir = GetNewDirection(pos, last_dirs, oldIndex);
 
     while (m_PauseTracking){}
 
     if (dir.magnitude()<0.0001f)
       return tractLength;
   }
   return tractLength;
 }
 
 
 float StreamlineTrackingFilter::CheckCurvature(DirectionContainer* fib, bool front)
 {
   if (fib->size()<8)
     return 0;
   float m_Distance = std::max(m_Parameters->GetMinVoxelSizeMm()*4, m_Parameters->GetStepSizeMm()*8);
   float dist = 0;
 
   std::vector< vnl_vector_fixed< float, 3 > > vectors;
   vnl_vector_fixed< float, 3 > meanV; meanV.fill(0);
   float dev = 0;
 
   if (front)
   {
     int c = 0;
     while(dist<m_Distance && c<static_cast<int>(fib->size())-1)
     {
       dist += m_Parameters->GetStepSizeMm();
       vnl_vector_fixed< float, 3 > v = fib->at(static_cast<unsigned int>(c));
       if (dot_product(v,meanV)<0)
         v = -v;
       vectors.push_back(v);
       meanV += v;
       c++;
     }
   }
   else
   {
     int c = static_cast<int>(fib->size())-1;
     while(dist<m_Distance && c>=0)
     {
       dist += m_Parameters->GetStepSizeMm();
       vnl_vector_fixed< float, 3 > v = fib->at(static_cast<unsigned int>(c));
       if (dot_product(v,meanV)<0)
         v = -v;
       vectors.push_back(v);
       meanV += v;
       c--;
     }
   }
   meanV.normalize();
 
   for (unsigned int c=0; c<vectors.size(); c++)
   {
     float angle = std::fabs(dot_product(meanV, vectors.at(c)));
     if (angle>1.0f)
       angle = 1.0;
     dev += acos(angle)*180.0f/static_cast<float>(itk::Math::pi);
   }
   if (vectors.size()>0)
     dev /= vectors.size();
 
   return dev;
 }
 
 std::shared_ptr<mitk::StreamlineTractographyParameters> StreamlineTrackingFilter::GetParameters() const
 {
   return m_Parameters;
 }
 
 void StreamlineTrackingFilter::SetParameters(std::shared_ptr< mitk::StreamlineTractographyParameters > Parameters)
 {
   m_Parameters = Parameters;
 }
 
 void StreamlineTrackingFilter::SetTrackingPriorHandler(mitk::TrackingDataHandler *TrackingPriorHandler)
 {
   m_TrackingPriorHandler = TrackingPriorHandler;
 }
 
 void StreamlineTrackingFilter::GetSeedPointsFromSeedImage()
 {
   MITK_INFO << "StreamlineTracking - Calculating seed points.";
   m_SeedPoints.clear();
 
   typedef ImageRegionConstIterator< ItkFloatImgType >     MaskIteratorType;
   MaskIteratorType    sit(m_SeedImage, m_SeedImage->GetLargestPossibleRegion());
   sit.GoToBegin();
 
   while (!sit.IsAtEnd())
   {
     if (sit.Value()>0)
     {
       ItkFloatImgType::IndexType index = sit.GetIndex();
       itk::ContinuousIndex<float, 3> start;
       start[0] = index[0];
       start[1] = index[1];
       start[2] = index[2];
       itk::Point<float> worldPos;
       m_SeedImage->TransformContinuousIndexToPhysicalPoint(start, worldPos);
 
       if ( mitk::imv::IsInsideMask<float>(worldPos, m_Parameters->m_InterpolateRoiImages, m_MaskInterpolator) )
       {
         m_SeedPoints.push_back(worldPos);
         for (unsigned int s = 1; s < m_Parameters->m_SeedsPerVoxel; s++)
         {
           start[0] = index[0] + static_cast<float>(m_TrackingHandler->GetRandDouble(-0.5, 0.5));
           start[1] = index[1] + static_cast<float>(m_TrackingHandler->GetRandDouble(-0.5, 0.5));
           start[2] = index[2] + static_cast<float>(m_TrackingHandler->GetRandDouble(-0.5, 0.5));
 
           itk::Point<float> worldPos;
           m_SeedImage->TransformContinuousIndexToPhysicalPoint(start, worldPos);
           m_SeedPoints.push_back(worldPos);
         }
       }
     }
     ++sit;
   }
   if (m_SeedPoints.empty())
     mitkThrow() << "No valid seed point in seed image! Is your seed image registered with the image you are tracking on?";
 }
 
 void StreamlineTrackingFilter::GenerateData()
 {
   this->BeforeTracking();
   if (!m_Parameters->m_FixRandomSeed)
-    std::random_shuffle(m_SeedPoints.begin(), m_SeedPoints.end());
+  {
+    std::random_device rd;
+    std::mt19937 g(rd());
+    std::shuffle(m_SeedPoints.begin(), m_SeedPoints.end(), g);
+  }
 
   m_CurrentTracts = 0;
   int num_seeds = static_cast<int>(m_SeedPoints.size());
   itk::Index<3> zeroIndex; zeroIndex.Fill(0);
   m_Progress = 0;
   int i = 0;
   int print_interval = num_seeds/100;
   if (print_interval<100)
     m_Verbose=false;
 
   unsigned int trials_per_seed = 1;
   if(m_Parameters->m_Mode==mitk::TrackingDataHandler::MODE::PROBABILISTIC)
     trials_per_seed = m_Parameters->m_TrialsPerSeed;
 
 #pragma omp parallel
   while (i<num_seeds && !m_StopTracking)
   {
     int temp_i = 0;
 #pragma omp critical
     {
       temp_i = i;
       i++;
     }
 
     if (temp_i>=num_seeds || m_StopTracking)
       continue;
     else if (m_Verbose && i%print_interval==0)
 #pragma omp critical
     {
       m_Progress += static_cast<unsigned int>(print_interval);
       std::cout << "                                                                                                     \r";
       if (m_Parameters->m_MaxNumFibers>0)
         std::cout << "Tried: " << m_Progress << "/" << num_seeds << " | Accepted: " << m_CurrentTracts << "/" << m_Parameters->m_MaxNumFibers << '\r';
       else
         std::cout << "Tried: " << m_Progress << "/" << num_seeds << " | Accepted: " << m_CurrentTracts << '\r';
       cout.flush();
     }
 
     const itk::Point<float> worldPos = m_SeedPoints.at(static_cast<unsigned int>(temp_i));
 
     for (unsigned int trials=0; trials<trials_per_seed; ++trials)
     {
       FiberType fib;
       DirectionContainer direction_container;
       float tractLength = 0;
       unsigned long counter = 0;
 
       // get starting direction
       vnl_vector_fixed<float,3> dir; dir.fill(0.0);
       std::deque< vnl_vector_fixed<float,3> > olddirs;
       dir = GetNewDirection(worldPos, olddirs, zeroIndex) * 0.5f;
 
       bool exclude = false;
       if (m_ExclusionRegions.IsNotNull() && mitk::imv::IsInsideMask<float>(worldPos, m_Parameters->m_InterpolateRoiImages, m_ExclusionInterpolator))
         exclude = true;
 
       bool success = false;
       if (dir.magnitude()>0.0001f && !exclude)
       {
         // forward tracking
         tractLength = FollowStreamline(worldPos, dir, &fib, &direction_container, 0, false, exclude);
         fib.push_front(worldPos);
 
         // backward tracking
         if (!exclude)
           tractLength = FollowStreamline(worldPos, -dir, &fib, &direction_container, tractLength, true, exclude);
 
         counter = fib.size();
 
         if (tractLength>=m_Parameters->m_MinTractLengthMm && counter>=2 && !exclude)
         {
 #pragma omp critical
           if ( IsValidFiber(&fib) )
           {
             if (!m_StopTracking)
             {
               if (!m_Parameters->m_OutputProbMap)
                 m_Tractogram.push_back(fib);
               else
                 FiberToProbmap(&fib);
               m_CurrentTracts++;
               success = true;
             }
             if (m_Parameters->m_MaxNumFibers > 0 && m_CurrentTracts>=static_cast<unsigned int>(m_Parameters->m_MaxNumFibers))
             {
               if (!m_StopTracking)
               {
                 std::cout << "                                                                                                     \r";
                 MITK_INFO << "Reconstructed maximum number of tracts (" << m_CurrentTracts << "). Stopping tractography.";
               }
               m_StopTracking = true;
             }
           }
         }
       }
 
       if (success || m_Parameters->m_Mode!=MODE::PROBABILISTIC)
         break;  // we only try one seed point multiple times if we use a probabilistic tracker and have not found a valid streamline yet
 
     }// trials per seed
 
   }// seed points
 
   this->AfterTracking();
 }
 
 bool StreamlineTrackingFilter::IsValidFiber(FiberType* fib)
 {
   if (m_Parameters->m_EpConstraints==EndpointConstraints::NONE)
   {
     return true;
   }
   else if (m_Parameters->m_EpConstraints==EndpointConstraints::EPS_IN_TARGET)
   {
     if (m_TargetImageSet)
     {
       if ( mitk::imv::IsInsideMask<float>(fib->front(), m_Parameters->m_InterpolateRoiImages, m_TargetInterpolator)
            && mitk::imv::IsInsideMask<float>(fib->back(), m_Parameters->m_InterpolateRoiImages, m_TargetInterpolator) )
         return true;
       return false;
     }
     else
       mitkThrow() << "No target image set but endpoint constraint EPS_IN_TARGET chosen!";
   }
   else if (m_Parameters->m_EpConstraints==EndpointConstraints::EPS_IN_TARGET_LABELDIFF)
   {
     if (m_TargetImageSet)
     {
       float v1 = mitk::imv::GetImageValue<float>(fib->front(), false, m_TargetInterpolator);
       float v2 = mitk::imv::GetImageValue<float>(fib->back(), false, m_TargetInterpolator);
       if ( v1>0.0f && v2>0.0f && v1!=v2  )
         return true;
       return false;
     }
     else
       mitkThrow() << "No target image set but endpoint constraint EPS_IN_TARGET_LABELDIFF chosen!";
   }
   else if (m_Parameters->m_EpConstraints==EndpointConstraints::EPS_IN_SEED_AND_TARGET)
   {
     if (m_TargetImageSet && m_SeedImageSet)
     {
       if ( mitk::imv::IsInsideMask<float>(fib->front(), m_Parameters->m_InterpolateRoiImages, m_SeedInterpolator)
            && mitk::imv::IsInsideMask<float>(fib->back(), m_Parameters->m_InterpolateRoiImages, m_TargetInterpolator) )
         return true;
       if ( mitk::imv::IsInsideMask<float>(fib->back(), m_Parameters->m_InterpolateRoiImages, m_SeedInterpolator)
            && mitk::imv::IsInsideMask<float>(fib->front(), m_Parameters->m_InterpolateRoiImages, m_TargetInterpolator) )
         return true;
       return false;
     }
     else
       mitkThrow() << "No target or seed image set but endpoint constraint EPS_IN_SEED_AND_TARGET chosen!";
   }
   else if (m_Parameters->m_EpConstraints==EndpointConstraints::MIN_ONE_EP_IN_TARGET)
   {
     if (m_TargetImageSet)
     {
       if ( mitk::imv::IsInsideMask<float>(fib->front(), m_Parameters->m_InterpolateRoiImages, m_TargetInterpolator)
            || mitk::imv::IsInsideMask<float>(fib->back(), m_Parameters->m_InterpolateRoiImages, m_TargetInterpolator) )
         return true;
       return false;
     }
     else
       mitkThrow() << "No target image set but endpoint constraint MIN_ONE_EP_IN_TARGET chosen!";
   }
   else if (m_Parameters->m_EpConstraints==EndpointConstraints::ONE_EP_IN_TARGET)
   {
     if (m_TargetImageSet)
     {
       if ( mitk::imv::IsInsideMask<float>(fib->front(), m_Parameters->m_InterpolateRoiImages, m_TargetInterpolator)
            && !mitk::imv::IsInsideMask<float>(fib->back(), m_Parameters->m_InterpolateRoiImages, m_TargetInterpolator) )
         return true;
       if ( !mitk::imv::IsInsideMask<float>(fib->back(), m_Parameters->m_InterpolateRoiImages, m_TargetInterpolator)
            && mitk::imv::IsInsideMask<float>(fib->front(), m_Parameters->m_InterpolateRoiImages, m_TargetInterpolator) )
         return true;
       return false;
     }
     else
       mitkThrow() << "No target image set but endpoint constraint ONE_EP_IN_TARGET chosen!";
   }
   else if (m_Parameters->m_EpConstraints==EndpointConstraints::NO_EP_IN_TARGET)
   {
     if (m_TargetImageSet)
     {
       if ( mitk::imv::IsInsideMask<float>(fib->front(), m_Parameters->m_InterpolateRoiImages, m_TargetInterpolator)
            || mitk::imv::IsInsideMask<float>(fib->back(), m_Parameters->m_InterpolateRoiImages, m_TargetInterpolator) )
         return false;
       return true;
     }
     else
       mitkThrow() << "No target image set but endpoint constraint NO_EP_IN_TARGET chosen!";
   }
 
   return true;
 }
 
 void StreamlineTrackingFilter::FiberToProbmap(FiberType* fib)
 {
   ItkDoubleImgType::IndexType last_idx; last_idx.Fill(0);
   for (auto p : *fib)
   {
     ItkDoubleImgType::IndexType idx;
     m_OutputProbabilityMap->TransformPhysicalPointToIndex(p, idx);
 
     if (idx != last_idx)
     {
       if (m_OutputProbabilityMap->GetLargestPossibleRegion().IsInside(idx))
         m_OutputProbabilityMap->SetPixel(idx, m_OutputProbabilityMap->GetPixel(idx)+1);
       last_idx = idx;
     }
   }
 }
 
 void StreamlineTrackingFilter::BuildFibers(bool check)
 {
   if (m_BuildFibersReady<omp_get_num_threads() && check)
     return;
 
   m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
   vtkSmartPointer<vtkCellArray> vNewLines = vtkSmartPointer<vtkCellArray>::New();
   vtkSmartPointer<vtkPoints> vNewPoints = vtkSmartPointer<vtkPoints>::New();
 
   for (unsigned int i=0; i<m_Tractogram.size(); i++)
   {
     vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
     FiberType fib = m_Tractogram.at(i);
     for (FiberType::iterator it = fib.begin(); it!=fib.end(); ++it)
     {
       vtkIdType id = vNewPoints->InsertNextPoint((*it).GetDataPointer());
       container->GetPointIds()->InsertNextId(id);
     }
     vNewLines->InsertNextCell(container);
   }
 
   if (check)
     for (int i=0; i<m_BuildFibersReady; i++)
       m_Tractogram.pop_back();
   m_BuildFibersReady = 0;
 
   m_FiberPolyData->SetPoints(vNewPoints);
   m_FiberPolyData->SetLines(vNewLines);
   m_BuildFibersFinished = true;
 }
 
 void StreamlineTrackingFilter::AfterTracking()
 {
   if (m_Verbose)
     std::cout << "                                                                                                     \r";
   if (!m_Parameters->m_OutputProbMap)
   {
     MITK_INFO << "Reconstructed " << m_Tractogram.size() << " fibers.";
     MITK_INFO << "Generating polydata ";
     BuildFibers(false);
   }
   else
   {
     itk::RescaleIntensityImageFilter< ItkDoubleImgType, ItkDoubleImgType >::Pointer filter = itk::RescaleIntensityImageFilter< ItkDoubleImgType, ItkDoubleImgType >::New();
     filter->SetInput(m_OutputProbabilityMap);
     filter->SetOutputMaximum(1.0);
     filter->SetOutputMinimum(0.0);
     filter->Update();
     m_OutputProbabilityMap = filter->GetOutput();
   }
   MITK_INFO << "done";
 
   m_EndTime = std::chrono::system_clock::now();
   std::chrono::hours   hh = std::chrono::duration_cast<std::chrono::hours>(m_EndTime - m_StartTime);
   std::chrono::minutes mm = std::chrono::duration_cast<std::chrono::minutes>(m_EndTime - m_StartTime);
   std::chrono::seconds ss = std::chrono::duration_cast<std::chrono::seconds>(m_EndTime - m_StartTime);
   mm %= 60;
   ss %= 60;
   MITK_INFO << "Tracking took " << hh.count() << "h, " << mm.count() << "m and " << ss.count() << "s";
 
   m_SeedPoints.clear();
 }
 
 void StreamlineTrackingFilter::SetDicomProperties(mitk::FiberBundle::Pointer fib)
 {
   std::string model_code_value = "-";
   std::string model_code_meaning = "-";
   std::string algo_code_value = "-";
   std::string algo_code_meaning = "-";
 
   if ( m_Parameters->m_Mode==MODE::DETERMINISTIC && dynamic_cast<mitk::TrackingHandlerTensor*>(m_TrackingHandler))
   {
     algo_code_value = "sup181_ee01";
     algo_code_meaning = "Deterministic";
 
     if (m_Parameters->m_F > 0.99 && m_Parameters->m_G < 0.01)
     {
       if (m_Parameters->m_InterpolateTractographyData)
       {
         algo_code_value = "sup181_ee08";
         algo_code_meaning = "Euler";
       }
       else
       {
         algo_code_value = "sup181_ee04";
         algo_code_meaning = "FACT";
       }
     }
     else if (m_Parameters->m_G > 0.99 && m_Parameters->m_F < 0.01)
     {
       algo_code_value = "sup181_ee06";
       algo_code_meaning = "TEND";
     }
   }
   else if (m_Parameters->m_Mode==MODE::DETERMINISTIC)
   {
     algo_code_value = "sup181_ee01";
     algo_code_meaning = "Deterministic";
   }
   else if (m_Parameters->m_Mode==MODE::PROBABILISTIC)
   {
     algo_code_value = "sup181_ee02";
     algo_code_meaning = "Probabilistic";
   }
 
   if (dynamic_cast<mitk::TrackingHandlerTensor*>(m_TrackingHandler) || (dynamic_cast<mitk::TrackingHandlerOdf*>(m_TrackingHandler) && dynamic_cast<mitk::TrackingHandlerOdf*>(m_TrackingHandler)->GetIsOdfFromTensor() ) )
   {
     if ( dynamic_cast<mitk::TrackingHandlerTensor*>(m_TrackingHandler) && dynamic_cast<mitk::TrackingHandlerTensor*>(m_TrackingHandler)->GetNumTensorImages()>1 )
     {
       model_code_value = "sup181_bb02";
       model_code_meaning = "Multi Tensor";
     }
     else
     {
       model_code_value = "sup181_bb01";
       model_code_meaning = "Single Tensor";
     }
   }
   else if (dynamic_cast<mitk::TrackingHandlerRandomForest<6, 28>*>(m_TrackingHandler) || dynamic_cast<mitk::TrackingHandlerRandomForest<6, 100>*>(m_TrackingHandler))
   {
     model_code_value = "sup181_bb03";
     model_code_meaning = "Model Free";
   }
   else if (dynamic_cast<mitk::TrackingHandlerOdf*>(m_TrackingHandler))
   {
     model_code_value = "-";
     model_code_meaning = "ODF";
   }
   else if (dynamic_cast<mitk::TrackingHandlerPeaks*>(m_TrackingHandler))
   {
     model_code_value = "-";
     model_code_meaning = "Peaks";
   }
 
   fib->SetProperty("DICOM.anatomy.value", mitk::StringProperty::New("T-A0095"));
   fib->SetProperty("DICOM.anatomy.meaning", mitk::StringProperty::New("White matter of brain and spinal cord"));
 
   fib->SetProperty("DICOM.algo_family_code.value", mitk::StringProperty::New(algo_code_value));
   fib->SetProperty("DICOM.algo_family_code.meaning", mitk::StringProperty::New(algo_code_meaning));
 
   fib->SetProperty("DICOM.model_code.value", mitk::StringProperty::New(model_code_value));
   fib->SetProperty("DICOM.model_code.meaning", mitk::StringProperty::New(model_code_meaning));
 }
 
 }