diff --git a/Modules/Core/src/Rendering/mitkPointSetVtkMapper2D.cpp b/Modules/Core/src/Rendering/mitkPointSetVtkMapper2D.cpp
index 833c828976..4db5bddb8b 100644
--- a/Modules/Core/src/Rendering/mitkPointSetVtkMapper2D.cpp
+++ b/Modules/Core/src/Rendering/mitkPointSetVtkMapper2D.cpp
@@ -1,713 +1,715 @@
 /*===================================================================
 
 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 "mitkPointSetVtkMapper2D.h"
 
 //mitk includes
 #include "mitkVtkPropRenderer.h"
 #include <mitkDataNode.h>
 #include <mitkProperties.h>
 #include <mitkPointSet.h>
 #include <mitkPlaneGeometry.h>
 
 //vtk includes
 #include <vtkActor.h>
 #include <vtkPropAssembly.h>
 #include <vtkPolyDataMapper.h>
 #include <vtkTransform.h>
 #include <vtkGlyph3D.h>
 #include <vtkTransformFilter.h>
 #include <vtkLine.h>
 #include <vtkGlyphSource2D.h>
 #include <vtkFloatArray.h>
 #include <vtkPointData.h>
 #include <vtkTextActor.h>
 #include <vtkTextProperty.h>
 #include <vtkCellArray.h>
 
 #include <stdlib.h>
 
 // constructor LocalStorage
 mitk::PointSetVtkMapper2D::LocalStorage::LocalStorage()
 {
   // points
   m_UnselectedPoints = vtkSmartPointer<vtkPoints>::New();
   m_SelectedPoints = vtkSmartPointer<vtkPoints>::New();
   m_ContourPoints = vtkSmartPointer<vtkPoints>::New();
 
   // scales
   m_UnselectedScales = vtkSmartPointer<vtkFloatArray>::New();
   m_SelectedScales = vtkSmartPointer<vtkFloatArray>::New();
 
   // distances
   m_DistancesBetweenPoints = vtkSmartPointer<vtkFloatArray>::New();
 
   // lines
   m_ContourLines = vtkSmartPointer<vtkCellArray>::New();
 
   // glyph source (provides the different shapes)
   m_UnselectedGlyphSource2D = vtkSmartPointer<vtkGlyphSource2D>::New();
   m_SelectedGlyphSource2D = vtkSmartPointer<vtkGlyphSource2D>::New();
 
   // glyphs
   m_UnselectedGlyph3D = vtkSmartPointer<vtkGlyph3D>::New();
   m_SelectedGlyph3D = vtkSmartPointer<vtkGlyph3D>::New();
 
   // polydata
   m_VtkUnselectedPointListPolyData = vtkSmartPointer<vtkPolyData>::New();
   m_VtkSelectedPointListPolyData = vtkSmartPointer <vtkPolyData>::New();
   m_VtkContourPolyData = vtkSmartPointer<vtkPolyData>::New();
 
   // actors
   m_UnselectedActor = vtkSmartPointer <vtkActor>::New();
   m_SelectedActor = vtkSmartPointer <vtkActor>::New();
   m_ContourActor = vtkSmartPointer <vtkActor>::New();
 
   // mappers
   m_VtkUnselectedPolyDataMapper = vtkSmartPointer<vtkPolyDataMapper>::New();
   m_VtkSelectedPolyDataMapper = vtkSmartPointer<vtkPolyDataMapper>::New();
   m_VtkContourPolyDataMapper = vtkSmartPointer<vtkPolyDataMapper>::New();
 
   // propassembly
   m_PropAssembly = vtkSmartPointer <vtkPropAssembly>::New();
 }
 // destructor LocalStorage
 mitk::PointSetVtkMapper2D::LocalStorage::~LocalStorage()
 {
 }
 
 // input for this mapper ( = point set)
 const mitk::PointSet* mitk::PointSetVtkMapper2D::GetInput() const
 {
   return static_cast<const mitk::PointSet * > ( GetDataNode()->GetData() );
 }
 
 // constructor PointSetVtkMapper2D
 mitk::PointSetVtkMapper2D::PointSetVtkMapper2D()
 : m_ShowContour(false),
 m_CloseContour(false),
 m_ShowPoints(true),
 m_ShowDistances(false),
 m_DistancesDecimalDigits(1),
 m_ShowAngles(false),
 m_ShowDistantLines(false),
 m_LineWidth(1),
 m_PointLineWidth(1),
 m_Point2DSize(6),
 m_IDShapeProperty(mitk::PointSetShapeProperty::CROSS),
 m_FillShape(false),
 m_DistanceToPlane(4.0f)
 {
 }
 
 // destructor
 mitk::PointSetVtkMapper2D::~PointSetVtkMapper2D()
 {
 }
 
 // reset mapper so that nothing is displayed e.g. toggle visiblity of the propassembly
 void mitk::PointSetVtkMapper2D::ResetMapper( BaseRenderer* renderer )
 {
   LocalStorage *ls = m_LSH.GetLocalStorage(renderer);
   ls->m_PropAssembly->VisibilityOff();
 }
 
 // returns propassembly
 vtkProp* mitk::PointSetVtkMapper2D::GetVtkProp(mitk::BaseRenderer * renderer)
 {
   LocalStorage *ls = m_LSH.GetLocalStorage(renderer);
   return ls->m_PropAssembly;
 }
 
 static bool makePerpendicularVector2D(const mitk::Vector2D& in, mitk::Vector2D& out)
 {
   // The dot product of orthogonal vectors is zero.
   // In two dimensions the slopes of perpendicular lines are negative reciprocals.
   if((fabs(in[0])>0) && ( (fabs(in[0])>fabs(in[1])) || (in[1] == 0) ) )
   {
     // negative reciprocal
     out[0]=-in[1]/in[0];
     out[1]=1;
     out.Normalize();
     return true;
   }
   else
     if(fabs(in[1])>0)
     {
       out[0]=1;
       // negative reciprocal
       out[1]=-in[0]/in[1];
       out.Normalize();
       return true;
     }
     else
       return false;
 }
 
 void mitk::PointSetVtkMapper2D::CreateVTKRenderObjects(mitk::BaseRenderer* renderer)
 {
   LocalStorage *ls = m_LSH.GetLocalStorage(renderer);
 
   unsigned i = 0;
 
   // The vtk text actors need to be removed manually from the propassembly
   // since the same vtk text actors are not overwriten within this function,
   // but new actors are added to the propassembly each time this function is executed.
   // Thus, the actors from the last call must be removed in the beginning.
   for(i=0; i< ls->m_VtkTextLabelActors.size(); i++)
   {
     if(ls->m_PropAssembly->GetParts()->IsItemPresent(ls->m_VtkTextLabelActors.at(i)))
       ls->m_PropAssembly->RemovePart(ls->m_VtkTextLabelActors.at(i));
   }
 
   for(i=0; i< ls->m_VtkTextDistanceActors.size(); i++)
   {
     if(ls->m_PropAssembly->GetParts()->IsItemPresent(ls->m_VtkTextDistanceActors.at(i)))
       ls->m_PropAssembly->RemovePart(ls->m_VtkTextDistanceActors.at(i));
   }
 
   for(i=0; i< ls->m_VtkTextAngleActors.size(); i++)
   {
     if(ls->m_PropAssembly->GetParts()->IsItemPresent(ls->m_VtkTextAngleActors.at(i)))
       ls->m_PropAssembly->RemovePart(ls->m_VtkTextAngleActors.at(i));
   }
 
   // initialize polydata here, otherwise we have update problems when
   // executing this function again
   ls->m_VtkUnselectedPointListPolyData = vtkSmartPointer<vtkPolyData>::New();
   ls->m_VtkSelectedPointListPolyData = vtkSmartPointer <vtkPolyData>::New();
   ls->m_VtkContourPolyData = vtkSmartPointer<vtkPolyData>::New();
 
   // get input point set and update the PointSet
   mitk::PointSet::Pointer input  = const_cast<mitk::PointSet*>(this->GetInput());
 
   // only update the input data, if the property tells us to
   bool update = true;
   this->GetDataNode()->GetBoolProperty("updateDataOnRender", update);
   if (update == true)
     input->Update();
 
   int timestep = this->GetTimestep();
   mitk::PointSet::DataType::Pointer itkPointSet = input->GetPointSet( timestep );
 
   if ( itkPointSet.GetPointer() == NULL)
   {
     ls->m_PropAssembly->VisibilityOff();
     return;
   }
 
   //iterator for point set
   mitk::PointSet::PointsContainer::Iterator pointsIter = itkPointSet->GetPoints()->Begin();
 
   // PointDataContainer has additional information to each point, e.g. whether
   // it is selected or not
   mitk::PointSet::PointDataContainer::Iterator pointDataIter;
   pointDataIter = itkPointSet->GetPointData()->Begin();
 
   //check if the list for the PointDataContainer is the same size as the PointsContainer.
   //If not, then the points were inserted manually and can not be visualized according to the PointData (selected/unselected)
   bool pointDataBroken = (itkPointSet->GetPointData()->Size() != itkPointSet->GetPoints()->Size());
 
   if( itkPointSet->GetPointData()->size() == 0 || pointDataBroken)
   {
     ls->m_PropAssembly->VisibilityOff();
     return;
   }
 
   ls->m_PropAssembly->VisibilityOn();
 
   // empty point sets, cellarrays, scalars
   ls->m_UnselectedPoints->Reset();
   ls->m_SelectedPoints->Reset();
 
   ls->m_ContourPoints->Reset();
   ls->m_ContourLines->Reset();
 
   ls->m_UnselectedScales->Reset();
   ls->m_SelectedScales->Reset();
 
   ls->m_DistancesBetweenPoints->Reset();
 
   ls->m_VtkTextLabelActors.clear();
   ls->m_VtkTextDistanceActors.clear();
   ls->m_VtkTextAngleActors.clear();
 
   ls->m_UnselectedScales->SetNumberOfComponents(3);
   ls->m_SelectedScales->SetNumberOfComponents(3);
 
   int NumberContourPoints = 0;
   bool pointsOnSameSideOfPlane = false;
 
   const int text2dDistance = 10;
 
   // initialize points with a random start value
 
   // current point in point set
   itk::Point<ScalarType> point = pointsIter->Value();
 
   mitk::Point3D p = point;              // currently visited point
   mitk::Point3D lastP = point;          // last visited point (predecessor in point set of "point")
   mitk::Vector3D vec;                   // p - lastP
   mitk::Vector3D lastVec;               // lastP - point before lastP
   vec.Fill(0.0);
   lastVec.Fill(0.0);
 
   mitk::Point3D projected_p = point;     // p projected on viewplane
 
   mitk::Point2D pt2d;
   pt2d[0] = point[0];                    // projected_p in display coordinates
   pt2d[1] = point[1];
   mitk::Point2D lastPt2d = pt2d;         // last projected_p in display coordinates (predecessor in point set of "pt2d")
   mitk::Point2D preLastPt2d = pt2d ;     // projected_p in display coordinates before lastPt2
 
   mitk::DisplayGeometry::Pointer displayGeometry = renderer->GetDisplayGeometry();
   const mitk::PlaneGeometry* geo2D = renderer->GetCurrentWorldPlaneGeometry();
 
   vtkLinearTransform* dataNodeTransform = input->GetGeometry()->GetVtkTransform();
 
   int count = 0;
 
   for (pointsIter=itkPointSet->GetPoints()->Begin();
     pointsIter!=itkPointSet->GetPoints()->End();
     pointsIter++)
   {
     lastP = p; // valid for number of points count > 0
     preLastPt2d = lastPt2d; // valid only for count > 1
     lastPt2d = pt2d;  // valid for number of points count > 0
 
     lastVec = vec;    // valid only for counter > 1
 
     // get current point in point set
     point = pointsIter->Value();
 
     // transform point
     {
       float vtkp[3];
       itk2vtk(point, vtkp);
       dataNodeTransform->TransformPoint(vtkp, vtkp);
       vtk2itk(vtkp,point);
     }
 
     p[0] = point[0];
     p[1] = point[1];
     p[2] = point[2];
 
     displayGeometry->Project(p, projected_p);
     displayGeometry->Map(projected_p, pt2d);
     displayGeometry->WorldToDisplay(pt2d, pt2d);
 
     vec = p-lastP;    // valid only for counter > 0
 
     // compute distance to current plane
-    float diff = geo2D->Distance(point);
-    diff = diff * diff;
+    float diff = fabs(geo2D->Distance(point));
+//    diff = diff * diff;
 
     //draw markers on slices a certain distance away from the points
     //location according to the tolerance threshold (m_DistanceToPlane)
     if(diff < m_DistanceToPlane)
     {
       // is point selected or not?
       if (pointDataIter->Value().selected)
       {
         ls->m_SelectedPoints->InsertNextPoint(point[0],point[1],point[2]);
         // point is scaled according to its distance to the plane
-        ls->m_SelectedScales->InsertNextTuple3(m_Point2DSize - (2*diff),0,0);
+        ls->m_SelectedScales->InsertNextTuple3((double)m_Point2DSize*(1-diff/m_DistanceToPlane),0,0);
       }
       else
       {
         ls->m_UnselectedPoints->InsertNextPoint(point[0],point[1],point[2]);
         // point is scaled according to its distance to the plane
-        ls->m_UnselectedScales->InsertNextTuple3(m_Point2DSize - (2*diff),0,0);
+        ls->m_UnselectedScales->InsertNextTuple3((double)m_Point2DSize*(1-diff/m_DistanceToPlane),0,0);
       }
 
       //---- LABEL -----//
       //paint label for each point if available
       if (dynamic_cast<mitk::StringProperty *>(this->GetDataNode()->GetProperty("label")) != NULL)
       {
         const char * pointLabel = dynamic_cast<mitk::StringProperty *>(
           this->GetDataNode()->GetProperty("label"))->GetValue();
         std::string l = pointLabel;
         if (input->GetSize()>1)
         {
           std::stringstream ss;
           ss << pointsIter->Index();
           l.append(ss.str());
         }
 
         ls->m_VtkTextActor = vtkSmartPointer<vtkTextActor>::New();
 
         ls->m_VtkTextActor->SetPosition(pt2d[0] + text2dDistance, pt2d[1] + text2dDistance);
         ls->m_VtkTextActor->SetInput(l.c_str());
         ls->m_VtkTextActor->GetTextProperty()->SetOpacity( 100 );
 
         float unselectedColor[4] = {1.0, 1.0, 0.0, 1.0};
 
         //check if there is a color property
         GetDataNode()->GetColor(unselectedColor);
 
         ls->m_VtkTextActor->GetTextProperty()->SetColor(unselectedColor[0], unselectedColor[1], unselectedColor[2]);
 
         ls->m_VtkTextLabelActors.push_back(ls->m_VtkTextActor);
       }
     }
 
     // draw contour, distance text and angle text in render window
 
     // lines between points, which intersect the current plane, are drawn
     if( m_ShowContour && count > 0 )
     {
       ScalarType distance =    displayGeometry->GetWorldGeometry()->SignedDistance(point);
       ScalarType lastDistance =    displayGeometry->GetWorldGeometry()->SignedDistance(lastP);
 
       pointsOnSameSideOfPlane = (distance * lastDistance) > 0.5;
 
       // Points must be on different side of plane in order to draw a contour.
       // If "show distant lines" is enabled this condition is disregarded.
       if ( !pointsOnSameSideOfPlane || m_ShowDistantLines)
       {
         vtkSmartPointer<vtkLine> line = vtkSmartPointer<vtkLine>::New();
 
         ls->m_ContourPoints->InsertNextPoint(lastP[0],lastP[1],lastP[2]);
         line->GetPointIds()->SetId(0, NumberContourPoints);
         NumberContourPoints++;
 
         ls->m_ContourPoints->InsertNextPoint(point[0], point[1], point[2]);
         line->GetPointIds()->SetId(1, NumberContourPoints);
         NumberContourPoints++;
 
         ls->m_ContourLines->InsertNextCell(line);
 
         if(m_ShowDistances) // calculate and print distance between adjacent points
         {
           float distancePoints = point.EuclideanDistanceTo(lastP);
 
           std::stringstream buffer;
           buffer<<std::fixed <<std::setprecision(m_DistancesDecimalDigits)<<distancePoints<<" mm";
 
           // compute desired display position of text
           Vector2D vec2d = pt2d-lastPt2d;
           makePerpendicularVector2D(vec2d, vec2d); // text is rendered within text2dDistance perpendicular to current line
           Vector2D pos2d = (lastPt2d.GetVectorFromOrigin() + pt2d.GetVectorFromOrigin() ) * 0.5 + vec2d * text2dDistance;
 
           ls->m_VtkTextActor = vtkSmartPointer<vtkTextActor>::New();
 
           ls->m_VtkTextActor->SetPosition(pos2d[0],pos2d[1]);
           ls->m_VtkTextActor->SetInput(buffer.str().c_str());
           ls->m_VtkTextActor->GetTextProperty()->SetColor(0.0, 1.0, 0.0);
 
           ls->m_VtkTextDistanceActors.push_back(ls->m_VtkTextActor);
         }
 
         if(m_ShowAngles && count > 1) // calculate and print angle between connected lines
         {
           std::stringstream buffer;
           buffer << angle(vec.GetVnlVector(), -lastVec.GetVnlVector())*180/vnl_math::pi << "°";
 
           //compute desired display position of text
           Vector2D vec2d = pt2d-lastPt2d;             // first arm enclosing the angle
           vec2d.Normalize();
           Vector2D lastVec2d = lastPt2d-preLastPt2d;  // second arm enclosing the angle
           lastVec2d.Normalize();
           vec2d=vec2d-lastVec2d;                      // vector connecting both arms
           vec2d.Normalize();
 
           // middle between two vectors that enclose the angle
           Vector2D pos2d = lastPt2d.GetVectorFromOrigin() + vec2d * text2dDistance * text2dDistance;
 
           ls->m_VtkTextActor = vtkSmartPointer<vtkTextActor>::New();
 
           ls->m_VtkTextActor->SetPosition(pos2d[0],pos2d[1]);
           ls->m_VtkTextActor->SetInput(buffer.str().c_str());
           ls->m_VtkTextActor->GetTextProperty()->SetColor(0.0, 1.0, 0.0);
 
           ls->m_VtkTextAngleActors.push_back(ls->m_VtkTextActor);
         }
       }
     }
 
     if(pointDataIter != itkPointSet->GetPointData()->End())
     {
       pointDataIter++;
       count++;
     }
   }
 
   // add each single text actor to the assembly
   for(i=0; i< ls->m_VtkTextLabelActors.size(); i++)
   {
     ls->m_PropAssembly->AddPart(ls->m_VtkTextLabelActors.at(i));
   }
 
   for(i=0; i< ls->m_VtkTextDistanceActors.size(); i++)
   {
     ls->m_PropAssembly->AddPart(ls->m_VtkTextDistanceActors.at(i));
   }
 
   for(i=0; i< ls->m_VtkTextAngleActors.size(); i++)
   {
     ls->m_PropAssembly->AddPart(ls->m_VtkTextAngleActors.at(i));
   }
 
   //---- CONTOUR -----//
 
   //create lines between the points which intersect the plane
   if (m_ShowContour)
   {
     // draw line between first and last point which is rendered
     if(m_CloseContour && NumberContourPoints > 1){
       vtkSmartPointer<vtkLine> closingLine = vtkSmartPointer<vtkLine>::New();
       closingLine->GetPointIds()->SetId(0, 0); // index of first point
       closingLine->GetPointIds()->SetId(1, NumberContourPoints-1); // index of last point
       ls->m_ContourLines->InsertNextCell(closingLine);
     }
 
     ls->m_VtkContourPolyData->SetPoints(ls->m_ContourPoints);
     ls->m_VtkContourPolyData->SetLines(ls->m_ContourLines);
 
     ls->m_VtkContourPolyDataMapper->SetInputData(ls->m_VtkContourPolyData);
     ls->m_ContourActor->SetMapper(ls->m_VtkContourPolyDataMapper);
     ls->m_ContourActor->GetProperty()->SetLineWidth(m_LineWidth);
 
     ls->m_PropAssembly->AddPart(ls->m_ContourActor);
   }
 
   // the point set must be transformed in order to obtain the appropriate glyph orientation
   // according to the current view
   vtkSmartPointer<vtkTransform> transform = vtkSmartPointer<vtkTransform>::New();
   vtkSmartPointer<vtkMatrix4x4> a,b = vtkSmartPointer<vtkMatrix4x4>::New();
 
   a = geo2D->GetVtkTransform()->GetMatrix();
   b->DeepCopy( a );
 
   // delete transformation from matrix, only take orientation
   b->SetElement(3,3,1);
   b->SetElement(2,3,0);
   b->SetElement(1,3,0);
   b->SetElement(0,3,0);
   b->SetElement(3,2,0);
   b->SetElement(3,1,0);
   b->SetElement(3,0,0);
 
   transform->SetMatrix(  b );
 
   //---- UNSELECTED POINTS  -----//
 
   // apply properties to glyph
   ls->m_UnselectedGlyphSource2D->SetGlyphType(m_IDShapeProperty);
 
   if(m_FillShape)
     ls->m_UnselectedGlyphSource2D->FilledOn();
   else
     ls->m_UnselectedGlyphSource2D->FilledOff();
 
+  ls->m_UnselectedGlyphSource2D->SetScale((double)m_Point2DSize/100);
+
   // apply transform
   vtkSmartPointer<vtkTransformFilter> transformFilterU = vtkSmartPointer<vtkTransformFilter>::New();
   transformFilterU->SetInputConnection(ls->m_UnselectedGlyphSource2D->GetOutputPort());
   transformFilterU->SetTransform(transform);
 
   ls->m_VtkUnselectedPointListPolyData->SetPoints(ls->m_UnselectedPoints);
   ls->m_VtkUnselectedPointListPolyData->GetPointData()->SetVectors(ls->m_UnselectedScales);
 
   // apply transform of current plane to glyphs
   ls->m_UnselectedGlyph3D->SetSourceConnection(transformFilterU->GetOutputPort());
   ls->m_UnselectedGlyph3D->SetInputData(ls->m_VtkUnselectedPointListPolyData);
   ls->m_UnselectedGlyph3D->SetScaleModeToScaleByVector();
   ls->m_UnselectedGlyph3D->SetVectorModeToUseVector();
 
   ls->m_VtkUnselectedPolyDataMapper->SetInputConnection(ls->m_UnselectedGlyph3D->GetOutputPort());
   ls->m_UnselectedActor->SetMapper(ls->m_VtkUnselectedPolyDataMapper);
   ls->m_UnselectedActor->GetProperty()->SetLineWidth(m_PointLineWidth);
 
   ls->m_PropAssembly->AddPart(ls->m_UnselectedActor);
 
   //---- SELECTED POINTS  -----//
 
   ls->m_SelectedGlyphSource2D->SetGlyphTypeToDiamond();
   ls->m_SelectedGlyphSource2D->CrossOn();
   ls->m_SelectedGlyphSource2D->FilledOff();
 
   // apply transform
   vtkSmartPointer<vtkTransformFilter> transformFilterS = vtkSmartPointer<vtkTransformFilter>::New();
   transformFilterS->SetInputConnection(ls->m_SelectedGlyphSource2D->GetOutputPort());
   transformFilterS->SetTransform(transform);
 
   ls->m_VtkSelectedPointListPolyData->SetPoints(ls->m_SelectedPoints);
   ls->m_VtkSelectedPointListPolyData->GetPointData()->SetVectors(ls->m_SelectedScales);
 
   // apply transform of current plane to glyphs
   ls->m_SelectedGlyph3D->SetSourceConnection(transformFilterS->GetOutputPort());
   ls->m_SelectedGlyph3D->SetInputData(ls->m_VtkSelectedPointListPolyData);
   ls->m_SelectedGlyph3D->SetScaleModeToScaleByVector();
   ls->m_SelectedGlyph3D->SetVectorModeToUseVector();
 
   ls->m_VtkSelectedPolyDataMapper->SetInputConnection(ls->m_SelectedGlyph3D->GetOutputPort());
   ls->m_SelectedActor->SetMapper(ls->m_VtkSelectedPolyDataMapper);
   ls->m_SelectedActor->GetProperty()->SetLineWidth(m_PointLineWidth);
 
   ls->m_PropAssembly->AddPart(ls->m_SelectedActor);
 }
 
 void mitk::PointSetVtkMapper2D::GenerateDataForRenderer( mitk::BaseRenderer *renderer )
 {
   const mitk::DataNode* node = GetDataNode();
   if( node == NULL )
     return;
 
   LocalStorage *ls = m_LSH.GetLocalStorage(renderer);
 
   // check whether the input data has been changed
   bool needGenerateData = ls->IsGenerateDataRequired( renderer, this, GetDataNode() );
 
   // toggle visibility
   bool visible = true;
   node->GetVisibility(visible, renderer, "visible");
   if(!visible)
   {
     ls->m_UnselectedActor->VisibilityOff();
     ls->m_SelectedActor->VisibilityOff();
     ls->m_ContourActor->VisibilityOff();
     ls->m_PropAssembly->VisibilityOff();
     return;
   }else{
     ls->m_PropAssembly->VisibilityOn();
   }
 
   node->GetBoolProperty("show contour",       m_ShowContour, renderer);
   node->GetBoolProperty("close contour",      m_CloseContour, renderer);
   node->GetBoolProperty("show points",        m_ShowPoints, renderer);
   node->GetBoolProperty("show distances",     m_ShowDistances, renderer);
   node->GetIntProperty("distance decimal digits",     m_DistancesDecimalDigits, renderer);
   node->GetBoolProperty("show angles",        m_ShowAngles, renderer);
   node->GetBoolProperty("show distant lines", m_ShowDistantLines, renderer);
   node->GetIntProperty("line width",          m_LineWidth, renderer);
   node->GetIntProperty("point line width",    m_PointLineWidth, renderer);
   node->GetIntProperty("point 2D size",       m_Point2DSize, renderer);
   node->GetBoolProperty("Pointset.2D.fill shape", m_FillShape, renderer);
   node->GetFloatProperty("Pointset.2D.distance to plane", m_DistanceToPlane, renderer );
 
   mitk::PointSetShapeProperty::Pointer shape = dynamic_cast<mitk::PointSetShapeProperty*>(this->GetDataNode()->GetProperty( "Pointset.2D.shape", renderer ));
   if(shape.IsNotNull())
   {
     m_IDShapeProperty = shape->GetPointSetShape();
   }
 
   //check for color props and use it for rendering of selected/unselected points and contour
   //due to different params in VTK (double/float) we have to convert
 
   float unselectedColor[4];
   double selectedColor[4]={1.0f,0.0f,0.0f,1.0f};    //red
   double contourColor[4]={1.0f,0.0f,0.0f,1.0f};     //red
 
   float opacity = 1.0;
 
   GetDataNode()->GetOpacity(opacity, renderer);
 
   // apply color and opacity
   if(m_ShowPoints)
   {
     ls->m_UnselectedActor->VisibilityOn();
     ls->m_SelectedActor->VisibilityOn();
 
     //check if there is a color property
     GetDataNode()->GetColor(unselectedColor);
 
     //get selected color property
     if (dynamic_cast<mitk::ColorProperty*>(this->GetDataNode()->GetPropertyList(renderer)->GetProperty("selectedcolor")) != NULL)
     {
       mitk::Color tmpColor = dynamic_cast<mitk::ColorProperty *>(this->GetDataNode()->GetPropertyList(renderer)->GetProperty("selectedcolor"))->GetValue();
       selectedColor[0] = tmpColor[0];
       selectedColor[1] = tmpColor[1];
       selectedColor[2] = tmpColor[2];
       selectedColor[3] = 1.0f; // alpha value
     }
     else if (dynamic_cast<mitk::ColorProperty*>(this->GetDataNode()->GetPropertyList(NULL)->GetProperty("selectedcolor")) != NULL)
     {
       mitk::Color tmpColor = dynamic_cast<mitk::ColorProperty *>(this->GetDataNode()->GetPropertyList(NULL)->GetProperty("selectedcolor"))->GetValue();
       selectedColor[0] = tmpColor[0];
       selectedColor[1] = tmpColor[1];
       selectedColor[2] = tmpColor[2];
       selectedColor[3] = 1.0f; // alpha value
     }
 
     ls->m_SelectedActor->GetProperty()->SetColor(selectedColor);
     ls->m_SelectedActor->GetProperty()->SetOpacity(opacity);
 
     ls->m_UnselectedActor->GetProperty()->SetColor(unselectedColor[0],unselectedColor[1],unselectedColor[2]);
     ls->m_UnselectedActor->GetProperty()->SetOpacity(opacity);
   }
   else
   {
     ls->m_UnselectedActor->VisibilityOff();
     ls-> m_SelectedActor->VisibilityOff();
   }
 
   if (m_ShowContour)
   {
     ls->m_ContourActor->VisibilityOn();
 
     //get contour color property
     if (dynamic_cast<mitk::ColorProperty*>(this->GetDataNode()->GetPropertyList(renderer)->GetProperty("contourcolor")) != NULL)
     {
       mitk::Color tmpColor = dynamic_cast<mitk::ColorProperty *>(this->GetDataNode()->GetPropertyList(renderer)->GetProperty("contourcolor"))->GetValue();
       contourColor[0] = tmpColor[0];
       contourColor[1] = tmpColor[1];
       contourColor[2] = tmpColor[2];
       contourColor[3] = 1.0f;
     }
     else if (dynamic_cast<mitk::ColorProperty*>(this->GetDataNode()->GetPropertyList(NULL)->GetProperty("contourcolor")) != NULL)
     {
       mitk::Color tmpColor = dynamic_cast<mitk::ColorProperty *>(this->GetDataNode()->GetPropertyList(NULL)->GetProperty("contourcolor"))->GetValue();
       contourColor[0] = tmpColor[0];
       contourColor[1] = tmpColor[1];
       contourColor[2] = tmpColor[2];
       contourColor[3] = 1.0f;
     }
 
     ls->m_ContourActor->GetProperty()->SetColor(contourColor);
     ls->m_ContourActor->GetProperty()->SetOpacity(opacity);
   }
   else
   {
     ls->m_ContourActor->VisibilityOff();
   }
 
   if(needGenerateData)
   {
     // create new vtk render objects (e.g. a circle for a point)
     this->CreateVTKRenderObjects(renderer);
   }
 }
 
 void mitk::PointSetVtkMapper2D::SetDefaultProperties(mitk::DataNode* node, mitk::BaseRenderer* renderer, bool overwrite)
 {
   node->AddProperty( "line width", mitk::IntProperty::New(2), renderer, overwrite );
   node->AddProperty( "point line width", mitk::IntProperty::New(1), renderer, overwrite );
   node->AddProperty( "point 2D size", mitk::IntProperty::New(6), renderer, overwrite );
   node->AddProperty( "show contour", mitk::BoolProperty::New(false), renderer, overwrite );
   node->AddProperty( "close contour", mitk::BoolProperty::New(false), renderer, overwrite );
   node->AddProperty( "show points", mitk::BoolProperty::New(true), renderer, overwrite );
   node->AddProperty( "show distances", mitk::BoolProperty::New(false), renderer, overwrite );
   node->AddProperty( "distance decimal digits", mitk::IntProperty::New(2), renderer, overwrite );
   node->AddProperty( "show angles", mitk::BoolProperty::New(false), renderer, overwrite );
   node->AddProperty( "show distant lines", mitk::BoolProperty::New(false), renderer, overwrite );
   node->AddProperty( "layer", mitk::IntProperty::New(1), renderer, overwrite );
   node->AddProperty( "Pointset.2D.fill shape", mitk::BoolProperty::New(false), renderer, overwrite); // fill or do not fill the glyph shape
   mitk::PointSetShapeProperty::Pointer pointsetShapeProperty = mitk::PointSetShapeProperty::New();
   node->AddProperty( "Pointset.2D.shape", pointsetShapeProperty, renderer, overwrite);
   node->AddProperty( "Pointset.2D.distance to plane", mitk::FloatProperty::New(4.0f), renderer, overwrite ); //show the point at a certain distance above/below the 2D imaging plane.
 
   Superclass::SetDefaultProperties(node, renderer, overwrite);
 }
diff --git a/Modules/DiffusionImaging/DiffusionIO/mitkFiberBundleVtkReader.cpp b/Modules/DiffusionImaging/DiffusionIO/mitkFiberBundleVtkReader.cpp
index e492d58937..afee3da252 100644
--- a/Modules/DiffusionImaging/DiffusionIO/mitkFiberBundleVtkReader.cpp
+++ b/Modules/DiffusionImaging/DiffusionIO/mitkFiberBundleVtkReader.cpp
@@ -1,277 +1,277 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 
 #include "mitkFiberBundleVtkReader.h"
 #include <itkMetaDataObject.h>
 #include <vtkPolyData.h>
 #include <vtkDataReader.h>
 #include <vtkPolyDataReader.h>
 #include <vtkMatrix4x4.h>
 #include <vtkPolyLine.h>
 #include <vtkCellArray.h>
 #include <vtkDataArray.h>
 #include <vtkFloatArray.h>
 #include <vtkCellData.h>
 #include <vtkPointData.h>
 #include <itksys/SystemTools.hxx>
 #include <tinyxml.h>
 #include <vtkCleanPolyData.h>
 #include <mitkTrackvis.h>
 #include <mitkCustomMimeType.h>
 #include <vtkXMLPolyDataReader.h>
 #include <vtkXMLDataReader.h>
 #include "mitkDiffusionIOMimeTypes.h"
 
 
 mitk::FiberBundleVtkReader::FiberBundleVtkReader()
     : mitk::AbstractFileReader( mitk::DiffusionIOMimeTypes::FIBERBUNDLE_VTK_MIMETYPE_NAME(), "VTK Fiber Bundle Reader" )
 {
     m_ServiceReg = this->RegisterService();
 }
 
 mitk::FiberBundleVtkReader::FiberBundleVtkReader(const FiberBundleVtkReader &other)
     :mitk::AbstractFileReader(other)
 {
 }
 
 mitk::FiberBundleVtkReader * mitk::FiberBundleVtkReader::Clone() const
 {
     return new FiberBundleVtkReader(*this);
 }
 
 
 std::vector<itk::SmartPointer<mitk::BaseData> > mitk::FiberBundleVtkReader::Read()
 {
 
     std::vector<itk::SmartPointer<mitk::BaseData> > result;
 
     const std::string& locale = "C";
     const std::string& currLocale = setlocale( LC_ALL, NULL );
     setlocale(LC_ALL, locale.c_str());
 
     std::string filename = this->GetInputLocation();
 
     std::string ext = itksys::SystemTools::GetFilenameLastExtension(filename);
     ext = itksys::SystemTools::LowerCase(ext);
 
     try
     {
         MITK_INFO << "Trying to load fiber file as VTK format.";
         vtkSmartPointer<vtkPolyDataReader> reader = vtkSmartPointer<vtkPolyDataReader>::New();
         reader->SetFileName( this->GetInputLocation().c_str() );
 
         if (reader->IsFilePolyData())
         {
             reader->Update();
 
             if ( reader->GetOutput() != NULL )
             {
                 vtkSmartPointer<vtkPolyData> fiberPolyData = reader->GetOutput();
                 FiberBundle::Pointer fiberBundle = FiberBundle::New(fiberPolyData);
 
                 vtkSmartPointer<vtkFloatArray> weights = vtkFloatArray::SafeDownCast(fiberPolyData->GetCellData()->GetArray("FIBER_WEIGHTS"));
                 if (weights!=NULL)
                 {
-//                                float weight=0;
-//                                for (int i=0; i<weights->GetSize(); i++)
-//                                    if (!mitk::Equal(weights->GetValue(i),weight,0.00001))
-//                                    {
-//                                        MITK_INFO << "Weight: " << weights->GetValue(i);
-//                                        weight = weights->GetValue(i);
-//                                    }
+                    float weight=0;
+                    for (int i=0; i<weights->GetSize(); i++)
+                        if (!mitk::Equal(weights->GetValue(i),weight,0.00001))
+                        {
+                            MITK_INFO << "Weight: " << weights->GetValue(i);
+                            weight = weights->GetValue(i);
+                        }
                     fiberBundle->SetFiberWeights(weights);
                 }
 
                 vtkSmartPointer<vtkUnsignedCharArray> fiberColors = vtkUnsignedCharArray::SafeDownCast(fiberPolyData->GetPointData()->GetArray("FIBER_COLORS"));
                 if (fiberColors!=NULL)
                     fiberBundle->SetFiberColors(fiberColors);
 
                 result.push_back(fiberBundle.GetPointer());
                 return result;
             }
         }
         else
             MITK_INFO << "File is not VTK format.";
     }
     catch(...)
     {
         throw;
     }
 
     try
     {
         MITK_INFO << "Trying to load fiber file as VTP format.";
         vtkSmartPointer<vtkXMLPolyDataReader> reader = vtkSmartPointer<vtkXMLPolyDataReader>::New();
         reader->SetFileName( this->GetInputLocation().c_str() );
 
         if ( reader->CanReadFile(this->GetInputLocation().c_str()) )
         {
             reader->Update();
 
             if ( reader->GetOutput() != NULL )
             {
                 vtkSmartPointer<vtkPolyData> fiberPolyData = reader->GetOutput();
                 FiberBundle::Pointer fiberBundle = FiberBundle::New(fiberPolyData);
 
                 vtkSmartPointer<vtkFloatArray> weights = vtkFloatArray::SafeDownCast(fiberPolyData->GetCellData()->GetArray("FIBER_WEIGHTS"));
 
                 if (weights!=NULL)
                 {
 //                                float weight=0;
 //                                for (int i=0; i<weights->GetSize(); i++)
 //                                    if (!mitk::Equal(weights->GetValue(i),weight,0.00001))
 //                                    {
 //                                        MITK_INFO << "Weight: " << weights->GetValue(i);
 //                                        weight = weights->GetValue(i);
 //                                    }
                     fiberBundle->SetFiberWeights(weights);
                 }
 
                 vtkSmartPointer<vtkUnsignedCharArray> fiberColors = vtkUnsignedCharArray::SafeDownCast(fiberPolyData->GetPointData()->GetArray("FIBER_COLORS"));
                 if (fiberColors!=NULL)
                     fiberBundle->SetFiberColors(fiberColors);
 
                 result.push_back(fiberBundle.GetPointer());
                 return result;
             }
         }
         else
             MITK_INFO << "File is not VTP format.";
     }
     catch(...)
     {
         throw;
     }
 
     try
     {
         MITK_INFO << "Trying to load fiber file as deprecated XML format.";
         vtkSmartPointer<vtkPolyData> fiberPolyData = vtkSmartPointer<vtkPolyData>::New();
         vtkSmartPointer<vtkCellArray> cellArray = vtkSmartPointer<vtkCellArray>::New();
         vtkSmartPointer<vtkPoints>    points = vtkSmartPointer<vtkPoints>::New();
         TiXmlDocument doc( this->GetInputLocation().c_str() );
         if(doc.LoadFile())
         {
             TiXmlHandle hDoc(&doc);
             TiXmlElement* pElem;
             TiXmlHandle hRoot(0);
             pElem = hDoc.FirstChildElement().Element();
             // save this for later
             hRoot = TiXmlHandle(pElem);
             pElem = hRoot.FirstChildElement("geometry").Element();
             // read geometry
             mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New();
             // read origin
             mitk::Point3D origin;
             double temp = 0;
             pElem->Attribute("origin_x", &temp);
             origin[0] = temp;
             pElem->Attribute("origin_y", &temp);
             origin[1] = temp;
             pElem->Attribute("origin_z", &temp);
             origin[2] = temp;
             geometry->SetOrigin(origin);
             // read spacing
             ScalarType spacing[3];
             pElem->Attribute("spacing_x", &temp);
             spacing[0] = temp;
             pElem->Attribute("spacing_y", &temp);
             spacing[1] = temp;
             pElem->Attribute("spacing_z", &temp);
             spacing[2] = temp;
             geometry->SetSpacing(spacing);
             // read transform
             vtkMatrix4x4* m = vtkMatrix4x4::New();
             pElem->Attribute("xx", &temp);
             m->SetElement(0,0,temp);
             pElem->Attribute("xy", &temp);
             m->SetElement(1,0,temp);
             pElem->Attribute("xz", &temp);
             m->SetElement(2,0,temp);
             pElem->Attribute("yx", &temp);
             m->SetElement(0,1,temp);
             pElem->Attribute("yy", &temp);
             m->SetElement(1,1,temp);
             pElem->Attribute("yz", &temp);
             m->SetElement(2,1,temp);
             pElem->Attribute("zx", &temp);
             m->SetElement(0,2,temp);
             pElem->Attribute("zy", &temp);
             m->SetElement(1,2,temp);
             pElem->Attribute("zz", &temp);
             m->SetElement(2,2,temp);
             m->SetElement(0,3,origin[0]);
             m->SetElement(1,3,origin[1]);
             m->SetElement(2,3,origin[2]);
             m->SetElement(3,3,1);
             geometry->SetIndexToWorldTransformByVtkMatrix(m);
             // read bounds
             float bounds[] = {0, 0, 0, 0, 0, 0};
             pElem->Attribute("size_x", &temp);
             bounds[1] = temp;
             pElem->Attribute("size_y", &temp);
             bounds[3] = temp;
             pElem->Attribute("size_z", &temp);
             bounds[5] = temp;
             geometry->SetFloatBounds(bounds);
             geometry->SetImageGeometry(true);
             pElem = hRoot.FirstChildElement("fiber_bundle").FirstChild().Element();
             for( ; pElem ; pElem=pElem->NextSiblingElement())
             {
                 TiXmlElement* pElem2 = pElem->FirstChildElement();
                 vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
                 for( ; pElem2; pElem2=pElem2->NextSiblingElement())
                 {
                     Point3D point;
                     pElem2->Attribute("pos_x", &temp);
                     point[0] = temp;
                     pElem2->Attribute("pos_y", &temp);
                     point[1] = temp;
                     pElem2->Attribute("pos_z", &temp);
                     point[2] = temp;
                     geometry->IndexToWorld(point, point);
                     vtkIdType id = points->InsertNextPoint(point.GetDataPointer());
                     container->GetPointIds()->InsertNextId(id);
                 }
                 cellArray->InsertNextCell(container);
             }
             fiberPolyData->SetPoints(points);
             fiberPolyData->SetLines(cellArray);
             vtkSmartPointer<vtkCleanPolyData> cleaner = vtkSmartPointer<vtkCleanPolyData>::New();
             cleaner->SetInputData(fiberPolyData);
             cleaner->Update();
             fiberPolyData = cleaner->GetOutput();
             FiberBundle::Pointer image = FiberBundle::New(fiberPolyData);
             result.push_back(image.GetPointer());
             return result;
         }
         else
         {
             MITK_INFO << "File is not deprectaed XML format.";
         }
 
         setlocale(LC_ALL, currLocale.c_str());
         MITK_INFO << "Fiber bundle read";
     }
     catch(...)
     {
         throw;
     }
 
     throw "Selected file is no vtk readable fiber format (binary or ascii vtk or vtp file).";
 
     return result;
 }
diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/MLTracking/itkMLBSTrackingFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/MLTracking/itkMLBSTrackingFilter.cpp
index da4a53b7ed..e7b01fcb32 100644
--- a/Modules/DiffusionImaging/FiberTracking/Algorithms/MLTracking/itkMLBSTrackingFilter.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/MLTracking/itkMLBSTrackingFilter.cpp
@@ -1,793 +1,812 @@
 /*===================================================================
 
 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 __itkMLBSTrackingFilter_txx
 #define __itkMLBSTrackingFilter_txx
 
 #include <time.h>
 #include <stdio.h>
 #include <stdlib.h>
 
 #include "itkMLBSTrackingFilter.h"
 #include <itkImageRegionConstIterator.h>
 #include <itkImageRegionConstIteratorWithIndex.h>
 #include <itkImageRegionIterator.h>
 #include <itkImageFileWriter.h>
 //#include <boost/thread/thread.hpp>
 
 #define _USE_MATH_DEFINES
 #include <math.h>
 
 namespace itk {
 
 template< int NumImageFeatures >
 MLBSTrackingFilter< NumImageFeatures >
 ::MLBSTrackingFilter()
     : m_FiberPolyData(NULL)
     , m_Points(NULL)
     , m_Cells(NULL)
     , m_AngularThreshold(0.7)
     , m_StepSize(0)
     , m_MaxLength(10000)
     , m_MinTractLength(20.0)
     , m_MaxTractLength(400.0)
     , m_SeedsPerVoxel(1)
-    , m_UseDirection(true)
     , m_NumberOfSamples(50)
     , m_SamplingDistance(-1)
     , m_SeedImage(NULL)
     , m_MaskImage(NULL)
     , m_DecisionForest(NULL)
     , m_StoppingRegions(NULL)
     , m_DemoMode(false)
     , m_PauseTracking(false)
     , m_AbortTracking(false)
     , m_RemoveWmEndFibers(false)
     , m_AposterioriCurvCheck(false)
     , m_AvoidStop(true)
+    , m_RandomSampling(false)
 {
     this->SetNumberOfRequiredInputs(1);
 }
 
 template< int NumImageFeatures >
 double MLBSTrackingFilter< NumImageFeatures >
 ::RoundToNearest(double num) {
     return (num > 0.0) ? floor(num + 0.5) : ceil(num - 0.5);
 }
 
 template< int NumImageFeatures >
 void MLBSTrackingFilter< NumImageFeatures >::BeforeThreadedGenerateData()
 {
     m_InputImage = const_cast<InputImageType*>(this->GetInput(0));
     PreprocessRawData();
 
     m_FiberPolyData = PolyDataType::New();
     m_Points = vtkSmartPointer< vtkPoints >::New();
     m_Cells = vtkSmartPointer< vtkCellArray >::New();
 
     m_ImageSize.resize(3);
     m_ImageSize[0] = m_FeatureImage->GetLargestPossibleRegion().GetSize()[0];
     m_ImageSize[1] = m_FeatureImage->GetLargestPossibleRegion().GetSize()[1];
     m_ImageSize[2] = m_FeatureImage->GetLargestPossibleRegion().GetSize()[2];
 
     m_ImageSpacing.resize(3);
     m_ImageSpacing[0] = m_FeatureImage->GetSpacing()[0];
     m_ImageSpacing[1] = m_FeatureImage->GetSpacing()[1];
     m_ImageSpacing[2] = m_FeatureImage->GetSpacing()[2];
 
     double minSpacing;
     if(m_ImageSpacing[0]<m_ImageSpacing[1] && m_ImageSpacing[0]<m_ImageSpacing[2])
         minSpacing = m_ImageSpacing[0];
     else if (m_ImageSpacing[1] < m_ImageSpacing[2])
         minSpacing = m_ImageSpacing[1];
     else
         minSpacing = m_ImageSpacing[2];
 
     m_StepSize *= minSpacing;
     if (m_StepSize<mitk::eps)
         m_StepSize = 0.5*minSpacing;
 
     m_SamplingDistance *= minSpacing;
     if (m_SamplingDistance<mitk::eps)
         m_SamplingDistance = minSpacing*0.5;
 
     m_PolyDataContainer.clear();
     for (unsigned int i=0; i<this->GetNumberOfThreads(); i++)
     {
         PolyDataType poly = PolyDataType::New();
         m_PolyDataContainer.push_back(poly);
     }
 
     m_NotWmImage = ItkDoubleImgType::New();
     m_NotWmImage->SetSpacing( m_FeatureImage->GetSpacing() );
     m_NotWmImage->SetOrigin( m_FeatureImage->GetOrigin() );
     m_NotWmImage->SetDirection( m_FeatureImage->GetDirection() );
     m_NotWmImage->SetRegions( m_FeatureImage->GetLargestPossibleRegion() );
     m_NotWmImage->Allocate();
     m_NotWmImage->FillBuffer(0);
 
     m_WmImage = ItkDoubleImgType::New();
     m_WmImage->SetSpacing( m_FeatureImage->GetSpacing() );
     m_WmImage->SetOrigin( m_FeatureImage->GetOrigin() );
     m_WmImage->SetDirection( m_FeatureImage->GetDirection() );
     m_WmImage->SetRegions( m_FeatureImage->GetLargestPossibleRegion() );
     m_WmImage->Allocate();
     m_WmImage->FillBuffer(0);
 
     m_AvoidStopImage = ItkDoubleImgType::New();
     m_AvoidStopImage->SetSpacing( m_FeatureImage->GetSpacing() );
     m_AvoidStopImage->SetOrigin( m_FeatureImage->GetOrigin() );
     m_AvoidStopImage->SetDirection( m_FeatureImage->GetDirection() );
     m_AvoidStopImage->SetRegions( m_FeatureImage->GetLargestPossibleRegion() );
     m_AvoidStopImage->Allocate();
     m_AvoidStopImage->FillBuffer(0);
 
     if (m_StoppingRegions.IsNull())
     {
         m_StoppingRegions = ItkUcharImgType::New();
         m_StoppingRegions->SetSpacing( m_FeatureImage->GetSpacing() );
         m_StoppingRegions->SetOrigin( m_FeatureImage->GetOrigin() );
         m_StoppingRegions->SetDirection( m_FeatureImage->GetDirection() );
         m_StoppingRegions->SetRegions( m_FeatureImage->GetLargestPossibleRegion() );
         m_StoppingRegions->Allocate();
         m_StoppingRegions->FillBuffer(0);
     }
 
     if (m_SeedImage.IsNull())
     {
         m_SeedImage = ItkUcharImgType::New();
         m_SeedImage->SetSpacing( m_FeatureImage->GetSpacing() );
         m_SeedImage->SetOrigin( m_FeatureImage->GetOrigin() );
         m_SeedImage->SetDirection( m_FeatureImage->GetDirection() );
         m_SeedImage->SetRegions( m_FeatureImage->GetLargestPossibleRegion() );
         m_SeedImage->Allocate();
         m_SeedImage->FillBuffer(1);
     }
 
     if (m_MaskImage.IsNull())
     {
         // initialize mask image
         m_MaskImage = ItkUcharImgType::New();
         m_MaskImage->SetSpacing( m_FeatureImage->GetSpacing() );
         m_MaskImage->SetOrigin( m_FeatureImage->GetOrigin() );
         m_MaskImage->SetDirection( m_FeatureImage->GetDirection() );
         m_MaskImage->SetRegions( m_FeatureImage->GetLargestPossibleRegion() );
         m_MaskImage->Allocate();
         m_MaskImage->FillBuffer(1);
     }
     else
         std::cout << "MLBSTrackingFilter: using mask image" << std::endl;
 
     if (m_AngularThreshold<0.0)
         m_AngularThreshold = 0.5*minSpacing;
 
     m_BuildFibersReady = 0;
     m_BuildFibersFinished = false;
     m_Threads = 0;
     m_Tractogram.clear();
+    m_SamplingPointset = mitk::PointSet::New();
+    m_AlternativePointset = mitk::PointSet::New();
 
     std::cout << "MLBSTrackingFilter: Angular threshold: " << m_AngularThreshold << std::endl;
     std::cout << "MLBSTrackingFilter: Stepsize: " << m_StepSize << " mm" << std::endl;
     std::cout << "MLBSTrackingFilter: Seeds per voxel: " << m_SeedsPerVoxel << std::endl;
     std::cout << "MLBSTrackingFilter: Max. sampling distance: " << m_SamplingDistance << " mm" << std::endl;
     std::cout << "MLBSTrackingFilter: Number of samples: " << m_NumberOfSamples << std::endl;
     std::cout << "MLBSTrackingFilter: Max. tract length: " << m_MaxTractLength << " mm" << std::endl;
     std::cout << "MLBSTrackingFilter: Min. tract length: " << m_MinTractLength << " mm" << std::endl;
     std::cout << "MLBSTrackingFilter: Starting streamline tracking using " << this->GetNumberOfThreads() << " threads." << std::endl;
 }
 
 template< int NumImageFeatures >
 void MLBSTrackingFilter< NumImageFeatures >::PreprocessRawData()
 {
     typedef itk::AnalyticalDiffusionQballReconstructionImageFilter<short,short,float,6, 2*NumImageFeatures> InterpolationFilterType;
 
     std::cout << "MLBSTrackingFilter: Spherical signal interpolation and sampling ..." << std::endl;
     typename InterpolationFilterType::Pointer filter = InterpolationFilterType::New();
     filter->SetGradientImage( m_GradientDirections, m_InputImage );
     filter->SetBValue( m_B_Value );
     filter->SetLambda(0.006);
     filter->SetNormalizationMethod(InterpolationFilterType::QBAR_RAW_SIGNAL);
     filter->Update();
     //      FeatureImageType::Pointer itkFeatureImage = qballfilter->GetCoefficientImage();
     //      featureImageVector.push_back(itkFeatureImage);
 
     std::cout << "MLBSTrackingFilter: Creating feature image ..." << std::endl;
     vnl_vector_fixed<double,3> ref; ref.fill(0); ref[0]=1;
     itk::OrientationDistributionFunction< double, NumImageFeatures*2 > odf;
     m_DirectionIndices.clear();
     for (unsigned int f=0; f<NumImageFeatures*2; f++)
     {
         if (dot_product(ref, odf.GetDirection(f))>0)    // only used directions on one hemisphere
             m_DirectionIndices.push_back(f);
     }
 
     m_FeatureImage = FeatureImageType::New();
     m_FeatureImage->SetSpacing(filter->GetOutput()->GetSpacing());
     m_FeatureImage->SetOrigin(filter->GetOutput()->GetOrigin());
     m_FeatureImage->SetDirection(filter->GetOutput()->GetDirection());
     m_FeatureImage->SetLargestPossibleRegion(filter->GetOutput()->GetLargestPossibleRegion());
     m_FeatureImage->SetBufferedRegion(filter->GetOutput()->GetLargestPossibleRegion());
     m_FeatureImage->SetRequestedRegion(filter->GetOutput()->GetLargestPossibleRegion());
     m_FeatureImage->Allocate();
 
     itk::ImageRegionIterator< typename InterpolationFilterType::OutputImageType > it(filter->GetOutput(), filter->GetOutput()->GetLargestPossibleRegion());
     while(!it.IsAtEnd())
     {
         typename FeatureImageType::PixelType pix;
         for (unsigned int f=0; f<NumImageFeatures; f++)
             pix[f] = it.Get()[m_DirectionIndices.at(f)];
         m_FeatureImage->SetPixel(it.GetIndex(), pix);
         ++it;
     }
 }
 
 template< int NumImageFeatures >
 void MLBSTrackingFilter< NumImageFeatures >::CalculateNewPosition(itk::Point<double, 3>& pos, vnl_vector_fixed<double, 3>& dir)
 {
     //    vnl_matrix_fixed< double, 3, 3 > rot = m_FeatureImage->GetDirection().GetTranspose();
     //    dir = rot*dir;
 
     dir *= m_StepSize;
     pos[0] += dir[0];
     pos[1] += dir[1];
     pos[2] += dir[2];
 }
 
 template< int NumImageFeatures >
 bool MLBSTrackingFilter< NumImageFeatures >
 ::IsValidPosition(itk::Point<double, 3> &pos)
 {
     typename FeatureImageType::IndexType idx;
     m_FeatureImage->TransformPhysicalPointToIndex(pos, idx);
     if (!m_FeatureImage->GetLargestPossibleRegion().IsInside(idx) || m_MaskImage->GetPixel(idx)==0)
         return false;
 
     return true;
 }
 
 template< int NumImageFeatures >
 typename MLBSTrackingFilter< NumImageFeatures >::FeatureImageType::PixelType MLBSTrackingFilter< NumImageFeatures >::GetImageValues(itk::Point<float, 3> itkP)
 {
     itk::Index<3> idx;
     itk::ContinuousIndex< double, 3> cIdx;
     m_FeatureImage->TransformPhysicalPointToIndex(itkP, idx);
     m_FeatureImage->TransformPhysicalPointToContinuousIndex(itkP, cIdx);
 
     typename FeatureImageType::PixelType pix; pix.Fill(0.0);
     if ( m_FeatureImage->GetLargestPossibleRegion().IsInside(idx) )
         pix = m_FeatureImage->GetPixel(idx);
     else
         return pix;
 
     double frac_x = cIdx[0] - idx[0];
     double frac_y = cIdx[1] - idx[1];
     double frac_z = cIdx[2] - idx[2];
     if (frac_x<0)
     {
         idx[0] -= 1;
         frac_x += 1;
     }
     if (frac_y<0)
     {
         idx[1] -= 1;
         frac_y += 1;
     }
     if (frac_z<0)
     {
         idx[2] -= 1;
         frac_z += 1;
     }
     frac_x = 1-frac_x;
     frac_y = 1-frac_y;
     frac_z = 1-frac_z;
 
     // int coordinates inside image?
     if (idx[0] >= 0 && idx[0] < m_FeatureImage->GetLargestPossibleRegion().GetSize(0)-1 &&
             idx[1] >= 0 && idx[1] < m_FeatureImage->GetLargestPossibleRegion().GetSize(1)-1 &&
             idx[2] >= 0 && idx[2] < m_FeatureImage->GetLargestPossibleRegion().GetSize(2)-1)
     {
         vnl_vector_fixed<double, 8> interpWeights;
         interpWeights[0] = (  frac_x)*(  frac_y)*(  frac_z);
         interpWeights[1] = (1-frac_x)*(  frac_y)*(  frac_z);
         interpWeights[2] = (  frac_x)*(1-frac_y)*(  frac_z);
         interpWeights[3] = (  frac_x)*(  frac_y)*(1-frac_z);
         interpWeights[4] = (1-frac_x)*(1-frac_y)*(  frac_z);
         interpWeights[5] = (  frac_x)*(1-frac_y)*(1-frac_z);
         interpWeights[6] = (1-frac_x)*(  frac_y)*(1-frac_z);
         interpWeights[7] = (1-frac_x)*(1-frac_y)*(1-frac_z);
 
         pix = m_FeatureImage->GetPixel(idx) * interpWeights[0];
         typename FeatureImageType::IndexType tmpIdx = idx; tmpIdx[0]++;
         pix +=  m_FeatureImage->GetPixel(tmpIdx) * interpWeights[1];
         tmpIdx = idx; tmpIdx[1]++;
         pix +=  m_FeatureImage->GetPixel(tmpIdx) * interpWeights[2];
         tmpIdx = idx; tmpIdx[2]++;
         pix +=  m_FeatureImage->GetPixel(tmpIdx) * interpWeights[3];
         tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++;
         pix +=  m_FeatureImage->GetPixel(tmpIdx) * interpWeights[4];
         tmpIdx = idx; tmpIdx[1]++; tmpIdx[2]++;
         pix +=  m_FeatureImage->GetPixel(tmpIdx) * interpWeights[5];
         tmpIdx = idx; tmpIdx[2]++; tmpIdx[0]++;
         pix +=  m_FeatureImage->GetPixel(tmpIdx) * interpWeights[6];
         tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++; tmpIdx[2]++;
         pix +=  m_FeatureImage->GetPixel(tmpIdx) * interpWeights[7];
     }
 
     return pix;
 }
 
 template< int NumImageFeatures >
 vnl_vector_fixed<double,3> MLBSTrackingFilter< NumImageFeatures >::Classify(itk::Point<double, 3>& pos, int& candidates, vnl_vector_fixed<double,3>& olddir, double angularThreshold, double& prob, bool avoidStop)
 {
     vnl_vector_fixed<double,3> direction; direction.fill(0);
 
-    vigra::MultiArray<2, double> featureData;
-    if(m_UseDirection)
-        featureData = vigra::MultiArray<2, double>( vigra::Shape2(1,NumImageFeatures+3) );
-    else
-        featureData = vigra::MultiArray<2, double>( vigra::Shape2(1,NumImageFeatures) );
+    vigra::MultiArray<2, double> featureData = vigra::MultiArray<2, double>( vigra::Shape2(1,NumImageFeatures+3) );
 
     typename FeatureImageType::PixelType featurePixel = GetImageValues(pos);
 
     // pixel values
     for (unsigned int f=0; f<NumImageFeatures; f++)
         featureData(0,f) = featurePixel[f];
 
     // direction features
     int c = 0;
-    if (m_UseDirection)
     {
         vnl_vector_fixed<double,3> ref; ref.fill(0); ref[0]=1;
         for (unsigned int f=NumImageFeatures; f<NumImageFeatures+3; f++)
         {
             if (dot_product(ref, olddir)<0)
                 featureData(0,f) = -olddir[c];
             else
                 featureData(0,f) = olddir[c];
             c++;
         }
     }
 
     vigra::MultiArray<2, double> probs(vigra::Shape2(1, m_DecisionForest->class_count()));
     m_DecisionForest->predictProbabilities(featureData, probs);
 
     double outProb = 0;
     prob = 0;
     candidates = 0; // directions with probability > 0
     for (int i=0; i<m_DecisionForest->class_count(); i++)
     {
         if (probs(0,i)>0)
         {
             int classLabel = 0;
             m_DecisionForest->ext_param_.to_classlabel(i, classLabel);
 
             if (classLabel<m_DirectionIndices.size())
             {
                 candidates++;
 
                 vnl_vector_fixed<double,3> d = m_ODF.GetDirection(m_DirectionIndices.at(classLabel));
                 double dot = dot_product(d, olddir);
 
                 if (olddir.magnitude()>0)
                 {
                     if (fabs(dot)>angularThreshold)
                     {
                         if (dot<0)
                             d *= -1;
                         dot = fabs(dot);
                         direction += probs(0,i)*dot*d;
                         prob += probs(0,i)*dot;
                     }
                 }
                 else
                 {
                     direction += probs(0,i)*d;
                     prob += probs(0,i);
                 }
             }
             else
                 outProb += probs(0,i);
         }
     }
 
     ItkDoubleImgType::IndexType idx;
     m_NotWmImage->TransformPhysicalPointToIndex(pos, idx);
     if (m_NotWmImage->GetLargestPossibleRegion().IsInside(idx))
     {
         m_NotWmImage->SetPixel(idx, m_NotWmImage->GetPixel(idx)+outProb);
         m_WmImage->SetPixel(idx, m_WmImage->GetPixel(idx)+prob);
     }
     if (outProb>prob && prob>0)
     {
         candidates = 0;
         prob = 0;
         direction.fill(0.0);
     }
     if (avoidStop && m_AvoidStopImage->GetLargestPossibleRegion().IsInside(idx) && candidates>0 && direction.magnitude()>0.001)
         m_AvoidStopImage->SetPixel(idx, m_AvoidStopImage->GetPixel(idx)+0.1);
 
     return direction;
 }
 
 
 template< int NumImageFeatures >
 double MLBSTrackingFilter< NumImageFeatures >::GetRandDouble(double min, double max)
 {
     return (double)(rand()%((int)(10000*(max-min))) + 10000*min)/10000;
 }
 
 template< int NumImageFeatures >
 vnl_vector_fixed<double,3> MLBSTrackingFilter< NumImageFeatures >::GetNewDirection(itk::Point<double, 3> &pos, vnl_vector_fixed<double, 3>& olddir)
 {
+    if (m_DemoMode)
+    {
+        m_SamplingPointset->Clear();
+        m_AlternativePointset->Clear();
+    }
     vnl_vector_fixed<double,3> direction; direction.fill(0);
 
     ItkUcharImgType::IndexType idx;
     m_StoppingRegions->TransformPhysicalPointToIndex(pos, idx);
     if (m_StoppingRegions->GetPixel(idx)>0)
         return direction;
 
     if (olddir.magnitude()>0)
         olddir.normalize();
 
     int candidates = 0; // number of directions with probability > 0
     double prob = 0;
     direction = Classify(pos, candidates, olddir, m_AngularThreshold, prob); // sample neighborhood
     direction *= prob;
 
 
     itk::OrientationDistributionFunction< double, 50 >  probeVecs;
 
+    itk::Point<double, 3> sample_pos;
+    int alternatives = 1;
     for (int i=0; i<m_NumberOfSamples; i++)
     {
         vnl_vector_fixed<double,3> d;
-//        probe[0] = GetRandDouble()*m_SamplingDistance;
-//        probe[1] = GetRandDouble()*m_SamplingDistance;
-//        probe[2] = GetRandDouble()*m_SamplingDistance;
 
-        d = probeVecs.GetDirection(i)*m_SamplingDistance;
+        if (m_RandomSampling)
+        {
+            d[0] = GetRandDouble();
+            d[1] = GetRandDouble();
+            d[2] = GetRandDouble();
+            d.normalize();
+            d *= GetRandDouble(0,m_SamplingDistance);
+        }
+        else
+        {
+            d = probeVecs.GetDirection(i)*m_SamplingDistance;
+        }
 
-        itk::Point<double, 3> sample_pos;
         sample_pos[0] = pos[0] + d[0];
         sample_pos[1] = pos[1] + d[1];
         sample_pos[2] = pos[2] + d[2];
+        if(m_DemoMode)
+            m_SamplingPointset->InsertPoint(i, sample_pos);
 
         candidates = 0;
         vnl_vector_fixed<double,3> tempDir = Classify(sample_pos, candidates, olddir, m_AngularThreshold, prob); // sample neighborhood
         if (candidates>0 && tempDir.magnitude()>0.001)
         {
             direction += tempDir*prob;
         }
         else if (m_AvoidStop && candidates==0 && olddir.magnitude()>0) // out of white matter
         {
             double dot = dot_product(d, olddir);
             if (dot >= 0.0) // 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]*2;
-            sample_pos[1] = pos[1] + d[1]*2;
-            sample_pos[2] = pos[2] + d[2]*2;
+            sample_pos[0] = pos[0] + d[0];
+            sample_pos[1] = pos[1] + d[1];
+            sample_pos[2] = pos[2] + d[2];
+            if(m_DemoMode)
+                m_AlternativePointset->InsertPoint(alternatives, sample_pos);
+            alternatives++;
             candidates = 0;
             vnl_vector_fixed<double,3> tempDir = Classify(sample_pos, candidates, olddir, m_AngularThreshold, prob, true); // sample neighborhood
 
             if (candidates>0 && tempDir.magnitude()>0.001)  // are we back in the white matter?
             {
                 direction += d;         // go into the direction of the white matter
                 direction += tempDir*prob;  // go into the direction of the white matter direction at this location
             }
         }
     }
 
     if (direction.magnitude()>0.001)
     {
         direction.normalize();
         olddir[0] = direction[0];
         olddir[1] = direction[1];
         olddir[2] = direction[2];
     }
     else
         direction.fill(0);
 
     return direction;
 }
 
 template< int NumImageFeatures >
 double MLBSTrackingFilter< NumImageFeatures >::FollowStreamline(ThreadIdType threadId, itk::Point<double, 3> pos, vnl_vector_fixed<double,3> dir, FiberType* fib, double tractLength, bool front)
 {
     vnl_vector_fixed<double,3> dirOld = dir;
     dirOld = dir;
 
     for (int step=0; step< m_MaxLength/2; step++)
     {
-        while (m_PauseTracking){}
-        if (m_DemoMode)
-        {
-            m_Mutex.Lock();
-            m_BuildFibersReady++;
-            m_Tractogram.push_back(*fib);
-            BuildFibers(true);
-            m_Stop = true;
-            m_Mutex.Unlock();
-            while (m_Stop){}
-        }
 
         // get new position
         CalculateNewPosition(pos, dir);
 
         // is new position inside of image and mask
         if (!IsValidPosition(pos) || m_AbortTracking)   // if not end streamline
         {
             return tractLength;
         }
         else    // if yes, add new point to streamline
         {
             tractLength +=  m_StepSize;
             if (front)
                 fib->push_front(pos);
             else
                 fib->push_back(pos);
 
             if (m_AposterioriCurvCheck)
             {
                 int curv = CheckCurvature(fib, front);  // TODO: Move into classification ???
                 if (curv>0)
                 {
                     tractLength -= m_StepSize*curv;
                     while (curv>0)
                     {
                         if (front)
                             fib->pop_front();
                         else
                             fib->pop_back();
                         curv--;
                     }
                     return tractLength;
                 }
             }
 
             if (tractLength>m_MaxTractLength)
                 return tractLength;
         }
+        if (m_DemoMode) // CHECK: warum sind die samplingpunkte der streamline in der visualisierung immer einen schritt voras?
+        {
+            m_Mutex.Lock();
+            m_BuildFibersReady++;
+            m_Tractogram.push_back(*fib);
+            BuildFibers(true);
+            m_Stop = true;
+            m_Mutex.Unlock();
 
+            while (m_Stop){
+            }
+        }
         dir = GetNewDirection(pos, dirOld);
 
+        while (m_PauseTracking){}
+
         if (dir.magnitude()<0.0001)
             return tractLength;
     }
     return tractLength;
 }
 
 template< int NumImageFeatures >
 int MLBSTrackingFilter<NumImageFeatures>::CheckCurvature(FiberType* fib, bool front)
 {
     double m_Distance = 5;
     if (fib->size()<3)
         return 0;
 
     double dist = 0;
     std::vector< vnl_vector_fixed< float, 3 > > vectors;
     vnl_vector_fixed< float, 3 > meanV; meanV.fill(0);
     double dev = 0;
 
     if (front)
     {
         int c=0;
         while(dist<m_Distance && c<fib->size()-1)
         {
             itk::Point<double> p1 = fib->at(c);
             itk::Point<double> p2 = fib->at(c+1);
 
             vnl_vector_fixed< float, 3 > v;
             v[0] = p2[0]-p1[0];
             v[1] = p2[1]-p1[1];
             v[2] = p2[2]-p1[2];
             dist += v.magnitude();
             v.normalize();
             vectors.push_back(v);
             if (c==0)
                 meanV += v;
             c++;
         }
     }
     else
     {
         int c=fib->size()-1;
         while(dist<m_Distance && c>0)
         {
             itk::Point<double> p1 = fib->at(c);
             itk::Point<double> p2 = fib->at(c-1);
 
             vnl_vector_fixed< float, 3 > v;
             v[0] = p2[0]-p1[0];
             v[1] = p2[1]-p1[1];
             v[2] = p2[2]-p1[2];
             dist += v.magnitude();
             v.normalize();
             vectors.push_back(v);
             if (c==fib->size()-1)
                 meanV += v;
             c--;
         }
     }
     meanV.normalize();
 
     for (int c=0; c<vectors.size(); c++)
     {
         double angle = dot_product(meanV, vectors.at(c));
         if (angle>1.0)
             angle = 1.0;
         if (angle<-1.0)
             angle = -1.0;
         dev += acos(angle)*180/M_PI;
     }
     if (vectors.size()>0)
         dev /= vectors.size();
 
     if (dev<30)
         return 0;
     else
         return vectors.size();
 }
 
 template< int NumImageFeatures >
 void MLBSTrackingFilter< NumImageFeatures >::ThreadedGenerateData(const InputImageRegionType &regionForThread, ThreadIdType threadId)
 {
     m_Mutex.Lock();
     m_Threads++;
     m_Mutex.Unlock();
     typedef ImageRegionConstIterator< ItkUcharImgType >     MaskIteratorType;
     MaskIteratorType    sit(m_SeedImage, regionForThread );
     MaskIteratorType    mit(m_MaskImage, regionForThread );
 
     sit.GoToBegin();
     mit.GoToBegin();
     itk::Point<double> worldPos;
     while( !sit.IsAtEnd() )
     {
         if (sit.Value()==0 || mit.Value()==0)
         {
             ++sit;
             ++mit;
             continue;
         }
 
         for (int s=0; s<m_SeedsPerVoxel; s++)
         {
             FiberType fib;
             double tractLength = 0;
             typename FeatureImageType::IndexType index = sit.GetIndex();
             itk::ContinuousIndex<double, 3> start;
             unsigned int counter = 0;
 
             if (m_SeedsPerVoxel>1)
             {
                 start[0] = index[0]+GetRandDouble(-0.5, 0.5);
                 start[1] = index[1]+GetRandDouble(-0.5, 0.5);
                 start[2] = index[2]+GetRandDouble(-0.5, 0.5);
             }
             else
             {
                 start[0] = index[0];
                 start[1] = index[1];
                 start[2] = index[2];
             }
 
             // get staring position
             m_SeedImage->TransformContinuousIndexToPhysicalPoint( start, worldPos );
 
             // get starting direction
             int candidates = 0;
             double prob = 0;
             vnl_vector_fixed<double,3> dirOld; dirOld.fill(0.0);
             vnl_vector_fixed<double,3> dir = Classify(worldPos, candidates, dirOld, 0, prob);
             if (dir.magnitude()<0.0001)
                 continue;
 
             // forward tracking
             tractLength = FollowStreamline(threadId, worldPos, dir, &fib, 0, false);
             fib.push_front(worldPos);
 
             if (m_RemoveWmEndFibers)
             {
                 itk::Point<double> check = fib.back();
                 dirOld.fill(0.0);
                 vnl_vector_fixed<double,3> check2 = GetNewDirection(check, dirOld);
                 if (check2.magnitude()>0.001)
                 {
                     MITK_INFO << "Detected WM ending. Discarding fiber.";
                     continue;
                 }
             }
 
             // backward tracking
             tractLength = FollowStreamline(threadId, worldPos, -dir, &fib, tractLength, true);
             counter = fib.size();
 
             if (m_RemoveWmEndFibers)
             {
                 itk::Point<double> check = fib.front();
                 dirOld.fill(0.0);
                 vnl_vector_fixed<double,3> check2 = GetNewDirection(check, dirOld);
                 if (check2.magnitude()>0.001)
                 {
                     MITK_INFO << "Detected WM ending. Discarding fiber.";
                     continue;
                 }
             }
 
             if (tractLength<m_MinTractLength || counter<2)
                 continue;
 
             m_Mutex.Lock();
             m_Tractogram.push_back(fib);
             m_Mutex.Unlock();
 
             if (m_AbortTracking)
                 break;
         }
         if (m_AbortTracking)
             break;
         ++sit;
         ++mit;
     }
     m_Threads--;
     std::cout << "Thread " << threadId << " finished tracking" << std::endl;
 }
 
 template< int NumImageFeatures >
 void MLBSTrackingFilter< NumImageFeatures >::BuildFibers(bool check)
 {
     if (m_BuildFibersReady<m_Threads && check)
         return;
 
     m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
     vtkSmartPointer<vtkCellArray> vNewLines = vtkSmartPointer<vtkCellArray>::New();
     vtkSmartPointer<vtkPoints> vNewPoints = vtkSmartPointer<vtkPoints>::New();
 
     for (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++)
+        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;
 }
 
 template< int NumImageFeatures >
 void MLBSTrackingFilter< NumImageFeatures >::AfterThreadedGenerateData()
 {
     MITK_INFO << "Generating polydata ";
     BuildFibers(false);
     MITK_INFO << "done";
 }
 
 }
 
 #endif // __itkDiffusionQballPrincipleDirectionsImageFilter_txx
diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/MLTracking/itkMLBSTrackingFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/MLTracking/itkMLBSTrackingFilter.h
index 5e1a62d0e2..14f612463e 100644
--- a/Modules/DiffusionImaging/FiberTracking/Algorithms/MLTracking/itkMLBSTrackingFilter.h
+++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/MLTracking/itkMLBSTrackingFilter.h
@@ -1,191 +1,194 @@
 /*===================================================================
 
 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.
 
 ===================================================================*/
 
 /*===================================================================
 
 This file is based heavily on a corresponding ITK filter.
 
 ===================================================================*/
 #ifndef __itkMLBSTrackingFilter_h_
 #define __itkMLBSTrackingFilter_h_
 
 #include <itkImageToImageFilter.h>
 #include <itkVectorContainer.h>
 #include <itkVectorImage.h>
 #include <vtkSmartPointer.h>
 #include <vtkPolyData.h>
 #include <vtkCellArray.h>
 #include <vtkPoints.h>
 #include <vtkPolyLine.h>
 #include <vtkCleanPolyData.h>
 #include <itkOrientationDistributionFunction.h>
 #include <itkMersenneTwisterRandomVariateGenerator.h>
 #include <itkAnalyticalDiffusionQballReconstructionImageFilter.h>
 #include <itkSimpleFastMutexLock.h>
 #include <mitkDiffusionPropertyHelper.h>
