diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.cpp
index b431eff166..08791bbaa9 100644
--- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.cpp
@@ -1,832 +1,554 @@
 #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_UseTrilinearInterpolation(false),
     m_MaxNumDirections(3),
-    m_Thres(0.5),
+    m_SizeThreshold(0.2),
     m_NumDirectionsImage(NULL)
 {
     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();
+    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 crossings image
-    m_CrossingsImage = ItkUcharImgType::New();
-    m_CrossingsImage->SetSpacing( spacing );
-    m_CrossingsImage->SetOrigin( origin );
-    m_CrossingsImage->SetDirection( direction );
-    m_CrossingsImage->SetRegions( imageRegion );
-    m_CrossingsImage->Allocate();
-    m_CrossingsImage->FillBuffer(0);
-
     // 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();
+    for (unsigned int i=0; i<m_MaxNumDirections; i++)
+    {
+        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);
+    }
+
     // 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->ResampleFibers(minSpacing/3);
+    m_FiberBundle->ResampleFibers(minSpacing/10);
 
     // iterate over all fibers
     vtkSmartPointer<vtkPolyData> fiberPolyData = m_FiberBundle->GetFiberPolyData();
-    vtkSmartPointer<vtkCellArray> vLines = fiberPolyData->GetLines();
-    vLines->InitTraversal();
     int numFibers = m_FiberBundle->GetNumFibers();
-    itk::TimeProbe clock;
     m_DirectionsContainer = ContainerType::New();
 
