diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp
index 925392df9a..6d1344e482 100644
--- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp
@@ -1,883 +1,890 @@
 /*===================================================================
 
 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 __itkStreamlineTrackingFilter_txx
 #define __itkStreamlineTrackingFilter_txx
 
 #include <time.h>
 #include <stdio.h>
 #include <stdlib.h>
 
 #include "itkStreamlineTrackingFilter.h"
 #include <itkImageRegionConstIterator.h>
 #include <itkImageRegionConstIteratorWithIndex.h>
 #include <itkImageRegionIterator.h>
 
 #define _USE_MATH_DEFINES
 #include <math.h>
 
 namespace itk {
 
-//#define QBALL_RECON_PI       M_PI
-
 template< class TTensorPixelType, class TPDPixelType>
 StreamlineTrackingFilter< TTensorPixelType,
 TPDPixelType>
-::StreamlineTrackingFilter():
-    m_FaThreshold(0.2),
-    m_StepSize(1),
-    m_MaxLength(10000),
-    m_MinTractLength(0.0),
-    m_SeedsPerVoxel(1),
-    m_F(1.0),
-    m_G(0.0),
-    m_Interpolate(true),
-    m_ResampleFibers(false)
+::StreamlineTrackingFilter()
+    : m_FiberPolyData(NULL)
+    , m_Points(NULL)
+    , m_Cells(NULL)
+    , m_FaImage(NULL)
+    , m_NumberOfInputs(1)
+    , m_FaThreshold(0.2)
+    , m_MinCurvatureRadius(0)
+    , m_StepSize(0)
+    , m_MaxLength(10000)
+    , m_MinTractLength(0.0)
+    , m_SeedsPerVoxel(1)
+    , m_F(1.0)
+    , m_G(0.0)
+    , m_Interpolate(true)
+    , m_PointPistance(0.0)
+    , m_ResampleFibers(false)
+    , m_SeedImage(NULL)
+    , m_MaskImage(NULL)
 {
     // At least 1 inputs is necessary for a vector image.
     // For images added one at a time we need at least six
     this->SetNumberOfRequiredInputs( 1 );
     this->SetNumberOfIndexedInputs(3);
 }
 
 template< class TTensorPixelType,
           class TPDPixelType>
 double StreamlineTrackingFilter< TTensorPixelType,
 TPDPixelType>
 ::RoundToNearest(double num) {
     return (num > 0.0) ? floor(num + 0.5) : ceil(num - 0.5);
 }
 
 template< class TTensorPixelType,
           class TPDPixelType>
 void StreamlineTrackingFilter< TTensorPixelType,
 TPDPixelType>
 ::BeforeThreadedGenerateData()
 {
     m_FiberPolyData = FiberPolyDataType::New();
-    m_Points = vtkPoints::New();
-    m_Cells = vtkCellArray::New();
+    m_Points = vtkSmartPointer< vtkPoints >::New();
+    m_Cells = vtkSmartPointer< vtkCellArray >::New();
 
     InputImageType* inputImage = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) );
     m_ImageSize.resize(3);
     m_ImageSize[0] = inputImage->GetLargestPossibleRegion().GetSize()[0];
     m_ImageSize[1] = inputImage->GetLargestPossibleRegion().GetSize()[1];
     m_ImageSize[2] = inputImage->GetLargestPossibleRegion().GetSize()[2];
 
     if (m_ImageSize[0]<3 || m_ImageSize[1]<3 || m_ImageSize[2]<3)
         m_Interpolate = false;
 
     m_ImageSpacing.resize(3);
     m_ImageSpacing[0] = inputImage->GetSpacing()[0];
     m_ImageSpacing[1] = inputImage->GetSpacing()[1];
     m_ImageSpacing[2] = inputImage->GetSpacing()[2];
 
     float 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];
     if (m_StepSize<0.1*minSpacing)
         m_StepSize = 0.1*minSpacing;
 
     if (m_ResampleFibers)
         m_PointPistance = 0.5*minSpacing;
 
     m_PolyDataContainer = itk::VectorContainer< int, FiberPolyDataType >::New();
     for (unsigned int i=0; i<this->GetNumberOfThreads(); i++)
     {
         FiberPolyDataType poly = FiberPolyDataType::New();
         m_PolyDataContainer->InsertElement(i, poly);
     }
 
     if (m_SeedImage.IsNull())
     {
         // initialize mask image
         m_SeedImage = ItkUcharImgType::New();
         m_SeedImage->SetSpacing( inputImage->GetSpacing() );
         m_SeedImage->SetOrigin( inputImage->GetOrigin() );
         m_SeedImage->SetDirection( inputImage->GetDirection() );
         m_SeedImage->SetRegions( inputImage->GetLargestPossibleRegion() );
         m_SeedImage->Allocate();
         m_SeedImage->FillBuffer(1);
     }
 
     if (m_MaskImage.IsNull())
     {
         // initialize mask image
         m_MaskImage = ItkUcharImgType::New();
         m_MaskImage->SetSpacing( inputImage->GetSpacing() );
         m_MaskImage->SetOrigin( inputImage->GetOrigin() );
         m_MaskImage->SetDirection( inputImage->GetDirection() );
         m_MaskImage->SetRegions( inputImage->GetLargestPossibleRegion() );
         m_MaskImage->Allocate();
         m_MaskImage->FillBuffer(1);
     }
 
     bool useUserFaImage = true;
     if (m_FaImage.IsNull())
     {
         m_FaImage = ItkFloatImgType::New();
         m_FaImage->SetSpacing( inputImage->GetSpacing() );
         m_FaImage->SetOrigin( inputImage->GetOrigin() );
         m_FaImage->SetDirection( inputImage->GetDirection() );
         m_FaImage->SetRegions( inputImage->GetLargestPossibleRegion() );
         m_FaImage->Allocate();
         m_FaImage->FillBuffer(0.0);
         useUserFaImage = false;
     }
 
     m_NumberOfInputs = 0;
     for (unsigned int i=0; i<this->GetNumberOfIndexedInputs(); i++)
     {
         if (this->ProcessObject::GetInput(i)==NULL)
             break;
 
         ItkPDImgType::Pointer pdImage = ItkPDImgType::New();
         pdImage->SetSpacing( inputImage->GetSpacing() );
         pdImage->SetOrigin( inputImage->GetOrigin() );
         pdImage->SetDirection( inputImage->GetDirection() );
         pdImage->SetRegions( inputImage->GetLargestPossibleRegion() );
         pdImage->Allocate();
         m_PdImage.push_back(pdImage);
 
         ItkFloatImgType::Pointer emaxImage = ItkFloatImgType::New();
         emaxImage->SetSpacing( inputImage->GetSpacing() );
         emaxImage->SetOrigin( inputImage->GetOrigin() );
         emaxImage->SetDirection( inputImage->GetDirection() );
         emaxImage->SetRegions( inputImage->GetLargestPossibleRegion() );
         emaxImage->Allocate();
         emaxImage->FillBuffer(1.0);
         m_EmaxImage.push_back(emaxImage);
 
         typename InputImageType::Pointer inputImage = static_cast< InputImageType * >( this->ProcessObject::GetInput(i) );
         m_InputImage.push_back(inputImage);
 
         m_NumberOfInputs++;
     }
     MITK_INFO << "Processing " << m_NumberOfInputs << " tensor files";
 
     typedef itk::DiffusionTensor3D<TTensorPixelType>    TensorType;
     typename TensorType::EigenValuesArrayType eigenvalues;
     typename TensorType::EigenVectorsMatrixType eigenvectors;
 
 
     for (int x=0; x<m_ImageSize[0]; x++)
         for (int y=0; y<m_ImageSize[1]; y++)
             for (int z=0; z<m_ImageSize[2]; z++)
             {
                 typename InputImageType::IndexType index;
                 index[0] = x; index[1] = y; index[2] = z;
                 for (int i=0; i<m_NumberOfInputs; i++)
                 {
                     typename InputImageType::PixelType tensor = m_InputImage.at(i)->GetPixel(index);
                     vnl_vector_fixed<double,3> dir;
                     tensor.ComputeEigenAnalysis(eigenvalues, eigenvectors);
                     dir[0] = eigenvectors(2, 0);
                     dir[1] = eigenvectors(2, 1);
                     dir[2] = eigenvectors(2, 2);
                     dir.normalize();
                     m_PdImage.at(i)->SetPixel(index, dir);
                     if (!useUserFaImage)
                         m_FaImage->SetPixel(index, m_FaImage->GetPixel(index)+tensor.GetFractionalAnisotropy());
                     m_EmaxImage.at(i)->SetPixel(index, 2/eigenvalues[2]);
                 }
                 if (!useUserFaImage)
                     m_FaImage->SetPixel(index, m_FaImage->GetPixel(index)/m_NumberOfInputs);
             }
 
     if (m_Interpolate)
         std::cout << "StreamlineTrackingFilter: using trilinear interpolation" << std::endl;
     else
     {
         if (m_MinCurvatureRadius<0.0)
             m_MinCurvatureRadius = 0.1*minSpacing;
 
         std::cout << "StreamlineTrackingFilter: using nearest neighbor interpolation" << std::endl;
     }
 
     if (m_MinCurvatureRadius<0.0)
         m_MinCurvatureRadius = 0.5*minSpacing;
 
     std::cout << "StreamlineTrackingFilter: Min. curvature radius: " << m_MinCurvatureRadius << std::endl;
     std::cout << "StreamlineTrackingFilter: FA threshold: " << m_FaThreshold << std::endl;
     std::cout << "StreamlineTrackingFilter: stepsize: " << m_StepSize << " mm" << std::endl;
     std::cout << "StreamlineTrackingFilter: f: " << m_F << std::endl;
     std::cout << "StreamlineTrackingFilter: g: " << m_G << std::endl;
     std::cout << "StreamlineTrackingFilter: starting streamline tracking using " << this->GetNumberOfThreads() << " threads." << std::endl;
 }
 
 template< class TTensorPixelType, class TPDPixelType>
 void StreamlineTrackingFilter< TTensorPixelType, TPDPixelType>
 ::CalculateNewPosition(itk::ContinuousIndex<double, 3>& pos, vnl_vector_fixed<double,3>& dir, typename InputImageType::IndexType& index)
 {
     vnl_matrix_fixed< double, 3, 3 > rot = m_InputImage.at(0)->GetDirection().GetTranspose();
     dir = rot*dir;
     if (true)
     {
         dir *= m_StepSize;
         pos[0] += dir[0]/m_ImageSpacing[0];
         pos[1] += dir[1]/m_ImageSpacing[1];
         pos[2] += dir[2]/m_ImageSpacing[2];
         index[0] = RoundToNearest(pos[0]);
         index[1] = RoundToNearest(pos[1]);
         index[2] = RoundToNearest(pos[2]);
     }
     else
     {
         dir[0] /= m_ImageSpacing[0];
         dir[1] /= m_ImageSpacing[1];
         dir[2] /= m_ImageSpacing[2];
 
         int smallest = 0;
         float x = 100000;
         if (dir[0]>0)
         {
             if (fabs(fabs(RoundToNearest(pos[0])-pos[0])-0.5)>mitk::eps)
                 x = fabs(pos[0]-RoundToNearest(pos[0])-0.5)/dir[0];
             else
                 x = fabs(pos[0]-std::ceil(pos[0])-0.5)/dir[0];
         }
         else if (dir[0]<0)
         {
             if (fabs(fabs(RoundToNearest(pos[0])-pos[0])-0.5)>mitk::eps)
                 x = -fabs(pos[0]-RoundToNearest(pos[0])+0.5)/dir[0];
             else
                 x = -fabs(pos[0]-std::floor(pos[0])+0.5)/dir[0];
         }
         float s = x;
 
         float y = 100000;
         if (dir[1]>0)
         {
             if (fabs(fabs(RoundToNearest(pos[1])-pos[1])-0.5)>mitk::eps)
                 y = fabs(pos[1]-RoundToNearest(pos[1])-0.5)/dir[1];
             else
                 y = fabs(pos[1]-std::ceil(pos[1])-0.5)/dir[1];
         }
         else if (dir[1]<0)
         {
             if (fabs(fabs(RoundToNearest(pos[1])-pos[1])-0.5)>mitk::eps)
                 y = -fabs(pos[1]-RoundToNearest(pos[1])+0.5)/dir[1];
             else
                 y = -fabs(pos[1]-std::floor(pos[1])+0.5)/dir[1];
         }
         if (s>y)
         {
             s=y;
             smallest = 1;
         }
 
         float z = 100000;
         if (dir[2]>0)
         {
             if (fabs(fabs(RoundToNearest(pos[2])-pos[2])-0.5)>mitk::eps)
                 z = fabs(pos[2]-RoundToNearest(pos[2])-0.5)/dir[2];
             else
                 z = fabs(pos[2]-std::ceil(pos[2])-0.5)/dir[2];
         }
         else if (dir[2]<0)
         {
             if (fabs(fabs(RoundToNearest(pos[2])-pos[2])-0.5)>mitk::eps)
                 z = -fabs(pos[2]-RoundToNearest(pos[2])+0.5)/dir[2];
             else
                 z = -fabs(pos[2]-std::floor(pos[2])+0.5)/dir[2];
         }
         if (s>z)
         {
             s=z;
             smallest = 2;
         }
 
         //        MITK_INFO << "---------------------------------------------";
         //        MITK_INFO << "s: " << s;
         //        MITK_INFO << "dir: " << dir;
         //        MITK_INFO << "old: " << pos[0] << ", " << pos[1] << ", " << pos[2];
 
         pos[0] += dir[0]*s;
         pos[1] += dir[1]*s;
         pos[2] += dir[2]*s;
 
         switch (smallest)
         {
         case 0:
             if (dir[0]<0)
                 index[0] = std::floor(pos[0]);
             else
                 index[0] = std::ceil(pos[0]);
             index[1] = RoundToNearest(pos[1]);
             index[2] = RoundToNearest(pos[2]);
             break;
 
         case 1:
             if (dir[1]<0)
                 index[1] = std::floor(pos[1]);
             else
                 index[1] = std::ceil(pos[1]);
             index[0] = RoundToNearest(pos[0]);
             index[2] = RoundToNearest(pos[2]);
             break;
 
         case 2:
             if (dir[2]<0)
                 index[2] = std::floor(pos[2]);
             else
                 index[2] = std::ceil(pos[2]);
             index[1] = RoundToNearest(pos[1]);
             index[0] = RoundToNearest(pos[0]);
         }
 
         //        float x = 100000;
         //        if (dir[0]>0)
         //            x = fabs(pos[0]-RoundToNearest(pos[0])-0.5)/dir[0];
         //        else if (dir[0]<0)
         //            x = -fabs(pos[0]-RoundToNearest(pos[0])+0.5)/dir[0];
         //        float s = x;
 
         //        float y = 100000;
         //        if (dir[1]>0)
         //            y = fabs(pos[1]-RoundToNearest(pos[1])-0.5)/dir[1];
         //        else if (dir[1]<0)
         //            y = -fabs(pos[1]-RoundToNearest(pos[1])+0.5)/dir[1];
         //        if (s>y)
         //            s=y;
 
         //        float z = 100000;
         //        if (dir[2]>0)
         //            z = fabs(pos[2]-RoundToNearest(pos[2])-0.5)/dir[2];
         //        else if (dir[2]<0)
         //            z = -fabs(pos[2]-RoundToNearest(pos[2])+0.5)/dir[2];
 
         //        if (s>z)
         //            s=z;
         //        s *= 1.001;
 
         //        pos[0] += dir[0]*s;
         //        pos[1] += dir[1]*s;
         //        pos[2] += dir[2]*s;
 
         //        index[0] = RoundToNearest(pos[0]);
         //        index[1] = RoundToNearest(pos[1]);
         //        index[2] = RoundToNearest(pos[2]);
 
         //        MITK_INFO << "new: " << pos[0] << ", " << pos[1] << ", " << pos[2];
     }
 }
 
 template< class TTensorPixelType, class TPDPixelType>
 bool StreamlineTrackingFilter< TTensorPixelType, TPDPixelType>
 ::IsValidPosition(itk::ContinuousIndex<double, 3>& pos, typename InputImageType::IndexType &index, vnl_vector_fixed< float, 8 >& interpWeights, int imageIdx)
 {
     if (!m_InputImage.at(imageIdx)->GetLargestPossibleRegion().IsInside(index) || m_MaskImage->GetPixel(index)==0)
         return false;
 
     if (m_Interpolate)
     {
         float frac_x = pos[0] - index[0];
         float frac_y = pos[1] - index[1];
         float frac_z = pos[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 || index[0] >= m_ImageSize[0]-1)
             return false;
         if (index[1] < 0 || index[1] >= m_ImageSize[1]-1)
             return false;
         if (index[2] < 0 || index[2] >= m_ImageSize[2]-1)
             return false;
 
         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);
 
         typename InputImageType::IndexType tmpIdx;
         float FA = m_FaImage->GetPixel(index) * interpWeights[0];
         tmpIdx = index; tmpIdx[0]++;
         FA +=  m_FaImage->GetPixel(tmpIdx) * interpWeights[1];
         tmpIdx = index; tmpIdx[1]++;
         FA +=  m_FaImage->GetPixel(tmpIdx) * interpWeights[2];
         tmpIdx = index; tmpIdx[2]++;
         FA +=  m_FaImage->GetPixel(tmpIdx) * interpWeights[3];
         tmpIdx = index; tmpIdx[0]++; tmpIdx[1]++;
         FA +=  m_FaImage->GetPixel(tmpIdx) * interpWeights[4];
         tmpIdx = index; tmpIdx[1]++; tmpIdx[2]++;
         FA +=  m_FaImage->GetPixel(tmpIdx) * interpWeights[5];
         tmpIdx = index; tmpIdx[2]++; tmpIdx[0]++;
         FA +=  m_FaImage->GetPixel(tmpIdx) * interpWeights[6];
         tmpIdx = index; tmpIdx[0]++; tmpIdx[1]++; tmpIdx[2]++;
         FA +=  m_FaImage->GetPixel(tmpIdx) * interpWeights[7];
 
         if (FA<m_FaThreshold)
             return false;
     }
     else if (m_FaImage->GetPixel(index)<m_FaThreshold)
         return false;
 
     return true;
 }
 
 template< class TTensorPixelType, class TPDPixelType>
 float StreamlineTrackingFilter< TTensorPixelType, TPDPixelType>
 ::FollowStreamline(itk::ContinuousIndex<double, 3> pos, int dirSign, vtkPoints* points, std::vector< vtkIdType >& ids, int imageIdx)
 {
     float tractLength = 0;
     typedef itk::DiffusionTensor3D<TTensorPixelType>    TensorType;
     typename TensorType::EigenValuesArrayType eigenvalues;
     typename TensorType::EigenVectorsMatrixType eigenvectors;
     vnl_vector_fixed< float, 8 > interpWeights;
 
     typename InputImageType::IndexType index, indexOld;
     indexOld[0] = -1; indexOld[1] = -1; indexOld[2] = -1;
     itk::Point<double> worldPos;
     float distance = 0;
     float distanceInVoxel = 0;
 
     // starting index and direction
     index[0] = RoundToNearest(pos[0]);
     index[1] = RoundToNearest(pos[1]);
     index[2] = RoundToNearest(pos[2]);
 
     vnl_vector_fixed<double,3> dir = m_PdImage.at(imageIdx)->GetPixel(index);
     dir *= dirSign;                     // reverse direction
     vnl_vector_fixed<double,3> dirOld = dir;
     if (dir.magnitude()<mitk::eps)
         return tractLength;
 
     for (int step=0; step< m_MaxLength/2; step++)
     {
         // get new position
         CalculateNewPosition(pos, dir, index);
         distance += m_StepSize;
         tractLength +=  m_StepSize;
         distanceInVoxel += m_StepSize;
 
         // is new position valid (inside image, above FA threshold etc.)
         if (!IsValidPosition(pos, index, interpWeights, imageIdx))   // if not add last point and end streamline
         {
             m_SeedImage->TransformContinuousIndexToPhysicalPoint( pos, worldPos );
             ids.push_back(points->InsertNextPoint(worldPos.GetDataPointer()));
             return tractLength;
         }
         else if (distance>=m_PointPistance)
         {
             m_SeedImage->TransformContinuousIndexToPhysicalPoint( pos, worldPos );
             ids.push_back(points->InsertNextPoint(worldPos.GetDataPointer()));
             distance = 0;
         }
 
         if (!m_Interpolate)                         // use nearest neighbour interpolation
         {
             if (indexOld!=index)                    // did we enter a new voxel? if yes, calculate new direction
             {
                 double minAngle = 0;
                 for (int img=0; img<m_NumberOfInputs; img++)
                 {
                     vnl_vector_fixed<double,3> newDir = m_PdImage.at(img)->GetPixel(index);   // get principal direction
                     if (newDir.magnitude()<mitk::eps)
                         continue;
 
                     typename InputImageType::PixelType tensor = m_InputImage.at(img)->GetPixel(index);
                     float scale = m_EmaxImage.at(img)->GetPixel(index);
                     newDir[0] = m_F*newDir[0] + (1-m_F)*( (1-m_G)*dirOld[0] + scale*m_G*(tensor[0]*dirOld[0] + tensor[1]*dirOld[1] + tensor[2]*dirOld[2]));
                     newDir[1] = m_F*newDir[1] + (1-m_F)*( (1-m_G)*dirOld[1] + scale*m_G*(tensor[1]*dirOld[0] + tensor[3]*dirOld[1] + tensor[4]*dirOld[2]));
                     newDir[2] = m_F*newDir[2] + (1-m_F)*( (1-m_G)*dirOld[2] + scale*m_G*(tensor[2]*dirOld[0] + tensor[4]*dirOld[1] + tensor[5]*dirOld[2]));
                     newDir.normalize();
 
                     float angle = dot_product(dirOld, newDir);
                     if (angle<0)
                     {
                         newDir *= -1;
                         angle *= -1;
                     }
 
                     if (angle>minAngle)
                     {
                         minAngle = angle;
                         dir = newDir;
                     }
                 }
 
                 //float r = m_StepSize/(2*std::asin(std::acos(minAngle)/2));
                 vnl_vector_fixed<double,3> v3 = dir+dirOld; v3 *= m_StepSize;
                 float a = m_StepSize;
                 float b = m_StepSize;
                 float c = v3.magnitude();
                 float r = a*b*c/std::sqrt((a+b+c)*(a+b-c)*(b+c-a)*(a-b+c)); // radius of triangle via Heron's formula (area of triangle)
 
                 if (r<m_MinCurvatureRadius)
                 {
                     return tractLength;
                 }
                 dirOld = dir;
                 indexOld = index;
                 distanceInVoxel = 0;
             }
             else
                 dir = dirOld;
         }
         else // use trilinear interpolation (weights calculated in IsValidPosition())
         {
             typename InputImageType::PixelType tensor;
 
             typename InputImageType::IndexType tmpIdx = index;
             typename InputImageType::PixelType tmpTensor;
 
             if (m_NumberOfInputs>1)
             {
                 double minAngle = 0;
                 for (int img=0; img<m_NumberOfInputs; img++)
                 {
                     float angle = dot_product(dirOld, m_PdImage.at(img)->GetPixel(tmpIdx));
                     if (fabs(angle)>minAngle)
                     {
                         minAngle = angle;
                         tmpTensor = m_InputImage.at(img)->GetPixel(tmpIdx);
                     }
                 }
                 tensor = tmpTensor * interpWeights[0];
 
                 minAngle = 0;
                 tmpIdx = index; tmpIdx[0]++;
                 for (int img=0; img<m_NumberOfInputs; img++)
                 {
                     float angle = dot_product(dirOld, m_PdImage.at(img)->GetPixel(tmpIdx));
                     if (fabs(angle)>minAngle)
                     {
                         minAngle = angle;
                         tmpTensor = m_InputImage.at(img)->GetPixel(tmpIdx);
                     }
                 }
                 tensor += tmpTensor * interpWeights[1];
 
                 minAngle = 0;
                 tmpIdx = index; tmpIdx[1]++;
                 for (int img=0; img<m_NumberOfInputs; img++)
                 {
                     float angle = dot_product(dirOld, m_PdImage.at(img)->GetPixel(tmpIdx));
                     if (fabs(angle)>minAngle)
                     {
                         minAngle = angle;
                         tmpTensor = m_InputImage.at(img)->GetPixel(tmpIdx);
                     }
                 }
                 tensor += tmpTensor * interpWeights[2];
 
                 minAngle = 0;
                 tmpIdx = index; tmpIdx[2]++;
                 for (int img=0; img<m_NumberOfInputs; img++)
                 {
                     float angle = dot_product(dirOld, m_PdImage.at(img)->GetPixel(tmpIdx));
                     if (fabs(angle)>minAngle)
                     {
                         minAngle = angle;
                         tmpTensor = m_InputImage.at(img)->GetPixel(tmpIdx);
                     }
                 }
                 tensor += tmpTensor * interpWeights[3];
 
                 minAngle = 0;
                 tmpIdx = index; tmpIdx[0]++; tmpIdx[1]++;
                 for (int img=0; img<m_NumberOfInputs; img++)
                 {
                     float angle = dot_product(dirOld, m_PdImage.at(img)->GetPixel(tmpIdx));
                     if (fabs(angle)>minAngle)
                     {
                         minAngle = angle;
                         tmpTensor = m_InputImage.at(img)->GetPixel(tmpIdx);
                     }
                 }
                 tensor += tmpTensor * interpWeights[4];
 
                 minAngle = 0;
                 tmpIdx = index; tmpIdx[1]++; tmpIdx[2]++;
                 for (int img=0; img<m_NumberOfInputs; img++)
                 {
                     float angle = dot_product(dirOld, m_PdImage.at(img)->GetPixel(tmpIdx));
                     if (fabs(angle)>minAngle)
                     {
                         minAngle = angle;
                         tmpTensor = m_InputImage.at(img)->GetPixel(tmpIdx);
                     }
                 }
                 tensor += tmpTensor * interpWeights[5];
 
                 minAngle = 0;
                 tmpIdx = index; tmpIdx[2]++; tmpIdx[0]++;
                 for (int img=0; img<m_NumberOfInputs; img++)
                 {
                     float angle = dot_product(dirOld, m_PdImage.at(img)->GetPixel(tmpIdx));
                     if (fabs(angle)>minAngle)
                     {
                         minAngle = angle;
                         tmpTensor = m_InputImage.at(img)->GetPixel(tmpIdx);
                     }
                 }
                 tensor += tmpTensor * interpWeights[6];
 
                 minAngle = 0;
                 tmpIdx = index; tmpIdx[0]++; tmpIdx[1]++; tmpIdx[2]++;
                 for (int img=0; img<m_NumberOfInputs; img++)
                 {
                     float angle = dot_product(dirOld, m_PdImage.at(img)->GetPixel(tmpIdx));
                     if (fabs(angle)>minAngle)
                     {
                         minAngle = angle;
                         tmpTensor = m_InputImage.at(img)->GetPixel(tmpIdx);
                     }
                 }
                 tensor += tmpTensor * interpWeights[7];
             }
             else
             {
                 tensor = m_InputImage.at(0)->GetPixel(index) * interpWeights[0];
                 typename InputImageType::IndexType tmpIdx = index; tmpIdx[0]++;
                 tensor +=  m_InputImage.at(0)->GetPixel(tmpIdx) * interpWeights[1];
                 tmpIdx = index; tmpIdx[1]++;
                 tensor +=  m_InputImage.at(0)->GetPixel(tmpIdx) * interpWeights[2];
                 tmpIdx = index; tmpIdx[2]++;
                 tensor +=  m_InputImage.at(0)->GetPixel(tmpIdx) * interpWeights[3];
                 tmpIdx = index; tmpIdx[0]++; tmpIdx[1]++;
                 tensor +=  m_InputImage.at(0)->GetPixel(tmpIdx) * interpWeights[4];
                 tmpIdx = index; tmpIdx[1]++; tmpIdx[2]++;
                 tensor +=  m_InputImage.at(0)->GetPixel(tmpIdx) * interpWeights[5];
                 tmpIdx = index; tmpIdx[2]++; tmpIdx[0]++;
                 tensor +=  m_InputImage.at(0)->GetPixel(tmpIdx) * interpWeights[6];
                 tmpIdx = index; tmpIdx[0]++; tmpIdx[1]++; tmpIdx[2]++;
                 tensor +=  m_InputImage.at(0)->GetPixel(tmpIdx) * interpWeights[7];
             }
 
             tensor.ComputeEigenAnalysis(eigenvalues, eigenvectors);
             dir[0] = eigenvectors(2, 0);
             dir[1] = eigenvectors(2, 1);
             dir[2] = eigenvectors(2, 2);
             if (dir.magnitude()<mitk::eps)
                 continue;
             dir.normalize();
 
             float scale = 2/eigenvalues[2];
             dir[0] = m_F*dir[0] + (1-m_F)*( (1-m_G)*dirOld[0] + scale*m_G*(tensor[0]*dirOld[0] + tensor[1]*dirOld[1] + tensor[2]*dirOld[2]));
             dir[1] = m_F*dir[1] + (1-m_F)*( (1-m_G)*dirOld[1] + scale*m_G*(tensor[1]*dirOld[0] + tensor[3]*dirOld[1] + tensor[4]*dirOld[2]));
             dir[2] = m_F*dir[2] + (1-m_F)*( (1-m_G)*dirOld[2] + scale*m_G*(tensor[2]*dirOld[0] + tensor[4]*dirOld[1] + tensor[5]*dirOld[2]));
             dir.normalize();
 
             float angle = dot_product(dirOld, dir);
             if (angle<0)
             {
                 dir *= -1;
                 angle *= -1;
             }
 
             vnl_vector_fixed<double,3> v3 = dir+dirOld; v3 *= m_StepSize;
             float a = m_StepSize;
             float b = m_StepSize;
             float c = v3.magnitude();
             float r = a*b*c/std::sqrt((a+b+c)*(a+b-c)*(b+c-a)*(a-b+c)); // radius of triangle via Heron's formula (area of triangle)
 
             if (r<m_MinCurvatureRadius)
                 return tractLength;
 
             dirOld = dir;
             indexOld = index;
         }
     }
     return tractLength;
 }
 
 template< class TTensorPixelType,
           class TPDPixelType>
 void StreamlineTrackingFilter< TTensorPixelType,
 TPDPixelType>
 ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
                        ThreadIdType threadId)
 {
     FiberPolyDataType poly = m_PolyDataContainer->GetElement(threadId);
     vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
     vtkSmartPointer<vtkCellArray> Cells = vtkSmartPointer<vtkCellArray>::New();
 
     typedef itk::DiffusionTensor3D<TTensorPixelType>        TensorType;
     typedef ImageRegionConstIterator< InputImageType >      InputIteratorType;
     typedef ImageRegionConstIterator< ItkUcharImgType >     MaskIteratorType;
     typedef ImageRegionConstIterator< ItkFloatImgType >     FloatIteratorType;
     typedef typename InputImageType::PixelType              InputTensorType;
 
     MaskIteratorType    sit(m_SeedImage, outputRegionForThread );
     FloatIteratorType   fit(m_FaImage, outputRegionForThread );
     MaskIteratorType    mit(m_MaskImage, outputRegionForThread );
 
     for (int img=0; img<m_NumberOfInputs; img++)
     {
         sit.GoToBegin();
         mit.GoToBegin();
         fit.GoToBegin();
         itk::Point<double> worldPos;
         while( !sit.IsAtEnd() )
         {
             if (sit.Value()==0 || fit.Value()<m_FaThreshold || mit.Value()==0)
             {
                 ++sit;
                 ++mit;
                 ++fit;
                 continue;
             }
 
             for (int s=0; s<m_SeedsPerVoxel; s++)
             {
                 vtkSmartPointer<vtkPolyLine> line = vtkSmartPointer<vtkPolyLine>::New();
                 std::vector< vtkIdType > pointIDs;
                 typename InputImageType::IndexType index = sit.GetIndex();
                 itk::ContinuousIndex<double, 3> start;
                 unsigned int counter = 0;
 
                 if (m_SeedsPerVoxel>1)
                 {
                     start[0] = index[0]+(double)(rand()%99-49)/100;
                     start[1] = index[1]+(double)(rand()%99-49)/100;
                     start[2] = index[2]+(double)(rand()%99-49)/100;
                 }
                 else
                 {
                     start[0] = index[0];
                     start[1] = index[1];
                     start[2] = index[2];
                 }
 
                 // forward tracking
                 float tractLength = FollowStreamline(start, 1, points, pointIDs, img);
 
                 // add ids to line
                 counter += pointIDs.size();
                 while (!pointIDs.empty())
                 {
                     line->GetPointIds()->InsertNextId(pointIDs.back());
                     pointIDs.pop_back();
                 }
 
                 // insert start point
                 m_SeedImage->TransformContinuousIndexToPhysicalPoint( start, worldPos );
                 line->GetPointIds()->InsertNextId(points->InsertNextPoint(worldPos.GetDataPointer()));
 
                 // backward tracking
                 tractLength += FollowStreamline(start, -1, points, pointIDs, img);
 
                 counter += pointIDs.size();
 
                 //MITK_INFO << "Tract length " << tractLength;
 
                 if (tractLength<m_MinTractLength || counter<2)
                     continue;
 
                 // add ids to line
                 for (unsigned int i=0; i<pointIDs.size(); i++)
                     line->GetPointIds()->InsertNextId(pointIDs.at(i));
 
                 Cells->InsertNextCell(line);
             }
             ++sit;
             ++mit;
             ++fit;
         }
     }
     poly->SetPoints(points);
     poly->SetLines(Cells);
 
     std::cout << "Thread " << threadId << " finished tracking" << std::endl;
 }
 
 template< class TTensorPixelType,
           class TPDPixelType>
 vtkSmartPointer< vtkPolyData > StreamlineTrackingFilter< TTensorPixelType,
 TPDPixelType>
 ::AddPolyData(FiberPolyDataType poly1, FiberPolyDataType poly2)
 {
     vtkSmartPointer<vtkPolyData> vNewPolyData = vtkSmartPointer<vtkPolyData>::New();
     vtkSmartPointer<vtkCellArray> vNewLines = poly1->GetLines();
     vtkSmartPointer<vtkPoints> vNewPoints = poly1->GetPoints();
 
     for( int i=0; i<poly2->GetNumberOfLines(); i++ )
     {
         vtkCell* cell = poly2->GetCell(i);
         int numPoints = cell->GetNumberOfPoints();
         vtkPoints* points = cell->GetPoints();
 
         vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
         for (int j=0; j<numPoints; j++)
         {
             double p[3];
             points->GetPoint(j, p);
 
             vtkIdType id = vNewPoints->InsertNextPoint(p);
             container->GetPointIds()->InsertNextId(id);
         }
         vNewLines->InsertNextCell(container);
     }
 
     // initialize polydata
     vNewPolyData->SetPoints(vNewPoints);
     vNewPolyData->SetLines(vNewLines);
 
     return vNewPolyData;
 }
 template< class TTensorPixelType,
           class TPDPixelType>
 void StreamlineTrackingFilter< TTensorPixelType,
 TPDPixelType>
 ::AfterThreadedGenerateData()
 {
     MITK_INFO << "Generating polydata ";
     m_FiberPolyData = m_PolyDataContainer->GetElement(0);
     for (unsigned int i=1; i<this->GetNumberOfThreads(); i++)
     {
         m_FiberPolyData = AddPolyData(m_FiberPolyData, m_PolyDataContainer->GetElement(i));
     }
     MITK_INFO << "done";
 }
 
 template< class TTensorPixelType,
           class TPDPixelType>
 void StreamlineTrackingFilter< TTensorPixelType,
 TPDPixelType>
 ::PrintSelf(std::ostream&, Indent) const
 {
 }
 
 }
 #endif // __itkDiffusionQballPrincipleDirectionsImageFilter_txx
diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp
index ed12ca35d6..07b33e4053 100755
--- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp
@@ -1,899 +1,901 @@
 /*===================================================================
 
 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 "itkTractsToDWIImageFilter.h"
 #include <boost/progress.hpp>
 #include <vtkSmartPointer.h>
 #include <vtkPolyData.h>
 #include <vtkCellArray.h>
 #include <vtkPoints.h>
 #include <vtkPolyLine.h>
 #include <itkImageRegionIteratorWithIndex.h>
 #include <itkResampleImageFilter.h>
 #include <itkNearestNeighborInterpolateImageFunction.h>
 #include <itkBSplineInterpolateImageFunction.h>
 #include <itkCastImageFilter.h>
 #include <itkImageFileWriter.h>
 #include <itkRescaleIntensityImageFilter.h>
 #include <itkWindowedSincInterpolateImageFunction.h>
 #include <itkResampleDwiImageFilter.h>
 #include <itkKspaceImageFilter.h>
 #include <itkDftImageFilter.h>
 #include <itkAddImageFilter.h>
 #include <itkConstantPadImageFilter.h>
 #include <itkCropImageFilter.h>
 #include <mitkAstroStickModel.h>
 #include <vtkTransform.h>
 #include <iostream>
 #include <fstream>
 #include <itkImageDuplicator.h>
 #include <itksys/SystemTools.hxx>
 #include <mitkIOUtil.h>
 #include <boost/lexical_cast.hpp>
 
 namespace itk
 {
 
 template< class PixelType >
 TractsToDWIImageFilter< PixelType >::TractsToDWIImageFilter()
-    : m_UseConstantRandSeed(false)
+    : m_FiberBundle(NULL)
+    , m_StatusText("")
+    , m_UseConstantRandSeed(false)
+    , m_RandGen(itk::Statistics::MersenneTwisterRandomVariateGenerator::New())
 {
-    m_RandGen = itk::Statistics::MersenneTwisterRandomVariateGenerator::New();
     m_RandGen->SetSeed();
 }
 
 template< class PixelType >
 TractsToDWIImageFilter< PixelType >::~TractsToDWIImageFilter()
 {
 
 }
 
 template< class PixelType >
 TractsToDWIImageFilter< PixelType >::DoubleDwiType::Pointer TractsToDWIImageFilter< PixelType >::DoKspaceStuff( std::vector< DoubleDwiType::Pointer >& images )
 {
     // create slice object
     ImageRegion<2> sliceRegion;
     sliceRegion.SetSize(0, m_UpsampledImageRegion.GetSize()[0]);
     sliceRegion.SetSize(1, m_UpsampledImageRegion.GetSize()[1]);
     Vector< double, 2 > sliceSpacing;
     sliceSpacing[0] = m_UpsampledSpacing[0];
     sliceSpacing[1] = m_UpsampledSpacing[1];
 
     // frequency map slice
     SliceType::Pointer fMapSlice = NULL;
     if (m_Parameters.m_FrequencyMap.IsNotNull())
     {
         fMapSlice = SliceType::New();
         ImageRegion<2> region;
         region.SetSize(0, m_UpsampledImageRegion.GetSize()[0]);
         region.SetSize(1, m_UpsampledImageRegion.GetSize()[1]);
         fMapSlice->SetLargestPossibleRegion( region );
         fMapSlice->SetBufferedRegion( region );
         fMapSlice->SetRequestedRegion( region );
         fMapSlice->Allocate();
         fMapSlice->FillBuffer(0.0);
     }
 
     DoubleDwiType::Pointer newImage = DoubleDwiType::New();
     newImage->SetSpacing( m_Parameters.m_ImageSpacing );
     newImage->SetOrigin( m_Parameters.m_ImageOrigin );
     newImage->SetDirection( m_Parameters.m_ImageDirection );
     newImage->SetLargestPossibleRegion( m_Parameters.m_ImageRegion );
     newImage->SetBufferedRegion( m_Parameters.m_ImageRegion );
     newImage->SetRequestedRegion( m_Parameters.m_ImageRegion );
     newImage->SetVectorLength( images.at(0)->GetVectorLength() );
     newImage->Allocate();
 
     std::vector< unsigned int > spikeVolume;
     for (unsigned int i=0; i<m_Parameters.m_Spikes; i++)
         spikeVolume.push_back(m_RandGen->GetIntegerVariate()%images.at(0)->GetVectorLength());
     std::sort (spikeVolume.begin(), spikeVolume.end());
     std::reverse (spikeVolume.begin(), spikeVolume.end());
 
     m_StatusText += "0%   10   20   30   40   50   60   70   80   90   100%\n";
     m_StatusText += "|----|----|----|----|----|----|----|----|----|----|\n*";
     unsigned long lastTick = 0;
 
     boost::progress_display disp(2*images.at(0)->GetVectorLength()*images.at(0)->GetLargestPossibleRegion().GetSize(2));
     for (unsigned int g=0; g<images.at(0)->GetVectorLength(); g++)
     {
         std::vector< unsigned int > spikeSlice;
         while (!spikeVolume.empty() && spikeVolume.back()==g)
         {
             spikeSlice.push_back(m_RandGen->GetIntegerVariate()%images.at(0)->GetLargestPossibleRegion().GetSize(2));
             spikeVolume.pop_back();
         }
         std::sort (spikeSlice.begin(), spikeSlice.end());
         std::reverse (spikeSlice.begin(), spikeSlice.end());
 
         for (unsigned int z=0; z<images.at(0)->GetLargestPossibleRegion().GetSize(2); z++)
         {
             std::vector< SliceType::Pointer > compartmentSlices;
             std::vector< double > t2Vector;
 
             for (unsigned int i=0; i<images.size(); i++)
             {
                 DiffusionSignalModel<double>* signalModel;
                 if (i<m_Parameters.m_FiberModelList.size())
                     signalModel = m_Parameters.m_FiberModelList.at(i);
                 else
                     signalModel = m_Parameters.m_NonFiberModelList.at(i-m_Parameters.m_FiberModelList.size());
 
                 SliceType::Pointer slice = SliceType::New();
                 slice->SetLargestPossibleRegion( sliceRegion );
                 slice->SetBufferedRegion( sliceRegion );
                 slice->SetRequestedRegion( sliceRegion );
                 slice->SetSpacing(sliceSpacing);
                 slice->Allocate();
                 slice->FillBuffer(0.0);
 
                 // extract slice from channel g
                 for (unsigned int y=0; y<images.at(0)->GetLargestPossibleRegion().GetSize(1); y++)
                     for (unsigned int x=0; x<images.at(0)->GetLargestPossibleRegion().GetSize(0); x++)
                     {
                         SliceType::IndexType index2D; index2D[0]=x; index2D[1]=y;
                         DoubleDwiType::IndexType index3D; index3D[0]=x; index3D[1]=y; index3D[2]=z;
 
                         slice->SetPixel(index2D, images.at(i)->GetPixel(index3D)[g]);
 
                         if (fMapSlice.IsNotNull() && i==0)
                             fMapSlice->SetPixel(index2D, m_Parameters.m_FrequencyMap->GetPixel(index3D));
                     }
 
                 compartmentSlices.push_back(slice);
                 t2Vector.push_back(signalModel->GetT2());
             }
 
             if (this->GetAbortGenerateData())
                 return NULL;
 
             // create k-sapce (inverse fourier transform slices)
             itk::Size<2> outSize; outSize.SetElement(0, m_Parameters.m_ImageRegion.GetSize(0)); outSize.SetElement(1, m_Parameters.m_ImageRegion.GetSize(1));
             itk::KspaceImageFilter< SliceType::PixelType >::Pointer idft = itk::KspaceImageFilter< SliceType::PixelType >::New();
             idft->SetCompartmentImages(compartmentSlices);
             idft->SetT2(t2Vector);
             idft->SetUseConstantRandSeed(m_UseConstantRandSeed);
             idft->SetParameters(m_Parameters);
             idft->SetZ((double)z-(double)images.at(0)->GetLargestPossibleRegion().GetSize(2)/2.0);
             idft->SetDiffusionGradientDirection(m_Parameters.GetGradientDirection(g));
             idft->SetFrequencyMapSlice(fMapSlice);
             idft->SetOutSize(outSize);
             int numSpikes = 0;
             while (!spikeSlice.empty() && spikeSlice.back()==z)
             {
                 numSpikes++;
                 spikeSlice.pop_back();
             }
             idft->SetSpikesPerSlice(numSpikes);
             idft->Update();
 
             ComplexSliceType::Pointer fSlice;
             fSlice = idft->GetOutput();
 
             ++disp;
             unsigned long newTick = 50*disp.count()/disp.expected_count();
             for (unsigned  long tick = 0; tick<(newTick-lastTick); tick++)
                 m_StatusText += "*";
             lastTick = newTick;
 
             // fourier transform slice
             SliceType::Pointer newSlice;
             itk::DftImageFilter< SliceType::PixelType >::Pointer dft = itk::DftImageFilter< SliceType::PixelType >::New();
             dft->SetInput(fSlice);
             dft->Update();
             newSlice = dft->GetOutput();
 
             // put slice back into channel g
             for (unsigned int y=0; y<fSlice->GetLargestPossibleRegion().GetSize(1); y++)
                 for (unsigned int x=0; x<fSlice->GetLargestPossibleRegion().GetSize(0); x++)
                 {
                     DoubleDwiType::IndexType index3D; index3D[0]=x; index3D[1]=y; index3D[2]=z;
                     SliceType::IndexType index2D; index2D[0]=x; index2D[1]=y;
 
                     DoubleDwiType::PixelType pix3D = newImage->GetPixel(index3D);
                     pix3D[g] = newSlice->GetPixel(index2D);
                     newImage->SetPixel(index3D, pix3D);
                 }
             ++disp;
             newTick = 50*disp.count()/disp.expected_count();
             for (unsigned long tick = 0; tick<(newTick-lastTick); tick++)
                 m_StatusText += "*";
             lastTick = newTick;
         }
     }
     m_StatusText += "\n\n";
     return newImage;
 }
 
 template< class PixelType >
 void TractsToDWIImageFilter< PixelType >::GenerateData()
 {
     m_TimeProbe.Start();
     m_StatusText = "Starting simulation\n";
 
     // check input data
     if (m_FiberBundle.IsNull())
         itkExceptionMacro("Input fiber bundle is NULL!");
 
     if (m_Parameters.m_DoDisablePartialVolume)
         while (m_Parameters.m_FiberModelList.size()>1)
             m_Parameters.m_FiberModelList.pop_back();
 
     if (m_Parameters.m_NonFiberModelList.empty())
         itkExceptionMacro("No diffusion model for non-fiber compartments defined!");
 
     int baselineIndex = m_Parameters.GetFirstBaselineIndex();
     if (baselineIndex<0)
         itkExceptionMacro("No baseline index found!");
 
     if (m_UseConstantRandSeed)  // always generate the same random numbers?
         m_RandGen->SetSeed(0);
     else
         m_RandGen->SetSeed();
 
     // initialize output dwi image
     ImageRegion<3> croppedRegion = m_Parameters.m_ImageRegion; croppedRegion.SetSize(1, croppedRegion.GetSize(1)*m_Parameters.m_CroppingFactor);
     itk::Point<double,3> shiftedOrigin = m_Parameters.m_ImageOrigin; shiftedOrigin[1] += (m_Parameters.m_ImageRegion.GetSize(1)-croppedRegion.GetSize(1))*m_Parameters.m_ImageSpacing[1]/2;
     typename OutputImageType::Pointer outImage = OutputImageType::New();
     outImage->SetSpacing( m_Parameters.m_ImageSpacing );
     outImage->SetOrigin( shiftedOrigin );
     outImage->SetDirection( m_Parameters.m_ImageDirection );
     outImage->SetLargestPossibleRegion( croppedRegion );
     outImage->SetBufferedRegion( croppedRegion );
     outImage->SetRequestedRegion( croppedRegion );
     outImage->SetVectorLength( m_Parameters.GetNumVolumes() );
     outImage->Allocate();
     typename OutputImageType::PixelType temp;
     temp.SetSize(m_Parameters.GetNumVolumes());
     temp.Fill(0.0);
     outImage->FillBuffer(temp);
 
     // ADJUST GEOMETRY FOR FURTHER PROCESSING
     // is input slize size a power of two?
     unsigned int x=m_Parameters.m_ImageRegion.GetSize(0); unsigned int y=m_Parameters.m_ImageRegion.GetSize(1);
     ItkDoubleImgType::SizeType pad; pad[0]=x%2; pad[1]=y%2; pad[2]=0;
     m_Parameters.m_ImageRegion.SetSize(0, x+pad[0]);
     m_Parameters.m_ImageRegion.SetSize(1, y+pad[1]);
     if (m_Parameters.m_FrequencyMap.IsNotNull() && (pad[0]>0 || pad[1]>0))
     {
         itk::ConstantPadImageFilter<ItkDoubleImgType, ItkDoubleImgType>::Pointer zeroPadder = itk::ConstantPadImageFilter<ItkDoubleImgType, ItkDoubleImgType>::New();
         zeroPadder->SetInput(m_Parameters.m_FrequencyMap);
         zeroPadder->SetConstant(0);
         zeroPadder->SetPadUpperBound(pad);
         zeroPadder->Update();
         m_Parameters.m_FrequencyMap = zeroPadder->GetOutput();
     }
     if (m_Parameters.m_MaskImage.IsNotNull() && (pad[0]>0 || pad[1]>0))
     {
         itk::ConstantPadImageFilter<ItkUcharImgType, ItkUcharImgType>::Pointer zeroPadder = itk::ConstantPadImageFilter<ItkUcharImgType, ItkUcharImgType>::New();
         zeroPadder->SetInput(m_Parameters.m_MaskImage);
         zeroPadder->SetConstant(0);
         zeroPadder->SetPadUpperBound(pad);
         zeroPadder->Update();
         m_Parameters.m_MaskImage = zeroPadder->GetOutput();
     }
 
     // Apply in-plane upsampling for Gibbs ringing artifact
     double upsampling = 1;
     if (m_Parameters.m_DoAddGibbsRinging)
         upsampling = 2;
     m_UpsampledSpacing = m_Parameters.m_ImageSpacing;
     m_UpsampledSpacing[0] /= upsampling;
     m_UpsampledSpacing[1] /= upsampling;
     m_UpsampledImageRegion = m_Parameters.m_ImageRegion;
     m_UpsampledImageRegion.SetSize(0, m_Parameters.m_ImageRegion.GetSize()[0]*upsampling);
     m_UpsampledImageRegion.SetSize(1, m_Parameters.m_ImageRegion.GetSize()[1]*upsampling);
     m_UpsampledOrigin = m_Parameters.m_ImageOrigin;
     m_UpsampledOrigin[0] -= m_Parameters.m_ImageSpacing[0]/2; m_UpsampledOrigin[0] += m_UpsampledSpacing[0]/2;
     m_UpsampledOrigin[1] -= m_Parameters.m_ImageSpacing[1]/2; m_UpsampledOrigin[1] += m_UpsampledSpacing[1]/2;
     m_UpsampledOrigin[2] -= m_Parameters.m_ImageSpacing[2]/2; m_UpsampledOrigin[2] += m_UpsampledSpacing[2]/2;
 
     // generate double images to store the individual compartment signals
     std::vector< DoubleDwiType::Pointer > compartments;
     for (unsigned int i=0; i<m_Parameters.m_FiberModelList.size()+m_Parameters.m_NonFiberModelList.size(); i++)
     {
         DoubleDwiType::Pointer doubleDwi = DoubleDwiType::New();
         doubleDwi->SetSpacing( m_UpsampledSpacing );
         doubleDwi->SetOrigin( m_UpsampledOrigin );
         doubleDwi->SetDirection( m_Parameters.m_ImageDirection );
         doubleDwi->SetLargestPossibleRegion( m_UpsampledImageRegion );
         doubleDwi->SetBufferedRegion( m_UpsampledImageRegion );
         doubleDwi->SetRequestedRegion( m_UpsampledImageRegion );
         doubleDwi->SetVectorLength( m_Parameters.GetNumVolumes() );
         doubleDwi->Allocate();
         DoubleDwiType::PixelType pix;
         pix.SetSize(m_Parameters.GetNumVolumes());
         pix.Fill(0.0);
         doubleDwi->FillBuffer(pix);
         compartments.push_back(doubleDwi);
     }
 
     // initialize volume fraction images
     m_VolumeFractions.clear();
     for (unsigned int i=0; i<m_Parameters.m_FiberModelList.size()+m_Parameters.m_NonFiberModelList.size(); i++)
     {
         ItkDoubleImgType::Pointer doubleImg = ItkDoubleImgType::New();
         doubleImg->SetSpacing( m_UpsampledSpacing );
         doubleImg->SetOrigin( m_UpsampledOrigin );
         doubleImg->SetDirection( m_Parameters.m_ImageDirection );
         doubleImg->SetLargestPossibleRegion( m_UpsampledImageRegion );
         doubleImg->SetBufferedRegion( m_UpsampledImageRegion );
         doubleImg->SetRequestedRegion( m_UpsampledImageRegion );
         doubleImg->Allocate();
         doubleImg->FillBuffer(0);
         m_VolumeFractions.push_back(doubleImg);
     }
 
     // resample mask image and frequency map to fit upsampled geometry
     if (m_Parameters.m_DoAddGibbsRinging)
     {
         if (m_Parameters.m_MaskImage.IsNotNull())
         {
             // rescale mask image (otherwise there are problems with the resampling)
             itk::RescaleIntensityImageFilter<ItkUcharImgType,ItkUcharImgType>::Pointer rescaler = itk::RescaleIntensityImageFilter<ItkUcharImgType,ItkUcharImgType>::New();
             rescaler->SetInput(0,m_Parameters.m_MaskImage);
             rescaler->SetOutputMaximum(100);
             rescaler->SetOutputMinimum(0);
             rescaler->Update();
 
             // resample mask image
             itk::ResampleImageFilter<ItkUcharImgType, ItkUcharImgType>::Pointer resampler = itk::ResampleImageFilter<ItkUcharImgType, ItkUcharImgType>::New();
             resampler->SetInput(rescaler->GetOutput());
             resampler->SetOutputParametersFromImage(m_Parameters.m_MaskImage);
             resampler->SetSize(m_UpsampledImageRegion.GetSize());
             resampler->SetOutputSpacing(m_UpsampledSpacing);
             resampler->SetOutputOrigin(m_UpsampledOrigin);
             resampler->Update();
             m_Parameters.m_MaskImage = resampler->GetOutput();
         }
         // resample frequency map
         if (m_Parameters.m_FrequencyMap.IsNotNull())
         {
             itk::ResampleImageFilter<ItkDoubleImgType, ItkDoubleImgType>::Pointer resampler = itk::ResampleImageFilter<ItkDoubleImgType, ItkDoubleImgType>::New();
             resampler->SetInput(m_Parameters.m_FrequencyMap);
             resampler->SetOutputParametersFromImage(m_Parameters.m_FrequencyMap);
             resampler->SetSize(m_UpsampledImageRegion.GetSize());
             resampler->SetOutputSpacing(m_UpsampledSpacing);
             resampler->SetOutputOrigin(m_UpsampledOrigin);
             resampler->Update();
             m_Parameters.m_FrequencyMap = resampler->GetOutput();
         }
     }
 
     // no input tissue mask is set -> create default
     bool maskImageSet = true;
     if (m_Parameters.m_MaskImage.IsNull())
     {
         m_StatusText += "No tissue mask set\n";
         MITK_INFO << "No tissue mask set";
         m_Parameters.m_MaskImage = ItkUcharImgType::New();
         m_Parameters.m_MaskImage->SetSpacing( m_UpsampledSpacing );
         m_Parameters.m_MaskImage->SetOrigin( m_UpsampledOrigin );
         m_Parameters.m_MaskImage->SetDirection( m_Parameters.m_ImageDirection );
         m_Parameters.m_MaskImage->SetLargestPossibleRegion( m_UpsampledImageRegion );
         m_Parameters.m_MaskImage->SetBufferedRegion( m_UpsampledImageRegion );
         m_Parameters.m_MaskImage->SetRequestedRegion( m_UpsampledImageRegion );
         m_Parameters.m_MaskImage->Allocate();
         m_Parameters.m_MaskImage->FillBuffer(1);
         maskImageSet = false;
     }
     else
     {
         m_StatusText += "Using tissue mask\n";
         MITK_INFO << "Using tissue mask";
     }
 
     m_Parameters.m_ImageRegion = croppedRegion;
     x=m_Parameters.m_ImageRegion.GetSize(0); y=m_Parameters.m_ImageRegion.GetSize(1);
     if ( x%2 == 1 )
         m_Parameters.m_ImageRegion.SetSize(0, x+1);
     if ( y%2 == 1 )
         m_Parameters.m_ImageRegion.SetSize(1, y+1);
 
     // resample fiber bundle for sufficient voxel coverage
     m_StatusText += "\n"+this->GetTime()+" > Resampling fibers ...\n";
     double segmentVolume = 0.0001;
     float minSpacing = 1;
     if(m_UpsampledSpacing[0]<m_UpsampledSpacing[1] && m_UpsampledSpacing[0]<m_UpsampledSpacing[2])
         minSpacing = m_UpsampledSpacing[0];
     else if (m_UpsampledSpacing[1] < m_UpsampledSpacing[2])
         minSpacing = m_UpsampledSpacing[1];
     else
         minSpacing = m_UpsampledSpacing[2];
     FiberBundleType fiberBundle = m_FiberBundle->GetDeepCopy();
     double volumeAccuracy = 10;
     fiberBundle->ResampleFibers(minSpacing/volumeAccuracy);
     double mmRadius = m_Parameters.m_AxonRadius/1000;
     if (mmRadius>0)
         segmentVolume = M_PI*mmRadius*mmRadius*minSpacing/volumeAccuracy;
 
     double maxVolume = 0;
     double voxelVolume = m_UpsampledSpacing[0]*m_UpsampledSpacing[1]*m_UpsampledSpacing[2];
 
     ofstream logFile;
     if (m_Parameters.m_DoAddMotion)
     {
         std::string fileName = "fiberfox_motion_0.log";
         int c = 1;
         while (itksys::SystemTools::FileExists(mitk::IOUtil::GetTempPath().append(fileName).c_str()))
         {
             fileName = "fiberfox_motion_";
             fileName += boost::lexical_cast<std::string>(c);
             fileName += ".log";
             c++;
         }
 
         logFile.open(mitk::IOUtil::GetTempPath().append(fileName).c_str());
         logFile << "0 rotation: 0,0,0; translation: 0,0,0\n";
 
         if (m_Parameters.m_DoRandomizeMotion)
         {
             m_StatusText += "Adding random motion artifacts:\n";
             m_StatusText += "Maximum rotation: +/-" + boost::lexical_cast<std::string>(m_Parameters.m_Rotation) + "°\n";
             m_StatusText += "Maximum translation: +/-" + boost::lexical_cast<std::string>(m_Parameters.m_Translation) + "mm\n";
         }
         else
         {
             m_StatusText += "Adding linear motion artifacts:\n";
             m_StatusText += "Maximum rotation: " + boost::lexical_cast<std::string>(m_Parameters.m_Rotation) + "°\n";
             m_StatusText += "Maximum translation: " + boost::lexical_cast<std::string>(m_Parameters.m_Translation) + "mm\n";
         }
         m_StatusText += "Motion logfile: " + mitk::IOUtil::GetTempPath().append(fileName) + "\n";
         MITK_INFO << "Adding motion artifacts";
         MITK_INFO << "Maximum rotation: " << m_Parameters.m_Rotation;
         MITK_INFO << "Maxmimum translation: " << m_Parameters.m_Translation;
     }
     maxVolume = 0;
 
     m_StatusText += "\n"+this->GetTime()+" > Generating signal of " + boost::lexical_cast<std::string>(m_Parameters.m_FiberModelList.size()) + " fiber compartments\n";
     MITK_INFO << "Generating signal of " << m_Parameters.m_FiberModelList.size() << " fiber compartments";
     int numFibers = m_FiberBundle->GetNumFibers();
     boost::progress_display disp(numFibers*m_Parameters.GetNumVolumes());
 
 
     // get transform for motion artifacts
     FiberBundleType fiberBundleTransformed = fiberBundle;
     VectorType rotation = m_Parameters.m_Rotation/m_Parameters.GetNumVolumes();
     VectorType translation = m_Parameters.m_Translation/m_Parameters.GetNumVolumes();
 
     // creat image to hold transformed mask (motion artifact)
     ItkUcharImgType::Pointer tempTissueMask = ItkUcharImgType::New();
     itk::ImageDuplicator<ItkUcharImgType>::Pointer duplicator = itk::ImageDuplicator<ItkUcharImgType>::New();
     duplicator->SetInputImage(m_Parameters.m_MaskImage);
     duplicator->Update();
     tempTissueMask = duplicator->GetOutput();
 
     // second upsampling needed for motion artifacts
     ImageRegion<3>                      upsampledImageRegion = m_UpsampledImageRegion;
     itk::Vector<double,3>               upsampledSpacing = m_UpsampledSpacing;
     upsampledSpacing[0] /= 4;
     upsampledSpacing[1] /= 4;
     upsampledSpacing[2] /= 4;
     upsampledImageRegion.SetSize(0, m_UpsampledImageRegion.GetSize()[0]*4);
     upsampledImageRegion.SetSize(1, m_UpsampledImageRegion.GetSize()[1]*4);
     upsampledImageRegion.SetSize(2, m_UpsampledImageRegion.GetSize()[2]*4);
     itk::Point<double,3> upsampledOrigin = m_UpsampledOrigin;
     upsampledOrigin[0] -= m_UpsampledSpacing[0]/2; upsampledOrigin[0] += upsampledSpacing[0]/2;
     upsampledOrigin[1] -= m_UpsampledSpacing[1]/2; upsampledOrigin[1] += upsampledSpacing[1]/2;
     upsampledOrigin[2] -= m_UpsampledSpacing[2]/2; upsampledOrigin[2] += upsampledSpacing[2]/2;
     ItkUcharImgType::Pointer upsampledTissueMask = ItkUcharImgType::New();
     itk::ResampleImageFilter<ItkUcharImgType, ItkUcharImgType>::Pointer upsampler = itk::ResampleImageFilter<ItkUcharImgType, ItkUcharImgType>::New();
     upsampler->SetInput(m_Parameters.m_MaskImage);
     upsampler->SetOutputParametersFromImage(m_Parameters.m_MaskImage);
     upsampler->SetSize(upsampledImageRegion.GetSize());
     upsampler->SetOutputSpacing(upsampledSpacing);
     upsampler->SetOutputOrigin(upsampledOrigin);
     itk::NearestNeighborInterpolateImageFunction<ItkUcharImgType>::Pointer nn_interpolator
             = itk::NearestNeighborInterpolateImageFunction<ItkUcharImgType>::New();
     upsampler->SetInterpolator(nn_interpolator);
     upsampler->Update();
     upsampledTissueMask = upsampler->GetOutput();
 
     m_StatusText += "0%   10   20   30   40   50   60   70   80   90   100%\n";
     m_StatusText += "|----|----|----|----|----|----|----|----|----|----|\n*";
     unsigned long lastTick = 0;
 
     for (unsigned int g=0; g<m_Parameters.GetNumVolumes(); g++)
     {
         vtkPolyData* fiberPolyData = fiberBundleTransformed->GetFiberPolyData();
 
         ItkDoubleImgType::Pointer intraAxonalVolumeImage = ItkDoubleImgType::New();
         intraAxonalVolumeImage->SetSpacing( m_UpsampledSpacing );
         intraAxonalVolumeImage->SetOrigin( m_UpsampledOrigin );
         intraAxonalVolumeImage->SetDirection( m_Parameters.m_ImageDirection );
         intraAxonalVolumeImage->SetLargestPossibleRegion( m_UpsampledImageRegion );
         intraAxonalVolumeImage->SetBufferedRegion( m_UpsampledImageRegion );
         intraAxonalVolumeImage->SetRequestedRegion( m_UpsampledImageRegion );
         intraAxonalVolumeImage->Allocate();
         intraAxonalVolumeImage->FillBuffer(0);
 
         // generate fiber signal (if there are any fiber models present)
         if (!m_Parameters.m_FiberModelList.empty())
             for( int i=0; i<numFibers; i++ )
             {
                 vtkCell* cell = fiberPolyData->GetCell(i);
                 int numPoints = cell->GetNumberOfPoints();
                 vtkPoints* points = cell->GetPoints();
 
                 if (numPoints<2)
                     continue;
 
                 for( int j=0; j<numPoints; j++)
                 {
                     if (this->GetAbortGenerateData())
                     {
                         m_StatusText += "\n"+this->GetTime()+" > Simulation aborted\n";
                         return;
                     }
 
                     double* temp = points->GetPoint(j);
                     itk::Point<float, 3> vertex = GetItkPoint(temp);
                     itk::Vector<double> v = GetItkVector(temp);
 
                     itk::Vector<double, 3> dir(3);
                     if (j<numPoints-1)
                         dir = GetItkVector(points->GetPoint(j+1))-v;
                     else
                         dir = v-GetItkVector(points->GetPoint(j-1));
 
                     if (dir.GetSquaredNorm()<0.0001 || dir[0]!=dir[0] || dir[1]!=dir[1] || dir[2]!=dir[2])
                         continue;
 
                     itk::Index<3> idx;
                     itk::ContinuousIndex<float, 3> contIndex;
                     tempTissueMask->TransformPhysicalPointToIndex(vertex, idx);
                     tempTissueMask->TransformPhysicalPointToContinuousIndex(vertex, contIndex);
 
                     if (!tempTissueMask->GetLargestPossibleRegion().IsInside(idx) || tempTissueMask->GetPixel(idx)<=0)
                         continue;
 
                     // generate signal for each fiber compartment
                     for (unsigned int k=0; k<m_Parameters.m_FiberModelList.size(); k++)
                     {
                         m_Parameters.m_FiberModelList[k]->SetFiberDirection(dir);
                         DoubleDwiType::PixelType pix = compartments.at(k)->GetPixel(idx);
                         pix[g] += segmentVolume*m_Parameters.m_FiberModelList[k]->SimulateMeasurement(g);
                         compartments.at(k)->SetPixel(idx, pix);
                     }
 
                     // update fiber volume image
                     double vol = intraAxonalVolumeImage->GetPixel(idx) + segmentVolume;
                     intraAxonalVolumeImage->SetPixel(idx, vol);
                     if (g==0 && vol>maxVolume)
                         maxVolume = vol;
                 }
 
                 // progress report
                 ++disp;
                 unsigned long newTick = 50*disp.count()/disp.expected_count();
                 for (unsigned int tick = 0; tick<(newTick-lastTick); tick++)
                     m_StatusText += "*";
                 lastTick = newTick;
             }
 
         // generate non-fiber signal
         ImageRegionIterator<ItkUcharImgType> it3(tempTissueMask, tempTissueMask->GetLargestPossibleRegion());
         double fact = 1;
         if (m_Parameters.m_AxonRadius<0.0001 || maxVolume>voxelVolume)
             fact = voxelVolume/maxVolume;
         while(!it3.IsAtEnd())
         {
             if (it3.Get()>0)
             {
                 DoubleDwiType::IndexType index = it3.GetIndex();
 
                 // get fiber volume fraction
                 double intraAxonalVolume = intraAxonalVolumeImage->GetPixel(index)*fact;
 
                 for (unsigned int i=0; i<m_Parameters.m_FiberModelList.size(); i++)
                 {
                     DoubleDwiType::PixelType pix = compartments.at(i)->GetPixel(index);
                     pix[g] *= fact;
                     compartments.at(i)->SetPixel(index, pix);
                 }
 
                 if (intraAxonalVolume>0.0001 && m_Parameters.m_DoDisablePartialVolume)  // only fiber in voxel
                 {
                     DoubleDwiType::PixelType pix = compartments.at(0)->GetPixel(index);
                     pix[g] *= voxelVolume/intraAxonalVolume;
                     compartments.at(0)->SetPixel(index, pix);
                     m_VolumeFractions.at(0)->SetPixel(index, 1);
                     for (unsigned int i=1; i<m_Parameters.m_FiberModelList.size(); i++)
                     {
                         DoubleDwiType::PixelType pix = compartments.at(i)->GetPixel(index);
                         pix[g] = 0;
                         compartments.at(i)->SetPixel(index, pix);
                     }
                 }
                 else
                 {
                     m_VolumeFractions.at(0)->SetPixel(index, intraAxonalVolume/voxelVolume);
 
                     double extraAxonalVolume = voxelVolume-intraAxonalVolume;    // non-fiber volume
                     double interAxonalVolume = 0;
                     if (m_Parameters.m_FiberModelList.size()>1)
                         interAxonalVolume = extraAxonalVolume * intraAxonalVolume/voxelVolume;   // inter-axonal fraction of non fiber compartment scales linearly with f
                     double other = extraAxonalVolume - interAxonalVolume;        // rest of compartment
                     double singleinter = interAxonalVolume/(m_Parameters.m_FiberModelList.size()-1);
 
                     // adjust non-fiber and intra-axonal signal
                     for (unsigned int i=1; i<m_Parameters.m_FiberModelList.size(); i++)
                     {
                         DoubleDwiType::PixelType pix = compartments.at(i)->GetPixel(index);
                         if (intraAxonalVolume>0)    // remove scaling by intra-axonal volume from inter-axonal compartment
                             pix[g] /= intraAxonalVolume;
                         pix[g] *= singleinter;
                         compartments.at(i)->SetPixel(index, pix);
                         m_VolumeFractions.at(i)->SetPixel(index, singleinter/voxelVolume);
                     }
                     for (unsigned int i=0; i<m_Parameters.m_NonFiberModelList.size(); i++)
                     {
                         DoubleDwiType::Pointer doubleDwi = compartments.at(i+m_Parameters.m_FiberModelList.size());
                         DoubleDwiType::PixelType pix = doubleDwi->GetPixel(index);
                         pix[g] += m_Parameters.m_NonFiberModelList[i]->SimulateMeasurement(g)*other*m_Parameters.m_NonFiberModelList[i]->GetWeight();
                         doubleDwi->SetPixel(index, pix);
                         m_VolumeFractions.at(i+m_Parameters.m_FiberModelList.size())->SetPixel(index, other/voxelVolume*m_Parameters.m_NonFiberModelList[i]->GetWeight());
                     }
                 }
             }
             ++it3;
         }
 
         // move fibers
         if (m_Parameters.m_DoAddMotion)
         {
             if (m_Parameters.m_DoRandomizeMotion)
             {
                 fiberBundleTransformed = fiberBundle->GetDeepCopy();
                 rotation[0] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_Rotation[0]*2)-m_Parameters.m_Rotation[0];
                 rotation[1] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_Rotation[1]*2)-m_Parameters.m_Rotation[1];
                 rotation[2] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_Rotation[2]*2)-m_Parameters.m_Rotation[2];
                 translation[0] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_Translation[0]*2)-m_Parameters.m_Translation[0];
                 translation[1] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_Translation[1]*2)-m_Parameters.m_Translation[1];
                 translation[2] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_Translation[2]*2)-m_Parameters.m_Translation[2];
             }
 
             // rotate mask image
             if (maskImageSet)
             {
                 ImageRegionIterator<ItkUcharImgType> maskIt(upsampledTissueMask, upsampledTissueMask->GetLargestPossibleRegion());
                 tempTissueMask->FillBuffer(0);
 
                 while(!maskIt.IsAtEnd())
                 {
                     if (maskIt.Get()<=0)
                     {
                         ++maskIt;
                         continue;
                     }
 
                     DoubleDwiType::IndexType index = maskIt.GetIndex();
                     itk::Point<double, 3> point;
                     upsampledTissueMask->TransformIndexToPhysicalPoint(index, point);
                     if (m_Parameters.m_DoRandomizeMotion)
                         point = fiberBundle->TransformPoint(point.GetVnlVector(), rotation[0],rotation[1],rotation[2],translation[0],translation[1],translation[2]);
                     else
                         point = fiberBundle->TransformPoint(point.GetVnlVector(), rotation[0]*(g+1),rotation[1]*(g+1),rotation[2]*(g+1),translation[0]*(g+1),translation[1]*(g+1),translation[2]*(g+1));
 
                     tempTissueMask->TransformPhysicalPointToIndex(point, index);
                     if (tempTissueMask->GetLargestPossibleRegion().IsInside(index))
                         tempTissueMask->SetPixel(index,100);
                     ++maskIt;
                 }
             }
 
             // rotate fibers
             if (logFile.is_open())
             {
                 logFile << g+1 << " rotation: " << rotation[0] << "," << rotation[1] << "," << rotation[2] << ";";
                 logFile << " translation: " << translation[0] << "," << translation[1] << "," << translation[2] << "\n";
             }
             fiberBundleTransformed->TransformFibers(rotation[0],rotation[1],rotation[2],translation[0],translation[1],translation[2]);
         }
     }
     if (logFile.is_open())
     {
         logFile << "DONE";
         logFile.close();
     }
     m_StatusText += "\n\n";
     if (this->GetAbortGenerateData())
     {
         m_StatusText += "\n"+this->GetTime()+" > Simulation aborted\n";
         return;
     }
 
     // do k-space stuff
     DoubleDwiType::Pointer doubleOutImage;
     if (m_Parameters.m_Spikes>0 || m_Parameters.m_FrequencyMap.IsNotNull() || m_Parameters.m_KspaceLineOffset>0 || m_Parameters.m_DoSimulateRelaxation || m_Parameters.m_EddyStrength>0 || m_Parameters.m_DoAddGibbsRinging || m_Parameters.m_CroppingFactor<1.0)
     {
         m_StatusText += this->GetTime()+" > Adjusting complex signal\n";
         MITK_INFO << "Adjusting complex signal:";
 
         if (m_Parameters.m_DoSimulateRelaxation)
             m_StatusText += "Simulating signal relaxation\n";
         if (m_Parameters.m_FrequencyMap.IsNotNull())
             m_StatusText += "Simulating distortions\n";
         if (m_Parameters.m_DoAddGibbsRinging)
             m_StatusText += "Simulating ringing artifacts\n";
         if (m_Parameters.m_EddyStrength>0)
             m_StatusText += "Simulating eddy currents\n";
         if (m_Parameters.m_Spikes>0)
             m_StatusText += "Simulating spikes\n";
         if (m_Parameters.m_CroppingFactor<1.0)
             m_StatusText += "Simulating aliasing artifacts\n";
         if (m_Parameters.m_KspaceLineOffset>0)
             m_StatusText += "Simulating ghosts\n";
 
         doubleOutImage = DoKspaceStuff(compartments);
         m_Parameters.m_SignalScale = 1;
     }
     else
     {
         m_StatusText += this->GetTime()+" > Summing compartments\n";
         MITK_INFO << "Summing compartments";
         doubleOutImage = compartments.at(0);
 
         for (unsigned int i=1; i<compartments.size(); i++)
         {
             itk::AddImageFilter< DoubleDwiType, DoubleDwiType, DoubleDwiType>::Pointer adder = itk::AddImageFilter< DoubleDwiType, DoubleDwiType, DoubleDwiType>::New();
             adder->SetInput1(doubleOutImage);
             adder->SetInput2(compartments.at(i));
             adder->Update();
             doubleOutImage = adder->GetOutput();
         }
     }
     if (this->GetAbortGenerateData())
     {
         m_StatusText += "\n"+this->GetTime()+" > Simulation aborted\n";
         return;
     }
 
     m_StatusText += this->GetTime()+" > Finalizing image\n";
     MITK_INFO << "Finalizing image";
     if (m_Parameters.m_SignalScale>1)
         m_StatusText += " Scaling signal\n";
     if (m_Parameters.m_NoiseModel!=NULL)
         m_StatusText += " Adding noise\n";
     unsigned int window = 0;
     unsigned int min = itk::NumericTraits<unsigned int>::max();
     ImageRegionIterator<OutputImageType> it4 (outImage, outImage->GetLargestPossibleRegion());
     DoubleDwiType::PixelType signal; signal.SetSize(m_Parameters.GetNumVolumes());
     boost::progress_display disp2(outImage->GetLargestPossibleRegion().GetNumberOfPixels());
 
     m_StatusText += "0%   10   20   30   40   50   60   70   80   90   100%\n";
     m_StatusText += "|----|----|----|----|----|----|----|----|----|----|\n*";
     lastTick = 0;
 
     while(!it4.IsAtEnd())
     {
         if (this->GetAbortGenerateData())
         {
             m_StatusText += "\n"+this->GetTime()+" > Simulation aborted\n";
             return;
         }
 
         ++disp2;
         unsigned long newTick = 50*disp2.count()/disp2.expected_count();
         for (unsigned long tick = 0; tick<(newTick-lastTick); tick++)
             m_StatusText += "*";
         lastTick = newTick;
 
         typename OutputImageType::IndexType index = it4.GetIndex();
         signal = doubleOutImage->GetPixel(index)*m_Parameters.m_SignalScale;
 
         if (m_Parameters.m_NoiseModel!=NULL)
         {
             DoubleDwiType::PixelType accu = signal; accu.Fill(0.0);
             for (unsigned int i=0; i<m_Parameters.m_Repetitions; i++)
             {
                 DoubleDwiType::PixelType temp = signal;
                 m_Parameters.m_NoiseModel->AddNoise(temp);
                 accu += temp;
             }
             signal = accu/m_Parameters.m_Repetitions;
         }
 
         for (unsigned int i=0; i<signal.Size(); i++)
         {
             if (signal[i]>0)
                 signal[i] = floor(signal[i]+0.5);
             else
                 signal[i] = ceil(signal[i]-0.5);
 
             if (!m_Parameters.IsBaselineIndex(i) && signal[i]>window)
                 window = signal[i];
             if (!m_Parameters.IsBaselineIndex(i) && signal[i]<min)
                 min = signal[i];
         }
         it4.Set(signal);
         ++it4;
     }
     window -= min;
     unsigned int level = window/2 + min;
     m_LevelWindow.SetLevelWindow(level, window);
     this->SetNthOutput(0, outImage);
 
     m_StatusText += "\n\n";
     m_StatusText += "Finished simulation\n";
     m_StatusText += "Simulation time: "+GetTime();
     m_TimeProbe.Stop();
 }
 
 template< class PixelType >
 itk::Point<float, 3> TractsToDWIImageFilter< PixelType >::GetItkPoint(double point[3])
 {
     itk::Point<float, 3> itkPoint;
     itkPoint[0] = point[0];
     itkPoint[1] = point[1];
     itkPoint[2] = point[2];
     return itkPoint;
 }
 
 template< class PixelType >
 itk::Vector<double, 3> TractsToDWIImageFilter< PixelType >::GetItkVector(double point[3])
 {
     itk::Vector<double, 3> itkVector;
     itkVector[0] = point[0];
     itkVector[1] = point[1];
     itkVector[2] = point[2];
     return itkVector;
 }
 
 template< class PixelType >
 vnl_vector_fixed<double, 3> TractsToDWIImageFilter< PixelType >::GetVnlVector(double point[3])
 {
     vnl_vector_fixed<double, 3> vnlVector;
     vnlVector[0] = point[0];
     vnlVector[1] = point[1];
     vnlVector[2] = point[2];
     return vnlVector;
 }
 
 template< class PixelType >
 vnl_vector_fixed<double, 3> TractsToDWIImageFilter< PixelType >::GetVnlVector(Vector<float,3>& vector)
 {
     vnl_vector_fixed<double, 3> vnlVector;
     vnlVector[0] = vector[0];
     vnlVector[1] = vector[1];
     vnlVector[2] = vector[2];
     return vnlVector;
 }
 
 template< class PixelType >
 double TractsToDWIImageFilter< PixelType >::RoundToNearest(double num) {
     return (num > 0.0) ? floor(num + 0.5) : ceil(num - 0.5);
 }
 
 template< class PixelType >
 std::string TractsToDWIImageFilter< PixelType >::GetTime()
 {
     m_TimeProbe.Stop();
     unsigned long total = RoundToNearest(m_TimeProbe.GetTotal());
     unsigned long hours = total/3600;
     unsigned long minutes = (total%3600)/60;
     unsigned long seconds = total%60;
     std::string out = "";
     out.append(boost::lexical_cast<std::string>(hours));
     out.append(":");
     out.append(boost::lexical_cast<std::string>(minutes));
     out.append(":");
     out.append(boost::lexical_cast<std::string>(seconds));
     m_TimeProbe.Start();
     return out;
 }
 
 }
diff --git a/Modules/DiffusionImaging/FiberTracking/Testing/mitkStreamlineTrackingTest.cpp b/Modules/DiffusionImaging/FiberTracking/Testing/mitkStreamlineTrackingTest.cpp
index af85732258..718788205f 100755
--- a/Modules/DiffusionImaging/FiberTracking/Testing/mitkStreamlineTrackingTest.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/Testing/mitkStreamlineTrackingTest.cpp
@@ -1,138 +1,141 @@
 /*===================================================================
 
 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 <mitkImageCast.h>
 #include <mitkTensorImage.h>
 #include <mitkIOUtil.h>
 #include <mitkBaseDataIOFactory.h>
 #include <mitkFiberBundleX.h>
 #include <itkStreamlineTrackingFilter.h>
 #include <itkDiffusionTensor3D.h>
 #include <mitkTestingMacros.h>
 #include <mitkFiberBundleXWriter.h>
 
 using namespace std;
 
 int mitkStreamlineTrackingTest(int argc, char* argv[])
 {
     MITK_TEST_BEGIN("mitkStreamlineTrackingTest");
 
     MITK_TEST_CONDITION_REQUIRED(argc>3,"check for input data")
 
     string dtiFileName = argv[1];
     string maskFileName = argv[2];
     string referenceFileName = argv[3];
 
     MITK_INFO << "DTI file: " << dtiFileName;
     MITK_INFO << "Mask file: " << maskFileName;
     MITK_INFO << "Reference fiber file: " << referenceFileName;
 
     float minFA = 0.05;
     float minCurv = -1;
     float stepSize = -1;
     float tendf = 1;
     float tendg = 0;
     float minLength = 20;
     int numSeeds = 1;
     bool interpolate = true;
 
     try
     {
         // load input image
         const std::string s1="", s2="";
         std::vector<mitk::BaseData::Pointer> infile = mitk::BaseDataIO::LoadBaseDataFromFile( dtiFileName, s1, s2, false );
 
         MITK_INFO << "Loading tensor image ...";
         typedef itk::Image< itk::DiffusionTensor3D<float>, 3 >    ItkTensorImage;
         mitk::TensorImage::Pointer mitkTensorImage = dynamic_cast<mitk::TensorImage*>(infile.at(0).GetPointer());
         ItkTensorImage::Pointer itk_dti = ItkTensorImage::New();
         mitk::CastToItkImage<ItkTensorImage>(mitkTensorImage, itk_dti);
 
         MITK_INFO << "Loading seed image ...";
         typedef itk::Image< unsigned char, 3 >    ItkUCharImageType;
         mitk::Image::Pointer mitkSeedImage = mitk::IOUtil::LoadImage(maskFileName);
 
         MITK_INFO << "Loading mask image ...";
         mitk::Image::Pointer mitkMaskImage = mitk::IOUtil::LoadImage(maskFileName);
 
         // instantiate tracker
         typedef itk::StreamlineTrackingFilter< float > FilterType;
         FilterType::Pointer filter = FilterType::New();
         filter->SetInput(itk_dti);
         filter->SetSeedsPerVoxel(numSeeds);
         filter->SetFaThreshold(minFA);
         filter->SetMinCurvatureRadius(minCurv);
         filter->SetStepSize(stepSize);
         filter->SetF(tendf);
         filter->SetG(tendg);
         filter->SetInterpolate(interpolate);
         filter->SetMinTractLength(minLength);
         filter->SetNumberOfThreads(1);
 
         if (mitkSeedImage.IsNotNull())
         {
             ItkUCharImageType::Pointer mask = ItkUCharImageType::New();
             mitk::CastToItkImage<ItkUCharImageType>(mitkSeedImage, mask);
             filter->SetSeedImage(mask);
         }
 
         if (mitkMaskImage.IsNotNull())
         {
             ItkUCharImageType::Pointer mask = ItkUCharImageType::New();
             mitk::CastToItkImage<ItkUCharImageType>(mitkMaskImage, mask);
             filter->SetMaskImage(mask);
         }
 
         filter->Update();
 
         vtkSmartPointer<vtkPolyData> fiberBundle = filter->GetFiberPolyData();
         mitk::FiberBundleX::Pointer fib1 = mitk::FiberBundleX::New(fiberBundle);
 
         infile = mitk::BaseDataIO::LoadBaseDataFromFile( referenceFileName, s1, s2, false );
         mitk::FiberBundleX::Pointer fib2 = dynamic_cast<mitk::FiberBundleX*>(infile.at(0).GetPointer());
         MITK_TEST_CONDITION_REQUIRED(fib2.IsNotNull(), "Check if reference tractogram is not null.");
-        if (!fib1->Equals(fib2))
+        bool ok = fib1->Equals(fib2);
+        if (!ok)
         {
             MITK_WARN << "TEST FAILED. TRACTOGRAMS ARE NOT EQUAL!";
             mitk::FiberBundleXWriter::Pointer writer = mitk::FiberBundleXWriter::New();
-            writer->SetFileName("testBundle.fib");
+            writer->SetFileName(mitk::IOUtil::GetTempPath()+"testBundle.fib");
             writer->SetInputFiberBundleX(fib1);
             writer->Update();
 
-            writer->SetFileName("refBundle.fib");
+            writer->SetFileName(mitk::IOUtil::GetTempPath()+"refBundle.fib");
             writer->SetInputFiberBundleX(fib2);
             writer->Update();
+
+            MITK_INFO << "OUTPUT: " << mitk::IOUtil::GetTempPath();
         }
-        //MITK_TEST_CONDITION_REQUIRED(fib1->Equals(fib2), "Check if tractograms are equal.");
+        MITK_TEST_CONDITION_REQUIRED(ok, "Check if tractograms 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();
 }