+#include <mitkPointSet.h>
 
 // classification includes
 #include <vigra/random_forest.hxx>
 #include <vigra/multi_array.hxx>
 #include <vigra/random_forest_hdf5_impex.hxx>
 
 namespace itk{
 
 /**
 * \brief Performes deterministic streamline tracking on the input tensor image.   */
 
 template< int NumImageFeatures=100 >
 class MLBSTrackingFilter : public ImageToImageFilter< VectorImage< short, 3 >, Image< double, 3 > >
 {
 
 public:
 
     typedef MLBSTrackingFilter Self;
     typedef SmartPointer<Self>                      Pointer;
     typedef SmartPointer<const Self>                ConstPointer;
     typedef ImageToImageFilter< VectorImage< short, 3 >, Image< double, 3 > > Superclass;
 
     typedef vigra::RandomForest<int>                        DecisionForestType;
     typedef typename Superclass::InputImageType             InputImageType;
     typedef typename Superclass::InputImageRegionType       InputImageRegionType;
     typedef Image< Vector< float, NumImageFeatures > , 3 >     FeatureImageType;
 
     /** Method for creation through the object factory. */
     itkFactorylessNewMacro(Self)
     itkCloneMacro(Self)
 
     /** Runtime information support. */
     itkTypeMacro(MLBSTrackingFilter, ImageToImageFilter)
 
     typedef itk::Image<unsigned char, 3>                ItkUcharImgType;
     typedef itk::Image<double, 3>                       ItkDoubleImgType;
     typedef itk::Image<float, 3>                        ItkFloatImgType;
     typedef vtkSmartPointer< vtkPolyData >              PolyDataType;
 
     typedef std::deque< itk::Point<double> > FiberType;
     typedef std::vector< FiberType > BundleType;
 
-    bool    m_PauseTracking;
+    volatile bool    m_PauseTracking;
     bool    m_AbortTracking;
     bool    m_BuildFibersFinished;
     int     m_BuildFibersReady;
-    bool    m_Stop;
+    volatile bool m_Stop;
+    mitk::PointSet::Pointer             m_SamplingPointset;
+    mitk::PointSet::Pointer             m_AlternativePointset;
 //    void RequestFibers(){ m_Stop=true; m_BuildFibersReady=0; m_BuildFibersFinished=false; }
 
     itkGetMacro( FiberPolyData, PolyDataType )          ///< Output fibers
     itkSetMacro( SeedImage, ItkUcharImgType::Pointer)   ///< Seeds are only placed inside of this mask.
     itkSetMacro( MaskImage, ItkUcharImgType::Pointer)   ///< Tracking is only performed inside of this mask image.
     itkSetMacro( SeedsPerVoxel, int)                    ///< One seed placed in the center of each voxel or multiple seeds randomly placed inside each voxel.
     itkSetMacro( StepSize, double)                      ///< Integration step size in mm
     itkSetMacro( MinTractLength, double )               ///< Shorter tracts are discarded.
     itkSetMacro( MaxTractLength, double )
     itkSetMacro( AngularThreshold, double )
-    itkSetMacro( UseDirection, bool )
     itkSetMacro( SamplingDistance, double )
     itkSetMacro( NumberOfSamples, int )
     itkSetMacro( StoppingRegions, ItkUcharImgType::Pointer)
     itkSetMacro( B_Value, float )
     itkSetMacro( GradientDirections, mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer )
     itkSetMacro( DemoMode, bool )
     itkSetMacro( RemoveWmEndFibers, bool )
     itkSetMacro( AposterioriCurvCheck, bool )
     itkSetMacro( AvoidStop, bool )
+    itkSetMacro( RandomSampling, bool )
 
     void SetDecisionForest( DecisionForestType* forest )
     {
         m_DecisionForest = forest;
     }
 
     itkGetMacro( WmImage, ItkDoubleImgType::Pointer )
     itkGetMacro( NotWmImage, ItkDoubleImgType::Pointer )
     itkGetMacro( AvoidStopImage, ItkDoubleImgType::Pointer )
 
     protected:
         MLBSTrackingFilter();
     ~MLBSTrackingFilter() {}
 
     void CalculateNewPosition(itk::Point<double, 3>& pos, vnl_vector_fixed<double,3>& dir);    ///< Calculate next integration step.
     double FollowStreamline(ThreadIdType threadId, itk::Point<double, 3> pos, vnl_vector_fixed<double,3> dir, FiberType* fib, double tractLength, bool front);       ///< Start streamline in one direction.
     bool IsValidPosition(itk::Point<double, 3>& pos);   ///< Are we outside of the mask image?
     vnl_vector_fixed<double,3> GetNewDirection(itk::Point<double, 3>& pos, vnl_vector_fixed<double,3>& olddir);
     vnl_vector_fixed<double,3> Classify(itk::Point<double, 3>& pos, int& candidates, vnl_vector_fixed<double,3>& olddir, double angularThreshold, double& prob, bool avoidStop=false);
 
     typename FeatureImageType::PixelType GetImageValues(itk::Point<float, 3> itkP);
     double GetRandDouble(double min=-1, double max=1);
     double RoundToNearest(double num);
 
     void BeforeThreadedGenerateData();
     void PreprocessRawData();
     void ThreadedGenerateData( const InputImageRegionType &outputRegionForThread, ThreadIdType threadId);
     void AfterThreadedGenerateData();
 
     PolyDataType                        m_FiberPolyData;
     vtkSmartPointer<vtkPoints>          m_Points;
     vtkSmartPointer<vtkCellArray>       m_Cells;
     BundleType                          m_Tractogram;
 
     double                              m_AngularThreshold;
     double                              m_StepSize;
     int                                 m_MaxLength;
     double                              m_MinTractLength;
     double                              m_MaxTractLength;
     int                                 m_SeedsPerVoxel;
-    bool                                m_UseDirection;
+    bool                                m_RandomSampling;
     double                              m_SamplingDistance;
     int                                 m_NumberOfSamples;
     std::vector< int >                  m_ImageSize;
     std::vector< double >               m_ImageSpacing;
 
     SimpleFastMutexLock                 m_Mutex;
     ItkUcharImgType::Pointer            m_StoppingRegions;
     ItkDoubleImgType::Pointer           m_WmImage;
     ItkDoubleImgType::Pointer           m_NotWmImage;
     ItkDoubleImgType::Pointer           m_AvoidStopImage;
     ItkUcharImgType::Pointer            m_SeedImage;
     ItkUcharImgType::Pointer            m_MaskImage;
     typename FeatureImageType::Pointer  m_FeatureImage;
     typename InputImageType::Pointer    m_InputImage;
     mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer m_GradientDirections;
     float                               m_B_Value;
 
     bool                                m_AposterioriCurvCheck;
     bool                                m_RemoveWmEndFibers;
     bool                                m_AvoidStop;
 
     int                                 m_Threads;
     bool                                m_DemoMode;
     void BuildFibers(bool check);
     int CheckCurvature(FiberType* fib, bool front);
 
     // decision forest
     DecisionForestType*                                                 m_DecisionForest;
     itk::OrientationDistributionFunction< double, NumImageFeatures*2 >  m_ODF;
     std::vector< int >                                                  m_DirectionIndices;
 
     std::vector< PolyDataType >         m_PolyDataContainer;
 
 private:
 
 };
 
 }
 
 #ifndef ITK_MANUAL_INSTANTIATION
 #include "itkMLBSTrackingFilter.cpp"
 #endif
 
 #endif //__itkMLBSTrackingFilter_h_
 
diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/MLTracking/mitkTrackingForestHandler.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/MLTracking/mitkTrackingForestHandler.cpp
index aa50256bdb..566b167309 100644
--- a/Modules/DiffusionImaging/FiberTracking/Algorithms/MLTracking/mitkTrackingForestHandler.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/MLTracking/mitkTrackingForestHandler.cpp
@@ -1,540 +1,477 @@
 /*===================================================================
 
 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 _TrackingForestHandler_cpp
 #define _TrackingForestHandler_cpp
 
 #include "mitkTrackingForestHandler.h"
 #include <itkTractDensityImageFilter.h>
 #include <mitkDiffusionPropertyHelper.h>
 
 namespace mitk
 {
 
 template< int NumberOfSignalFeatures >
 TrackingForestHandler< NumberOfSignalFeatures >::TrackingForestHandler()
     : m_GrayMatterSamplesPerVoxel(50)
     , m_StepSize(-1)
-    , m_UsePreviousDirection(true)
     , m_NumTrees(30)
     , m_MaxTreeDepth(50)
     , m_SampleFraction(1.0)
 {
 
 }
 
 template< int NumberOfSignalFeatures >
 TrackingForestHandler< NumberOfSignalFeatures >::~TrackingForestHandler()
 {
 
 }
 
 template< int NumberOfSignalFeatures >
 typename TrackingForestHandler< NumberOfSignalFeatures >::InterpolatedRawImageType::PixelType TrackingForestHandler< NumberOfSignalFeatures >::GetImageValues(itk::Point<float, 3> itkP, typename InterpolatedRawImageType::Pointer image)
 {
     itk::Index<3> idx;
     itk::ContinuousIndex< double, 3> cIdx;
     image->TransformPhysicalPointToIndex(itkP, idx);
     image->TransformPhysicalPointToContinuousIndex(itkP, cIdx);
 
     typename InterpolatedRawImageType::PixelType pix; pix.Fill(0.0);
     if ( image->GetLargestPossibleRegion().IsInside(idx) )
         pix = image->GetPixel(idx);
     else
         return pix;
 
     double frac_x = cIdx[0] - idx[0];
     double frac_y = cIdx[1] - idx[1];
     double frac_z = cIdx[2] - idx[2];
     if (frac_x<0)
     {
         idx[0] -= 1;
         frac_x += 1;
     }
     if (frac_y<0)
     {
         idx[1] -= 1;
         frac_y += 1;
     }
     if (frac_z<0)
     {
         idx[2] -= 1;
         frac_z += 1;
     }
     frac_x = 1-frac_x;
     frac_y = 1-frac_y;
     frac_z = 1-frac_z;
 
     // int coordinates inside image?
     if (idx[0] >= 0 && idx[0] < image->GetLargestPossibleRegion().GetSize(0)-1 &&
             idx[1] >= 0 && idx[1] < image->GetLargestPossibleRegion().GetSize(1)-1 &&
             idx[2] >= 0 && idx[2] < image->GetLargestPossibleRegion().GetSize(2)-1)
     {
         vnl_vector_fixed<double, 8> interpWeights;
         interpWeights[0] = (  frac_x)*(  frac_y)*(  frac_z);
         interpWeights[1] = (1-frac_x)*(  frac_y)*(  frac_z);
         interpWeights[2] = (  frac_x)*(1-frac_y)*(  frac_z);
         interpWeights[3] = (  frac_x)*(  frac_y)*(1-frac_z);
         interpWeights[4] = (1-frac_x)*(1-frac_y)*(  frac_z);
         interpWeights[5] = (  frac_x)*(1-frac_y)*(1-frac_z);
         interpWeights[6] = (1-frac_x)*(  frac_y)*(1-frac_z);
         interpWeights[7] = (1-frac_x)*(1-frac_y)*(1-frac_z);
 
 
         pix = image->GetPixel(idx) * interpWeights[0];
         typename InterpolatedRawImageType::IndexType tmpIdx = idx; tmpIdx[0]++;
         pix +=  image->GetPixel(tmpIdx) * interpWeights[1];
         tmpIdx = idx; tmpIdx[1]++;
         pix +=  image->GetPixel(tmpIdx) * interpWeights[2];
         tmpIdx = idx; tmpIdx[2]++;
         pix +=  image->GetPixel(tmpIdx) * interpWeights[3];
         tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++;
         pix +=  image->GetPixel(tmpIdx) * interpWeights[4];
         tmpIdx = idx; tmpIdx[1]++; tmpIdx[2]++;
         pix +=  image->GetPixel(tmpIdx) * interpWeights[5];
         tmpIdx = idx; tmpIdx[2]++; tmpIdx[0]++;
         pix +=  image->GetPixel(tmpIdx) * interpWeights[6];
         tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++; tmpIdx[2]++;
         pix +=  image->GetPixel(tmpIdx) * interpWeights[7];
     }
 
     return pix;
 }
 
 template< int NumberOfSignalFeatures >
 void TrackingForestHandler< NumberOfSignalFeatures >::InputDataValidForTracking()
 {
     if (m_RawData.empty())
         mitkThrow() << "No diffusion-weighted images set!";
 }
 
 template< int NumberOfSignalFeatures >
 void TrackingForestHandler< NumberOfSignalFeatures >::StartTraining()
 {
     InputDataValidForTraining();
     PreprocessInputData();
     CalculateFeatures();
     TrainForest();
 }
 
 template< int NumberOfSignalFeatures >
 void TrackingForestHandler< NumberOfSignalFeatures >::InputDataValidForTraining()
 {
     if (m_RawData.empty())
         mitkThrow() << "No diffusion-weighted images set!";
     if (m_Tractograms.empty())
         mitkThrow() << "No tractograms set!";
     if (m_RawData.size()!=m_Tractograms.size())
         mitkThrow() << "Unequal number of diffusion-weighted images and tractograms detected!";
 }
 
 template< int NumberOfSignalFeatures >
 void TrackingForestHandler< NumberOfSignalFeatures >::PreprocessInputData()
 {
     typedef itk::AnalyticalDiffusionQballReconstructionImageFilter<short,short,float,6, 2*NumberOfSignalFeatures> InterpolationFilterType;
 
     MITK_INFO << "Spherical signal interpolation and sampling ...";
     for (unsigned int i=0; i<m_RawData.size(); i++)
     {
         typename InterpolationFilterType::Pointer qballfilter = InterpolationFilterType::New();
         qballfilter->SetGradientImage( mitk::DiffusionPropertyHelper::GetGradientContainer(m_RawData.at(i)), mitk::DiffusionPropertyHelper::GetItkVectorImage(m_RawData.at(i)) );
         qballfilter->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(m_RawData.at(i)));
         qballfilter->SetLambda(0.006);
         qballfilter->SetNormalizationMethod(InterpolationFilterType::QBAR_RAW_SIGNAL);
         qballfilter->Update();
         //      FeatureImageType::Pointer itkFeatureImage = qballfilter->GetCoefficientImage();
         //      featureImageVector.push_back(itkFeatureImage);
         m_InterpolatedRawImages.push_back(qballfilter->GetOutput());
 
         if (i>=m_MaskImages.size())
         {
             ItkUcharImgType::Pointer newMask = ItkUcharImgType::New();
             newMask->SetSpacing( m_InterpolatedRawImages.at(i)->GetSpacing() );
             newMask->SetOrigin( m_InterpolatedRawImages.at(i)->GetOrigin() );
             newMask->SetDirection( m_InterpolatedRawImages.at(i)->GetDirection() );
             newMask->SetLargestPossibleRegion( m_InterpolatedRawImages.at(i)->GetLargestPossibleRegion() );
             newMask->SetBufferedRegion( m_InterpolatedRawImages.at(i)->GetLargestPossibleRegion() );
             newMask->SetRequestedRegion( m_InterpolatedRawImages.at(i)->GetLargestPossibleRegion() );
             newMask->Allocate();
             newMask->FillBuffer(1);
             m_MaskImages.push_back(newMask);
         }
     }
 
     MITK_INFO << "Resampling fibers and calculating number of samples ...";
     m_NumberOfSamples = 0;
     for (unsigned int t=0; t<m_Tractograms.size(); t++)
     {
         ItkUcharImgType::Pointer mask = m_MaskImages.at(t);
         ItkUcharImgType::Pointer wmmask;
         if (t<m_WhiteMatterImages.size())
             wmmask = m_WhiteMatterImages.at(t);
         else
         {
             MITK_INFO << "Calculating fiber envelope";
             itk::TractDensityImageFilter< ItkUcharImgType >::Pointer env = itk::TractDensityImageFilter< ItkUcharImgType >::New();
             env->SetFiberBundle(m_Tractograms.at(t));
             env->SetInputImage(mask);
             env->SetBinaryOutput(true);
             env->SetUseImageGeometry(true);
             env->Update();
             wmmask = env->GetOutput();
             m_WhiteMatterImages.push_back(wmmask);
         }
 
         itk::ImageRegionConstIterator<ItkUcharImgType> it(wmmask, wmmask->GetLargestPossibleRegion());
         int OUTOFWM = 0;
         while(!it.IsAtEnd())
         {
             if (it.Get()==0 && mask->GetPixel(it.GetIndex())>0)
                 OUTOFWM++;
             ++it;
         }
         m_NumberOfSamples += m_GrayMatterSamplesPerVoxel*OUTOFWM;
         MITK_INFO << "Samples outside of WM: " << m_NumberOfSamples;
 
         if (m_StepSize<0)
         {
             typename InterpolatedRawImageType::Pointer image = m_InterpolatedRawImages.at(t);
             float minSpacing = 1;
             if(image->GetSpacing()[0]<image->GetSpacing()[1] && image->GetSpacing()[0]<image->GetSpacing()[2])
                 minSpacing = image->GetSpacing()[0];
             else if (image->GetSpacing()[1] < image->GetSpacing()[2])
                 minSpacing = image->GetSpacing()[1];
             else
                 minSpacing = image->GetSpacing()[2];
             m_StepSize = minSpacing*0.5;
         }
 
         m_Tractograms.at(t)->ResampleSpline(m_StepSize);
         m_NumberOfSamples += m_Tractograms.at(t)->GetNumberOfPoints();
         m_NumberOfSamples -= 2*m_Tractograms.at(t)->GetNumFibers();
     }
     MITK_INFO << "Number of samples: " << m_NumberOfSamples;
 }
 
 template< int NumberOfSignalFeatures >
 void TrackingForestHandler< NumberOfSignalFeatures >::CalculateFeatures()
 {
     vnl_vector_fixed<double,3> ref; ref.fill(0); ref[0]=1;
     itk::OrientationDistributionFunction< double, 2*NumberOfSignalFeatures > directions;
     std::vector< int > directionIndices;
     for (unsigned int f=0; f<2*NumberOfSignalFeatures; f++)
     {
         if (dot_product(ref, directions.GetDirection(f))>0)
             directionIndices.push_back(f);
     }
 
-    int numDirectionFeatures = 0;
-    if (m_UsePreviousDirection)
-        numDirectionFeatures = 3;
+    int numDirectionFeatures = 3;
 
     m_FeatureData.reshape( vigra::Shape2(m_NumberOfSamples, NumberOfSignalFeatures+numDirectionFeatures) );
     m_LabelData.reshape( vigra::Shape2(m_NumberOfSamples,1) );
     MITK_INFO << "Number of features: " << m_FeatureData.shape(1);
 
     itk::Statistics::MersenneTwisterRandomVariateGenerator::Pointer m_RandGen = itk::Statistics::MersenneTwisterRandomVariateGenerator::New();
     m_RandGen->SetSeed();
     MITK_INFO <<  "Creating training data ...";
     int sampleCounter = 0;
     for (unsigned int t=0; t<m_Tractograms.size(); t++)
     {
         typename InterpolatedRawImageType::Pointer image = m_InterpolatedRawImages.at(t);
         ItkUcharImgType::Pointer wmMask = m_WhiteMatterImages.at(t);
         ItkUcharImgType::Pointer mask;
         if (t<m_MaskImages.size())
             mask = m_MaskImages.at(t);
 
         // non-white matter samples
         itk::ImageRegionConstIterator<ItkUcharImgType> it(wmMask, wmMask->GetLargestPossibleRegion());
         while(!it.IsAtEnd())
         {
             if (it.Get()==0 && (mask.IsNull() || (mask.IsNotNull() && mask->GetPixel(it.GetIndex())>0)))
             {
                 typename InterpolatedRawImageType::PixelType pix = image->GetPixel(it.GetIndex());
-                if (m_UsePreviousDirection)
-                {
+
                     // null direction
                     for (unsigned int f=0; f<NumberOfSignalFeatures; f++)
                     {
                         m_FeatureData(sampleCounter,f) = pix[directionIndices[f]];
                         if(m_FeatureData(sampleCounter,f)!=m_FeatureData(sampleCounter,f))
                             m_FeatureData(sampleCounter,f) = 0;
                     }
                     m_FeatureData(sampleCounter,NumberOfSignalFeatures) = 0;
                     m_FeatureData(sampleCounter,NumberOfSignalFeatures+1) = 0;
                     m_FeatureData(sampleCounter,NumberOfSignalFeatures+2) = 0;
                     m_LabelData(sampleCounter,0) = directionIndices.size();
                     sampleCounter++;
 
                     // random directions
                     for (int i=1; i<m_GrayMatterSamplesPerVoxel; i++)
                     {
                         for (unsigned int f=0; f<NumberOfSignalFeatures; f++)
                         {
                             m_FeatureData(sampleCounter,f) = pix[directionIndices[f]];
                             if(m_FeatureData(sampleCounter,f)!=m_FeatureData(sampleCounter,f))
                                 m_FeatureData(sampleCounter,f) = 0;
                         }
                         int c=0;
                         vnl_vector_fixed<double,3> probe;
                         probe[0] = m_RandGen->GetVariate()*2-1;
                         probe[1] = m_RandGen->GetVariate()*2-1;
                         probe[2] = m_RandGen->GetVariate()*2-1;
                         probe.normalize();
                         if (dot_product(ref, probe)<0)
                             probe *= -1;
                         for (unsigned int f=NumberOfSignalFeatures; f<NumberOfSignalFeatures+3; f++)
                         {
                             m_FeatureData(sampleCounter,f) = probe[c];
                             c++;
                         }
                         m_LabelData(sampleCounter,0) = directionIndices.size();
                         sampleCounter++;
                     }
-                }
-                else
-                {
-                    itk::ContinuousIndex<double, 3> idx;
-                    idx[0] = it.GetIndex()[0];
-                    idx[1] = it.GetIndex()[1];
-                    idx[2] = it.GetIndex()[2];
-                    itk::Point<float, 3> itkP1;
-                    image->TransformContinuousIndexToPhysicalPoint(idx, itkP1);
-                    typename InterpolatedRawImageType::PixelType pix = GetImageValues(itkP1, image);;
-                    for (unsigned int f=0; f<NumberOfSignalFeatures; f++)
-                    {
-                        m_FeatureData(sampleCounter,f) = pix[directionIndices[f]];
-                        if(m_FeatureData(sampleCounter,f)!=m_FeatureData(sampleCounter,f))
-                            m_FeatureData(sampleCounter,f) = 0;
-                    }
-                    m_LabelData(sampleCounter,0) = directionIndices.size();
-                    sampleCounter++;
-
-                    for (int i=1; i<m_GrayMatterSamplesPerVoxel; i++)
-                    {
-                        itk::ContinuousIndex<double, 3> idx;
-                        idx[0] = it.GetIndex()[0] + m_RandGen->GetVariate()-0.5;
-                        idx[1] = it.GetIndex()[1] + m_RandGen->GetVariate()-0.5;
-                        idx[2] = it.GetIndex()[2] + m_RandGen->GetVariate()-0.5;
-                        itk::Point<float, 3> itkP1;
-                        image->TransformContinuousIndexToPhysicalPoint(idx, itkP1);
-                        typename InterpolatedRawImageType::PixelType pix = GetImageValues(itkP1, image);;
-                        for (unsigned int f=0; f<NumberOfSignalFeatures; f++)
-                        {
-                            m_FeatureData(sampleCounter,f) = pix[directionIndices[f]];
-                            if(m_FeatureData(sampleCounter,f)!=m_FeatureData(sampleCounter,f))
-                                m_FeatureData(sampleCounter,f) = 0;
-                        }
-                        m_LabelData(sampleCounter,0) = directionIndices.size();
-                        sampleCounter++;
-                    }
-                }
             }
             ++it;
         }
 
         // white matter samples
         mitk::FiberBundle::Pointer fib = m_Tractograms.at(t);
         vtkSmartPointer< vtkPolyData > polyData = fib->GetFiberPolyData();
         for (int i=0; i<fib->GetNumFibers(); i++)
         {
             vtkCell* cell = polyData->GetCell(i);
             int numPoints = cell->GetNumberOfPoints();
             vtkPoints* points = cell->GetPoints();
 
             vnl_vector_fixed<double,3> dirOld; dirOld.fill(0.0);
 
             for (int j=0; j<numPoints-1; j++)
             {
                 // calculate direction
                 double* p1 = points->GetPoint(j);
                 itk::Point<float, 3> itkP1;
                 itkP1[0] = p1[0]; itkP1[1] = p1[1]; itkP1[2] = p1[2];
 
                 vnl_vector_fixed<double,3> dir; dir.fill(0.0);
 
                 itk::Point<float, 3> itkP2;
                 double* p2 = points->GetPoint(j+1);
                 itkP2[0] = p2[0]; itkP2[1] = p2[1]; itkP2[2] = p2[2];
                 dir[0]=itkP2[0]-itkP1[0];
                 dir[1]=itkP2[1]-itkP1[1];
                 dir[2]=itkP2[2]-itkP1[2];
 
                 if (dir.magnitude()<0.0001)
                 {
                     MITK_INFO << "streamline error!";
                     continue;
                 }
                 dir.normalize();
                 if (dir[0]!=dir[0] || dir[1]!=dir[1] || dir[2]!=dir[2])
                 {
                     MITK_INFO << "ERROR: NaN direction!";
                     continue;
                 }
 
                 if (j==0)
                 {
                     dirOld = dir;
                     continue;
                 }
 
                 // get voxel values
                 typename InterpolatedRawImageType::PixelType pix = GetImageValues(itkP1, image);
                 for (unsigned int f=0; f<NumberOfSignalFeatures; f++)
                     m_FeatureData(sampleCounter,f) = pix[directionIndices[f]];
 
                 // direction training features
                 int c = 0;
                 if (dot_product(ref, dirOld)<0)
                     dirOld *= -1;
 
-                if (m_UsePreviousDirection)
-                    for (unsigned int f=NumberOfSignalFeatures; f<NumberOfSignalFeatures+3; f++)
-                    {
-                        m_FeatureData(sampleCounter,f) = dirOld[c];
-                        c++;
-                    }
+                for (unsigned int f=NumberOfSignalFeatures; f<NumberOfSignalFeatures+3; f++)
+                {
+                    m_FeatureData(sampleCounter,f) = dirOld[c];
+                    c++;
+                }
 
                 // set label values
                 double angle = 0;
                 double m = dir.magnitude();
                 if (m>0.0001)
                 {
                     int label = 0;
                     for (unsigned int f=0; f<NumberOfSignalFeatures; f++)
                     {
                         double a = fabs(dot_product(dir, directions.GetDirection(directionIndices[f])));
                         if (a>angle)
                         {
                             m_LabelData(sampleCounter,0) = f;
                             angle = a;
                             label = f;
                         }
                     }
                 }
 
                 dirOld = dir;
                 sampleCounter++;
             }
         }
     }
 }
 
 template< int NumberOfSignalFeatures >
 void TrackingForestHandler< NumberOfSignalFeatures >::TrainForest()
 {
+
     MITK_INFO << "Maximum tree depths: " << m_MaxTreeDepth;
     MITK_INFO << "Sample fraction per tree: " << m_SampleFraction;
     MITK_INFO << "Number of trees: " << m_NumTrees;
-    bool random_split = false;
-    vigra::rf::visitors::OOB_Error oob_v;
-    MITK_INFO << "Create Split Function";
-    //   typedef ThresholdSplit<int, vigra::BestGiniOfColumn<vigra::GiniCriterion>, vigra::ClassificationTag> DefaultSplitType;
-
-    m_Forest.set_options().use_stratification(vigra::RF_NONE); // How the data should be made equal
-    m_Forest.set_options().sample_with_replacement(true); // if sampled with replacement or not
-    m_Forest.set_options().samples_per_tree(m_SampleFraction); // Fraction of samples that are used to train a tree
-    m_Forest.set_options().tree_count(1); // Number of trees that are calculated;
-    m_Forest.set_options().min_split_node_size(5); // Minimum number of datapoints that must be in a node
-    //    rf.set_options().features_per_node(10);
-
-    m_Forest.learn(m_FeatureData, m_LabelData, vigra::rf::visitors::create_visitor(oob_v));
-
-    // Prepare parallel VariableImportance Calculation
-    int numMod =  m_FeatureData.shape(1);
-    const int numClass = 2 + 2;
-
-    float** varImp = new float*[numMod];
-
-    for(int i = 0; i < numMod; i++)
-        varImp[i] = new float[numClass];
-
-    for (int i = 0; i < numMod; ++i)
-        for (int j = 0; j < numClass; ++j)
-            varImp[i][j] = 0.0;
 
+    std::vector< vigra::RandomForest<int> > trees;
+    int count = 0;
 #pragma omp parallel for
-    for (int i = 0; i < m_NumTrees - 1; ++i)
+    for (int i = 0; i < m_NumTrees; ++i)
     {
         vigra::RandomForest<int> lrf;
         vigra::rf::visitors::OOB_Error loob_v;
 
         lrf.set_options().use_stratification(vigra::RF_NONE); // How the data should be made equal
         lrf.set_options().sample_with_replacement(true); // if sampled with replacement or not
         lrf.set_options().samples_per_tree(m_SampleFraction); // Fraction of samples that are used to train a tree
         lrf.set_options().tree_count(1); // Number of trees that are calculated;
         lrf.set_options().min_split_node_size(5); // Minimum number of datapoints that must be in a node
-        //        lrf.set_options().features_per_node(10);
-
-        vigra::rf::visitors::VariableImportanceVisitor lvariableImportance;
-        lrf.learn(m_FeatureData, m_LabelData, vigra::rf::visitors::create_visitor(loob_v));
+        lrf.ext_param_.max_tree_depth = m_MaxTreeDepth;
 
+        lrf.learn(m_FeatureData, m_LabelData);
 #pragma omp critical
         {
-            m_Forest.trees_.push_back(lrf.trees_[0]);
+            count++;
+            MITK_INFO << "Tree " << count << " finished training.";
+            trees.push_back(lrf);
         }
     }
 
+    for (int i = 1; i < m_NumTrees; ++i)
+        trees.at(0).trees_.push_back(trees.at(i).trees_[0]);
+
+    m_Forest = trees.at(0);
     m_Forest.options_.tree_count_ = m_NumTrees;
     MITK_INFO << "Training finsihed";
-    MITK_INFO << "The out-of-bag error is: " << oob_v.oob_breiman << std::endl;
 }
 
 template< int NumberOfSignalFeatures >
 void TrackingForestHandler< NumberOfSignalFeatures >::SaveForest(std::string forestFile)
 {
     MITK_INFO << "Saving forest to " << forestFile;
     vigra::rf_export_HDF5( m_Forest, forestFile, "" );
 }
 
 template< int NumberOfSignalFeatures >
 void TrackingForestHandler< NumberOfSignalFeatures >::LoadForest(std::string forestFile)
 {
     MITK_INFO << "Loading forest from " << forestFile;
     vigra::rf_import_HDF5(m_Forest, forestFile);
 }
 
 //// superclass implementations
 //template< int NumberOfSignalFeatures >
 //void TrackingForestHandler< NumberOfSignalFeatures >::UpdateOutputInformation()
 //{
 
 //}
 //template< int NumberOfSignalFeatures >
 //void TrackingForestHandler< NumberOfSignalFeatures >::SetRequestedRegionToLargestPossibleRegion()
 //{
 
 //}
 //template< int NumberOfSignalFeatures >
 //bool TrackingForestHandler< NumberOfSignalFeatures >::RequestedRegionIsOutsideOfTheBufferedRegion()
 //{
 //    return false;
 //}
 //template< int NumberOfSignalFeatures >
 //bool TrackingForestHandler< NumberOfSignalFeatures >::VerifyRequestedRegion()
 //{
 //    return true;
 //}
 //template< int NumberOfSignalFeatures >
 //void TrackingForestHandler< NumberOfSignalFeatures >::SetRequestedRegion(const itk::DataObject* )
 //{
 
 //}
 
 }
 
 #endif
diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/MLTracking/mitkTrackingForestHandler.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/MLTracking/mitkTrackingForestHandler.h
index 4aecb1ac41..5c34f6a0a8 100644
--- a/Modules/DiffusionImaging/FiberTracking/Algorithms/MLTracking/mitkTrackingForestHandler.h
+++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/MLTracking/mitkTrackingForestHandler.h
@@ -1,115 +1,113 @@
 /*===================================================================
 
 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 _TrackingForestHandler
 #define _TrackingForestHandler
 
 #include "mitkBaseData.h"
 
 #include <mitkFiberBundle.h>
 #include <itkAnalyticalDiffusionQballReconstructionImageFilter.h>
 #include <itkOrientationDistributionFunction.h>
 #include <itkMersenneTwisterRandomVariateGenerator.h>
 
 // classification includes
 //#include "RegressionForestClasses.hxx"
 #undef DIFFERENCE
 #define VIGRA_STATIC_LIB
 #include <vigra/random_forest.hxx>
 #include <vigra/multi_array.hxx>
 #include <vigra/random_forest_hdf5_impex.hxx>
 
 #define _USE_MATH_DEFINES
 #include <math.h>
 
 namespace mitk
 {
 
 /**
 * \brief  */
 
 template< int NumberOfSignalFeatures=100 >
 class TrackingForestHandler
 {
 
 public:
 
     TrackingForestHandler();
     ~TrackingForestHandler();
 
     typedef itk::Image<unsigned char, 3>    ItkUcharImgType;
     typedef itk::Image< itk::Vector< float, NumberOfSignalFeatures*2 > , 3 > InterpolatedRawImageType;
 
     void SetRawData( std::vector< Image::Pointer > images ){ m_RawData = images; }
     void SetTractograms( std::vector< FiberBundle::Pointer > tractograms )
     {
         m_Tractograms.clear();
         for (int i=0; i<tractograms.size(); i++)
         {
             m_Tractograms.push_back(tractograms.at(i)->GetDeepCopy());
         }
     }
     void SetMaskImages( std::vector< ItkUcharImgType::Pointer > images ){ m_MaskImages = images; }
     void SetWhiteMatterImages( std::vector< ItkUcharImgType::Pointer > images ){ m_WhiteMatterImages = images; }
 
     void StartTraining();
     void SaveForest(std::string forestFile);
     void LoadForest(std::string forestFile);
 
     void SetNumTrees(int num){ m_NumTrees = num; }
     void SetMaxTreeDepth(int depth){ m_MaxTreeDepth = depth; }
-    void SetUsePreviousDirection(bool use){ m_UsePreviousDirection = use; }
     void SetStepSize(double step){ m_StepSize = step; }
     void SetGrayMatterSamplesPerVoxel(int samples){ m_GrayMatterSamplesPerVoxel = samples; }
     void SetSampleFraction(double fraction){ m_SampleFraction = fraction; }
     vigra::RandomForest<int> GetForest(){ return m_Forest; }
 
 protected:
 
     void InputDataValidForTracking();
     void InputDataValidForTraining();
     void PreprocessInputData();
     void CalculateFeatures();
     void TrainForest();
 
     int         m_GrayMatterSamplesPerVoxel;
     double      m_StepSize;
-    bool        m_UsePreviousDirection;
     int         m_NumTrees;
     int         m_MaxTreeDepth;
     double      m_SampleFraction;
 
     std::vector< Image::Pointer > m_RawData;
     std::vector< FiberBundle::Pointer > m_Tractograms;
     std::vector< ItkUcharImgType::Pointer > m_MaskImages;
     std::vector< ItkUcharImgType::Pointer > m_WhiteMatterImages;
     std::vector< ItkUcharImgType::Pointer > m_SeedImages;
     std::vector< ItkUcharImgType::Pointer > m_StopImages;
 
     int m_NumberOfSamples;
     vigra::RandomForest<int> m_Forest;
     vigra::MultiArray<2, double> m_FeatureData;
     vigra::MultiArray<2, double> m_LabelData;
     std::vector< typename InterpolatedRawImageType::Pointer > m_InterpolatedRawImages;
 
     typename InterpolatedRawImageType::PixelType GetImageValues(itk::Point<float, 3> itkP, typename InterpolatedRawImageType::Pointer image);
 };
 
 }
 
 #include "mitkTrackingForestHandler.cpp"
 
 #endif
diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.cpp
index 4d434011cb..d844388d70 100644
--- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.cpp
@@ -1,565 +1,575 @@
 #include "itkTractsToVectorImageFilter.h"
 
 // VTK
 #include <vtkPolyLine.h>
 #include <vtkCellArray.h>
 #include <vtkCellData.h>
 
 // ITK
 #include <itkTimeProbe.h>
 #include <itkImageRegionIterator.h>
 
 // misc
 #define _USE_MATH_DEFINES
 #include <math.h>
 #include <boost/progress.hpp>
 
 
 namespace itk{
 
 static bool CompareVectorLengths(const vnl_vector_fixed< double, 3 >& v1, const vnl_vector_fixed< double, 3 >& v2)
 {
     return (v1.magnitude()>v2.magnitude());
 }
 
 template< class PixelType >
 TractsToVectorImageFilter< PixelType >::TractsToVectorImageFilter():
     m_AngularThreshold(0.7),
     m_Epsilon(0.999),
     m_MaskImage(NULL),
     m_NormalizeVectors(false),
     m_UseWorkingCopy(true),
     m_MaxNumDirections(3),
     m_SizeThreshold(0.3),
     m_NumDirectionsImage(NULL),
     m_CreateDirectionImages(true)
 {
     this->SetNumberOfRequiredOutputs(1);
 }
 
 
 template< class PixelType >
 TractsToVectorImageFilter< PixelType >::~TractsToVectorImageFilter()
 {
 }
 
 
 template< class PixelType >
 vnl_vector_fixed<double, 3> TractsToVectorImageFilter< PixelType >::GetVnlVector(double point[])
 {
     vnl_vector_fixed<double, 3> vnlVector;
     vnlVector[0] = point[0];
     vnlVector[1] = point[1];
     vnlVector[2] = point[2];
     return vnlVector;
 }
 
 
 template< class PixelType >
 itk::Point<double, 3> TractsToVectorImageFilter< PixelType >::GetItkPoint(double point[])
 {
     itk::Point<double, 3> itkPoint;
     itkPoint[0] = point[0];
     itkPoint[1] = point[1];
     itkPoint[2] = point[2];
     return itkPoint;
 }
 
 template< class PixelType >
 void TractsToVectorImageFilter< PixelType >::GenerateData()
 {
     mitk::BaseGeometry::Pointer geometry = m_FiberBundle->GetGeometry();
 
     // calculate new image parameters
     itk::Vector<double> spacing;
     itk::Point<double> origin;
     itk::Matrix<double, 3, 3> direction;
     ImageRegion<3> imageRegion;
     if (!m_MaskImage.IsNull())
     {
         spacing = m_MaskImage->GetSpacing();
         imageRegion = m_MaskImage->GetLargestPossibleRegion();
         origin = m_MaskImage->GetOrigin();
         direction = m_MaskImage->GetDirection();
     }
     else
     {
         spacing = geometry->GetSpacing();
         origin = geometry->GetOrigin();
         mitk::BaseGeometry::BoundsArrayType bounds = geometry->GetBounds();
         origin[0] += bounds.GetElement(0);
         origin[1] += bounds.GetElement(2);
         origin[2] += bounds.GetElement(4);
 
         for (int i=0; i<3; i++)
             for (int j=0; j<3; j++)
                 direction[j][i] = geometry->GetMatrixColumn(i)[j];
         imageRegion.SetSize(0, geometry->GetExtent(0));
         imageRegion.SetSize(1, geometry->GetExtent(1));
         imageRegion.SetSize(2, geometry->GetExtent(2));
 
 
         m_MaskImage = ItkUcharImgType::New();
         m_MaskImage->SetSpacing( spacing );
         m_MaskImage->SetOrigin( origin );
         m_MaskImage->SetDirection( direction );
         m_MaskImage->SetRegions( imageRegion );
         m_MaskImage->Allocate();
         m_MaskImage->FillBuffer(1);
     }
     OutputImageType::RegionType::SizeType outImageSize = imageRegion.GetSize();
     m_OutImageSpacing = m_MaskImage->GetSpacing();
     m_ClusteredDirectionsContainer = ContainerType::New();
 
     // initialize num directions image
     m_NumDirectionsImage = ItkUcharImgType::New();
     m_NumDirectionsImage->SetSpacing( spacing );
     m_NumDirectionsImage->SetOrigin( origin );
     m_NumDirectionsImage->SetDirection( direction );
     m_NumDirectionsImage->SetRegions( imageRegion );
     m_NumDirectionsImage->Allocate();
     m_NumDirectionsImage->FillBuffer(0);
 
     // initialize direction images
     m_DirectionImageContainer = DirectionImageContainerType::New();
 
     // resample fiber bundle
     double minSpacing = 1;
     if(m_OutImageSpacing[0]<m_OutImageSpacing[1] && m_OutImageSpacing[0]<m_OutImageSpacing[2])
         minSpacing = m_OutImageSpacing[0];
     else if (m_OutImageSpacing[1] < m_OutImageSpacing[2])
         minSpacing = m_OutImageSpacing[1];
     else
         minSpacing = m_OutImageSpacing[2];
 
     if (m_UseWorkingCopy)
         m_FiberBundle = m_FiberBundle->GetDeepCopy();
 
     // resample fiber bundle for sufficient voxel coverage
     m_FiberBundle->ResampleSpline(minSpacing/10);
 
     // iterate over all fibers
     vtkSmartPointer<vtkPolyData> fiberPolyData = m_FiberBundle->GetFiberPolyData();
     int numFibers = m_FiberBundle->GetNumFibers();
     m_DirectionsContainer = ContainerType::New();
 
+    VectorContainer< unsigned int, std::vector< double > >::Pointer peakLengths = VectorContainer< unsigned int, std::vector< double > >::New();
+
     MITK_INFO << "Generating directions from tractogram";
     boost::progress_display disp(numFibers);
     for( int i=0; i<numFibers; i++ )
     {
         ++disp;
         vtkCell* cell = fiberPolyData->GetCell(i);
         int numPoints = cell->GetNumberOfPoints();
         vtkPoints* points = cell->GetPoints();
         if (numPoints<2)
             continue;
 
         vnl_vector_fixed<double, 3> dir;
         itk::Point<double, 3> worldPos;
         vnl_vector<double> v;
 
+
+        float fiberWeight = m_FiberBundle->GetFiberWeight(i);
+
         for( int j=0; j<numPoints-1; j++)
         {
             // get current position along fiber in world coordinates
             double* temp = points->GetPoint(j);
             worldPos = GetItkPoint(temp);
             itk::Index<3> index;
             m_MaskImage->TransformPhysicalPointToIndex(worldPos, index);
             if (!m_MaskImage->GetLargestPossibleRegion().IsInside(index) || m_MaskImage->GetPixel(index)==0)
                 continue;
 
             // get fiber tangent direction at this position
             v = GetVnlVector(temp);
             dir = GetVnlVector(points->GetPoint(j+1))-v;
             if (dir.is_zero())
                 continue;
             dir.normalize();
 
             // add direction to container
             unsigned int idx = index[0] + outImageSize[0]*(index[1] + outImageSize[1]*index[2]);
             DirectionContainerType::Pointer dirCont;
             if (m_DirectionsContainer->IndexExists(idx))
             {
+                peakLengths->ElementAt(idx).push_back(fiberWeight);
+
                 dirCont = m_DirectionsContainer->GetElement(idx);
                 if (dirCont.IsNull())
                 {
                     dirCont = DirectionContainerType::New();
                     dirCont->push_back(dir);
                     m_DirectionsContainer->InsertElement(idx, dirCont);
                 }
                 else
                     dirCont->push_back(dir);
             }
             else
             {
                 dirCont = DirectionContainerType::New();
                 dirCont->push_back(dir);
                 m_DirectionsContainer->InsertElement(idx, dirCont);
+
+                std::vector< double > lengths; lengths.push_back(fiberWeight);
+                peakLengths->InsertElement(idx, lengths);
             }
         }
     }
 
     vtkSmartPointer<vtkCellArray> m_VtkCellArray = vtkSmartPointer<vtkCellArray>::New();
     vtkSmartPointer<vtkPoints>    m_VtkPoints = vtkSmartPointer<vtkPoints>::New();
 
     itk::ImageRegionIterator<ItkUcharImgType> dirIt(m_NumDirectionsImage, m_NumDirectionsImage->GetLargestPossibleRegion());
 
     MITK_INFO << "Clustering directions";
     boost::progress_display disp2(outImageSize[0]*outImageSize[1]*outImageSize[2]);
     while(!dirIt.IsAtEnd())
     {
         ++disp2;
         OutputImageType::IndexType index = dirIt.GetIndex();
         int idx = index[0]+(index[1]+index[2]*outImageSize[1])*outImageSize[0];
 
         if (!m_DirectionsContainer->IndexExists(idx))
         {
             ++dirIt;
             continue;
         }
         DirectionContainerType::Pointer dirCont = m_DirectionsContainer->GetElement(idx);
         if (dirCont.IsNull() || dirCont->empty())
         {
             ++dirIt;
             continue;
         }
 
-        std::vector< double > lengths; lengths.resize(dirCont->size(), 1);  // all peaks have size 1
+//        std::vector< double > lengths; lengths.resize(dirCont->size(), 1);  // all peaks have size 1
         DirectionContainerType::Pointer directions;
         if (m_MaxNumDirections>0)
         {
-            directions = FastClustering(dirCont, lengths);
+            directions = FastClustering(dirCont, peakLengths->GetElement(idx));
             std::sort( directions->begin(), directions->end(), CompareVectorLengths );
         }
         else
             directions = dirCont;
 
         unsigned int numDir = directions->size();
         if (m_MaxNumDirections>0 && numDir>m_MaxNumDirections)
             numDir = m_MaxNumDirections;
 
         int count = 0;
         for (unsigned int i=0; i<numDir; i++)
         {
             vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
             itk::ContinuousIndex<double, 3> center;
             center[0] = index[0];
             center[1] = index[1];
             center[2] = index[2];
             itk::Point<double> worldCenter;
             m_MaskImage->TransformContinuousIndexToPhysicalPoint( center, worldCenter );
             DirectionType dir = directions->at(i);
 
             if (dir.magnitude()<m_SizeThreshold)
                 continue;
             if (m_NormalizeVectors)
                 dir.normalize();
             count++;
 
             if (m_CreateDirectionImages && i<10)
             {
                 if (i==m_DirectionImageContainer->size())
                 {
                     ItkDirectionImageType::Pointer directionImage = ItkDirectionImageType::New();
                     directionImage->SetSpacing( spacing );
                     directionImage->SetOrigin( origin );
                     directionImage->SetDirection( direction );
                     directionImage->SetRegions( imageRegion );
                     directionImage->Allocate();
                     Vector< float, 3 > nullVec; nullVec.Fill(0.0);
                     directionImage->FillBuffer(nullVec);
                     m_DirectionImageContainer->InsertElement(i, directionImage);
                 }
 
                 // set direction image pixel
                 ItkDirectionImageType::Pointer directionImage = m_DirectionImageContainer->GetElement(i);
                 Vector< float, 3 > pixel;
                 pixel.SetElement(0, dir[0]);
                 pixel.SetElement(1, dir[1]);
                 pixel.SetElement(2, dir[2]);
                 directionImage->SetPixel(index, pixel);
             }
 
             // add direction to vector field (with spacing compensation)
             itk::Point<double> worldStart;
             worldStart[0] = worldCenter[0]-dir[0]/2*minSpacing;
             worldStart[1] = worldCenter[1]-dir[1]/2*minSpacing;
             worldStart[2] = worldCenter[2]-dir[2]/2*minSpacing;
             vtkIdType id = m_VtkPoints->InsertNextPoint(worldStart.GetDataPointer());
             container->GetPointIds()->InsertNextId(id);
             itk::Point<double> worldEnd;
             worldEnd[0] = worldCenter[0]+dir[0]/2*minSpacing;
             worldEnd[1] = worldCenter[1]+dir[1]/2*minSpacing;
             worldEnd[2] = worldCenter[2]+dir[2]/2*minSpacing;
             id = m_VtkPoints->InsertNextPoint(worldEnd.GetDataPointer());
             container->GetPointIds()->InsertNextId(id);
             m_VtkCellArray->InsertNextCell(container);
         }
         dirIt.Set(count);
         ++dirIt;
     }
 
     vtkSmartPointer<vtkPolyData> directionsPolyData = vtkSmartPointer<vtkPolyData>::New();
     directionsPolyData->SetPoints(m_VtkPoints);
     directionsPolyData->SetLines(m_VtkCellArray);
     m_OutputFiberBundle = mitk::FiberBundle::New(directionsPolyData);
 }
 
 
 template< class PixelType >
 TractsToVectorImageFilter< PixelType >::DirectionContainerType::Pointer TractsToVectorImageFilter< PixelType >::FastClustering(DirectionContainerType::Pointer inDirs, std::vector< double > lengths)
 {
     DirectionContainerType::Pointer outDirs = DirectionContainerType::New();
     if (inDirs->size()<2)
         return inDirs;
 
     DirectionType oldMean, currentMean;
     std::vector< int > touched;
 
     // initialize
     touched.resize(inDirs->size(), 0);
     bool free = true;
     currentMean = inDirs->at(0);  // initialize first seed
     currentMean.normalize();
     double length = lengths.at(0);
     touched[0] = 1;
     std::vector< double > newLengths;
     bool meanChanged = false;
 
     double max = 0;
     while (free)
     {
         oldMean.fill(0.0);
 
         // start mean-shift clustering
         double angle = 0;
 
         while (fabs(dot_product(currentMean, oldMean))<0.99)
         {
             oldMean = currentMean;
             currentMean.fill(0.0);
             for (unsigned int i=0; i<inDirs->size(); i++)
             {
                 angle = dot_product(oldMean, inDirs->at(i));
                 if (angle>=m_AngularThreshold)
                 {
                     currentMean += inDirs->at(i);
                     if (meanChanged)
                         length += lengths.at(i);
                     touched[i] = 1;
                     meanChanged = true;
                 }
                 else if (-angle>=m_AngularThreshold)
                 {
                     currentMean -= inDirs->at(i);
                     if (meanChanged)
                         length += lengths.at(i);
                     touched[i] = 1;
                     meanChanged = true;
                 }
             }
             if(!meanChanged)
                 currentMean = oldMean;
             else
                 currentMean.normalize();
         }
 
         // found stable mean
         outDirs->push_back(currentMean);
         newLengths.push_back(length);
         if (length>max)
             max = length;
 
         // find next unused seed
         free = false;
         for (unsigned int i=0; i<touched.size(); i++)
             if (touched[i]==0)
             {
                 currentMean = inDirs->at(i);
                 free = true;
                 meanChanged = false;
                 length = lengths.at(i);
                 touched[i] = 1;
                 break;
             }
     }
 
     if (inDirs->size()==outDirs->size())
     {
         if (max>0)
             for (unsigned int i=0; i<outDirs->size(); i++)
                 outDirs->SetElement(i, outDirs->at(i)*newLengths.at(i)/max);
         return outDirs;
     }
     else
         return FastClustering(outDirs, newLengths);
 }
 
 
 //template< class PixelType >
 //std::vector< DirectionType > TractsToVectorImageFilter< PixelType >::Clustering(std::vector< DirectionType >& inDirs)
 //{
 //    std::vector< DirectionType > outDirs;
 //    if (inDirs.empty())
 //        return outDirs;
 //    DirectionType oldMean, currentMean, workingMean;
 
 //    std::vector< DirectionType > normalizedDirs;
 //    std::vector< int > touched;
 //    for (std::size_t i=0; i<inDirs.size(); i++)
 //    {
 //        normalizedDirs.push_back(inDirs[i]);
 //        normalizedDirs.back().normalize();
 //    }
 
 //    // initialize
 //    double max = 0.0;
 //    touched.resize(inDirs.size(), 0);
 //    for (std::size_t j=0; j<inDirs.size(); j++)
 //    {
 //        currentMean = inDirs[j];
 //        oldMean.fill(0.0);
 
 //        // start mean-shift clustering
 //        double angle = 0.0;
 //        int counter = 0;
 //        while ((currentMean-oldMean).magnitude()>0.0001)
 //        {
 //            counter = 0;
 //            oldMean = currentMean;
 //            workingMean = oldMean;
 //            workingMean.normalize();
 //            currentMean.fill(0.0);
 //            for (std::size_t i=0; i<normalizedDirs.size(); i++)
 //            {
 //                angle = dot_product(workingMean, normalizedDirs[i]);
 //                if (angle>=m_AngularThreshold)
 //                {
 //                    currentMean += inDirs[i];
 //                    counter++;
 //                }
 //                else if (-angle>=m_AngularThreshold)
 //                {
 //                    currentMean -= inDirs[i];
 //                    counter++;
 //                }
 //            }
 //        }
 
 //        // found stable mean
 //        if (counter>0)
 //        {
 //            bool add = true;
 //            DirectionType normMean = currentMean;
 //            normMean.normalize();
 //            for (std::size_t i=0; i<outDirs.size(); i++)
 //            {
 //                DirectionType dir = outDirs[i];
 //                dir.normalize();
 //                if ((normMean-dir).magnitude()<=0.0001)
 //                {
 //                    add = false;
 //                    break;
 //                }
 //            }
 
 //            currentMean /= counter;
 //            if (add)
 //            {
 //                double mag = currentMean.magnitude();
 //                if (mag>0)
 //                {
 //                    if (mag>max)
 //                        max = mag;
 
 //                    outDirs.push_back(currentMean);
 //                }
 //            }
 //        }
 //    }
 
 //    if (m_NormalizeVectors)
 //        for (std::size_t i=0; i<outDirs.size(); i++)
 //            outDirs[i].normalize();
 //    else if (max>0)
 //        for (std::size_t i=0; i<outDirs.size(); i++)
 //            outDirs[i] /= max;
 
 //    if (inDirs.size()==outDirs.size())
 //        return outDirs;
 //    else
 //        return FastClustering(outDirs);
 //}
 
 
 //template< class PixelType >
 //TractsToVectorImageFilter< PixelType >::DirectionContainerType::Pointer TractsToVectorImageFilter< PixelType >::MeanShiftClustering(DirectionContainerType::Pointer dirCont)
 //{
 //    DirectionContainerType::Pointer container = DirectionContainerType::New();
 
 //    double max = 0;
 //    for (DirectionContainerType::ConstIterator it = dirCont->Begin(); it!=dirCont->End(); ++it)
 //    {
 //        vnl_vector_fixed<double, 3> mean = ClusterStep(dirCont, it.Value());
 
 //        if (mean.is_zero())
 //            continue;
 //        bool addMean = true;
 
 //        for (DirectionContainerType::ConstIterator it2 = container->Begin(); it2!=container->End(); ++it2)
 //        {
 //            vnl_vector_fixed<double, 3> dir = it2.Value();
 //            double angle = fabs(dot_product(mean, dir)/(mean.magnitude()*dir.magnitude()));
 //            if (angle>=m_Epsilon)
 //            {
 //                addMean = false;
 //                break;
 //            }
 //        }
 
 //        if (addMean)
 //        {
 //            if (m_NormalizeVectors)
 //                mean.normalize();
 //            else if (mean.magnitude()>max)
 //                max = mean.magnitude();
 //            container->InsertElement(container->Size(), mean);
 //        }
 //    }
 
 //    // max normalize voxel directions
 //    if (max>0 && !m_NormalizeVectors)
 //        for (std::size_t i=0; i<container->Size(); i++)
 //            container->ElementAt(i) /= max;
 
 //    if (container->Size()<dirCont->Size())
 //        return MeanShiftClustering(container);
 //    else
 //        return container;
 //}
 
 
 //template< class PixelType >
 //vnl_vector_fixed<double, 3> TractsToVectorImageFilter< PixelType >::ClusterStep(DirectionContainerType::Pointer dirCont, vnl_vector_fixed<double, 3> currentMean)
 //{
 //    vnl_vector_fixed<double, 3> newMean; newMean.fill(0);
 
 //    for (DirectionContainerType::ConstIterator it = dirCont->Begin(); it!=dirCont->End(); ++it)
 //    {
 //        vnl_vector_fixed<double, 3> dir = it.Value();
 //        double angle = dot_product(currentMean, dir)/(currentMean.magnitude()*dir.magnitude());
 //        if (angle>=m_AngularThreshold)
 //            newMean += dir;
 //        else if (-angle>=m_AngularThreshold)
 //            newMean -= dir;
 //    }
 
 //    if (fabs(dot_product(currentMean, newMean)/(currentMean.magnitude()*newMean.magnitude()))>=m_Epsilon || newMean.is_zero())
 //        return newMean;
 //    else
 //        return ClusterStep(dirCont, newMean);
 //}
 }
 
 
 