-    if (m_UseTrilinearInterpolation)
-        MITK_INFO << "Generating directions from tractogram (trilinear interpolation)";
-    else
-        MITK_INFO << "Generating directions from tractogram";
-
+    MITK_INFO << "Generating directions from tractogram";
     boost::progress_display disp(numFibers);
     for( int i=0; i<numFibers; i++ )
     {
         ++disp;
-        clock.Start();
-        vtkIdType   numPoints(0);
-        vtkIdType*  points(NULL);
-        vLines->GetNextCell ( numPoints, points );
+        vtkCell* cell = fiberPolyData->GetCell(i);
+        int numPoints = cell->GetNumberOfPoints();
+        vtkPoints* points = cell->GetPoints();
         if (numPoints<2)
             continue;
 
-        itk::Index<3> index; index.Fill(0);
-        itk::ContinuousIndex<double, 3> contIndex;
-        vnl_vector_fixed<double, 3> dir, wDir;
+        vnl_vector_fixed<double, 3> dir;
         itk::Point<double, 3> worldPos;
         vnl_vector<double> v;
+
         for( int j=0; j<numPoints-1; j++)
         {
-            double* temp = fiberPolyData->GetPoint(points[j]);
+            // get current position along fiber in world coordinates
+            double* temp = points->GetPoint(j);
             worldPos = GetItkPoint(temp);
-            v = GetVnlVector(temp);
-
-            dir = GetVnlVector(fiberPolyData->GetPoint(points[j+1]))-v;
-            dir.normalize();
-
+            itk::Index<3> index;
             m_MaskImage->TransformPhysicalPointToIndex(worldPos, index);
-            m_MaskImage->TransformPhysicalPointToContinuousIndex(worldPos, contIndex);
-
-            if (m_MaskImage->GetPixel(index)==0)
+            if (!m_MaskImage->GetLargestPossibleRegion().IsInside(index) || m_MaskImage->GetPixel(index)==0)
                 continue;
 
-            if (!m_UseTrilinearInterpolation)
-            {
-                if (index[0] < 0 || (unsigned long)index[0] >= outImageSize[0])
-                    continue;
-                if (index[1] < 0 || (unsigned long)index[1] >= outImageSize[1])
-                    continue;
-                if (index[2] < 0 || (unsigned long)index[2] >= outImageSize[2])
-                    continue;
-
-                unsigned int idx = index[0] + outImageSize[0]*(index[1] + outImageSize[1]*index[2]);
-                DirectionContainerType::Pointer dirCont = DirectionContainerType::New();
-                if (m_DirectionsContainer->IndexExists(idx))
-                {
-                    dirCont = m_DirectionsContainer->GetElement(idx);
-                    if (dirCont.IsNull())
-                    {
-                        dirCont = DirectionContainerType::New();
-                        dirCont->InsertElement(0, dir);
-                        m_DirectionsContainer->InsertElement(idx, dirCont);
-                    }
-                    else
-                        dirCont->InsertElement(dirCont->Size(), dir);
-                }
-                else
-                {
-                    dirCont->InsertElement(0, dir);
-                    m_DirectionsContainer->InsertElement(idx, dirCont);
-                }
-
-                continue;
-            }
-
-            double frac_x = contIndex[0] - index[0];
-            double frac_y = contIndex[1] - index[1];
-            double frac_z = contIndex[2] - index[2];
-
-            if (frac_x<0)
-            {
-                index[0] -= 1;
-                frac_x += 1;
-            }
-            if (frac_y<0)
-            {
-                index[1] -= 1;
-                frac_y += 1;
-            }
-            if (frac_z<0)
-            {
-                index[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 (index[0] < 0 || (unsigned long)index[0] >= outImageSize[0]-1)
-                continue;
-            if (index[1] < 0 || (unsigned long)index[1] >= outImageSize[1]-1)
-                continue;
-            if (index[2] < 0 || (unsigned long)index[2] >= outImageSize[2]-1)
+            // 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;
-            int idx;
-            wDir = dir;
-            double weight = (  frac_x)*(  frac_y)*(  frac_z);
-            if (weight>m_Thres)
-            {
-                wDir *= weight;
-                idx = index[0]   + outImageSize[0]*(index[1]  + outImageSize[1]*index[2]  );
-                dirCont = DirectionContainerType::New();
-                if (m_DirectionsContainer->IndexExists(idx))
-                {
-                    dirCont = m_DirectionsContainer->GetElement(idx);
-                    if (dirCont.IsNull())
-                    {
-                        dirCont = DirectionContainerType::New();
-                        dirCont->InsertElement(0, wDir);
-                        m_DirectionsContainer->InsertElement(idx, dirCont);
-                    }
-                    else
-                        dirCont->InsertElement(dirCont->Size(), wDir);
-                }
-                else
-                {
-                    dirCont->InsertElement(0, wDir);
-                    m_DirectionsContainer->InsertElement(idx, dirCont);
-                }
-            }
-
-            wDir = dir;
-            weight = (  frac_x)*(1-frac_y)*(  frac_z);
-            if (weight>m_Thres)
-            {
-                wDir *= weight;
-                idx = index[0]   + outImageSize[0]*(index[1]+1+ outImageSize[1]*index[2]  );
-                dirCont = DirectionContainerType::New();
-                if (m_DirectionsContainer->IndexExists(idx))
-                {
-                    dirCont = m_DirectionsContainer->GetElement(idx);
-                    if (dirCont.IsNull())
-                    {
-                        dirCont = DirectionContainerType::New();
-                        dirCont->InsertElement(0, wDir);
-                        m_DirectionsContainer->InsertElement(idx, dirCont);
-                    }
-                    else
-                        dirCont->InsertElement(dirCont->Size(), wDir);
-                }
-                else
-                {
-                    dirCont->InsertElement(0, wDir);
-                    m_DirectionsContainer->InsertElement(idx, dirCont);
-                }
-            }
-
-            wDir = dir;
-            weight = (  frac_x)*(  frac_y)*(1-frac_z);
-            if (weight>m_Thres)
+            if (m_DirectionsContainer->IndexExists(idx))
             {
-                wDir *= weight;
-                idx = index[0]   + outImageSize[0]*(index[1]  + outImageSize[1]*index[2]+outImageSize[1]);
-                dirCont = DirectionContainerType::New();
-                if (m_DirectionsContainer->IndexExists(idx))
-                {
-                    dirCont = m_DirectionsContainer->GetElement(idx);
-                    if (dirCont.IsNull())
-                    {
-                        dirCont = DirectionContainerType::New();
-                        dirCont->InsertElement(0, wDir);
-                        m_DirectionsContainer->InsertElement(idx, dirCont);
-                    }
-                    else
-                        dirCont->InsertElement(dirCont->Size(), wDir);
-                }
-                else
+                dirCont = m_DirectionsContainer->GetElement(idx);
+                if (dirCont.IsNull())
                 {
-                    dirCont->InsertElement(0, wDir);
+                    dirCont = DirectionContainerType::New();
+                    dirCont->push_back(dir);
                     m_DirectionsContainer->InsertElement(idx, dirCont);
                 }
-            }
-
-            wDir = dir;
-            weight = (  frac_x)*(1-frac_y)*(1-frac_z);
-            if (weight>m_Thres)
-            {
-                wDir *= weight;
-                idx = index[0]   + outImageSize[0]*(index[1]+1+ outImageSize[1]*index[2]+outImageSize[1]);
-                dirCont = DirectionContainerType::New();
-                if (m_DirectionsContainer->IndexExists(idx))
-                {
-                    dirCont = m_DirectionsContainer->GetElement(idx);
-                    if (dirCont.IsNull())
-                    {
-                        dirCont = DirectionContainerType::New();
-                        dirCont->InsertElement(0, wDir);
-                        m_DirectionsContainer->InsertElement(idx, dirCont);
-                    }
-                    else
-                        dirCont->InsertElement(dirCont->Size(), wDir);
-                }
                 else
-                {
-                    dirCont->InsertElement(0, wDir);
-                    m_DirectionsContainer->InsertElement(idx, dirCont);
-                }
+                    dirCont->push_back(dir);
             }
-
-            wDir = dir;
-            weight = (1-frac_x)*(  frac_y)*(  frac_z);
-            if (weight>m_Thres)
+            else
             {
-                wDir *= weight;
-                idx = index[0]+1 + outImageSize[0]*(index[1]  + outImageSize[1]*index[2]  );
                 dirCont = DirectionContainerType::New();
-                if (m_DirectionsContainer->IndexExists(idx))
-                {
-                    dirCont = m_DirectionsContainer->GetElement(idx);
-                    if (dirCont.IsNull())
-                    {
-                        dirCont = DirectionContainerType::New();
-                        dirCont->InsertElement(0, wDir);
-                        m_DirectionsContainer->InsertElement(idx, dirCont);
-                    }
-                    else
-                        dirCont->InsertElement(dirCont->Size(), wDir);
-                }
-                else
-                {
-                    dirCont->InsertElement(0, wDir);
-                    m_DirectionsContainer->InsertElement(idx, dirCont);
-                }
-            }
-
-            wDir = dir;
-            weight = (1-frac_x)*(  frac_y)*(1-frac_z);
-            if (weight>m_Thres)
-            {
-                wDir *= weight;
-                idx = index[0]+1 + outImageSize[0]*(index[1]  + outImageSize[1]*index[2]+outImageSize[1]);
-                dirCont = DirectionContainerType::New();
-                if (m_DirectionsContainer->IndexExists(idx))
-                {
-                    dirCont = m_DirectionsContainer->GetElement(idx);
-                    if (dirCont.IsNull())
-                    {
-                        dirCont = DirectionContainerType::New();
-                        dirCont->InsertElement(0, wDir);
-                        m_DirectionsContainer->InsertElement(idx, dirCont);
-                    }
-                    else
-                        dirCont->InsertElement(dirCont->Size(), wDir);
-                }
-                else
-                {
-                    dirCont->InsertElement(0, wDir);
-                    m_DirectionsContainer->InsertElement(idx, dirCont);
-                }
-            }
-
-            wDir = dir;
-            weight = (1-frac_x)*(1-frac_y)*(  frac_z);
-            if (weight>m_Thres)
-            {
-                wDir *= weight;
-                idx = index[0]+1 + outImageSize[0]*(index[1]+1+ outImageSize[1]*index[2]  );
-                dirCont = DirectionContainerType::New();
-                if (m_DirectionsContainer->IndexExists(idx))
-                {
-                    dirCont = m_DirectionsContainer->GetElement(idx);
-                    if (dirCont.IsNull())
-                    {
-                        dirCont = DirectionContainerType::New();
-                        dirCont->InsertElement(0, wDir);
-                        m_DirectionsContainer->InsertElement(idx, dirCont);
-                    }
-                    else
-                        dirCont->InsertElement(dirCont->Size(), wDir);
-                }
-                else
-                {
-                    dirCont->InsertElement(0, wDir);
-                    m_DirectionsContainer->InsertElement(idx, dirCont);
-                }
-            }
-
-            wDir = dir;
-            weight = (1-frac_x)*(1-frac_y)*(1-frac_z);
-            if (weight>m_Thres)
-            {
-                wDir *= weight;
-                idx = index[0]+1 + outImageSize[0]*(index[1]+1+ outImageSize[1]*index[2]+outImageSize[1]);
-                dirCont = DirectionContainerType::New();
-                if (m_DirectionsContainer->IndexExists(idx))
-                {
-                    dirCont = m_DirectionsContainer->GetElement(idx);
-                    if (dirCont.IsNull())
-                    {
-                        dirCont = DirectionContainerType::New();
-                        dirCont->InsertElement(0, wDir);
-                        m_DirectionsContainer->InsertElement(idx, dirCont);
-                    }
-                    else
-                        dirCont->InsertElement(dirCont->Size(), wDir);
-                }
-                else
-                {
-                    dirCont->InsertElement(0, wDir);
-                    m_DirectionsContainer->InsertElement(idx, dirCont);
-                }
+                dirCont->push_back(dir);
+                m_DirectionsContainer->InsertElement(idx, dirCont);
             }
         }
-        clock.Stop();
     }
 
     vtkSmartPointer<vtkCellArray> m_VtkCellArray = vtkSmartPointer<vtkCellArray>::New();
     vtkSmartPointer<vtkPoints>    m_VtkPoints = vtkSmartPointer<vtkPoints>::New();
 
     itk::ImageRegionIterator<ItkUcharImgType> dirIt(m_NumDirectionsImage, m_NumDirectionsImage->GetLargestPossibleRegion());
-    itk::ImageRegionIterator<ItkUcharImgType> crossIt(m_CrossingsImage, m_CrossingsImage->GetLargestPossibleRegion());
-
-    m_DirectionImageContainer = DirectionImageContainerType::New();
-    unsigned int maxNumDirections = 0;
 
     MITK_INFO << "Clustering directions";
     boost::progress_display disp2(outImageSize[0]*outImageSize[1]*outImageSize[2]);
-    for(crossIt.GoToBegin(); !crossIt.IsAtEnd(); ++crossIt)
+    while(!dirIt.IsAtEnd())
     {
         ++disp2;
-        OutputImageType::IndexType index = crossIt.GetIndex();
+        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() || index[0] < 0 || (unsigned long)index[0] >= outImageSize[0] || index[1] < 0 || (unsigned long)index[1] >= outImageSize[1] || index[2] < 0 || (unsigned long)index[2] >= outImageSize[2])
+        if (dirCont.IsNull() || dirCont->empty())
         {
             ++dirIt;
             continue;
         }
 
-        std::vector< DirectionType > directions;
-
-        for (unsigned int i=0; i<dirCont->Size(); i++)
-            if (dirCont->ElementAt(i).magnitude()>0.0001)
-                directions.push_back(dirCont->ElementAt(i));
-
-        if (!directions.empty())
-            directions = FastClustering(directions);
-
-        std::sort( directions.begin(), directions.end(), CompareVectorLengths );
-
-        if ( directions.size() > maxNumDirections )
-        {
-            for (unsigned int i=maxNumDirections; i<std::min<unsigned int>(directions.size(), m_MaxNumDirections); i++)
-            {
-                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);
-            }
-            maxNumDirections = std::min<unsigned int>(directions.size(), m_MaxNumDirections);
-        }
+        std::vector< double > lengths; lengths.resize(dirCont->size(), 1);  // all peaks have size 1
+        DirectionContainerType::Pointer directions = FastClustering(dirCont, lengths);
+        std::sort( directions->begin(), directions->end(), CompareVectorLengths );
 
-        unsigned int numDir = directions.size();
+        unsigned int numDir = directions->size();
         if (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);
+            DirectionType dir = directions->at(i);
+
+            if (dir.magnitude()<m_SizeThreshold)
+                continue;
+            count++;
 
             // 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(numDir);
+        dirIt.Set(count);
         ++dirIt;
     }
 
     vtkSmartPointer<vtkPolyData> directionsPolyData = vtkSmartPointer<vtkPolyData>::New();
     directionsPolyData->SetPoints(m_VtkPoints);
     directionsPolyData->SetLines(m_VtkCellArray);
     m_OutputFiberBundle = mitk::FiberBundleX::New(directionsPolyData);
 }
 
 
 template< class PixelType >
