diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkAstroStickModel.cpp b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkAstroStickModel.cpp
index ee2b139e23..6ea65bcbaa 100644
--- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkAstroStickModel.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkAstroStickModel.cpp
@@ -1,127 +1,129 @@
 /*===================================================================
 
 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 <vnl/vnl_cross.h>
 #include <vnl/vnl_quaternion.h>
 #include <mitkAstroStickModel.h>
 
 using namespace mitk;
 
 template< class ScalarType >
 AstroStickModel< ScalarType >::AstroStickModel()
     : m_BValue(1000)
     , m_Diffusivity(0.001)
     , m_NumSticks(42)
     , m_RandomizeSticks(false)
 {
     this->m_RandGen = itk::Statistics::MersenneTwisterRandomVariateGenerator::New();
     this->m_RandGen->SetSeed();
 
     vnl_matrix_fixed<double,3,42>* sticks = itk::PointShell<42, vnl_matrix_fixed<double, 3, 42> >::DistributePointShell();
     for (unsigned int i=0; i<m_NumSticks; i++)
     {
         GradientType stick;
         stick[0] = sticks->get(0,i); stick[1] = sticks->get(1,i); stick[2] = sticks->get(2,i);
         stick.Normalize();
         m_Sticks.push_back(stick);
     }
 }
 
 template< class ScalarType >
 AstroStickModel< ScalarType >::~AstroStickModel()
 {
 
 }
 
 template< class ScalarType >
 ScalarType AstroStickModel< ScalarType >::SimulateMeasurement(unsigned int dir)
 {
     ScalarType signal = 0;
 
     if (dir>=this->m_GradientList.size())
         return signal;
 
     ScalarType b = -m_BValue*m_Diffusivity;
 
     if (m_RandomizeSticks)  // random number of sticks
         m_NumSticks = 30 + this->m_RandGen->GetIntegerVariate()%31;
 
     GradientType g = this->m_GradientList[dir];
-    ScalarType bVal = g.GetNorm(); bVal *= bVal;
+    float norm = g.GetNorm();
+    ScalarType bVal = norm*norm;
 
     if (bVal>0.0001)    // is weighted direction
     {
         for (unsigned int j=0; j<m_NumSticks; j++)
         {
             ScalarType dot = 0;
             if(m_RandomizeSticks)
-                dot = GetRandomDirection()*g;
+                dot = GetRandomDirection()*g/norm;
             else
-                dot = m_Sticks[j]*g;
+                dot = m_Sticks[j]*g/norm;
             signal += std::exp( (double)(b*bVal*dot*dot) );
         }
         signal /= m_NumSticks;
     }
     else    // is baseline direction
         signal = 1;
 
     return signal;
 }
 
 template< class ScalarType >
 typename AstroStickModel< ScalarType >::GradientType AstroStickModel< ScalarType >::GetRandomDirection()
 {
     GradientType vec;
     vec[0] = this->m_RandGen->GetNormalVariate();
     vec[1] = this->m_RandGen->GetNormalVariate();
     vec[2] = this->m_RandGen->GetNormalVariate();
     vec.Normalize();
     return vec;
 }
 
 template< class ScalarType >
 typename AstroStickModel< ScalarType >::PixelType AstroStickModel< ScalarType >::SimulateMeasurement()
 {
     PixelType signal;
     signal.SetSize(this->m_GradientList.size());
     ScalarType b = -m_BValue*m_Diffusivity;
 
     if (m_RandomizeSticks)
         m_NumSticks = 30 + this->m_RandGen->GetIntegerVariate()%31;
 
     for( unsigned int i=0; i<this->m_GradientList.size(); i++)
     {
         GradientType g = this->m_GradientList[i];
-        ScalarType bVal = g.GetNorm(); bVal *= bVal;
+        float norm = g.GetNorm();
+        ScalarType bVal = norm*norm;
 
         if (bVal>0.0001)
         {
             for (unsigned int j=0; j<m_NumSticks; j++)
             {
                 ScalarType dot = 0;
                 if(m_RandomizeSticks)
-                    dot = GetRandomDirection()*g;
+                    dot = GetRandomDirection()*g/norm;
                 else
-                    dot = m_Sticks[j]*g;
+                    dot = m_Sticks[j]*g/norm;
                 signal[i] += std::exp( (double)(b*bVal*dot*dot) );
             }
             signal[i] /= m_NumSticks;
         }
         else
             signal[i] = 1;
     }
 
     return signal;
 }
diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkStickModel.cpp b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkStickModel.cpp
index 29b8cdbb29..13ec41e63a 100644
--- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkStickModel.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkStickModel.cpp
@@ -1,82 +1,84 @@
 /*===================================================================
 
 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 <vnl/vnl_cross.h>
 #include <vnl/vnl_quaternion.h>
 #include <mitkStickModel.h>
 
 using namespace mitk;
 
 template< class ScalarType >
 StickModel< ScalarType >::StickModel()
     : m_Diffusivity(0.001)
     , m_BValue(1000)
 {
 
 }
 
 template< class ScalarType >
 StickModel< ScalarType >::~StickModel()
 {
 
 }
 
 template< class ScalarType >
 ScalarType StickModel< ScalarType >::SimulateMeasurement(unsigned int dir)
 {
     ScalarType signal = 0;
 
     if (dir>=this->m_GradientList.size())
         return signal;
 
     this->m_FiberDirection.Normalize();
 
     GradientType g = this->m_GradientList[dir];
-    ScalarType bVal = g.GetNorm(); bVal *= bVal;
+    float norm = g.GetNorm();
+    ScalarType bVal = norm*norm;
 
     if (bVal>0.0001)
     {
-        ScalarType dot = this->m_FiberDirection*g;
+        ScalarType dot = this->m_FiberDirection*g/norm;
         signal = std::exp( -m_BValue * bVal * m_Diffusivity*dot*dot );
     }
     else
         signal = 1;
 
     return signal;
 }
 
 template< class ScalarType >
 typename StickModel< ScalarType >::PixelType StickModel< ScalarType >::SimulateMeasurement()
 {
     this->m_FiberDirection.Normalize();
     PixelType signal;
     signal.SetSize(this->m_GradientList.size());
 
     for( unsigned int i=0; i<this->m_GradientList.size(); i++)
     {
         GradientType g = this->m_GradientList[i];
-        ScalarType bVal = g.GetNorm(); bVal *= bVal;
+        float norm = g.GetNorm();
+        ScalarType bVal = norm*norm;
 
         if (bVal>0.0001)
         {
-            ScalarType dot = this->m_FiberDirection*g;
+            ScalarType dot = this->m_FiberDirection*g/norm;
             signal[i] = std::exp( -m_BValue * bVal * m_Diffusivity*dot*dot );
         }
         else
             signal[i] = 1;
     }
 
     return signal;
 }
diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkTensorModel.cpp b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkTensorModel.cpp
index 610d39661e..57a54176bc 100644
--- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkTensorModel.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkTensorModel.cpp
@@ -1,129 +1,130 @@
 /*===================================================================
 
 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 <vnl/vnl_cross.h>
 #include <vnl/vnl_quaternion.h>
 #include <mitkTensorModel.h>
 
 using namespace mitk;
 
 template< class ScalarType >
 TensorModel< ScalarType >::TensorModel()
     : m_BValue(1000)
 {
     m_KernelDirection[0]=1; m_KernelDirection[1]=0; m_KernelDirection[2]=0;
     m_KernelTensorMatrix.fill(0.0);
     m_KernelTensorMatrix[0][0] = 0.002;
     m_KernelTensorMatrix[1][1] = 0.0005;
     m_KernelTensorMatrix[2][2] = 0.0005;
 }
 
 template< class ScalarType >
 TensorModel< ScalarType >::~TensorModel()
 {
 
 }
 
 template< class ScalarType >
 ScalarType TensorModel< ScalarType >::SimulateMeasurement(unsigned int dir)
 {
     ScalarType signal = 0;
 
     if (dir>=this->m_GradientList.size())
         return signal;
 
     ItkTensorType tensor; tensor.Fill(0.0);
     this->m_FiberDirection.Normalize();
     vnl_vector_fixed<double, 3> axis = itk::CrossProduct(m_KernelDirection, this->m_FiberDirection).GetVnlVector(); axis.normalize();
     vnl_quaternion<double> rotation(axis, acos(m_KernelDirection*this->m_FiberDirection));
     rotation.normalize();
     vnl_matrix_fixed<double, 3, 3> matrix = rotation.rotation_matrix_transpose();
 
     vnl_matrix_fixed<double, 3, 3> tensorMatrix = matrix.transpose()*m_KernelTensorMatrix*matrix;
     tensor[0] = tensorMatrix[0][0]; tensor[1] = tensorMatrix[0][1]; tensor[2] = tensorMatrix[0][2];
     tensor[3] = tensorMatrix[1][1]; tensor[4] = tensorMatrix[1][2]; tensor[5] = tensorMatrix[2][2];
 
     GradientType g = this->m_GradientList[dir];
-    ScalarType bVal = g.GetNorm(); bVal *= bVal;
+    float norm = g.GetNorm();
+    ScalarType bVal = norm*norm;
 
     if (bVal>0.0001)
     {
         itk::DiffusionTensor3D< ScalarType > S;
-        S[0] = g[0]*g[0];
-        S[1] = g[1]*g[0];
-        S[2] = g[2]*g[0];
-        S[3] = g[1]*g[1];
-        S[4] = g[2]*g[1];
-        S[5] = g[2]*g[2];
+        S[0] = g[0]*g[0]/bVal;
+        S[1] = g[1]*g[0]/bVal;
+        S[2] = g[2]*g[0]/bVal;
+        S[3] = g[1]*g[1]/bVal;
+        S[4] = g[2]*g[1]/bVal;
+        S[5] = g[2]*g[2]/bVal;
 
         ScalarType D = tensor[0]*S[0] + tensor[1]*S[1] + tensor[2]*S[2] +
                        tensor[1]*S[1] + tensor[3]*S[3] + tensor[4]*S[4] +
                        tensor[2]*S[2] + tensor[4]*S[4] + tensor[5]*S[5];
 
         // check for corrupted tensor and generate signal
         if (D>=0)
             signal = std::exp ( -m_BValue * bVal * D );
     }
     else
         signal = 1;
 
     return signal;
 }
 
 template< class ScalarType >
 typename TensorModel< ScalarType >::PixelType TensorModel< ScalarType >::SimulateMeasurement()
 {
     PixelType signal; signal.SetSize(this->m_GradientList.size()); signal.Fill(0.0);
 
     ItkTensorType tensor; tensor.Fill(0.0);
     this->m_FiberDirection.Normalize();
     vnl_vector_fixed<double, 3> axis = itk::CrossProduct(m_KernelDirection, this->m_FiberDirection).GetVnlVector(); axis.normalize();
     vnl_quaternion<double> rotation(axis, acos(m_KernelDirection*this->m_FiberDirection));
     rotation.normalize();
     vnl_matrix_fixed<double, 3, 3> matrix = rotation.rotation_matrix_transpose();
 
     vnl_matrix_fixed<double, 3, 3> tensorMatrix = matrix.transpose()*m_KernelTensorMatrix*matrix;
     tensor[0] = tensorMatrix[0][0]; tensor[1] = tensorMatrix[0][1]; tensor[2] = tensorMatrix[0][2];
     tensor[3] = tensorMatrix[1][1]; tensor[4] = tensorMatrix[1][2]; tensor[5] = tensorMatrix[2][2];
 
     for( unsigned int i=0; i<this->m_GradientList.size(); i++)
     {
         GradientType g = this->m_GradientList[i];
         ScalarType bVal = g.GetNorm(); bVal *= bVal;
 
         if (bVal>0.0001)
         {
             itk::DiffusionTensor3D< ScalarType > S;
-            S[0] = g[0]*g[0];
-            S[1] = g[1]*g[0];
-            S[2] = g[2]*g[0];
-            S[3] = g[1]*g[1];
-            S[4] = g[2]*g[1];
-            S[5] = g[2]*g[2];
+            S[0] = g[0]*g[0]/bVal;
+            S[1] = g[1]*g[0]/bVal;
+            S[2] = g[2]*g[0]/bVal;
+            S[3] = g[1]*g[1]/bVal;
+            S[4] = g[2]*g[1]/bVal;
+            S[5] = g[2]*g[2]/bVal;
 
             ScalarType D = tensor[0]*S[0] + tensor[1]*S[1] + tensor[2]*S[2] +
                            tensor[1]*S[1] + tensor[3]*S[3] + tensor[4]*S[4] +
                            tensor[2]*S[2] + tensor[4]*S[4] + tensor[5]*S[5];
 
             // check for corrupted tensor and generate signal
             if (D>=0)
                 signal[i] = std::exp ( -m_BValue * bVal * D );
         }
         else
             signal[i] = 1;
     }
 
     return signal;
 }