diff --git a/Modules/DiffusionImaging/MiniApps/DFTracking.cpp b/Modules/DiffusionImaging/MiniApps/DFTracking.cpp
index 93ab11c30f..cedda2c683 100755
--- a/Modules/DiffusionImaging/MiniApps/DFTracking.cpp
+++ b/Modules/DiffusionImaging/MiniApps/DFTracking.cpp
@@ -1,199 +1,193 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 
 #include <mitkBaseData.h>
 #include <mitkImageCast.h>
 #include <mitkImageToItk.h>
 #include <metaCommand.h>
 #include "mitkCommandLineParser.h"
 #include <usAny.h>
 #include <itkImageFileWriter.h>
 #include <mitkIOUtil.h>
 //#include <boost/lexical_cast.hpp>
 #include <iostream>
 #include <fstream>
 #include <itksys/SystemTools.hxx>
 #include <mitkCoreObjectFactory.h>
 
 #include <mitkFiberBundle.h>
 #include <itkDwiNormilzationFilter.h>
 //#include <itkAnalyticalDiffusionQballReconstructionImageFilter.h>
 #include <itkMLBSTrackingFilter.h>
 #include <vtkCell.h>
 #include <itkOrientationDistributionFunction.h>
 #include <itkMersenneTwisterRandomVariateGenerator.h>
 
 #define _USE_MATH_DEFINES
 #include <math.h>
 
 const int numOdfSamples = 200;
 typedef itk::Image< itk::Vector< float, numOdfSamples > , 3 > SampledShImageType;
 
 int main(int argc, char* argv[])
 {
     mitkCommandLineParser parser;
 
     parser.setTitle("Machine Learning Based Streamline Tractography");
     parser.setCategory("Fiber Tracking and Processing Methods");
     parser.setDescription("");
     parser.setContributor("MBI");
 
     parser.setArgumentPrefix("--", "-");
     parser.addArgument("image", "i", mitkCommandLineParser::String, "DWIs:", "input diffusion-weighted image", us::Any(), false);
     parser.addArgument("forest", "f", mitkCommandLineParser::String, "Forest:", "input forest", us::Any(), false);
     parser.addArgument("out", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output fiberbundle", us::Any(), false);
 
     parser.addArgument("stop", "st", mitkCommandLineParser::String, "Stop image:", "stop image", us::Any());
     parser.addArgument("mask", "m", mitkCommandLineParser::String, "Mask image:", "mask image", us::Any());
     parser.addArgument("seed", "s", mitkCommandLineParser::String, "Seed image:", "seed image", us::Any());
 
     parser.addArgument("athres", "a", mitkCommandLineParser::Float, "Angular threshold:", "angular threshold (in radians)", us::Any());
     parser.addArgument("stepsize", "se", mitkCommandLineParser::Float, "Stepsize:", "stepsize", us::Any());
     parser.addArgument("samples", "ns", mitkCommandLineParser::Int, "Samples:", "samples", us::Any());
     parser.addArgument("samplingdist", "sd", mitkCommandLineParser::Float, "Sampling distance:", "sampling distance (in voxels)", us::Any());
     parser.addArgument("seeds", "nse", mitkCommandLineParser::Int, "Seeds per voxel:", "seeds per voxel", us::Any());
 
-    parser.addArgument("usedirection", "ud", mitkCommandLineParser::Bool, "Use previous direction:", "use previous direction as feature", us::Any());
     parser.addArgument("verbose", "v", mitkCommandLineParser::Bool, "Verbose:", "output additional images", us::Any());
 
     map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
     if (parsedArgs.size()==0)
         return EXIT_FAILURE;
 
     string imageFile = us::any_cast<string>(parsedArgs["image"]);
     string forestFile = us::any_cast<string>(parsedArgs["forest"]);
     string outFile = us::any_cast<string>(parsedArgs["out"]);
 
     string maskFile = "";
     if (parsedArgs.count("mask"))
         maskFile = us::any_cast<string>(parsedArgs["mask"]);
 
     string seedFile = "";
     if (parsedArgs.count("seed"))
         seedFile = us::any_cast<string>(parsedArgs["seed"]);
 
     string stopFile = "";
     if (parsedArgs.count("stop"))
         stopFile = us::any_cast<string>(parsedArgs["stop"]);
 
     float stepsize = -1;
     if (parsedArgs.count("stepsize"))
         stepsize = us::any_cast<float>(parsedArgs["stepsize"]);
 
     float athres = 0.7;
     if (parsedArgs.count("athres"))
         athres = us::any_cast<float>(parsedArgs["athres"]);
 
     float samplingdist = 0.25;
     if (parsedArgs.count("samplingdist"))
         samplingdist = us::any_cast<float>(parsedArgs["samplingdist"]);
 
-    bool useDirection = false;
-    if (parsedArgs.count("usedirection"))
-        useDirection = true;
-
     bool verbose = false;
     if (parsedArgs.count("verbose"))
         verbose = true;
 
     int samples = 10;
     if (parsedArgs.count("samples"))
         samples = us::any_cast<int>(parsedArgs["samples"]);
 
     int seeds = 1;
     if (parsedArgs.count("seeds"))
         seeds = us::any_cast<int>(parsedArgs["seeds"]);
 
     typedef itk::Image<unsigned char, 3> ItkUcharImgType;
 
     MITK_INFO << "loading diffusion-weighted image";
     mitk::Image::Pointer dwi = dynamic_cast<mitk::Image*>(mitk::IOUtil::LoadImage(imageFile).GetPointer());
 
     ItkUcharImgType::Pointer mask;
     if (!maskFile.empty())
     {
         MITK_INFO << "loading mask image";
         mitk::Image::Pointer img = dynamic_cast<mitk::Image*>(mitk::IOUtil::LoadImage(maskFile).GetPointer());
         mask = ItkUcharImgType::New();
         mitk::CastToItkImage(img, mask);
     }
 
     ItkUcharImgType::Pointer seed;
     if (!seedFile.empty())
     {
         MITK_INFO << "loading seed image";
         mitk::Image::Pointer img = dynamic_cast<mitk::Image*>(mitk::IOUtil::LoadImage(seedFile).GetPointer());
         seed = ItkUcharImgType::New();
         mitk::CastToItkImage(img, seed);
     }
 
     ItkUcharImgType::Pointer stop;
     if (!stopFile.empty())
     {
         MITK_INFO << "loading stop image";
         mitk::Image::Pointer img = dynamic_cast<mitk::Image*>(mitk::IOUtil::LoadImage(stopFile).GetPointer());
         stop = ItkUcharImgType::New();
         mitk::CastToItkImage(img, stop);
     }
 
     MITK_INFO << "loading forest";
     vigra::RandomForest<int> rf;
     vigra::rf_import_HDF5(rf, forestFile);
 
     typedef itk::MLBSTrackingFilter<100> TrackerType;
     TrackerType::Pointer tracker = TrackerType::New();
     tracker->SetInput(0, mitk::DiffusionPropertyHelper::GetItkVectorImage(dwi));
     tracker->SetGradientDirections( mitk::DiffusionPropertyHelper::GetGradientContainer(dwi) );
     tracker->SetB_Value( mitk::DiffusionPropertyHelper::GetReferenceBValue(dwi) );
     tracker->SetMaskImage(mask);
     tracker->SetSeedImage(seed);
     tracker->SetStoppingRegions(stop);
     tracker->SetSeedsPerVoxel(seeds);
-    tracker->SetUseDirection(useDirection);
     tracker->SetStepSize(stepsize);
     tracker->SetAngularThreshold(athres);
     tracker->SetDecisionForest(&rf);
     tracker->SetSamplingDistance(samplingdist);
     tracker->SetNumberOfSamples(samples);
     //tracker->SetAvoidStop(false);
     tracker->SetAposterioriCurvCheck(false);
     tracker->SetRemoveWmEndFibers(false);
     tracker->Update();
     vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData();
     mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly);
 
     mitk::IOUtil::SaveBaseData(outFib, outFile);
 
     if (verbose)
     {
         MITK_INFO << "Writing images...";
         string outName = itksys::SystemTools::GetFilenamePath(outFile)+"/"+itksys::SystemTools::GetFilenameWithoutLastExtension(outFile);
         itk::ImageFileWriter< TrackerType::ItkDoubleImgType >::Pointer writer = itk::ImageFileWriter< TrackerType::ItkDoubleImgType >::New();
         writer->SetFileName(outName+"_WhiteMatter.nrrd");
         writer->SetInput(tracker->GetWmImage());
         writer->Update();
 
         writer->SetFileName(outName+"_NotWhiteMatter.nrrd");
         writer->SetInput(tracker->GetNotWmImage());
         writer->Update();
 
         writer->SetFileName(outName+"_AvoidStop.nrrd");
         writer->SetInput(tracker->GetAvoidStopImage());
         writer->Update();
     }
 
     return EXIT_SUCCESS;
 }