-std::vector< vnl_vector_fixed< double, 3 > > TractsToVectorImageFilter< PixelType >::FastClustering(std::vector< vnl_vector_fixed< double, 3 > >& inDirs)
+typename TractsToVectorImageFilter< PixelType >::DirectionContainerType::Pointer TractsToVectorImageFilter< PixelType >::FastClustering(DirectionContainerType::Pointer inDirs, std::vector< double > lengths)
 {
-    std::vector< vnl_vector_fixed< double, 3 > > outDirs;
-    if (inDirs.empty())
-        return outDirs;
-    vnl_vector_fixed< double, 3 > oldMean, currentMean, workingMean;
+    DirectionContainerType::Pointer outDirs = DirectionContainerType::New();
+    if (inDirs->size()<2)
+        return inDirs;
 
-    std::vector< vnl_vector_fixed< double, 3 > > normalizedDirs;
+    DirectionType oldMean, currentMean;
     std::vector< int > touched;
-    for (unsigned int i=0; i<inDirs.size(); i++)
-    {
-        normalizedDirs.push_back(inDirs[i]);
-        normalizedDirs.back().normalize();
-    }
 
     // initialize
-    double max = 0.0;
-    touched.resize(inDirs.size(), 0);
+    touched.resize(inDirs->size(), 0);
     bool free = true;
-    currentMean = inDirs[0];  // initialize first seed
+    currentMean = inDirs->at(0);  // initialize first seed
+    double length = lengths.at(0);
+    std::vector< double > newLengths;
+    bool meanChanged = false;
+
+    double max = 0;
     while (free)
     {
         oldMean.fill(0.0);
 
         // start mean-shift clustering
-        double angle = 0.0;
+        double angle = 0;
         int counter = 0;
-        while ((currentMean-oldMean).magnitude()>0.0001)
+
+        while (dot_product(currentMean, oldMean)<0.99)
         {
             counter = 0;
             oldMean = currentMean;
-            workingMean = oldMean;
-            workingMean.normalize();
             currentMean.fill(0.0);
-            for (unsigned int i=0; i<normalizedDirs.size(); i++)
+            for (unsigned int i=0; i<inDirs->size(); i++)
             {
-                angle = dot_product(workingMean, normalizedDirs[i]);
+                angle = dot_product(oldMean, inDirs->at(i));
                 if (angle>=m_AngularThreshold)
                 {
-                    currentMean += inDirs[i];
+                    currentMean += inDirs->at(i);
+                    if (meanChanged)
+                        length += lengths.at(i);
                     touched[i] = 1;
                     counter++;
+                    meanChanged = true;
                 }
                 else if (-angle>=m_AngularThreshold)
                 {
-                    currentMean -= inDirs[i];
+                    currentMean -= inDirs->at(i);
+                    if (meanChanged)
+                        length += lengths.at(i);
                     touched[i] = 1;
                     counter++;
+                    meanChanged = true;
                 }
             }
+            if (currentMean.magnitude()>0.001)
+                currentMean.normalize();
         }
 
         // found stable mean
         if (counter>0)
         {
-            currentMean /= counter;
-            double mag = currentMean.magnitude();
-
-            if (mag>0)
-            {
-                if (mag>max)
-                    max = mag;
-
-                outDirs.push_back(currentMean);
-            }
+            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[i];
+                currentMean = inDirs->at(i);
                 free = true;
+                meanChanged = false;
+                length = lengths.at(i);
+                break;
             }
     }
 
-    if (m_NormalizeVectors)
-        for (unsigned int i=0; i<outDirs.size(); i++)
-            outDirs[i].normalize();
-    else if (max>0)
-        for (unsigned int i=0; i<outDirs.size(); i++)
-            outDirs[i] /= max;
-
-    if (inDirs.size()==outDirs.size())
-        return outDirs;
-    else
-        return FastClustering(outDirs);
-}
-
-
-template< class PixelType >
-std::vector< vnl_vector_fixed< double, 3 > > TractsToVectorImageFilter< PixelType >::Clustering(std::vector< vnl_vector_fixed< double, 3 > >& inDirs)
-{
-    std::vector< vnl_vector_fixed< double, 3 > > outDirs;
-    if (inDirs.empty())
-        return outDirs;
-    vnl_vector_fixed< double, 3 > oldMean, currentMean, workingMean;
-
-    std::vector< vnl_vector_fixed< double, 3 > > 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++)
+    if (inDirs->size()==outDirs->size())
     {
-        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;
-            vnl_vector_fixed< double, 3 > normMean = currentMean;
-            normMean.normalize();
-            for (std::size_t i=0; i<outDirs.size(); i++)
-            {
-                vnl_vector_fixed< double, 3 > 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())
+        if (!m_NormalizeVectors && 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);
-}
-
-
-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;
+        return FastClustering(outDirs, newLengths);
 }
 
 