diff --git a/Modules/DiffusionImaging/MiniApps/DFTraining.cpp b/Modules/DiffusionImaging/MiniApps/DFTraining.cpp
index d4c2bc3d40..03d9f9107e 100755
--- a/Modules/DiffusionImaging/MiniApps/DFTraining.cpp
+++ b/Modules/DiffusionImaging/MiniApps/DFTraining.cpp
@@ -1,529 +1,483 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 
 #include <mitkBaseData.h>
 #include <mitkImageCast.h>
 #include <mitkImageToItk.h>
 #include <metaCommand.h>
 #include "mitkCommandLineParser.h"
 #include <usAny.h>
 #include <itkImageFileWriter.h>
 #include <mitkIOUtil.h>
 #include <iostream>
 #include <fstream>
 #include <itksys/SystemTools.hxx>
 #include <mitkCoreObjectFactory.h>
 
 #include <mitkFiberBundle.h>
 #include <itkDwiNormilzationFilter.h>
 #include <itkAnalyticalDiffusionQballReconstructionImageFilter.h>
 #include <vtkCell.h>
 #include <itkOrientationDistributionFunction.h>
 #include <itkMersenneTwisterRandomVariateGenerator.h>
 
 #include <vigra/random_forest.hxx>
 #include <vigra/multi_array.hxx>
 #include <vigra/random_forest_hdf5_impex.hxx>
 
 #define _USE_MATH_DEFINES
 #include <math.h>
 
 const int numOdfSamples = 200;  // ODF is sampled in 200 directions but actuyll only 100 are used (symmetric)
 typedef itk::Image< itk::Vector< float, numOdfSamples > , 3 > SampledShImageType;
 
 void TrainForest( vigra::RandomForest<int> &rf, vigra::MultiArray<2, double> &labelData, vigra::MultiArray<2, double> &featureData, int numTrees, int max_tree_depth, double sample_fraction )
 {
     MITK_INFO << "Maximum tree depths: " << max_tree_depth;
     MITK_INFO << "Sample fraction per tree: " << sample_fraction;
     MITK_INFO << "Number of trees: " << numTrees;
     vigra::rf::visitors::OOB_Error oob_v;
 
 //    rf.set_options().use_stratification(vigra::RF_NONE); // How the data should be made equal
 //    rf.set_options().sample_with_replacement(true); // if sampled with replacement or not
 //    rf.set_options().samples_per_tree(sample_fraction); // Fraction of samples that are used to train a tree
 //    rf.set_options().tree_count(1); // Number of trees that are calculated;
 //    rf.set_options().min_split_node_size(5); // Minimum number of datapoints that must be in a node
 //    rf.ext_param_.max_tree_depth = max_tree_depth;
 //    //    rf.set_options().features_per_node(10);
 //    rf.learn(featureData, labelData, vigra::rf::visitors::create_visitor(oob_v));
 
     std::vector< vigra::RandomForest<int> > trees;
     int count = 0;
 #pragma omp parallel for
     for (int i = 0; i < numTrees; ++i)
     {
         vigra::RandomForest<int> lrf;
         vigra::rf::visitors::OOB_Error loob_v;
 
         lrf.set_options().use_stratification(vigra::RF_NONE); // How the data should be made equal
         lrf.set_options().sample_with_replacement(true); // if sampled with replacement or not
         lrf.set_options().samples_per_tree(sample_fraction); // Fraction of samples that are used to train a tree
         lrf.set_options().tree_count(1); // Number of trees that are calculated;
         lrf.set_options().min_split_node_size(5); // Minimum number of datapoints that must be in a node
         lrf.ext_param_.max_tree_depth = max_tree_depth;
         //        lrf.set_options().features_per_node(10);
 
         lrf.learn(featureData, labelData);//, vigra::rf::visitors::create_visitor(loob_v));
 #pragma omp critical
         {
             count++;
             MITK_INFO << "Tree " << count << " finished training.";
             trees.push_back(lrf);
             //rf.trees_.push_back(lrf.trees_[0]);
         }
     }
 
     for (int i = 1; i < numTrees; ++i)
         trees.at(0).trees_.push_back(trees.at(i).trees_[0]);
 
     rf = trees.at(0);
     rf.options_.tree_count_ = numTrees;
     MITK_INFO << "Training finsihed";
     //MITK_INFO << "The out-of-bag error is: " << oob_v.oob_breiman << std::endl;
 }
 
 SampledShImageType::PixelType GetImageValues(itk::Point<float, 3> itkP, SampledShImageType::Pointer image)
 {
     itk::Index<3> idx;
     itk::ContinuousIndex< double, 3> cIdx;
     image->TransformPhysicalPointToIndex(itkP, idx);
     image->TransformPhysicalPointToContinuousIndex(itkP, cIdx);
 
     SampledShImageType::PixelType pix; pix.Fill(0.0);
     if ( image->GetLargestPossibleRegion().IsInside(idx) )
         pix = image->GetPixel(idx);
     else
         return pix;
 
     double frac_x = cIdx[0] - idx[0];
     double frac_y = cIdx[1] - idx[1];
     double frac_z = cIdx[2] - idx[2];
     if (frac_x<0)
     {
         idx[0] -= 1;
         frac_x += 1;
     }
     if (frac_y<0)
     {
         idx[1] -= 1;
         frac_y += 1;
     }
     if (frac_z<0)
     {
         idx[2] -= 1;
         frac_z += 1;
     }
     frac_x = 1-frac_x;
     frac_y = 1-frac_y;
     frac_z = 1-frac_z;
 
     // int coordinates inside image?
     if (idx[0] >= 0 && idx[0] < image->GetLargestPossibleRegion().GetSize(0)-1 &&
             idx[1] >= 0 && idx[1] < image->GetLargestPossibleRegion().GetSize(1)-1 &&
             idx[2] >= 0 && idx[2] < image->GetLargestPossibleRegion().GetSize(2)-1)
     {
         vnl_vector_fixed<double, 8> interpWeights;
         interpWeights[0] = (  frac_x)*(  frac_y)*(  frac_z);
         interpWeights[1] = (1-frac_x)*(  frac_y)*(  frac_z);
         interpWeights[2] = (  frac_x)*(1-frac_y)*(  frac_z);
         interpWeights[3] = (  frac_x)*(  frac_y)*(1-frac_z);
         interpWeights[4] = (1-frac_x)*(1-frac_y)*(  frac_z);
         interpWeights[5] = (  frac_x)*(1-frac_y)*(1-frac_z);
         interpWeights[6] = (1-frac_x)*(  frac_y)*(1-frac_z);
         interpWeights[7] = (1-frac_x)*(1-frac_y)*(1-frac_z);
 
 
         pix = image->GetPixel(idx) * interpWeights[0];
         SampledShImageType::IndexType tmpIdx = idx; tmpIdx[0]++;
         pix +=  image->GetPixel(tmpIdx) * interpWeights[1];
         tmpIdx = idx; tmpIdx[1]++;
         pix +=  image->GetPixel(tmpIdx) * interpWeights[2];
         tmpIdx = idx; tmpIdx[2]++;
         pix +=  image->GetPixel(tmpIdx) * interpWeights[3];
         tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++;
         pix +=  image->GetPixel(tmpIdx) * interpWeights[4];
         tmpIdx = idx; tmpIdx[1]++; tmpIdx[2]++;
         pix +=  image->GetPixel(tmpIdx) * interpWeights[5];
         tmpIdx = idx; tmpIdx[2]++; tmpIdx[0]++;
         pix +=  image->GetPixel(tmpIdx) * interpWeights[6];
         tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++; tmpIdx[2]++;
         pix +=  image->GetPixel(tmpIdx) * interpWeights[7];
     }
 
     return pix;
 }
 
 int main(int argc, char* argv[])
 {
     MITK_INFO << "DFTraining";
     mitkCommandLineParser parser;
 
     parser.setTitle("Machine Learning Based Streamline Tractography");
     parser.setCategory("Fiber Tracking and Processing Methods");
     parser.setDescription("");
     parser.setContributor("MBI");
 
     parser.setArgumentPrefix("--", "-");
     parser.addArgument("images", "i", mitkCommandLineParser::StringList, "DWIs:", "input diffusion-weighted images", us::Any(), false);
     parser.addArgument("wmmasks", "w", mitkCommandLineParser::StringList, "WM-Masks:", "white matter mask images", us::Any(), false);
     parser.addArgument("tractograms", "t", mitkCommandLineParser::StringList, "Tractograms:", "input tractograms (.fib, vtk ascii file format)", us::Any(), false);
     parser.addArgument("masks", "m", mitkCommandLineParser::StringList, "Masks:", "mask images", us::Any());
     parser.addArgument("forest", "f", mitkCommandLineParser::OutputFile, "Forest:", "output forest", us::Any(), false);
 
     parser.addArgument("stepsize", "s", mitkCommandLineParser::Float, "Stepsize:", "stepsize", us::Any());
     parser.addArgument("gmsamples", "g", mitkCommandLineParser::Int, "Number of gray matter samples per voxel:", "Number of gray matter samples per voxel", us::Any());
     parser.addArgument("numtrees", "n", mitkCommandLineParser::Int, "Number of trees:", "number of trees", us::Any());
     parser.addArgument("max_tree_depth", "d", mitkCommandLineParser::Int, "Max. tree depth:", "maximum tree depth", us::Any());
     parser.addArgument("sample_fraction", "sf", mitkCommandLineParser::Float, "Sample fraction:", "fraction of samples used per tree", us::Any());
-    parser.addArgument("usedirection", "ud", mitkCommandLineParser::Bool, "bla:", "bla", us::Any());
 
     map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
     if (parsedArgs.size()==0)
         return EXIT_FAILURE;
 
     mitkCommandLineParser::StringContainerType imageFiles = us::any_cast<mitkCommandLineParser::StringContainerType>(parsedArgs["images"]);
     mitkCommandLineParser::StringContainerType wmMaskFiles = us::any_cast<mitkCommandLineParser::StringContainerType>(parsedArgs["wmmasks"]);
 
     mitkCommandLineParser::StringContainerType maskFiles;
     if (parsedArgs.count("masks"))
         maskFiles = us::any_cast<mitkCommandLineParser::StringContainerType>(parsedArgs["masks"]);
 
     string forestFile = us::any_cast<string>(parsedArgs["forest"]);
 
     mitkCommandLineParser::StringContainerType tractogramFiles;
     if (parsedArgs.count("tractograms"))
         tractogramFiles = us::any_cast<mitkCommandLineParser::StringContainerType>(parsedArgs["tractograms"]);
 
     int numTrees = 30;
     if (parsedArgs.count("numtrees"))
         numTrees = us::any_cast<int>(parsedArgs["numtrees"]);
 
     int gmsamples = 50;
     if (parsedArgs.count("gmsamples"))
         gmsamples = us::any_cast<int>(parsedArgs["gmsamples"]);
 
     float stepsize = -1;
     if (parsedArgs.count("stepsize"))
         stepsize = us::any_cast<float>(parsedArgs["stepsize"]);
 
     int max_tree_depth = 50;
     if (parsedArgs.count("max_tree_depth"))
         max_tree_depth = us::any_cast<int>(parsedArgs["max_tree_depth"]);
 
     double sample_fraction = 1.0;
     if (parsedArgs.count("sample_fraction"))
         sample_fraction = us::any_cast<float>(parsedArgs["sample_fraction"]);
 
     // load DWI images
     if (imageFiles.size()<tractogramFiles.size())
     {
         MITK_INFO << "Missing training images!";
         return EXIT_FAILURE;
     }
-    bool useDirection = false;
-    if (parsedArgs.count("usedirection"))
-        useDirection = true;
 
     const int order = 6;
     typedef itk::AnalyticalDiffusionQballReconstructionImageFilter<short,short,float,order,numOdfSamples> QballFilterType;
 
     MITK_INFO << "loading diffusion-weighted images and reconstructing feature images";
     std::vector< SampledShImageType::Pointer > sampledShImages;
     for (unsigned int i=0; i<imageFiles.size(); i++)
     {
         mitk::Image::Pointer dwi = dynamic_cast<mitk::Image*>(mitk::IOUtil::LoadImage(imageFiles.at(i)).GetPointer());
 
         QballFilterType::Pointer qballfilter = QballFilterType::New();
         qballfilter->SetGradientImage( mitk::DiffusionPropertyHelper::GetGradientContainer(dwi), mitk::DiffusionPropertyHelper::GetItkVectorImage(dwi) );
         qballfilter->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(dwi));
         qballfilter->SetLambda(0.006);
         qballfilter->SetNormalizationMethod(QballFilterType::QBAR_RAW_SIGNAL);
         qballfilter->Update();
         //            FeatureImageType::Pointer itkFeatureImage = qballfilter->GetCoefficientImage();
         //            featureImageVector.push_back(itkFeatureImage);
         sampledShImages.push_back(qballfilter->GetOutput());
     }
 
 
     typedef itk::Image<unsigned char, 3> ItkUcharImgType;
     std::vector< ItkUcharImgType::Pointer > maskImageVector;
     std::vector< ItkUcharImgType::Pointer > wmMaskImageVector;
 
     MITK_INFO << "loading tractograms";
     int numSamples = 0;
     std::vector< mitk::FiberBundle::Pointer > tractograms;
     for (unsigned int t=0; t<tractogramFiles.size(); t++)
     {
         ItkUcharImgType::Pointer mask;
         if (t<maskFiles.size())
         {
             mitk::Image::Pointer img = dynamic_cast<mitk::Image*>(mitk::IOUtil::LoadImage(maskFiles.at(t)).GetPointer());
             mask = ItkUcharImgType::New();
             mitk::CastToItkImage(img, mask);
             maskImageVector.push_back(mask);
         }
         mitk::Image::Pointer img2 = dynamic_cast<mitk::Image*>(mitk::IOUtil::LoadImage(wmMaskFiles.at(t)).GetPointer());
         ItkUcharImgType::Pointer wmmask = ItkUcharImgType::New();
         mitk::CastToItkImage(img2, wmmask);
         wmMaskImageVector.push_back(wmmask);
 
         itk::ImageRegionConstIterator<ItkUcharImgType> it(wmmask, wmmask->GetLargestPossibleRegion());
         int OUTOFWM = 0;    // count voxels outside of the white matter mask
         while(!it.IsAtEnd())
         {
             if (it.Get()==0)
                 if (mask.IsNull() || (mask.IsNotNull() && mask->GetPixel(it.GetIndex())>0))
                     OUTOFWM++;
             ++it;
         }
         numSamples += gmsamples*OUTOFWM;    // for each of the non-white matter voxels we add a certain number of sampling points. these sampling points are used to tell the classifier where to recognize non-WM tissue
 
         MITK_INFO << "Samples outside of WM: " << numSamples << " (" << gmsamples << " per non-WM voxel)";
 
         // load and resample training tractograms
         mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(mitk::IOUtil::Load(tractogramFiles.at(t)).at(0).GetPointer());
         if (stepsize<0)
         {
             SampledShImageType::Pointer image = sampledShImages.at(t);
             float minSpacing = 1;
             if(image->GetSpacing()[0]<image->GetSpacing()[1] && image->GetSpacing()[0]<image->GetSpacing()[2])
                 minSpacing = image->GetSpacing()[0];
             else if (image->GetSpacing()[1] < image->GetSpacing()[2])
                 minSpacing = image->GetSpacing()[1];
             else
                 minSpacing = image->GetSpacing()[2];
             stepsize = minSpacing*0.5;
         }
         fib->ResampleSpline(stepsize);
         tractograms.push_back(fib);
         numSamples += fib->GetNumberOfPoints(); // each point of the fiber gives us a training direction
         numSamples -= 2*fib->GetNumFibers();    // we don't use the first and last point because there we do not have a previous direction, which is needed as feature
     }
     MITK_INFO << "Number of samples: " << numSamples;
 
     // get ODF directions and number of features
     vnl_vector_fixed<double,3> ref; ref.fill(0); ref[0]=1;
     itk::OrientationDistributionFunction< double, numOdfSamples > odf;
     std::vector< int > directionIndices;
     for (unsigned int f=0; f<numOdfSamples; f++)
     {
         if (dot_product(ref, odf.GetDirection(f))>0)    // we only use directions on one hemisphere (symmetric)
             directionIndices.push_back(f);  // remember indices that are on the desired hemisphere
     }
     const int numSignalFeatures = numOdfSamples/2;
-    int numDirectionFeatures = 0;
-    if (useDirection)
-        numDirectionFeatures = 3;
+    int numDirectionFeatures = 3;
 
     vigra::MultiArray<2, double> featureData( vigra::Shape2(numSamples,numSignalFeatures+numDirectionFeatures) );
     MITK_INFO << "Number of features: " << featureData.shape(1);
     vigra::MultiArray<2, double> labelData( vigra::Shape2(numSamples,1) );
 
     itk::Statistics::MersenneTwisterRandomVariateGenerator::Pointer m_RandGen = itk::Statistics::MersenneTwisterRandomVariateGenerator::New();
     m_RandGen->SetSeed();
     MITK_INFO <<  "Creating training data from tractograms and feature images";
     int sampleCounter = 0;
     for (unsigned int t=0; t<tractograms.size(); t++)
     {
         SampledShImageType::Pointer image = sampledShImages.at(t);  // raw data
         ItkUcharImgType::Pointer wmMask = wmMaskImageVector.at(t);  // white matter mask
         ItkUcharImgType::Pointer mask;                              // mask image. only train inside mask (if set)
         if (t<maskImageVector.size())
             mask = maskImageVector.at(t);
 
         // get training data for non-WM tissue
         itk::ImageRegionConstIterator<ItkUcharImgType> it(wmMask, wmMask->GetLargestPossibleRegion());
         while(!it.IsAtEnd())
         {
             if (it.Get()==0 && (mask.IsNull() || (mask.IsNotNull() && mask->GetPixel(it.GetIndex())>0)))
             {
                 SampledShImageType::PixelType pix = image->GetPixel(it.GetIndex());
-                if (useDirection)
-                {
-                    // null direction
-                    for (unsigned int f=0; f<numSignalFeatures; f++)
-                    {
-                        featureData(sampleCounter,f) = pix[directionIndices[f]];
-                        if(featureData(sampleCounter,f)!=featureData(sampleCounter,f))
-                            featureData(sampleCounter,f) = 0;
-                    }
-                    featureData(sampleCounter,numSignalFeatures) = 0;
-                    featureData(sampleCounter,numSignalFeatures+1) = 0;
-                    featureData(sampleCounter,numSignalFeatures+2) = 0;
-                    labelData(sampleCounter,0) = directionIndices.size();
-                    sampleCounter++;
 
-                    // we use random "previous" directions for the non-WM samples
-                    for (int i=1; i<gmsamples; i++)
-                    {
-                        for (unsigned int f=0; f<numSignalFeatures; f++)
-                        {
-                            featureData(sampleCounter,f) = pix[directionIndices[f]];
-                            if(featureData(sampleCounter,f)!=featureData(sampleCounter,f))
-                                featureData(sampleCounter,f) = 0;
-                        }
-                        int c=0;
-                        vnl_vector_fixed<double,3> probe;
-                        probe[0] = m_RandGen->GetVariate()*2-1;
-                        probe[1] = m_RandGen->GetVariate()*2-1;
-                        probe[2] = m_RandGen->GetVariate()*2-1;
-                        probe.normalize();
-                        if (dot_product(ref, probe)<0)
-                            probe *= -1;
-                        for (unsigned int f=numSignalFeatures; f<numSignalFeatures+3; f++)
-                        {
-                            featureData(sampleCounter,f) = probe[c];
-                            c++;
-                        }
-                        labelData(sampleCounter,0) = directionIndices.size();
-                        sampleCounter++;
-                    }
+                // null direction
+                for (unsigned int f=0; f<numSignalFeatures; f++)
+                {
+                    featureData(sampleCounter,f) = pix[directionIndices[f]];
+                    if(featureData(sampleCounter,f)!=featureData(sampleCounter,f))
+                        featureData(sampleCounter,f) = 0;
                 }
-                else
+                featureData(sampleCounter,numSignalFeatures) = 0;
+                featureData(sampleCounter,numSignalFeatures+1) = 0;
+                featureData(sampleCounter,numSignalFeatures+2) = 0;
+                labelData(sampleCounter,0) = directionIndices.size();
+                sampleCounter++;
+
+                // we use random "previous" directions for the non-WM samples
+                for (int i=1; i<gmsamples; i++)
                 {
-                    itk::ContinuousIndex<double, 3> idx;
-                    idx[0] = it.GetIndex()[0];
-                    idx[1] = it.GetIndex()[1];
-                    idx[2] = it.GetIndex()[2];
-                    itk::Point<float, 3> itkP1;
-                    image->TransformContinuousIndexToPhysicalPoint(idx, itkP1);
-                    SampledShImageType::PixelType pix = GetImageValues(itkP1, image);;
                     for (unsigned int f=0; f<numSignalFeatures; f++)
                     {
                         featureData(sampleCounter,f) = pix[directionIndices[f]];
                         if(featureData(sampleCounter,f)!=featureData(sampleCounter,f))
                             featureData(sampleCounter,f) = 0;
                     }
-                    labelData(sampleCounter,0) = directionIndices.size();
-                    sampleCounter++;
-
-                    for (int i=1; i<gmsamples; i++)
+                    int c=0;
+                    vnl_vector_fixed<double,3> probe;
+                    probe[0] = m_RandGen->GetVariate()*2-1;
+                    probe[1] = m_RandGen->GetVariate()*2-1;
+                    probe[2] = m_RandGen->GetVariate()*2-1;
+                    probe.normalize();
+                    if (dot_product(ref, probe)<0)
+                        probe *= -1;
+                    for (unsigned int f=numSignalFeatures; f<numSignalFeatures+3; f++)
                     {
-                        itk::ContinuousIndex<double, 3> idx;
-                        idx[0] = it.GetIndex()[0] + m_RandGen->GetVariate()-0.5;
-                        idx[1] = it.GetIndex()[1] + m_RandGen->GetVariate()-0.5;
-                        idx[2] = it.GetIndex()[2] + m_RandGen->GetVariate()-0.5;
-                        itk::Point<float, 3> itkP1;
-                        image->TransformContinuousIndexToPhysicalPoint(idx, itkP1);
-                        SampledShImageType::PixelType pix = GetImageValues(itkP1, image);;
-                        for (unsigned int f=0; f<numSignalFeatures; f++)
-                        {
-                            featureData(sampleCounter,f) = pix[directionIndices[f]];
-                            if(featureData(sampleCounter,f)!=featureData(sampleCounter,f))
-                                featureData(sampleCounter,f) = 0;
-                        }
-                        labelData(sampleCounter,0) = directionIndices.size();
-                        sampleCounter++;
+                        featureData(sampleCounter,f) = probe[c];
+                        c++;
                     }
+                    labelData(sampleCounter,0) = directionIndices.size();
+                    sampleCounter++;
                 }
             }
             ++it;
         }
 
         mitk::FiberBundle::Pointer fib = tractograms.at(t);
         vtkSmartPointer< vtkPolyData > polyData = fib->GetFiberPolyData();
         for (int i=0; i<fib->GetNumFibers(); i++)
         {
             vtkCell* cell = polyData->GetCell(i);
             int numPoints = cell->GetNumberOfPoints();
             vtkPoints* points = cell->GetPoints();
 
             vnl_vector_fixed<double,3> dirOld; dirOld.fill(0.0);
 
             for (int j=0; j<numPoints-1; j++)
             {
                 // calculate direction
                 double* p1 = points->GetPoint(j);
                 itk::Point<float, 3> itkP1;
                 itkP1[0] = p1[0]; itkP1[1] = p1[1]; itkP1[2] = p1[2];
 
                 vnl_vector_fixed<double,3> dir; dir.fill(0.0);
 
                 itk::Point<float, 3> itkP2;
                 double* p2 = points->GetPoint(j+1);
                 itkP2[0] = p2[0]; itkP2[1] = p2[1]; itkP2[2] = p2[2];
                 dir[0]=itkP2[0]-itkP1[0];
                 dir[1]=itkP2[1]-itkP1[1];
                 dir[2]=itkP2[2]-itkP1[2];
 
                 if (dir.magnitude()<0.0001)
                 {
                     MITK_INFO << "streamline error!";
                     continue;
                 }
                 dir.normalize();
                 if (dir[0]!=dir[0] || dir[1]!=dir[1] || dir[2]!=dir[2])
                 {
                     MITK_INFO << "ERROR: NaN direction!";
                     continue;
                 }
 
                 if (j==0)
                 {
                     dirOld = dir;
                     continue;
                 }
 
                 // get voxel values
                 SampledShImageType::PixelType pix = GetImageValues(itkP1, image);
                 for (unsigned int f=0; f<numSignalFeatures; f++)
                     featureData(sampleCounter,f) = pix[directionIndices[f]];
 
                 // direction training features
                 int c = 0;
                 if (dot_product(ref, dirOld)<0)
                     dirOld *= -1;
 
-                if (useDirection)
-                    for (unsigned int f=numSignalFeatures; f<numSignalFeatures+3; f++)
-                    {
-                        featureData(sampleCounter,f) = dirOld[c];
-                        c++;
-                    }
+                for (unsigned int f=numSignalFeatures; f<numSignalFeatures+3; f++)
+                {
+                    featureData(sampleCounter,f) = dirOld[c];
+                    c++;
+                }
 
                 // set label values
                 double angle = 0;
                 double m = dir.magnitude();
                 if (m>0.0001)
                 {
                     int label = 0;
                     for (unsigned int f=0; f<numSignalFeatures; f++)
                     {
                         double a = fabs(dot_product(dir, odf.GetDirection(directionIndices[f])));
                         if (a>angle)
                         {
                             labelData(sampleCounter,0) = f;
                             angle = a;
                             label = f;
                         }
                     }
                 }
 
                 dirOld = dir;
                 sampleCounter++;
             }
         }
     }
 
     MITK_INFO << "Training forest";
     vigra::RandomForest<int> rf;
     TrainForest( rf, labelData, featureData, numTrees, max_tree_depth, sample_fraction );
     MITK_INFO << "Writing forest";
     vigra::rf_export_HDF5( rf, forestFile, "" );
     MITK_INFO << "Finished training";
 
     return EXIT_SUCCESS;
 }
diff --git a/Modules/DiffusionImaging/MiniApps/LocalDirectionalFiberPlausibility.cpp b/Modules/DiffusionImaging/MiniApps/LocalDirectionalFiberPlausibility.cpp
index 326838cef7..2154df23e4 100755
--- a/Modules/DiffusionImaging/MiniApps/LocalDirectionalFiberPlausibility.cpp
+++ b/Modules/DiffusionImaging/MiniApps/LocalDirectionalFiberPlausibility.cpp
@@ -1,300 +1,312 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 
 #include <mitkBaseData.h>
 #include <mitkImageCast.h>
 #include <mitkImageToItk.h>
 #include <itkEvaluateDirectionImagesFilter.h>
 #include <metaCommand.h>
 #include "mitkCommandLineParser.h"
 #include <itkTractsToVectorImageFilter.h>
 #include <usAny.h>
 #include <itkImageFileWriter.h>
 #include <mitkIOUtil.h>
 #include <boost/lexical_cast.hpp>
 #include <iostream>
 #include <fstream>
 #include <itksys/SystemTools.hxx>
 #include <mitkCoreObjectFactory.h>
 
 #define _USE_MATH_DEFINES
 #include <math.h>
 
 int main(int argc, char* argv[])
 {
     mitkCommandLineParser parser;
 
     parser.setTitle("Local Directional Fiber Plausibility");
     parser.setCategory("Fiber Tracking and Processing Methods");
     parser.setDescription(" ");
     parser.setContributor("MBI");
 
     parser.setArgumentPrefix("--", "-");
     parser.addArgument("input", "i", mitkCommandLineParser::InputFile, "Input:", "input tractogram (.fib, vtk ascii file format)", us::Any(), false);
     parser.addArgument("reference", "r", mitkCommandLineParser::StringList, "Reference images:", "reference direction images", us::Any(), false);
     parser.addArgument("out", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output root", us::Any(), false);
     parser.addArgument("mask", "m", mitkCommandLineParser::StringList, "Masks:", "mask images");
     parser.addArgument("athresh", "a", mitkCommandLineParser::Float, "Angular threshold:", "angular threshold in degrees. closer fiber directions are regarded as one direction and clustered together.", 25, true);
+    parser.addArgument("sthresh", "s", mitkCommandLineParser::Float, "Size threshold:", "Relative peak size threshold per voxel.", 0.0, true);
+    parser.addArgument("maxdirs", "md", mitkCommandLineParser::Int, "Max. Clusters:", "Maximum number of fiber clusters.", 0, true);
     parser.addArgument("verbose", "v", mitkCommandLineParser::Bool, "Verbose:", "output optional and intermediate calculation results");
     parser.addArgument("ignore", "n", mitkCommandLineParser::Bool, "Ignore:", "don't increase error for missing or too many directions");
     parser.addArgument("fileID", "id", mitkCommandLineParser::String, "ID:", "optional ID field");
 
     map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
     if (parsedArgs.size()==0)
         return EXIT_FAILURE;
 
     mitkCommandLineParser::StringContainerType referenceImages = us::any_cast<mitkCommandLineParser::StringContainerType>(parsedArgs["reference"]);
     mitkCommandLineParser::StringContainerType maskImages;
     if (parsedArgs.count("mask"))
         maskImages = us::any_cast<mitkCommandLineParser::StringContainerType>(parsedArgs["mask"]);
 
     string fibFile = us::any_cast<string>(parsedArgs["input"]);
 
     float angularThreshold = 25;
     if (parsedArgs.count("athresh"))
         angularThreshold = us::any_cast<float>(parsedArgs["athresh"]);
 
+    float sizeThreshold = 0;
+    if (parsedArgs.count("sthresh"))
+        sizeThreshold = us::any_cast<float>(parsedArgs["sthresh"]);
+
+    int maxDirs = 0;
+    if (parsedArgs.count("maxdirs"))
+        maxDirs = us::any_cast<int>(parsedArgs["maxdirs"]);
+
     string outRoot = us::any_cast<string>(parsedArgs["out"]);
 
     bool verbose = false;
     if (parsedArgs.count("verbose"))
         verbose = us::any_cast<bool>(parsedArgs["verbose"]);
 
     bool ignore = false;
     if (parsedArgs.count("ignore"))
         ignore = us::any_cast<bool>(parsedArgs["ignore"]);
 
     string fileID = "";
     if (parsedArgs.count("fileID"))
         fileID = us::any_cast<string>(parsedArgs["fileID"]);
 
 
     try
     {
         typedef itk::Image<unsigned char, 3>                                    ItkUcharImgType;
         typedef itk::Image< itk::Vector< float, 3>, 3 >                         ItkDirectionImage3DType;
         typedef itk::VectorContainer< unsigned int, ItkDirectionImage3DType::Pointer >   ItkDirectionImageContainerType;
         typedef itk::EvaluateDirectionImagesFilter< float >                     EvaluationFilterType;
 
         // load fiber bundle
         mitk::FiberBundle::Pointer inputTractogram = dynamic_cast<mitk::FiberBundle*>(mitk::IOUtil::LoadDataNode(fibFile)->GetData());
 
         // load reference directions
         ItkDirectionImageContainerType::Pointer referenceImageContainer = ItkDirectionImageContainerType::New();
         for (unsigned int i=0; i<referenceImages.size(); i++)
         {
             try
             {
                 mitk::Image::Pointer img = dynamic_cast<mitk::Image*>(mitk::IOUtil::LoadDataNode(referenceImages.at(i))->GetData());
                 typedef mitk::ImageToItk< ItkDirectionImage3DType > CasterType;
                 CasterType::Pointer caster = CasterType::New();
                 caster->SetInput(img);
                 caster->Update();
                 ItkDirectionImage3DType::Pointer itkImg = caster->GetOutput();
                 referenceImageContainer->InsertElement(referenceImageContainer->Size(),itkImg);
             }
             catch(...){ std::cout << "could not load: " << referenceImages.at(i); }
         }
 
         ItkUcharImgType::Pointer itkMaskImage = ItkUcharImgType::New();
         ItkDirectionImage3DType::Pointer dirImg = referenceImageContainer->GetElement(0);
         itkMaskImage->SetSpacing( dirImg->GetSpacing() );
         itkMaskImage->SetOrigin( dirImg->GetOrigin() );
         itkMaskImage->SetDirection( dirImg->GetDirection() );
         itkMaskImage->SetLargestPossibleRegion( dirImg->GetLargestPossibleRegion() );
         itkMaskImage->SetBufferedRegion( dirImg->GetLargestPossibleRegion() );
         itkMaskImage->SetRequestedRegion( dirImg->GetLargestPossibleRegion() );
         itkMaskImage->Allocate();
         itkMaskImage->FillBuffer(1);
 
         // extract directions from fiber bundle
         itk::TractsToVectorImageFilter<float>::Pointer fOdfFilter = itk::TractsToVectorImageFilter<float>::New();
         fOdfFilter->SetFiberBundle(inputTractogram);
         fOdfFilter->SetMaskImage(itkMaskImage);
         fOdfFilter->SetAngularThreshold(cos(angularThreshold*M_PI/180));
         fOdfFilter->SetNormalizeVectors(true);
         fOdfFilter->SetUseWorkingCopy(false);
+        fOdfFilter->SetSizeThreshold(sizeThreshold);
+        fOdfFilter->SetMaxNumDirections(maxDirs);
         fOdfFilter->Update();
         ItkDirectionImageContainerType::Pointer directionImageContainer = fOdfFilter->GetDirectionImageContainer();
 
         if (verbose)
         {
             // write vector field
             mitk::FiberBundle::Pointer directions = fOdfFilter->GetOutputFiberBundle();
 
             string outfilename = outRoot;
             outfilename.append("_VECTOR_FIELD.fib");
 
             mitk::IOUtil::SaveBaseData(directions.GetPointer(), outfilename );
 
             // write direction images
             for (unsigned int i=0; i<directionImageContainer->Size(); i++)
             {
                 itk::TractsToVectorImageFilter<float>::ItkDirectionImageType::Pointer itkImg = directionImageContainer->GetElement(i);
                 typedef itk::ImageFileWriter< itk::TractsToVectorImageFilter<float>::ItkDirectionImageType > WriterType;
                 WriterType::Pointer writer = WriterType::New();
 
                 string outfilename = outRoot;
                 outfilename.append("_DIRECTION_");
                 outfilename.append(boost::lexical_cast<string>(i));
                 outfilename.append(".nrrd");
 
                 writer->SetFileName(outfilename.c_str());
                 writer->SetInput(itkImg);
                 writer->Update();
             }
 
             // write num direction image
             {
                 ItkUcharImgType::Pointer numDirImage = fOdfFilter->GetNumDirectionsImage();
                 typedef itk::ImageFileWriter< ItkUcharImgType > WriterType;
                 WriterType::Pointer writer = WriterType::New();
 
                 string outfilename = outRoot;
                 outfilename.append("_NUM_DIRECTIONS.nrrd");
 
                 writer->SetFileName(outfilename.c_str());
                 writer->SetInput(numDirImage);
                 writer->Update();
             }
         }
 
         string logFile = outRoot;
         logFile.append("_ANGULAR_ERROR.csv");
         ofstream file;
         file.open (logFile.c_str());
 
         if (maskImages.size()>0)
         {
             for (unsigned int i=0; i<maskImages.size(); i++)
             {
                 mitk::Image::Pointer mitkMaskImage = dynamic_cast<mitk::Image*>(mitk::IOUtil::LoadDataNode(maskImages.at(i))->GetData());
                 mitk::CastToItkImage(mitkMaskImage, itkMaskImage);
 
                 // evaluate directions
                 EvaluationFilterType::Pointer evaluationFilter = EvaluationFilterType::New();
                 evaluationFilter->SetImageSet(directionImageContainer);
                 evaluationFilter->SetReferenceImageSet(referenceImageContainer);
                 evaluationFilter->SetMaskImage(itkMaskImage);
                 evaluationFilter->SetIgnoreMissingDirections(ignore);
                 evaluationFilter->Update();
 
                 if (verbose)
                 {
                     EvaluationFilterType::OutputImageType::Pointer angularErrorImage = evaluationFilter->GetOutput(0);
                     typedef itk::ImageFileWriter< EvaluationFilterType::OutputImageType > WriterType;
                     WriterType::Pointer writer = WriterType::New();
 
                     string outfilename = outRoot;
                     outfilename.append("_ERROR_IMAGE.nrrd");
 
                     writer->SetFileName(outfilename.c_str());
                     writer->SetInput(angularErrorImage);
                     writer->Update();
                 }
 
                 string maskFileName = itksys::SystemTools::GetFilenameWithoutExtension(maskImages.at(i));
                 unsigned found = maskFileName.find_last_of("_");
 
                 string sens = itksys::SystemTools::GetFilenameWithoutLastExtension(fibFile);
                 if (!fileID.empty())
                     sens = fileID;
                 sens.append(",");
 
                 sens.append(maskFileName.substr(found+1));
                 sens.append(",");
 
                 sens.append(boost::lexical_cast<string>(evaluationFilter->GetMeanAngularError()));
                 sens.append(",");
 
                 sens.append(boost::lexical_cast<string>(evaluationFilter->GetMedianAngularError()));
                 sens.append(",");
 
                 sens.append(boost::lexical_cast<string>(evaluationFilter->GetMaxAngularError()));
                 sens.append(",");
 
                 sens.append(boost::lexical_cast<string>(evaluationFilter->GetMinAngularError()));
                 sens.append(",");
 
                 sens.append(boost::lexical_cast<string>(std::sqrt(evaluationFilter->GetVarAngularError())));
                 sens.append(";\n");
                 file << sens;
             }
         }
         else
         {
             // evaluate directions
             EvaluationFilterType::Pointer evaluationFilter = EvaluationFilterType::New();
             evaluationFilter->SetImageSet(directionImageContainer);
             evaluationFilter->SetReferenceImageSet(referenceImageContainer);
             evaluationFilter->SetMaskImage(itkMaskImage);
             evaluationFilter->SetIgnoreMissingDirections(ignore);
             evaluationFilter->Update();
 
             if (verbose)
             {
                 EvaluationFilterType::OutputImageType::Pointer angularErrorImage = evaluationFilter->GetOutput(0);
                 typedef itk::ImageFileWriter< EvaluationFilterType::OutputImageType > WriterType;
                 WriterType::Pointer writer = WriterType::New();
 
                 string outfilename = outRoot;
                 outfilename.append("_ERROR_IMAGE.nrrd");
 
                 writer->SetFileName(outfilename.c_str());
                 writer->SetInput(angularErrorImage);
                 writer->Update();
             }
 
             string sens = itksys::SystemTools::GetFilenameWithoutLastExtension(fibFile);
             if (!fileID.empty())
                 sens = fileID;
             sens.append(",");
 
             sens.append(boost::lexical_cast<string>(evaluationFilter->GetMeanAngularError()));
             sens.append(",");
 
             sens.append(boost::lexical_cast<string>(evaluationFilter->GetMedianAngularError()));
             sens.append(",");
 
             sens.append(boost::lexical_cast<string>(evaluationFilter->GetMaxAngularError()));
             sens.append(",");
 
             sens.append(boost::lexical_cast<string>(evaluationFilter->GetMinAngularError()));
             sens.append(",");
 
             sens.append(boost::lexical_cast<string>(std::sqrt(evaluationFilter->GetVarAngularError())));
             sens.append(";\n");
             file << sens;
         }
         file.close();
     }
     catch (itk::ExceptionObject e)
     {
         std::cout << e;
         return EXIT_FAILURE;
     }
     catch (std::exception e)
     {
         std::cout << e.what();
         return EXIT_FAILURE;
     }
     catch (...)
     {
         std::cout << "ERROR!?!";
         return EXIT_FAILURE;
     }
     return EXIT_SUCCESS;
 }
diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkMLBTView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkMLBTView.cpp
index 992e0c577e..82ffd5cfe0 100644
--- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkMLBTView.cpp
+++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkMLBTView.cpp
@@ -1,329 +1,364 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 
 
 // Blueberry
 #include <berryISelectionService.h>
 #include <berryIWorkbenchWindow.h>
 
 // Qmitk
 #include "QmitkMLBTView.h"
 #include "QmitkStdMultiWidget.h"
 
 // Qt
 #include <QMessageBox>
 #include <QFileDialog>
 #include <QString>
 
 #include "mitkNodePredicateDataType.h"
 #include <mitkNodePredicateProperty.h>
 #include <mitkNodePredicateAnd.h>
 #include <mitkNodePredicateNot.h>
 #include <mitkNodePredicateOr.h>
 #include <mitkImageCast.h>
+#include <mitkEnumerationProperty.h>
+#include <mitkPointSetShapeProperty.h>
+#include <mitkPointSet.h>
 
 #define _USE_MATH_DEFINES
 #include <math.h>
 
 const std::string QmitkMLBTView::VIEW_ID = "org.mitk.views.mlbtview";
 using namespace berry;
 
 QmitkMLBTView::QmitkMLBTView()
     : QmitkFunctionality()
     , m_Controls( 0 )
     , m_MultiWidget( NULL )
 {
     m_TrackingTimer = new QTimer(this);
 }
 
 // Destructor
 QmitkMLBTView::~QmitkMLBTView()
 {
 
 }
 
 void QmitkMLBTView::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::QmitkMLBTViewControls;
         m_Controls->setupUi( parent );
 
         connect( m_Controls->m_StartTrainingButton, SIGNAL ( clicked() ), this, SLOT( StartTrainingThread() ) );
         connect( &m_TrainingWatcher, SIGNAL ( finished() ), this, SLOT( OnTrainingThreadStop() ) );
         connect( m_Controls->m_StartTrackingButton, SIGNAL ( clicked() ), this, SLOT( StartTrackingThread() ) );
         connect( &m_TrackingWatcher, SIGNAL ( finished() ), this, SLOT( OnTrackingThreadStop() ) );
         connect( m_Controls->m_SaveForestButton, SIGNAL ( clicked() ), this, SLOT( SaveForest() ) );
         connect( m_Controls->m_LoadForestButton, SIGNAL ( clicked() ), this, SLOT( LoadForest() ) );
         connect( m_TrackingTimer, SIGNAL(timeout()), this, SLOT(BuildFibers()) );
         connect( m_Controls->m_TimerIntervalBox, SIGNAL(valueChanged(int)), this, SLOT( ChangeTimerInterval(int) ));
         connect( m_Controls->m_DemoModeBox, SIGNAL(stateChanged(int)), this, SLOT( ToggleDemoMode(int) ));
         connect( m_Controls->m_PauseTrackingButton, SIGNAL ( clicked() ), this, SLOT( PauseTracking() ) );
         connect( m_Controls->m_AbortTrackingButton, SIGNAL ( clicked() ), this, SLOT( AbortTracking() ) );
 
         int numThread = itk::MultiThreader::GetGlobalDefaultNumberOfThreads();
         m_Controls->m_NumberOfThreadsBox->setMaximum(numThread);
         m_Controls->m_NumberOfThreadsBox->setValue(numThread);
 
         m_Controls->m_TrackingMaskImageBox->SetDataStorage(this->GetDataStorage());
         m_Controls->m_TrackingSeedImageBox->SetDataStorage(this->GetDataStorage());
         m_Controls->m_TrackingStopImageBox->SetDataStorage(this->GetDataStorage());
         m_Controls->m_TrackingRawImageBox->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 isQbi = mitk::NodePredicateDataType::New("QBallImage");
         mitk::NodePredicateOr::Pointer isDiffusionImage = mitk::NodePredicateOr::New(isDwi, isDti);
         isDiffusionImage = mitk::NodePredicateOr::New(isDiffusionImage, isQbi);
         mitk::NodePredicateNot::Pointer noDiffusionImage = mitk::NodePredicateNot::New(isDiffusionImage);
         mitk::NodePredicateAnd::Pointer finalPredicate = mitk::NodePredicateAnd::New(isMitkImage, noDiffusionImage);
         mitk::NodePredicateProperty::Pointer isBinaryPredicate = mitk::NodePredicateProperty::New("binary", mitk::BoolProperty::New(true));
         finalPredicate = mitk::NodePredicateAnd::New(finalPredicate, isBinaryPredicate);
         m_Controls->m_TrackingMaskImageBox->SetPredicate(finalPredicate);
         m_Controls->m_TrackingSeedImageBox->SetPredicate(finalPredicate);
         m_Controls->m_TrackingStopImageBox->SetPredicate(finalPredicate);
-        m_Controls->m_TrackingRawImageBox->SetPredicate(isDwi);
+        m_Controls->m_TrackingRawImageBox->SetPredicate(isMitkImage);
     }
 }
 
 void QmitkMLBTView::AbortTracking()
 {
     if (tracker.IsNotNull())
     {
         tracker->m_AbortTracking = true;
     }
 }
 
 void QmitkMLBTView::PauseTracking()
 {
     if (tracker.IsNotNull())
     {
         tracker->m_PauseTracking = !tracker->m_PauseTracking;
     }
 }
 
 void QmitkMLBTView::ChangeTimerInterval(int value)
 {
     m_TrackingTimer->setInterval(value);
 }
 
 void QmitkMLBTView::ToggleDemoMode(int state)
 {
     if (tracker.IsNotNull())
     {
         tracker->SetDemoMode(m_Controls->m_DemoModeBox->isChecked());
         tracker->m_Stop = false;
     }
 }
 
 void QmitkMLBTView::BuildFibers()
 {
     if (m_Controls->m_DemoModeBox->isChecked() && tracker.IsNotNull() && tracker->m_BuildFibersFinished)
     {
         vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData();
         mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly);
-        outFib->SetFiberColors(255,0,0);
+        outFib->SetFiberColors(255,255,255);
         m_TractogramNode->SetData(outFib);
+
+        m_SamplingPointsNode->SetData(tracker->m_SamplingPointset);
+        m_AlternativePointsNode->SetData(tracker->m_AlternativePointset);
+
         mitk::RenderingManager::GetInstance()->RequestUpdateAll();
         tracker->m_BuildFibersFinished = false;
-        tracker->m_Stop = false;
         tracker->m_BuildFibersReady = 0;
+        tracker->m_Stop = false;
     }
 }
 
 void QmitkMLBTView::LoadForest()
 {
 
     QString filename = QFileDialog::getOpenFileName(0, tr("Load Forest"), QDir::currentPath(), tr("HDF5 random forest file (*.rf)") );
     if(filename.isEmpty() || filename.isNull())
         return;
 
     m_ForestHandler.LoadForest( filename.toStdString() );
 }
 
 void QmitkMLBTView::StartTrackingThread()
 {
     m_TractogramNode = mitk::DataNode::New();
     m_TractogramNode->SetName("MLBT Result");
     m_TractogramNode->SetProperty("Fiber2DSliceThickness", mitk::FloatProperty::New(20));
     m_TractogramNode->SetProperty("Fiber2DfadeEFX", mitk::BoolProperty::New(false));
     m_TractogramNode->SetProperty("LineWidth", mitk::IntProperty::New(2));
     m_TractogramNode->SetProperty("color",mitk::ColorProperty::New(0, 1, 1));
     this->GetDataStorage()->Add(m_TractogramNode);
 
+    m_SamplingPointsNode = mitk::DataNode::New();
+    m_SamplingPointsNode->SetName("SamplingPoints");
+    m_SamplingPointsNode->SetProperty("pointsize", mitk::FloatProperty::New(0.2));
+    m_SamplingPointsNode->SetProperty("color", mitk::ColorProperty::New(1,1,1));
+    mitk::PointSetShapeProperty::Pointer bla = mitk::PointSetShapeProperty::New(); bla->SetValue(8);
+    m_SamplingPointsNode->SetProperty("Pointset.2D.shape", bla);
+    m_SamplingPointsNode->SetProperty("Pointset.2D.distance to plane", mitk::FloatProperty::New(1.5));
+    m_SamplingPointsNode->SetProperty("point 2D size", mitk::IntProperty::New(3));
+    m_SamplingPointsNode->SetProperty("Pointset.2D.fill shape", mitk::BoolProperty::New(true));
+    this->GetDataStorage()->Add(m_SamplingPointsNode);
+
+    m_AlternativePointsNode = mitk::DataNode::New();
+    m_AlternativePointsNode->SetName("AlternativePoints");
+    m_AlternativePointsNode->SetProperty("pointsize", mitk::FloatProperty::New(0.2));
+    m_AlternativePointsNode->SetProperty("color", mitk::ColorProperty::New(1,0,0));
+    m_AlternativePointsNode->SetProperty("Pointset.2D.shape", bla);
+    m_AlternativePointsNode->SetProperty("Pointset.2D.distance to plane", mitk::FloatProperty::New(1.5));
+    m_AlternativePointsNode->SetProperty("point 2D size", mitk::IntProperty::New(3));
+    m_AlternativePointsNode->SetProperty("Pointset.2D.fill shape", mitk::BoolProperty::New(true));
+    this->GetDataStorage()->Add(m_AlternativePointsNode);
+
     QFuture<void> future = QtConcurrent::run( this, &QmitkMLBTView::StartTracking );
     m_TrackingWatcher.setFuture(future);
     m_TrackingThreadIsRunning = true;
     m_Controls->m_StartTrackingButton->setEnabled(false);
-    m_TrackingTimer->start(10);
+    m_TrackingTimer->start(m_Controls->m_TimerIntervalBox->value());
 }
 
 void QmitkMLBTView::OnTrackingThreadStop()
 {
     m_TrackingThreadIsRunning = false;
     m_Controls->m_StartTrackingButton->setEnabled(true);
 
-
     vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData();
     mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly);
-    outFib->SetFiberColors(255,0,0);
+    outFib->SetFiberColors(255,255,255);
 //    mitk::DataNode::Pointer node = mitk::DataNode::New();
     m_TractogramNode->SetData(outFib);
+    m_SamplingPointsNode->SetData(tracker->m_SamplingPointset);
+    m_AlternativePointsNode->SetData(tracker->m_AlternativePointset);
 
     tracker = NULL;
     m_TrackingTimer->stop();
+
+    mitk::RenderingManager::GetInstance()->RequestUpdateAll();
 }
 
 void QmitkMLBTView::StartTracking()
 {
     if ( m_Controls->m_TrackingRawImageBox->GetSelectedNode().IsNull() )
         return;
 
     mitk::Image::Pointer dwi = dynamic_cast<mitk::Image*>(m_Controls->m_TrackingRawImageBox->GetSelectedNode()->GetData());
     tracker = TrackerType::New();
     tracker->SetNumberOfThreads(m_Controls->m_NumberOfThreadsBox->value());
     tracker->SetInput(0,  mitk::DiffusionPropertyHelper::GetItkVectorImage(dwi) );
     tracker->SetGradientDirections( mitk::DiffusionPropertyHelper::GetGradientContainer(dwi) );
     tracker->SetB_Value( mitk::DiffusionPropertyHelper::GetReferenceBValue(dwi) );
     tracker->SetDemoMode(m_Controls->m_DemoModeBox->isChecked());
+    if (m_Controls->m_DemoModeBox->isChecked())
+        tracker->SetNumberOfThreads(1);
     if (m_Controls->m_TrackingUseMaskImageBox->isChecked() && m_Controls->m_TrackingMaskImageBox->GetSelectedNode().IsNotNull())
     {
         mitk::Image::Pointer mask = dynamic_cast<mitk::Image*>(m_Controls->m_TrackingMaskImageBox->GetSelectedNode()->GetData());
         ItkUcharImgType::Pointer itkMask = ItkUcharImgType::New();
         mitk::CastToItkImage(mask, itkMask);
         tracker->SetMaskImage(itkMask);
     }
     if (m_Controls->m_TrackingUseSeedImageBox->isChecked() && m_Controls->m_TrackingSeedImageBox->GetSelectedNode().IsNotNull())
     {
         mitk::Image::Pointer img = dynamic_cast<mitk::Image*>(m_Controls->m_TrackingSeedImageBox->GetSelectedNode()->GetData());
         ItkUcharImgType::Pointer itkImg = ItkUcharImgType::New();
         mitk::CastToItkImage(img, itkImg);
         tracker->SetSeedImage(itkImg);
     }
     if (m_Controls->m_TrackingUseStopImageBox->isChecked() && m_Controls->m_TrackingStopImageBox->GetSelectedNode().IsNotNull())
     {
         mitk::Image::Pointer img = dynamic_cast<mitk::Image*>(m_Controls->m_TrackingStopImageBox->GetSelectedNode()->GetData());
         ItkUcharImgType::Pointer itkImg = ItkUcharImgType::New();
         mitk::CastToItkImage(img, itkImg);
         tracker->SetStoppingRegions(itkImg);
     }
     tracker->SetSeedsPerVoxel(m_Controls->m_NumberOfSeedsBox->value());
-    tracker->SetUseDirection(true);
     tracker->SetStepSize(m_Controls->m_TrackingStepSizeBox->value());
     tracker->SetAngularThreshold(cos((float)m_Controls->m_AngularThresholdBox->value()*M_PI/180));
     tracker->SetMinTractLength(m_Controls->m_MinLengthBox->value());
     tracker->SetMaxTractLength(m_Controls->m_MaxLengthBox->value());
+    tracker->SetAposterioriCurvCheck(m_Controls->m_Curvcheck2->isChecked());
+    tracker->SetRemoveWmEndFibers(false);
+    tracker->SetAvoidStop(m_Controls->m_AvoidStop->isChecked());
 
     vigra::RandomForest<int> forest = m_ForestHandler.GetForest();
     tracker->SetDecisionForest(&forest);
     tracker->SetSamplingDistance(m_Controls->m_SamplingDistanceBox->value());
     tracker->SetNumberOfSamples(m_Controls->m_NumSamplesBox->value());
+    tracker->SetRandomSampling(m_Controls->m_RandomSampling->isChecked());
     tracker->Update();
 //    vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData();
 //    mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly);
 //    outFib->SetColorCoding(mitk::FiberBundle::COLORCODING_CUSTOM);
 //    mitk::DataNode::Pointer node = mitk::DataNode::New();
 //    m_TractogramNode->SetData(outFib);
 //    node->SetData(outFib);
 //    node->SetName("MLBT Result");
 //    this->GetDataStorage()->Add(node);
 //    mitk::RenderingManager::GetInstance()->RequestUpdateAll();
 }
 
 void QmitkMLBTView::SaveForest()
 {
     QString filename = QFileDialog::getSaveFileName(0, tr("Save Forest"), QDir::currentPath()+"/forest.rf", tr("HDF5 random forest file (*.rf)") );
 
     if(filename.isEmpty() || filename.isNull())
         return;
     if(!filename.endsWith(".rf"))
         filename += ".rf";
 
     m_ForestHandler.SaveForest( filename.toStdString() );
 }
 
 void QmitkMLBTView::StartTrainingThread()
 {
     QFuture<void> future = QtConcurrent::run( this, &QmitkMLBTView::StartTraining );
     m_TrainingWatcher.setFuture(future);
     m_Controls->m_StartTrainingButton->setEnabled(false);
     m_Controls->m_SaveForestButton->setEnabled(false);
     m_Controls->m_LoadForestButton->setEnabled(false);
 }
 
 void QmitkMLBTView::OnTrainingThreadStop()
 {
     m_Controls->m_StartTrainingButton->setEnabled(true);
     m_Controls->m_SaveForestButton->setEnabled(true);
     m_Controls->m_LoadForestButton->setEnabled(true);
 }
 
 void QmitkMLBTView::StartTraining()
 {
     m_ForestHandler.SetRawData(m_SelectedDiffImages);
     m_ForestHandler.SetTractograms(m_SelectedFB);
     m_ForestHandler.SetNumTrees(m_Controls->m_NumTreesBox->value());
     m_ForestHandler.SetMaxTreeDepth(m_Controls->m_MaxDepthBox->value());
     m_ForestHandler.SetGrayMatterSamplesPerVoxel(m_Controls->m_GmSamplingBox->value());
     m_ForestHandler.SetSampleFraction(m_Controls->m_SampleFractionBox->value());
     m_ForestHandler.SetStepSize(m_Controls->m_TrainingStepSizeBox->value());
-    m_ForestHandler.SetUsePreviousDirection(m_Controls->m_TrackingUsePreviousDirectionBox->isChecked());
     m_ForestHandler.StartTraining();
 }
 
 void QmitkMLBTView::StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget)
 {
     m_MultiWidget = &stdMultiWidget;
 }
 
 
 void QmitkMLBTView::StdMultiWidgetNotAvailable()
 {
     m_MultiWidget = NULL;
 }
 
 void QmitkMLBTView::OnSelectionChanged( std::vector<mitk::DataNode*> nodes )
 {
     if ( !this->IsVisible() )
     {
         // do nothing if nobody wants to see me :-(
         return;
     }
 
     m_SelectedFB.clear();
     m_SelectedDiffImages.clear();
     m_MaskImages.clear();
     m_WhiteMatterImages.clear();
 
     for( std::vector<mitk::DataNode*>::iterator it = nodes.begin(); it != nodes.end(); ++it )
     {
         mitk::DataNode::Pointer node = *it;
         if ( dynamic_cast<mitk::FiberBundle*>(node->GetData()) )
             m_SelectedFB.push_back( dynamic_cast<mitk::FiberBundle*>(node->GetData()) );
         else if (mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(node))
             m_SelectedDiffImages.push_back( dynamic_cast<mitk::Image*>(node->GetData()) );
     }
 }
 
 void QmitkMLBTView::Activated()
 {
 
 }
 
 
diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkMLBTView.h b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkMLBTView.h
index 9b260c22db..06e706996f 100644
--- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkMLBTView.h
+++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkMLBTView.h
@@ -1,105 +1,107 @@
 /*===================================================================
 
 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 <QmitkFunctionality.h>
 #include "ui_QmitkMLBTViewControls.h"
 
 #ifndef Q_MOC_RUN
 #include "mitkDataStorage.h"
 #include "mitkDataStorageSelection.h"
 #include <mitkTrackingForestHandler.h>
 #include <itkMLBSTrackingFilter.h>
 #endif
 
 #include <QtConcurrentRun>
 #include <QFutureWatcher>
 #include <QTimer>
 
 /*!
 \brief
 */
 
 // Forward Qt class declarations
 
 
 class QmitkMLBTView : public QmitkFunctionality
 {
 
 
   // this is needed for all Qt objects that should have a Qt meta-object
   // (everything that derives from QObject and wants to have signal/slots)
   Q_OBJECT
 
 public:
 
   static const std::string VIEW_ID;
 
   typedef itk::Image<unsigned char, 3>      ItkUcharImgType;
   typedef itk::MLBSTrackingFilter<100>      TrackerType;
 
   QmitkMLBTView();
   virtual ~QmitkMLBTView();
 
   virtual void CreateQtPartControl(QWidget *parent);
 
   virtual void StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget);
   virtual void StdMultiWidgetNotAvailable();
   virtual void Activated();
 
   protected slots:
 
   void StartTrackingThread();
   void OnTrackingThreadStop();
   void StartTrainingThread();
   void OnTrainingThreadStop();
   void SaveForest();
   void LoadForest();
   void BuildFibers();
   void ChangeTimerInterval(int value);
   void ToggleDemoMode(int state);
   void PauseTracking();
   void AbortTracking();
 
 protected:
 
   void StartTracking();
   void StartTraining();
 
   /// \brief called by QmitkFunctionality when DataManager's selection has changed
   virtual void OnSelectionChanged( std::vector<mitk::DataNode*> nodes );
 
   Ui::QmitkMLBTViewControls* m_Controls;
   QmitkStdMultiWidget* m_MultiWidget;
 
   mitk::TrackingForestHandler<> m_ForestHandler;
   std::vector< mitk::Image::Pointer > m_SelectedDiffImages;
   std::vector< mitk::FiberBundle::Pointer > m_SelectedFB;
   std::vector< ItkUcharImgType::Pointer > m_MaskImages;
   std::vector< ItkUcharImgType::Pointer > m_WhiteMatterImages;
 
   QFutureWatcher<void> m_TrainingWatcher;
   QFutureWatcher<void> m_TrackingWatcher;
   bool m_TrackingThreadIsRunning;
   TrackerType::Pointer tracker;
   QTimer*   m_TrackingTimer;
   mitk::DataNode::Pointer m_TractogramNode;
+  mitk::DataNode::Pointer m_SamplingPointsNode;
+  mitk::DataNode::Pointer m_AlternativePointsNode;
 
 private:
 
  };
 
 
 
diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkMLBTViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkMLBTViewControls.ui
index cf19c52391..a641f8762c 100644
--- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkMLBTViewControls.ui
+++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkMLBTViewControls.ui
@@ -1,620 +1,640 @@
 <?xml version="1.0" encoding="UTF-8"?>
 <ui version="4.0">
  <class>QmitkMLBTViewControls</class>
  <widget class="QWidget" name="QmitkMLBTViewControls">
   <property name="geometry">
    <rect>
     <x>0</x>
     <y>0</y>
     <width>359</width>
     <height>1127</height>
    </rect>
   </property>
   <property name="windowTitle">
    <string>Form</string>
   </property>
   <layout class="QVBoxLayout" name="verticalLayout">
    <property name="leftMargin">
     <number>0</number>
    </property>
    <property name="topMargin">
     <number>0</number>
    </property>
    <property name="rightMargin">
     <number>0</number>
    </property>
    <property name="bottomMargin">
     <number>0</number>
    </property>
    <item>
     <widget class="QToolBox" name="toolBox">
      <property name="currentIndex">
       <number>1</number>
      </property>
      <widget class="QWidget" name="page">
       <property name="geometry">
        <rect>
         <x>0</x>
         <y>0</y>
-        <width>359</width>
-        <height>1065</height>
+        <width>290</width>
+        <height>376</height>
        </rect>
       </property>
       <attribute name="label">
        <string>Training</string>
       </attribute>
       <layout class="QGridLayout" name="gridLayout">
        <item row="4" column="0">
         <widget class="QCommandLinkButton" name="m_LoadForestButton">
          <property name="text">
           <string>Load Forest</string>
          </property>
         </widget>
        </item>
        <item row="5" column="0">
         <spacer name="verticalSpacer_2">
          <property name="orientation">
           <enum>Qt::Vertical</enum>
          </property>
          <property name="sizeHint" stdset="0">
           <size>
            <width>20</width>
            <height>40</height>
           </size>
          </property>
         </spacer>
        </item>
        <item row="0" column="0">
         <widget class="QFrame" name="frame">
          <property name="frameShape">
           <enum>QFrame::NoFrame</enum>
          </property>
          <property name="frameShadow">
           <enum>QFrame::Raised</enum>
          </property>
          <layout class="QGridLayout" name="gridLayout_3">
           <property name="leftMargin">
            <number>0</number>
           </property>
           <property name="topMargin">
            <number>0</number>
           </property>
           <property name="rightMargin">
            <number>0</number>
           </property>
           <property name="bottomMargin">
            <number>0</number>
           </property>
           <item row="4" column="0">
            <widget class="QLabel" name="label_5">
             <property name="text">
              <string>Step size:</string>
             </property>
            </widget>
           </item>
           <item row="0" column="0">
            <widget class="QLabel" name="label">
             <property name="text">
              <string>Num. Trees:</string>
             </property>
            </widget>
           </item>
           <item row="1" column="0">
            <widget class="QLabel" name="label_2">
             <property name="text">
              <string>Max. Depth:</string>
             </property>
            </widget>
           </item>
           <item row="4" column="1">
            <widget class="QDoubleSpinBox" name="m_TrainingStepSizeBox">
             <property name="decimals">
              <number>3</number>
             </property>
             <property name="minimum">
              <double>-1.000000000000000</double>
             </property>
             <property name="maximum">
              <double>999.000000000000000</double>
             </property>
             <property name="singleStep">
              <double>0.100000000000000</double>
             </property>
             <property name="value">
              <double>-1.000000000000000</double>
             </property>
            </widget>
           </item>
           <item row="0" column="1">
            <widget class="QSpinBox" name="m_NumTreesBox">
             <property name="minimum">
              <number>1</number>
             </property>
             <property name="maximum">
              <number>999999999</number>
             </property>
            </widget>
           </item>
           <item row="3" column="0">
            <widget class="QLabel" name="label_4">
             <property name="text">
              <string>Sample Fraction:</string>
             </property>
            </widget>
           </item>
           <item row="1" column="1">
            <widget class="QSpinBox" name="m_MaxDepthBox">
             <property name="minimum">
              <number>1</number>
             </property>
             <property name="maximum">
              <number>999999999</number>
             </property>
             <property name="value">
              <number>10</number>
             </property>
            </widget>
           </item>
           <item row="2" column="1">
            <widget class="QSpinBox" name="m_GmSamplingBox">
             <property name="minimum">
              <number>1</number>
             </property>
             <property name="maximum">
              <number>999999999</number>
             </property>
             <property name="value">
              <number>10</number>
             </property>
            </widget>
           </item>
           <item row="3" column="1">
            <widget class="QDoubleSpinBox" name="m_SampleFractionBox">
             <property name="decimals">
              <number>3</number>
             </property>
             <property name="maximum">
              <double>1.000000000000000</double>
             </property>
             <property name="singleStep">
              <double>0.100000000000000</double>
             </property>
             <property name="value">
              <double>1.000000000000000</double>
             </property>
            </widget>
           </item>
           <item row="2" column="0">
            <widget class="QLabel" name="label_3">
             <property name="text">
              <string>GM Sampling Points:</string>
             </property>
            </widget>
           </item>
           <item row="5" column="0">
            <widget class="QLabel" name="label_6">
             <property name="text">
              <string>Use Previous Direction:</string>
             </property>
            </widget>
           </item>
           <item row="5" column="1">
            <widget class="QCheckBox" name="m_TrainingUsePreviousDirectionBox">
             <property name="text">
              <string/>
             </property>
             <property name="checked">
              <bool>true</bool>
             </property>
            </widget>
           </item>
          </layout>
         </widget>
        </item>
        <item row="2" column="0">
         <widget class="QCommandLinkButton" name="m_StartTrainingButton">
          <property name="text">
           <string>Start Training</string>
          </property>
         </widget>
        </item>
        <item row="3" column="0">
         <widget class="QCommandLinkButton" name="m_SaveForestButton">
          <property name="text">
           <string>Save Forest</string>
          </property>
         </widget>
        </item>
       </layout>
      </widget>
      <widget class="QWidget" name="page_2">
       <property name="geometry">
        <rect>
         <x>0</x>
         <y>0</y>
         <width>359</width>
         <height>1065</height>
        </rect>
       </property>
       <attribute name="label">
        <string>Tractography</string>
       </attribute>
       <layout class="QGridLayout" name="gridLayout_2">
-       <item row="2" column="0">
-        <widget class="QCheckBox" name="m_TrackingUsePreviousDirectionBox">
+       <item row="3" column="0">
+        <widget class="QCheckBox" name="m_AvoidStop">
          <property name="text">
-          <string>Use Previous Direction</string>
+          <string>Avoid premature termination</string>
          </property>
          <property name="checked">
           <bool>true</bool>
          </property>
         </widget>
        </item>
+       <item row="12" column="0">
+        <spacer name="verticalSpacer_3">
+         <property name="orientation">
+          <enum>Qt::Vertical</enum>
+         </property>
+         <property name="sizeHint" stdset="0">
+          <size>
+           <width>20</width>
+           <height>40</height>
+          </size>
+         </property>
+        </spacer>
+       </item>
        <item row="4" column="0">
-        <widget class="QFrame" name="frame_5">
+        <widget class="QCheckBox" name="m_Curvcheck2">
+         <property name="text">
+          <string>Secondary curvature check</string>
+         </property>
+         <property name="checked">
+          <bool>true</bool>
+         </property>
+        </widget>
+       </item>
+       <item row="10" column="0">
+        <widget class="QFrame" name="frame_4">
          <property name="frameShape">
           <enum>QFrame::NoFrame</enum>
          </property>
          <property name="frameShadow">
           <enum>QFrame::Raised</enum>
          </property>
-         <layout class="QGridLayout" name="gridLayout_7">
+         <layout class="QGridLayout" name="gridLayout_6">
           <property name="leftMargin">
            <number>0</number>
           </property>
           <property name="topMargin">
            <number>0</number>
           </property>
           <property name="rightMargin">
            <number>0</number>
           </property>
           <property name="bottomMargin">
            <number>0</number>
           </property>
-          <item row="0" column="1">
-           <widget class="QToolButton" name="m_PauseTrackingButton">
-            <property name="text">
-             <string>...</string>
-            </property>
-            <property name="icon">
-             <iconset resource="../../../../../../../debug-binary/MBI-build/Plugins/org.mbi.gui.qt.diffusionimaginginternal/org_mbi_gui_qt_diffusionimaginginternal_cached.qrc">
-              <normaloff>:/org.mbi.gui.qt.diffusionimaginginternal/resources/Media-playback-pause.svg</normaloff>:/org.mbi.gui.qt.diffusionimaginginternal/resources/Media-playback-pause.svg</iconset>
-            </property>
-           </widget>
-          </item>
           <item row="0" column="0">
-           <widget class="QToolButton" name="m_StartTrackingButton">
+           <widget class="QCheckBox" name="m_DemoModeBox">
             <property name="text">
-             <string>...</string>
-            </property>
-            <property name="icon">
-             <iconset resource="../../../../../../../debug-binary/MBI-build/Plugins/org.mbi.gui.qt.diffusionimaginginternal/org_mbi_gui_qt_diffusionimaginginternal_cached.qrc">
-              <normaloff>:/org.mbi.gui.qt.diffusionimaginginternal/resources/Media-playback-start.svg</normaloff>:/org.mbi.gui.qt.diffusionimaginginternal/resources/Media-playback-start.svg</iconset>
+             <string>Demo Mode</string>
             </property>
            </widget>
           </item>
-          <item row="0" column="2">
-           <widget class="QToolButton" name="m_AbortTrackingButton">
-            <property name="text">
-             <string>...</string>
+          <item row="0" column="1">
+           <widget class="QSpinBox" name="m_TimerIntervalBox">
+            <property name="minimum">
+             <number>1</number>
             </property>
-            <property name="icon">
-             <iconset resource="../../../../../../../debug-binary/MBI-build/Plugins/org.mbi.gui.qt.diffusionimaginginternal/org_mbi_gui_qt_diffusionimaginginternal_cached.qrc">
-              <normaloff>:/org.mbi.gui.qt.diffusionimaginginternal/resources/Media-playback-stop.svg</normaloff>:/org.mbi.gui.qt.diffusionimaginginternal/resources/Media-playback-stop.svg</iconset>
+            <property name="maximum">
+             <number>1000</number>
+            </property>
+            <property name="value">
+             <number>10</number>
             </property>
            </widget>
           </item>
          </layout>
         </widget>
        </item>
-       <item row="8" column="0">
-        <widget class="QFrame" name="frame_4">
+       <item row="5" column="0">
+        <widget class="QFrame" name="frame_3">
          <property name="frameShape">
           <enum>QFrame::NoFrame</enum>
          </property>
          <property name="frameShadow">
           <enum>QFrame::Raised</enum>
          </property>
-         <layout class="QGridLayout" name="gridLayout_6">
+         <layout class="QGridLayout" name="gridLayout_5">
           <property name="leftMargin">
            <number>0</number>
           </property>
           <property name="topMargin">
            <number>0</number>
           </property>
           <property name="rightMargin">
            <number>0</number>
           </property>
           <property name="bottomMargin">
            <number>0</number>
           </property>
-          <item row="0" column="0">
-           <widget class="QCheckBox" name="m_DemoModeBox">
+          <item row="6" column="0">
+           <widget class="QCheckBox" name="m_TrackingUseStopImageBox">
             <property name="text">
-             <string>Demo Mode</string>
+             <string>Use Stop Image:</string>
             </property>
            </widget>
           </item>
-          <item row="0" column="1">
-           <widget class="QSpinBox" name="m_TimerIntervalBox">
-            <property name="minimum">
-             <number>1</number>
+          <item row="4" column="0">
+           <widget class="QCheckBox" name="m_TrackingUseSeedImageBox">
+            <property name="text">
+             <string>Use Seed Image:</string>
             </property>
-            <property name="maximum">
-             <number>1000</number>
+            <property name="checked">
+             <bool>false</bool>
             </property>
-            <property name="value">
-             <number>10</number>
+           </widget>
+          </item>
+          <item row="1" column="0">
+           <widget class="QCheckBox" name="m_TrackingUseMaskImageBox">
+            <property name="text">
+             <string>Use Mask Image:</string>
+            </property>
+            <property name="checked">
+             <bool>true</bool>
             </property>
            </widget>
           </item>
+          <item row="1" column="1">
+           <widget class="QmitkDataStorageComboBox" name="m_TrackingMaskImageBox"/>
+          </item>
+          <item row="4" column="1">
+           <widget class="QmitkDataStorageComboBox" name="m_TrackingSeedImageBox"/>
+          </item>
+          <item row="6" column="1">
+           <widget class="QmitkDataStorageComboBox" name="m_TrackingStopImageBox"/>
+          </item>
          </layout>
         </widget>
        </item>
        <item row="1" column="0">
         <widget class="QFrame" name="frame_2">
          <property name="frameShape">
           <enum>QFrame::NoFrame</enum>
          </property>
          <property name="frameShadow">
           <enum>QFrame::Raised</enum>
          </property>
          <layout class="QGridLayout" name="gridLayout_4">
           <property name="leftMargin">
            <number>0</number>
           </property>
           <property name="topMargin">
            <number>0</number>
           </property>
           <property name="rightMargin">
            <number>0</number>
           </property>
           <property name="bottomMargin">
            <number>0</number>
           </property>
           <item row="2" column="1">
            <widget class="QDoubleSpinBox" name="m_TrackingStepSizeBox">
             <property name="value">
              <double>0.500000000000000</double>
             </property>
            </widget>
           </item>
           <item row="7" column="0">
            <widget class="QLabel" name="label_14">
             <property name="text">
              <string>Max. Length</string>
             </property>
            </widget>
           </item>
           <item row="4" column="1">
            <widget class="QSpinBox" name="m_NumSamplesBox">
             <property name="maximum">
              <number>999999999</number>
             </property>
             <property name="value">
              <number>50</number>
             </property>
            </widget>
           </item>
           <item row="4" column="0">
            <widget class="QLabel" name="label_9">
             <property name="text">
              <string>Num. Samples:</string>
             </property>
            </widget>
           </item>
           <item row="0" column="0">
            <widget class="QLabel" name="label_12">
             <property name="text">
              <string>Input DWI:</string>
             </property>
            </widget>
           </item>
           <item row="2" column="0">
            <widget class="QLabel" name="label_8">
             <property name="text">
              <string>Step Size:</string>
             </property>
            </widget>
           </item>
           <item row="5" column="0">
            <widget class="QLabel" name="label_10">
             <property name="text">
              <string>Sampling Distance:</string>
             </property>
            </widget>
           </item>
           <item row="1" column="1">
            <widget class="QSpinBox" name="m_NumberOfSeedsBox">
             <property name="minimum">
              <number>1</number>
             </property>
             <property name="maximum">
              <number>999</number>
             </property>
            </widget>
           </item>
           <item row="6" column="0">
            <widget class="QLabel" name="label_13">
             <property name="text">
              <string>Min. Length</string>
             </property>
            </widget>
           </item>
           <item row="5" column="1">
            <widget class="QDoubleSpinBox" name="m_SamplingDistanceBox">
             <property name="singleStep">
              <double>0.100000000000000</double>
             </property>
             <property name="value">
              <double>0.500000000000000</double>
             </property>
            </widget>
           </item>
           <item row="6" column="1">
            <widget class="QDoubleSpinBox" name="m_MinLengthBox">
             <property name="maximum">
              <double>999999999.000000000000000</double>
             </property>
             <property name="singleStep">
              <double>1.000000000000000</double>
             </property>
             <property name="value">
              <double>20.000000000000000</double>
             </property>
            </widget>
           </item>
           <item row="3" column="0">
            <widget class="QLabel" name="label_11">
             <property name="text">
              <string>Angular Threshold:</string>
             </property>
            </widget>
           </item>
           <item row="1" column="0">
            <widget class="QLabel" name="label_7">
             <property name="text">
              <string>Num. Seeds:</string>
             </property>
            </widget>
           </item>
           <item row="7" column="1">
            <widget class="QDoubleSpinBox" name="m_MaxLengthBox">
             <property name="maximum">
              <double>999999999.000000000000000</double>
             </property>
             <property name="singleStep">
              <double>1.000000000000000</double>
             </property>
             <property name="value">
              <double>400.000000000000000</double>
             </property>
            </widget>
           </item>
           <item row="3" column="1">
            <widget class="QDoubleSpinBox" name="m_AngularThresholdBox">
             <property name="decimals">
              <number>1</number>
             </property>
             <property name="maximum">
              <double>90.000000000000000</double>
             </property>
             <property name="value">
              <double>45.000000000000000</double>
             </property>
            </widget>
           </item>
           <item row="0" column="1">
            <widget class="QmitkDataStorageComboBox" name="m_TrackingRawImageBox"/>
           </item>
           <item row="8" column="0">
            <widget class="QLabel" name="label_15">
             <property name="text">
              <string>Num. Threads:</string>
             </property>
            </widget>
           </item>
           <item row="8" column="1">
            <widget class="QSpinBox" name="m_NumberOfThreadsBox">
             <property name="minimum">
              <number>1</number>
             </property>
             <property name="maximum">
              <number>30</number>
             </property>
            </widget>
           </item>
          </layout>
         </widget>
        </item>
-       <item row="10" column="0">
-        <spacer name="verticalSpacer_3">
-         <property name="orientation">
-          <enum>Qt::Vertical</enum>
-         </property>
-         <property name="sizeHint" stdset="0">
-          <size>
-           <width>20</width>
-           <height>40</height>
-          </size>
-         </property>
-        </spacer>
-       </item>
-       <item row="3" column="0">
-        <widget class="QFrame" name="frame_3">
+       <item row="6" column="0">
+        <widget class="QFrame" name="frame_5">
          <property name="frameShape">
           <enum>QFrame::NoFrame</enum>
          </property>
          <property name="frameShadow">
           <enum>QFrame::Raised</enum>
          </property>
-         <layout class="QGridLayout" name="gridLayout_5">
+         <layout class="QGridLayout" name="gridLayout_7">
           <property name="leftMargin">
            <number>0</number>
           </property>
           <property name="topMargin">
            <number>0</number>
           </property>
           <property name="rightMargin">
            <number>0</number>
           </property>
           <property name="bottomMargin">
            <number>0</number>
           </property>
-          <item row="6" column="0">
-           <widget class="QCheckBox" name="m_TrackingUseStopImageBox">
+          <item row="0" column="1">
+           <widget class="QToolButton" name="m_PauseTrackingButton">
             <property name="text">
-             <string>Use Stop Image:</string>
+             <string>...</string>
+            </property>
+            <property name="icon">
+             <iconset resource="../../../org.mitk.gui.qt.ext/resources/org_mitk_icons.qrc">
+              <normaloff>:/org_mitk_icons/icons/tango/scalable/actions/media-playback-pause.svg</normaloff>:/org_mitk_icons/icons/tango/scalable/actions/media-playback-pause.svg</iconset>
             </property>
            </widget>
           </item>
-          <item row="4" column="0">
-           <widget class="QCheckBox" name="m_TrackingUseSeedImageBox">
+          <item row="0" column="0">
+           <widget class="QToolButton" name="m_StartTrackingButton">
             <property name="text">
-             <string>Use Seed Image:</string>
+             <string>...</string>
             </property>
-            <property name="checked">
-             <bool>false</bool>
+            <property name="icon">
+             <iconset resource="../../../org.mitk.gui.qt.ext/resources/org_mitk_icons.qrc">
+              <normaloff>:/org_mitk_icons/icons/tango/scalable/actions/media-playback-start.svg</normaloff>:/org_mitk_icons/icons/tango/scalable/actions/media-playback-start.svg</iconset>
             </property>
            </widget>
           </item>
-          <item row="1" column="0">
-           <widget class="QCheckBox" name="m_TrackingUseMaskImageBox">
+          <item row="0" column="2">
+           <widget class="QToolButton" name="m_AbortTrackingButton">
             <property name="text">
-             <string>Use Mask Image:</string>
+             <string>...</string>
             </property>
-            <property name="checked">
-             <bool>true</bool>
+            <property name="icon">
+             <iconset resource="../../../org.mitk.gui.qt.ext/resources/org_mitk_icons.qrc">
+              <normaloff>:/org_mitk_icons/icons/tango/scalable/actions/media-playback-stop.svg</normaloff>:/org_mitk_icons/icons/tango/scalable/actions/media-playback-stop.svg</iconset>
             </property>
            </widget>
           </item>
-          <item row="1" column="1">
-           <widget class="QmitkDataStorageComboBox" name="m_TrackingMaskImageBox"/>
-          </item>
-          <item row="4" column="1">
-           <widget class="QmitkDataStorageComboBox" name="m_TrackingSeedImageBox"/>
-          </item>
-          <item row="6" column="1">
-           <widget class="QmitkDataStorageComboBox" name="m_TrackingStopImageBox"/>
-          </item>
          </layout>
         </widget>
        </item>
+       <item row="2" column="0">
+        <widget class="QCheckBox" name="m_RandomSampling">
+         <property name="text">
+          <string>Random sampling</string>
+         </property>
+         <property name="checked">
+          <bool>false</bool>
+         </property>
+        </widget>
+       </item>
       </layout>
      </widget>
     </widget>
    </item>
   </layout>
  </widget>
  <customwidgets>
   <customwidget>
    <class>QmitkDataStorageComboBox</class>
    <extends>QComboBox</extends>
    <header location="global">QmitkDataStorageComboBox.h</header>
   </customwidget>
  </customwidgets>
  <resources>
-  <include location="../../../../../../../debug-binary/MBI-build/Plugins/org.mbi.gui.qt.diffusionimaginginternal/org_mbi_gui_qt_diffusionimaginginternal_cached.qrc"/>
+  <include location="../../../org.mitk.gui.qt.ext/resources/org_mitk_icons.qrc"/>
  </resources>
  <connections/>
 </ui>