-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);
-}
+//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/FiberTracking/Algorithms/itkTractsToVectorImageFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.h
index 17c38b9d6e..21d931e14c 100644
--- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.h
+++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.h
@@ -1,109 +1,108 @@
 #ifndef __itkTractsToVectorImageFilter_h__
 #define __itkTractsToVectorImageFilter_h__
 
 // MITK
 #include <mitkFiberBundleX.h>
 
 // ITK
 #include <itkImageSource.h>
 #include <itkVectorImage.h>
 
 // VTK
 #include <vtkSmartPointer.h>
 #include <vtkPolyData.h>
 #include <vtkCellArray.h>
 #include <vtkPoints.h>
 #include <vtkPolyLine.h>
 
 using namespace mitk;
 
 namespace itk{
 
 /**
 * \brief Extracts the voxel-wise main directions of the input fiber bundle.   */
 
 template< class PixelType >
 class TractsToVectorImageFilter : public ImageSource< VectorImage< float, 3 > >
 {
 
 public:
     typedef TractsToVectorImageFilter Self;
     typedef ProcessObject Superclass;
     typedef SmartPointer< Self > Pointer;
     typedef SmartPointer< const Self > ConstPointer;
 
     typedef itk::Vector<float,3>                    OutputVectorType;
     typedef itk::Image<OutputVectorType, 3>         OutputImageType;
     typedef std::vector< OutputImageType::Pointer > OutputImageContainerType;
 
-    typedef vnl_vector_fixed< double, 3 >                                         DirectionType;
-    typedef VectorContainer< unsigned int, DirectionType >                       DirectionContainerType;
-    typedef VectorContainer< unsigned int, DirectionContainerType::Pointer >     ContainerType;
-    typedef Image< Vector< float, 3 >, 3>                                        ItkDirectionImageType;
-    typedef VectorContainer< unsigned int, ItkDirectionImageType::Pointer >      DirectionImageContainerType;
-
-    typedef itk::Image<unsigned char, 3>  ItkUcharImgType;
+    typedef vnl_vector_fixed< double, 3 >                                       DirectionType;
+    typedef VectorContainer< unsigned int, DirectionType >                      DirectionContainerType;
+    typedef VectorContainer< unsigned int, DirectionContainerType::Pointer >    ContainerType;
+    typedef Image< Vector< float, 3 >, 3>                                       ItkDirectionImageType;
+    typedef VectorContainer< unsigned int, ItkDirectionImageType::Pointer >     DirectionImageContainerType;
+    typedef itk::Image<unsigned char, 3>                                        ItkUcharImgType;
+    typedef itk::Image<double, 3>                                               ItkDoubleImgType;
 
     itkFactorylessNewMacro(Self)
     itkCloneMacro(Self)
     itkTypeMacro( TractsToVectorImageFilter, ImageSource )
 
+    itkSetMacro( SizeThreshold, float)
+    itkGetMacro( SizeThreshold, float)
     itkSetMacro( AngularThreshold, float)                               ///< cluster directions that are closer together than the specified threshold
     itkGetMacro( AngularThreshold, float)                               ///< cluster directions that are closer together than the specified threshold
     itkSetMacro( NormalizeVectors, bool)                                ///< Normalize vectors to length 1
     itkGetMacro( NormalizeVectors, bool)                                ///< Normalize vectors to length 1
     itkSetMacro( UseWorkingCopy, bool)                                  ///< Do not modify input fiber bundle. Use a copy.
     itkGetMacro( UseWorkingCopy, bool)                                  ///< Do not modify input fiber bundle. Use a copy.
-    itkSetMacro( MaxNumDirections, unsigned long)                                 ///< If more directions are extracted, only the largest are kept.
-    itkGetMacro( MaxNumDirections, unsigned long)                                 ///< If more directions are extracted, only the largest are kept.
+    itkSetMacro( MaxNumDirections, unsigned long)                       ///< If more directions are extracted, only the largest are kept.
+    itkGetMacro( MaxNumDirections, unsigned long)                       ///< If more directions are extracted, only the largest are kept.
     itkSetMacro( MaskImage, ItkUcharImgType::Pointer)                   ///< only process voxels inside mask
     itkSetMacro( FiberBundle, FiberBundleX::Pointer)                    ///< input fiber bundle
     itkGetMacro( ClusteredDirectionsContainer, ContainerType::Pointer)  ///< output directions
-    itkGetMacro( NumDirectionsImage, ItkUcharImgType::Pointer)          ///< nimber of directions per voxel
-    itkGetMacro( CrossingsImage, ItkUcharImgType::Pointer)              ///< mask voxels containing crossings
+    itkGetMacro( NumDirectionsImage, ItkUcharImgType::Pointer)          ///< number of directions per voxel
     itkGetMacro( OutputFiberBundle, FiberBundleX::Pointer)              ///< vector field for visualization purposes
     itkGetMacro( DirectionImageContainer, DirectionImageContainerType::Pointer) ///< output directions
 
     void GenerateData();
 
 protected:
 
-    std::vector< DirectionType > Clustering(std::vector< DirectionType >& inDirs);
-    std::vector< DirectionType > FastClustering(std::vector< DirectionType >& inDirs);  ///< cluster fiber directions
-    DirectionContainerType::Pointer MeanShiftClustering(DirectionContainerType::Pointer dirCont);
-    vnl_vector_fixed<double, 3> ClusterStep(DirectionContainerType::Pointer dirCont, vnl_vector_fixed<double, 3> currentMean);
+    DirectionContainerType::Pointer FastClustering(DirectionContainerType::Pointer inDirs, std::vector< double > lengths);  ///< cluster fiber directions
+//    std::vector< DirectionType > Clustering(std::vector< DirectionType >& inDirs);
+//    DirectionContainerType::Pointer MeanShiftClustering(DirectionContainerType::Pointer dirCont);
+//    vnl_vector_fixed<double, 3> ClusterStep(DirectionContainerType::Pointer dirCont, vnl_vector_fixed<double, 3> currentMean);
 
     vnl_vector_fixed<double, 3> GetVnlVector(double point[3]);
     itk::Point<double, 3> GetItkPoint(double point[3]);
 
 
     TractsToVectorImageFilter();
     virtual ~TractsToVectorImageFilter();
 
     FiberBundleX::Pointer               m_FiberBundle;                      ///< input fiber bundle
     float                               m_AngularThreshold;                 ///< cluster directions that are closer together than the specified threshold
     float                               m_Epsilon;                          ///< epsilon for vector equality check
     ItkUcharImgType::Pointer            m_MaskImage;                        ///< only voxels inside the binary mask are processed
     bool                                m_NormalizeVectors;                 ///< normalize vectors to length 1
-    itk::Vector<float>                 m_OutImageSpacing;                  ///< spacing of output image
+    itk::Vector<float>                  m_OutImageSpacing;                  ///< spacing of output image
     ContainerType::Pointer              m_DirectionsContainer;              ///< container for fiber directions
     bool                                m_UseWorkingCopy;                   ///< do not modify input fiber bundle but work on copy
-    bool                                m_UseTrilinearInterpolation;        ///< trilinearly interpolate between neighbouring voxels
-    unsigned long                                 m_MaxNumDirections;                 ///< if more directions per voxel are extracted, only the largest are kept
-    float                               m_Thres;                            ///< distance threshold for trilinear interpolation
+    unsigned long                       m_MaxNumDirections;                 ///< if more directions per voxel are extracted, only the largest are kept
+    float                               m_SizeThreshold;
 
     // output datastructures
     ContainerType::Pointer                  m_ClusteredDirectionsContainer; ///< contains direction vectors for each voxel
     ItkUcharImgType::Pointer                m_NumDirectionsImage;           ///< shows number of fibers per voxel
-    ItkUcharImgType::Pointer                m_CrossingsImage;               ///< shows voxels containing more than one fiber
     DirectionImageContainerType::Pointer    m_DirectionImageContainer;      ///< contains images that contain the output directions
     FiberBundleX::Pointer                   m_OutputFiberBundle;            ///< vector field for visualization purposes
 };
 
 }
 
 #ifndef ITK_MANUAL_INSTANTIATION
 #include "itkTractsToVectorImageFilter.cpp"
 #endif
 
 #endif // __itkTractsToVectorImageFilter_h__
diff --git a/Modules/DiffusionImaging/FiberTracking/Testing/mitkLocalFiberPlausibilityTest.cpp b/Modules/DiffusionImaging/FiberTracking/Testing/mitkLocalFiberPlausibilityTest.cpp
index 2869dd74ca..43e0ede842 100755
--- a/Modules/DiffusionImaging/FiberTracking/Testing/mitkLocalFiberPlausibilityTest.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/Testing/mitkLocalFiberPlausibilityTest.cpp
@@ -1,173 +1,174 @@
 /*===================================================================
 
 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 <mitkBaseDataIOFactory.h>
 #include <mitkBaseData.h>
 #include <mitkImageCast.h>
 #include <mitkImageToItk.h>
 #include <itkEvaluateDirectionImagesFilter.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 <mitkTestingMacros.h>
 #include <mitkCompareImageDataFilter.h>
 #include <mitkFiberBundleXWriter.h>
 
 #define _USE_MATH_DEFINES
 #include <math.h>
 
 using namespace std;
 
 int mitkLocalFiberPlausibilityTest(int argc, char* argv[])
 {
     MITK_TEST_BEGIN("mitkLocalFiberPlausibilityTest");
     MITK_TEST_CONDITION_REQUIRED(argc==8,"check for input data")
 
             string fibFile = argv[1];
     vector< string > referenceImages;
     referenceImages.push_back(argv[2]);
     referenceImages.push_back(argv[3]);
     string LDFP_ERROR_IMAGE = argv[4];
     string LDFP_NUM_DIRECTIONS = argv[5];
     string LDFP_VECTOR_FIELD = argv[6];
     string LDFP_ERROR_IMAGE_IGNORE = argv[7];
 
     float angularThreshold = 25;
 
     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::FiberBundleX::Pointer inputTractogram = dynamic_cast<mitk::FiberBundleX*>(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(...){ MITK_INFO << "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->SetSizeThreshold(0.0);
         fOdfFilter->SetUseWorkingCopy(false);
         fOdfFilter->SetNumberOfThreads(1);
         fOdfFilter->Update();
         ItkDirectionImageContainerType::Pointer directionImageContainer = fOdfFilter->GetDirectionImageContainer();
 
         // Get directions and num directions image
         ItkUcharImgType::Pointer numDirImage = fOdfFilter->GetNumDirectionsImage();
         mitk::Image::Pointer mitkNumDirImage = mitk::Image::New();
         mitkNumDirImage->InitializeByItk( numDirImage.GetPointer() );
         mitkNumDirImage->SetVolume( numDirImage->GetBufferPointer() );
         mitk::FiberBundleX::Pointer testDirections = fOdfFilter->GetOutputFiberBundle();
 
         // evaluate directions with missing directions
         EvaluationFilterType::Pointer evaluationFilter = EvaluationFilterType::New();
         evaluationFilter->SetImageSet(directionImageContainer);
         evaluationFilter->SetReferenceImageSet(referenceImageContainer);
         evaluationFilter->SetMaskImage(itkMaskImage);
         evaluationFilter->SetIgnoreMissingDirections(false);
         evaluationFilter->Update();
 
         EvaluationFilterType::OutputImageType::Pointer angularErrorImage = evaluationFilter->GetOutput(0);
         mitk::Image::Pointer mitkAngularErrorImage = mitk::Image::New();
         mitkAngularErrorImage->InitializeByItk( angularErrorImage.GetPointer() );
         mitkAngularErrorImage->SetVolume( angularErrorImage->GetBufferPointer() );
 
         // evaluate directions without missing directions
         evaluationFilter->SetIgnoreMissingDirections(true);
         evaluationFilter->Update();
 
         EvaluationFilterType::OutputImageType::Pointer angularErrorImageIgnore = evaluationFilter->GetOutput(0);
         mitk::Image::Pointer mitkAngularErrorImageIgnore = mitk::Image::New();
         mitkAngularErrorImageIgnore->InitializeByItk( angularErrorImageIgnore.GetPointer() );
         mitkAngularErrorImageIgnore->SetVolume( angularErrorImageIgnore->GetBufferPointer() );
 
         mitk::Image::Pointer gtAngularErrorImageIgnore = dynamic_cast<mitk::Image*>(mitk::IOUtil::LoadDataNode(LDFP_ERROR_IMAGE_IGNORE)->GetData());
         mitk::Image::Pointer gtAngularErrorImage = dynamic_cast<mitk::Image*>(mitk::IOUtil::LoadDataNode(LDFP_ERROR_IMAGE)->GetData());
         mitk::Image::Pointer gtNumTestDirImage = dynamic_cast<mitk::Image*>(mitk::IOUtil::LoadDataNode(LDFP_NUM_DIRECTIONS)->GetData());
         mitk::FiberBundleX::Pointer gtTestDirections = dynamic_cast<mitk::FiberBundleX*>(mitk::IOUtil::LoadDataNode(LDFP_VECTOR_FIELD)->GetData());
 
         if (!testDirections->Equals(gtTestDirections))
         {
             MITK_INFO << "SAVING FILES TO " << mitk::IOUtil::GetTempPath();
             std::string out1 = mitk::IOUtil::GetTempPath().append("test.fib");
             std::string out2 = mitk::IOUtil::GetTempPath().append("reference.fib");
 
             mitk::FiberBundleXWriter::Pointer fibWriter = mitk::FiberBundleXWriter::New();
             fibWriter->SetFileName(out1.c_str());
             fibWriter->DoWrite(testDirections.GetPointer());
 
             fibWriter->SetFileName(out2.c_str());
             fibWriter->DoWrite(gtTestDirections.GetPointer());
         }
 
         MITK_TEST_CONDITION_REQUIRED(mitk::Equal(gtAngularErrorImageIgnore, mitkAngularErrorImageIgnore, 0.01, true), "Check if error images are equal (ignored missing directions).");
         MITK_TEST_CONDITION_REQUIRED(mitk::Equal(gtAngularErrorImage, mitkAngularErrorImage, 0.01, true), "Check if error images are equal.");
         MITK_TEST_CONDITION_REQUIRED(testDirections->Equals(gtTestDirections), "Check if vector fields are equal.");
         MITK_TEST_CONDITION_REQUIRED(mitk::Equal(gtNumTestDirImage, mitkNumDirImage, 0.0001, true), "Check if num direction images are equal.");
     }
     catch (itk::ExceptionObject e)
     {
         MITK_INFO << e;
         return EXIT_FAILURE;
     }
     catch (std::exception e)
     {
         MITK_INFO << e.what();
         return EXIT_FAILURE;
     }
     catch (...)
     {
         MITK_INFO << "ERROR!?!";
         return EXIT_FAILURE;
     }
     MITK_TEST_END();
 }
diff --git a/Modules/DiffusionImaging/MiniApps/FiberDirectionExtraction.cpp b/Modules/DiffusionImaging/MiniApps/FiberDirectionExtraction.cpp
index c6c0c2a327..57adf8e206 100755
--- a/Modules/DiffusionImaging/MiniApps/FiberDirectionExtraction.cpp
+++ b/Modules/DiffusionImaging/MiniApps/FiberDirectionExtraction.cpp
@@ -1,163 +1,179 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 
 #include "MiniAppManager.h"
 #include <mitkBaseDataIOFactory.h>
 #include <mitkBaseData.h>
 #include <mitkImageCast.h>
 #include <mitkImageToItk.h>
 #include <metaCommand.h>
 #include "ctkCommandLineParser.h"
 #include <usAny.h>
 #include <itkImageFileWriter.h>
 #include <mitkIOUtil.h>
 #include <boost/lexical_cast.hpp>
 #include <itkEvaluateDirectionImagesFilter.h>
 #include <itkTractsToVectorImageFilter.h>
 #include <mitkCoreObjectFactory.h>
 
 #define _USE_MATH_DEFINES
 #include <math.h>
 
 int FiberDirectionExtraction(int argc, char* argv[])
 {
     ctkCommandLineParser parser;
     parser.setArgumentPrefix("--", "-");
     parser.addArgument("input", "i", ctkCommandLineParser::String, "input tractogram (.fib, vtk ascii file format)", us::Any(), false);
     parser.addArgument("out", "o", ctkCommandLineParser::String, "output root", us::Any(), false);
     parser.addArgument("mask", "m", ctkCommandLineParser::String, "mask image");
     parser.addArgument("athresh", "a", ctkCommandLineParser::Float, "angular threshold in degrees. closer fiber directions are regarded as one direction and clustered together.", 25, true);
+    parser.addArgument("peakthresh", "t", ctkCommandLineParser::Float, "peak size threshold relative to largest peak in voxel", 0.2, true);
     parser.addArgument("verbose", "v", ctkCommandLineParser::Bool, "output optional and intermediate calculation results");
+    parser.addArgument("numdirs", "d", ctkCommandLineParser::Int, "maximum number of fibers per voxel", 3, true);
+    parser.addArgument("normalize", "n", ctkCommandLineParser::Bool, "normalize vectors");
 
     map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
     if (parsedArgs.size()==0)
         return EXIT_FAILURE;
 
     string fibFile = us::any_cast<string>(parsedArgs["input"]);
 
     string maskImage("");
     if (parsedArgs.count("mask"))
         maskImage = us::any_cast<string>(parsedArgs["mask"]);
 
+    float peakThreshold = 0.2;
+    if (parsedArgs.count("peakthresh"))
+        peakThreshold = us::any_cast<float>(parsedArgs["peakthresh"]);
+
     float angularThreshold = 25;
     if (parsedArgs.count("athresh"))
         angularThreshold = us::any_cast<float>(parsedArgs["athresh"]);
 
     string outRoot = us::any_cast<string>(parsedArgs["out"]);
 
     bool verbose = false;
     if (parsedArgs.count("verbose"))
         verbose = us::any_cast<bool>(parsedArgs["verbose"]);
 
+    int maxNumDirs = 3;
+    if (parsedArgs.count("numdirs"))
+        maxNumDirs = us::any_cast<int>(parsedArgs["numdirs"]);
+
+    bool normalize = false;
+    if (parsedArgs.count("normalize"))
+        normalize = us::any_cast<bool>(parsedArgs["normalize"]);
 
     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;
 
         // load fiber bundle
         mitk::FiberBundleX::Pointer inputTractogram = dynamic_cast<mitk::FiberBundleX*>(mitk::IOUtil::LoadDataNode(fibFile)->GetData());
 
         // load/create mask image
         ItkUcharImgType::Pointer itkMaskImage = NULL;
         if (maskImage.compare("")!=0)
         {
             MITK_INFO << "Using mask image";
             itkMaskImage = ItkUcharImgType::New();
             mitk::Image::Pointer mitkMaskImage = dynamic_cast<mitk::Image*>(mitk::IOUtil::LoadDataNode(maskImage)->GetData());
             mitk::CastToItkImage<ItkUcharImgType>(mitkMaskImage, itkMaskImage);
         }
 
         // 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->SetNormalizeVectors(normalize);
         fOdfFilter->SetUseWorkingCopy(false);
+        fOdfFilter->SetSizeThreshold(peakThreshold);
+        fOdfFilter->SetMaxNumDirections(maxNumDirs);
         fOdfFilter->Update();
         ItkDirectionImageContainerType::Pointer directionImageContainer = fOdfFilter->GetDirectionImageContainer();
 
         // 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");
-
-            MITK_INFO << "writing " << outfilename;
-            writer->SetFileName(outfilename.c_str());
-            writer->SetInput(itkImg);
-            writer->Update();
-        }
+//        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");
+
+//            MITK_INFO << "writing " << outfilename;
+//            writer->SetFileName(outfilename.c_str());
+//            writer->SetInput(itkImg);
+//            writer->Update();
+//        }
 
         if (verbose)
         {
             // write vector field
             mitk::FiberBundleX::Pointer directions = fOdfFilter->GetOutputFiberBundle();
             mitk::CoreObjectFactory::FileWriterList fileWriters = mitk::CoreObjectFactory::GetInstance()->GetFileWriters();
             for (mitk::CoreObjectFactory::FileWriterList::iterator it = fileWriters.begin() ; it != fileWriters.end() ; ++it)
             {
                 if ( (*it)->CanWriteBaseDataType(directions.GetPointer()) ) {
                     string outfilename = outRoot;
                     outfilename.append("_VECTOR_FIELD.fib");
                     (*it)->SetFileName( outfilename.c_str() );
                     (*it)->DoWrite( directions.GetPointer() );
                 }
             }
 
             // 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");
 
                 MITK_INFO << "writing " << outfilename;
                 writer->SetFileName(outfilename.c_str());
                 writer->SetInput(numDirImage);
                 writer->Update();
             }
         }
 
         MITK_INFO << "DONE";
     }
     catch (itk::ExceptionObject e)
     {
         MITK_INFO << e;
         return EXIT_FAILURE;
     }
     catch (std::exception e)
     {
         MITK_INFO << e.what();
         return EXIT_FAILURE;
     }
     catch (...)
     {
         MITK_INFO << "ERROR!?!";
         return EXIT_FAILURE;
     }
     return EXIT_SUCCESS;
 }
 RegisterDiffusionMiniApp(FiberDirectionExtraction);