diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkOrientationDistributionFunction.txx b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkOrientationDistributionFunction.txx
index 6c3a17ed2c..8002f2d6db 100644
--- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkOrientationDistributionFunction.txx
+++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkOrientationDistributionFunction.txx
@@ -1,1220 +1,1220 @@
 /*===================================================================
 
 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 _itkOrientationDistributionFunction_txx
 #define _itkOrientationDistributionFunction_txx
 
 #include <stdio.h>
 #include <stdlib.h>
 #include <time.h>
 #include <algorithm>
 #include <vector>
 #include <vnl/algo/vnl_matrix_inverse.h>
 #include "itkPointShell.h"
 
 #define _USE_MATH_DEFINES
 #include <math.h>
 
 namespace itk
 {
 
   template<class T, unsigned int N>
   vtkPolyData* itk::OrientationDistributionFunction<T,N>::m_BaseMesh = NULL;
 
   template<class T, unsigned int N>
   double itk::OrientationDistributionFunction<T,N>::m_MaxChordLength = -1.0;
 
   template<class T, unsigned int N>
   vnl_matrix_fixed<double, 3, N>* itk::OrientationDistributionFunction<T,N>::m_Directions
     = itk::PointShell<N, vnl_matrix_fixed<double, 3, N> >::DistributePointShell();
 
   template<class T, unsigned int N>
   std::vector< std::vector<int>* >* itk::OrientationDistributionFunction<T,N>::m_NeighborIdxs = NULL;
 
   template<class T, unsigned int N>
   std::vector< std::vector<int>* >* itk::OrientationDistributionFunction<T,N>::m_AngularRangeIdxs = NULL;
 
   template<class T, unsigned int N>
   std::vector<int>* itk::OrientationDistributionFunction<T,N>::m_HalfSphereIdxs = NULL;
 
   template<class T, unsigned int N>
   itk::SimpleFastMutexLock itk::OrientationDistributionFunction<T,N>::m_MutexBaseMesh;
   template<class T, unsigned int N>
   itk::SimpleFastMutexLock itk::OrientationDistributionFunction<T,N>::m_MutexHalfSphereIdxs;
   template<class T, unsigned int N>
   itk::SimpleFastMutexLock itk::OrientationDistributionFunction<T,N>::m_MutexNeighbors;
   template<class T, unsigned int N>
   itk::SimpleFastMutexLock itk::OrientationDistributionFunction<T,N>::m_MutexAngularRange;
 
 
 #define ODF_PI       M_PI
 
   /**
   * Assignment Operator
   */
   template<class T, unsigned int NOdfDirections>
   OrientationDistributionFunction<T, NOdfDirections>&
     OrientationDistributionFunction<T, NOdfDirections>
     ::operator= (const Self& r)
   {
     BaseArray::operator=(r);
     return *this;
   }
 
   /**
   * Assignment Operator from a scalar constant
   */
   template<class T, unsigned int NOdfDirections>
   OrientationDistributionFunction<T, NOdfDirections>&
     OrientationDistributionFunction<T, NOdfDirections>
     ::operator= (const ComponentType & r)
   {
     BaseArray::operator=(&r);
     return *this;
   }
 
   /**
   * Assigment from a plain array
   */
   template<class T, unsigned int NOdfDirections>
   OrientationDistributionFunction<T, NOdfDirections>&
     OrientationDistributionFunction<T, NOdfDirections>
     ::operator= (const ComponentArrayType r )
   {
     BaseArray::operator=(r);
     return *this;
   }
 
   /**
   * Returns a temporary copy of a vector
   */
   template<class T, unsigned int NOdfDirections>
   OrientationDistributionFunction<T, NOdfDirections>
     OrientationDistributionFunction<T, NOdfDirections>
     ::operator+(const Self & r) const
   {
     Self result;
     for( unsigned int i=0; i<InternalDimension; i++)
     {
       result[i] = (*this)[i] + r[i];
     }
     return result;
   }
 
   /**
   * Returns a temporary copy of a vector
   */
   template<class T, unsigned int NOdfDirections>
   OrientationDistributionFunction<T, NOdfDirections>
     OrientationDistributionFunction<T, NOdfDirections>
     ::operator-(const Self & r) const
   {
     Self result;
     for( unsigned int i=0; i<InternalDimension; i++)
     {
       result[i] = (*this)[i] - r[i];
     }
     return result;
   }
 
   /**
   * Performs addition in place
   */
   template<class T, unsigned int NOdfDirections>
   const OrientationDistributionFunction<T, NOdfDirections> &
     OrientationDistributionFunction<T, NOdfDirections>
     ::operator+=(const Self & r)
   {
     for( unsigned int i=0; i<InternalDimension; i++)
     {
       (*this)[i] += r[i];
     }
     return *this;
   }
 
   /**
   * Performs subtraction in place
   */
   template<class T, unsigned int NOdfDirections>
   const OrientationDistributionFunction<T, NOdfDirections> &
     OrientationDistributionFunction<T, NOdfDirections>
     ::operator-=(const Self & r)
   {
     for( unsigned int i=0; i<InternalDimension; i++)
     {
       (*this)[i] -= r[i];
     }
     return *this;
   }
 
   /**
   * Performs multiplication by a scalar, in place
   */
   template<class T, unsigned int NOdfDirections>
   const OrientationDistributionFunction<T, NOdfDirections> &
     OrientationDistributionFunction<T, NOdfDirections>
     ::operator*=(const RealValueType & r)
   {
     for( unsigned int i=0; i<InternalDimension; i++)
     {
       (*this)[i] *= r;
     }
     return *this;
   }
 
   /**
   * Performs division by a scalar, in place
   */
   template<class T, unsigned int NOdfDirections>
   const OrientationDistributionFunction<T, NOdfDirections> &
     OrientationDistributionFunction<T, NOdfDirections>
     ::operator/=(const RealValueType & r)
   {
     for( unsigned int i=0; i<InternalDimension; i++)
     {
       (*this)[i] /= r;
     }
     return *this;
   }
 
   /**
   * Performs multiplication with a scalar
   */
   template<class T, unsigned int NOdfDirections>
   OrientationDistributionFunction<T, NOdfDirections>
     OrientationDistributionFunction<T, NOdfDirections>
     ::operator*(const RealValueType & r) const
   {
     Self result;
     for( unsigned int i=0; i<InternalDimension; i++)
     {
       result[i] = (*this)[i] * r;
     }
     return result;
   }
 
   /**
   * Performs division by a scalar
   */
   template<class T, unsigned int NOdfDirections>
   OrientationDistributionFunction<T, NOdfDirections>
     OrientationDistributionFunction<T, NOdfDirections>
     ::operator/(const RealValueType & r) const
   {
     Self result;
     for( unsigned int i=0; i<InternalDimension; i++)
     {
       result[i] = (*this)[i] / r;
     }
     return result;
   }
 
   /**
   * Matrix notation access to elements
   */
   template<class T, unsigned int NOdfDirections>
   const typename OrientationDistributionFunction<T, NOdfDirections>::ValueType &
     OrientationDistributionFunction<T, NOdfDirections>
     ::operator()(unsigned int row, unsigned int col) const
   {
     unsigned int k;
 
     if( row < col )
     {
       k = row * InternalDimension + col - row * ( row + 1 ) / 2;
     }
     else
     {
       k = col * InternalDimension + row - col * ( col + 1 ) / 2;
     }
 
 
     if( k >= InternalDimension )
     {
       k = 0;
     }
 
     return (*this)[k];
   }
 
   /**
   * Matrix notation access to elements
   */
   template<class T, unsigned int NOdfDirections>
   typename OrientationDistributionFunction<T, NOdfDirections>::ValueType &
     OrientationDistributionFunction<T, NOdfDirections>
     ::operator()(unsigned int row, unsigned int col)
   {
     unsigned int k;
 
     if( row < col )
     {
       k = row * InternalDimension + col - row * ( row + 1 ) / 2;
     }
     else
     {
       k = col * InternalDimension + row - col * ( col + 1 ) / 2;
     }
 
 
     if( k >= InternalDimension )
     {
       k = 0;
     }
 
     return (*this)[k];
   }
 
   /**
   * Set the Tensor to an Identity.
   * Set ones in the diagonal and zeroes every where else.
   */
   template<class T, unsigned int NOdfDirections>
   void
     OrientationDistributionFunction<T, NOdfDirections>
     ::SetIsotropic()
   {
     this->Fill(NumericTraits< T >::One / NOdfDirections);
   }
 
   /**
   * Set the Tensor to an Identity.
   * Set ones in the diagonal and zeroes every where else.
   */
   template<class T, unsigned int NOdfDirections>
   void
     OrientationDistributionFunction<T, NOdfDirections>
     ::InitFromTensor(itk::DiffusionTensor3D<T> tensor)
   {
     for(unsigned int i=0; i<NOdfDirections; i++)
     {
       /*
       *               | t0  t1  t2 |    g0
       *  g0 g1 g2  *  | t1  t3  t4 | *  g1
       *               | t2  t4  t5 |    g2
       *
       * =   g0 * (t0g0*t1g1*t2g2)
       *   + g1 * (t1g0+t3g1+t4g2)
       *   + g2 * (t2g0+t4g1+t5g2)
       */
       T g0 = (*m_Directions)(0,i);
       T g1 = (*m_Directions)(1,i);
       T g2 = (*m_Directions)(2,i);
       T t0 = tensor[0];
       T t1 = tensor[1];
       T t2 = tensor[2];
       T t3 = tensor[3];
       T t4 = tensor[4];
       T t5 = tensor[5];
       (*this)[i] = g0 * (t0*g0+t1*g1+t2*g2)
         + g1 * (t1*g0+t3*g1+t4*g2)
         + g2 * (t2*g0+t4*g1+t5*g2);
 
       if ((*this)[i]<0 || (*this)[i]!=(*this)[i])
         (*this)[i] = 0;
     }
   }
 
   /**
   * L2-Normalization
   */
   template<class T, unsigned int NOdfDirections>
   void
     OrientationDistributionFunction<T, NOdfDirections>
     ::L2Normalize()
   {
     T sum = 0;
     for( unsigned int i=0; i<InternalDimension; i++)
     {
       sum += (*this)[i]*(*this)[i];
     }
     sum = std::sqrt(sum);
     for( unsigned int i=0; i<InternalDimension; i++)
     {
       (*this)[i] = (*this)[i] / sum;
     }
   }
 
   /**
   * Normalization to PDF
   */
   template<class T, unsigned int NOdfDirections>
   void
     OrientationDistributionFunction<T, NOdfDirections>
     ::Normalize()
   {
     T sum = 0;
     for( unsigned int i=0; i<InternalDimension; i++)
     {
       sum += (*this)[i];
     }
     if (sum>0)
     {
       for( unsigned int i=0; i<InternalDimension; i++)
       {
         (*this)[i] = (*this)[i] / sum;
       }
     }
   }
 
   /**
   * Min/Max-Normalization
   */
   template<class T, unsigned int NOdfDirections>
   OrientationDistributionFunction<T, NOdfDirections>
     OrientationDistributionFunction<T, NOdfDirections>
     ::MinMaxNormalize() const
   {
     T max = NumericTraits<T>::NonpositiveMin();
     T min = NumericTraits<T>::max();
     for( unsigned int i=0; i<InternalDimension; i++)
     {
       max = (*this)[i] > max ? (*this)[i] : max;
       min = (*this)[i] < min ? (*this)[i] : min;
     }
     Self retval;
     for( unsigned int i=0; i<InternalDimension; i++)
     {
       retval[i] = ((*this)[i] - min) / (max - min);
     }
     return retval;
   }
 
   /**
   * Max-Normalization
   */
   template<class T, unsigned int NOdfDirections>
   OrientationDistributionFunction<T, NOdfDirections>
     OrientationDistributionFunction<T, NOdfDirections>
     ::MaxNormalize() const
   {
     T max = NumericTraits<T>::NonpositiveMin();
     for( unsigned int i=0; i<InternalDimension; i++)
     {
       max = (*this)[i] > max ? (*this)[i] : max;
     }
     Self retval;
     for( unsigned int i=0; i<InternalDimension; i++)
     {
       retval[i] = (*this)[i] / max;
     }
     return retval;
   }
 
   template<class T, unsigned int NOdfDirections>
   T
     OrientationDistributionFunction<T, NOdfDirections>
     ::GetMaxValue() const
   {
     T max = NumericTraits<T>::NonpositiveMin();
     for( unsigned int i=0; i<InternalDimension; i++)
     {
       if((*this)[i] >= max )
       {
         max = (*this)[i];
       }
     }
     return max;
   }
 
   template<class T, unsigned int NOdfDirections>
   T
     OrientationDistributionFunction<T, NOdfDirections>
     ::GetMinValue() const
   {
     T min = NumericTraits<T>::max();
     for( unsigned int i=0; i<InternalDimension; i++)
     {
       if((*this)[i] >= min )
       {
         min = (*this)[i];
       }
     }
     return min;
   }
 
   template<class T, unsigned int NOdfDirections>
   T
     OrientationDistributionFunction<T, NOdfDirections>
     ::GetMeanValue() const
   {
     T sum = 0;
     for( unsigned int i=0; i<InternalDimension; i++)
       sum += (*this)[i];
     return sum/InternalDimension;
   }
 
   template<class T, unsigned int NOdfDirections>
   double
     OrientationDistributionFunction<T, NOdfDirections>
     ::GetMaxChordLength()
   {
     if(m_MaxChordLength<0.0)
     {
       ComputeBaseMesh();
       double max_dist = -1;
       vtkPoints* points = m_BaseMesh->GetPoints();
       for(int i=0; i<NOdfDirections; i++)
       {
         double p[3];
         points->GetPoint(i,p);
         std::vector<int> neighbors = GetNeighbors(i);
-        for(int j=0; j<neighbors.size(); j++)
+        for(std::size_type j=0; j<neighbors.size(); j++)
         {
           double n[3];
           points->GetPoint(neighbors[j],n);
           double d = sqrt(
             (p[0]-n[0])*(p[0]-n[0]) +
             (p[1]-n[1])*(p[1]-n[1]) +
             (p[2]-n[2])*(p[2]-n[2]));
           max_dist = d>max_dist ? d : max_dist;
         }
       }
       m_MaxChordLength = max_dist;
     }
     return m_MaxChordLength;
   }
 
   template<class T, unsigned int NOdfDirections>
   void
     OrientationDistributionFunction<T, NOdfDirections>
     ::ComputeBaseMesh()
   {
 
     m_MutexBaseMesh.Lock();
     if(m_BaseMesh == NULL)
     {
 
       vtkPoints* points = vtkPoints::New();
       for(unsigned int j=0; j<NOdfDirections; j++){
         double x = (*m_Directions)(0,j);
         double y = (*m_Directions)(1,j);
         double z = (*m_Directions)(2,j);
         double az = atan2(y,x);
         double elev = atan2(z,sqrt(x*x+y*y));
         double r = sqrt(x*x+y*y+z*z);
         points->InsertNextPoint(az,elev,r);
       }
 
       vtkPolyData* polydata = vtkPolyData::New();
       polydata->SetPoints( points );
       vtkDelaunay2D *delaunay = vtkDelaunay2D::New();
       delaunay->SetInputData( polydata );
       delaunay->Update();
 
       vtkCellArray* vtkpolys = delaunay->GetOutput()->GetPolys();
       vtkCellArray* vtknewpolys = vtkCellArray::New();
       vtkIdType npts; vtkIdType *pts;
       while(vtkpolys->GetNextCell(npts,pts))
       {
         bool insert = true;
         for(int i=0; i<npts; i++)
         {
           double *tmpPoint = points->GetPoint(pts[i]);
           double az = tmpPoint[0];
           double elev = tmpPoint[1];
           if((abs(az)>ODF_PI-0.5) || (abs(elev)>ODF_PI/2-0.5))
             insert = false;
         }
         if(insert)
           vtknewpolys->InsertNextCell(npts, pts);
       }
 
       vtkPoints* points2 = vtkPoints::New();
       for(unsigned int j=0; j<NOdfDirections; j++){
         double x = -(*m_Directions)(0,j);
         double y = -(*m_Directions)(2,j);
         double z = -(*m_Directions)(1,j);
         double az = atan2(y,x);
         double elev = atan2(z,sqrt(x*x+y*y));
         double r = sqrt(x*x+y*y+z*z);
         points2->InsertNextPoint(az,elev,r);
       }
 
       vtkPolyData* polydata2 = vtkPolyData::New();
       polydata2->SetPoints( points2 );
       vtkDelaunay2D *delaunay2 = vtkDelaunay2D::New();
       delaunay2->SetInputData( polydata2 );
       delaunay2->Update();
 
       vtkpolys = delaunay2->GetOutput()->GetPolys();
       while(vtkpolys->GetNextCell(npts,pts))
       {
         bool insert = true;
         for(int i=0; i<npts; i++)
         {
           double *tmpPoint = points2->GetPoint(pts[i]);
           double az = tmpPoint[0];
           double elev = tmpPoint[1];
           if((abs(az)>ODF_PI-0.5) || (abs(elev)>ODF_PI/2-0.5))
             insert = false;
         }
         if(insert)
           vtknewpolys->InsertNextCell(npts, pts);
       }
 
       polydata->SetPolys(vtknewpolys);
 
       for (vtkIdType p = 0; p < NOdfDirections; p++)
       {
         points->SetPoint(p,m_Directions->get_column(p).data_block());
       }
       polydata->SetPoints( points );
 
       m_BaseMesh = polydata;
     }
     m_MutexBaseMesh.Unlock();
   }
 
   /**
   * Extract the principle diffusion direction
   */
   template<class T, unsigned int NOdfDirections>
   int
     OrientationDistributionFunction<T, NOdfDirections>
     ::GetPrincipleDiffusionDirection() const
   {
     T max = NumericTraits<T>::NonpositiveMin();
     int maxidx = -1;
     for( unsigned int i=0; i<InternalDimension; i++)
     {
       if((*this)[i] >= max )
       {
         max = (*this)[i];
         maxidx = i;
       }
     }
 
     return maxidx;
   }
 
   template<class T, unsigned int NOdfDirections>
   std::vector<int>
     OrientationDistributionFunction<T, NOdfDirections>
     ::GetNeighbors(int idx)
   {
     ComputeBaseMesh();
 
     m_MutexNeighbors.Lock();
     if(m_NeighborIdxs == NULL)
     {
       m_NeighborIdxs = new std::vector< std::vector<int>* >();
       vtkCellArray* polys = m_BaseMesh->GetPolys();
 
       for(unsigned int i=0; i<NOdfDirections; i++)
       {
         std::vector<int> *idxs = new std::vector<int>();
         polys->InitTraversal();
         vtkIdType npts; vtkIdType *pts;
         while(polys->GetNextCell(npts,pts))
         {
           if( pts[0] == i )
           {
             idxs->push_back(pts[1]);
             idxs->push_back(pts[2]);
           }
           else if( pts[1] == i )
           {
             idxs->push_back(pts[0]);
             idxs->push_back(pts[2]);
           }
           else if( pts[2] == i )
           {
             idxs->push_back(pts[0]);
             idxs->push_back(pts[1]);
           }
         }
         std::sort(idxs->begin(), idxs->end());
         std::vector< int >::iterator endLocation;
         endLocation = std::unique( idxs->begin(), idxs->end() );
         idxs->erase(endLocation, idxs->end());
         m_NeighborIdxs->push_back(idxs);
       }
     }
     m_MutexNeighbors.Unlock();
 
     return *m_NeighborIdxs->at(idx);
   }
 
   /**
   * Extract the n-th diffusion direction
   */
   template<class T, unsigned int NOdfDirections>
   int
     OrientationDistributionFunction<T, NOdfDirections>
     ::GetNthDiffusionDirection(int n, vnl_vector_fixed<double,3> rndVec) const
   {
 
     if( n == 0 )
       return GetPrincipleDiffusionDirection();
 
     m_MutexHalfSphereIdxs.Lock();
     if( !m_HalfSphereIdxs )
     {
       m_HalfSphereIdxs = new std::vector<int>();
       for( unsigned int i=0; i<InternalDimension; i++)
       {
         if(dot_product(m_Directions->get_column(i),rndVec) > 0.0)
         {
           m_HalfSphereIdxs->push_back(i);
         }
       }
     }
     m_MutexHalfSphereIdxs.Unlock();
 
     // collect indices of directions
     // that are local maxima
     std::vector<int> localMaxima;
     std::vector<int>::iterator it;
     for( it=m_HalfSphereIdxs->begin();
       it!=m_HalfSphereIdxs->end();
       it++)
     {
       std::vector<int> nbs = GetNeighbors(*it);
       std::vector<int>::iterator it2;
       bool max = true;
       for(it2 = nbs.begin();
         it2 != nbs.end();
         it2++)
       {
         if((*this)[*it2] > (*this)[*it])
         {
           max = false;
           break;
         }
       }
       if(max)
         localMaxima.push_back(*it);
     }
 
     // delete n highest local maxima from list
     // and return remaining highest
     int maxidx = -1;
     std::vector<int>::iterator itMax;
     for( int i=0; i<=n; i++ )
     {
       maxidx = -1;
       T max = NumericTraits<T>::NonpositiveMin();
       for(it = localMaxima.begin();
         it != localMaxima.end();
         it++)
       {
         if((*this)[*it]>max)
         {
           max = (*this)[*it];
           maxidx = *it;
         }
       }
       it = find(localMaxima.begin(), localMaxima.end(), maxidx);
       if(it!=localMaxima.end())
         localMaxima.erase(it);
     }
 
     return maxidx;
   }
 
   template < typename TComponent, unsigned int NOdfDirections >
   vnl_vector_fixed<double,3> itk::OrientationDistributionFunction<TComponent, NOdfDirections>
     ::GetDirection( int i )
   {
     return m_Directions->get_column(i);
   }
 
   /**
   * Interpolate a position between sampled directions
   */
   template<class T, unsigned int NOdfDirections>
   T
     OrientationDistributionFunction<T, NOdfDirections>
     ::GetInterpolatedComponent(vnl_vector_fixed<double,3> dir, InterpolationMethods method) const
   {
 
     ComputeBaseMesh();
     double retval = -1.0;
 
     switch(method)
     {
     case ODF_NEAREST_NEIGHBOR_INTERP:
       {
 
         vtkPoints* points = m_BaseMesh->GetPoints();
         double current_min = NumericTraits<double>::max();
         int current_min_idx = -1;
         for(int i=0; i<NOdfDirections; i++)
         {
           vnl_vector_fixed<double,3> P(points->GetPoint(i));
           double dist = (dir-P).two_norm();
           current_min_idx = dist<current_min ? i : current_min_idx;
           current_min = dist<current_min ? dist : current_min;
         }
         retval = this->GetNthComponent(current_min_idx);
         break;
 
       }
     case ODF_TRILINEAR_BARYCENTRIC_INTERP:
       {
 
         double maxChordLength = GetMaxChordLength();
         vtkCellArray* polys = m_BaseMesh->GetPolys();
         vtkPoints* points = m_BaseMesh->GetPoints();
         vtkIdType npts; vtkIdType *pts;
         double current_min = NumericTraits<double>::max();
         polys->InitTraversal();
         while(polys->GetNextCell(npts,pts))
         {
           vnl_vector_fixed<double,3> A(points->GetPoint(pts[0]));
           vnl_vector_fixed<double,3> B(points->GetPoint(pts[1]));
           vnl_vector_fixed<double,3> C(points->GetPoint(pts[2]));
 
           vnl_vector_fixed<double,3> d1;
           d1.put(0,(dir-A).two_norm());
           d1.put(1,(dir-B).two_norm());
           d1.put(2,(dir-C).two_norm());
           double maxval = d1.max_value();
           if(maxval>maxChordLength)
           {
             continue;
           }
 
           // Compute vectors
           vnl_vector_fixed<double,3> v0 = C - A;
           vnl_vector_fixed<double,3> v1 = B - A;
 
           // Project direction to plane ABC
           vnl_vector_fixed<double,3> v6 = dir;
           vnl_vector_fixed<double,3> cross = vnl_cross_3d(v0, v1);
           cross = cross.normalize();
           vtkPlane::ProjectPoint(v6.data_block(),A.data_block(),cross.data_block(),v6.data_block());
           v6 = v6-A;
 
           // Calculate barycentric coords
           vnl_matrix_fixed<double,3,2> mat;
           mat.set_column(0, v0);
           mat.set_column(1, v1);
           vnl_matrix_inverse<double> inv(mat);
           vnl_matrix_fixed<double,2,3> inver = inv.pinverse();
           vnl_vector<double> uv = inv.pinverse()*v6;
 
           // Check if point is in triangle
           double eps = 0.01;
           if( (uv(0) >= 0-eps) && (uv(1) >= 0-eps) && (uv(0) + uv(1) <= 1+eps) )
           {
             // check if minimum angle is the max so far
             if(d1.two_norm() < current_min)
             {
               current_min = d1.two_norm();
               vnl_vector<double> barycentricCoords(3);
               barycentricCoords[2] = uv[0]<0 ? 0 : (uv[0]>1?1:uv[0]);
               barycentricCoords[1] = uv[1]<0 ? 0 : (uv[1]>1?1:uv[1]);
               barycentricCoords[0] = 1-(barycentricCoords[1]+barycentricCoords[2]);
               retval =  barycentricCoords[0]*this->GetNthComponent(pts[0]) +
                 barycentricCoords[1]*this->GetNthComponent(pts[1]) +
                 barycentricCoords[2]*this->GetNthComponent(pts[2]);
             }
           }
         }
 
         break;
       }
     case ODF_SPHERICAL_GAUSSIAN_BASIS_FUNCTIONS:
       {
         double maxChordLength = GetMaxChordLength();
         double sigma = asin(maxChordLength/2);
 
         // this is the contribution of each kernel to each sampling point on the
         // equator
         vnl_vector<double> contrib;
         contrib.set_size(NOdfDirections);
 
         vtkPoints* points = m_BaseMesh->GetPoints();
         double sum = 0;
         for(int i=0; i<NOdfDirections; i++)
         {
           vnl_vector_fixed<double,3> P(points->GetPoint(i));
           double stv =  dir[0]*P[0]
           + dir[1]*P[1]
           + dir[2]*P[2];
           stv = (stv<-1.0) ? -1.0 : ( (stv>1.0) ? 1.0 : stv);
           double x = acos(stv);
           contrib[i] = (1.0/(sigma*sqrt(2.0*ODF_PI)))
             *exp((-x*x)/(2*sigma*sigma));
           sum += contrib[i];
         }
 
         retval = 0;
         for(int i=0; i<NOdfDirections; i++)
         {
           retval += (contrib[i] / sum)*this->GetNthComponent(i);
         }
         break;
       }
 
     }
 
     if(retval==-1)
     {
       std::cout << "Interpolation failed" << std::endl;
       return 0;
     }
 
     return retval;
 
   }
 
   /**
   * Calculate Generalized Fractional Anisotropy
   */
   template<class T, unsigned int NOdfDirections>
   T
     OrientationDistributionFunction<T, NOdfDirections>
     ::GetGeneralizedFractionalAnisotropy() const
   {
     double mean = 0;
     double std = 0;
     double rms = 0;
     for( unsigned int i=0; i<InternalDimension; i++)
     {
       T val = (*this)[i];
       mean += val;
     }
     mean /= NOdfDirections;
     for( unsigned int i=0; i<InternalDimension; i++)
     {
       T val = (*this)[i];
       std += (val - mean) * (val - mean);
       rms += val*val;
     }
     std *= NOdfDirections;
     rms *= NOdfDirections - 1;
 
     if(rms == 0)
     {
       return 0;
     }
     else
     {
       return sqrt(std/rms);
     }
   }
 
 
   template < typename T, unsigned int N>
   T itk::OrientationDistributionFunction<T, N>
     ::GetGeneralizedGFA( int k, int p ) const
   {
     double mean = 0;
     double std = 0;
     double rms = 0;
     double max = NumericTraits<double>::NonpositiveMin();
     for( unsigned int i=0; i<InternalDimension; i++)
     {
       double val = (double)(*this)[i];
       mean += pow(val,(double)p);
       max = val > max ? val : max;
     }
     max = pow(max,(double)p);
     mean /= N;
     for( unsigned int i=0; i<InternalDimension; i++)
     {
       double val = (double)(*this)[i];
       std += (pow(val,(double)p) - mean) * (pow(val,(double)p) - mean);
       if(k>0)
       {
         rms += pow(val,(double)(p*k));
       }
     }
     std /= N - 1;
     std = sqrt(std);
 
     if(k>0)
     {
       rms /= N;
       rms = pow(rms,(double)(1.0/k));
     }
     else if(k<0) // lim k->inf gives us the maximum
     {
       rms = max;
     }
     else // k==0 undefined, we define zeros root from 1 as 1
     {
       rms = 1;
     }
 
     if(rms == 0)
     {
       return 0;
     }
     else
     {
       return (T)(std/rms);
     }
   }
 
   /**
   * Calculate Nematic Order Parameter
   */
   template < typename T, unsigned int N >
   T itk::OrientationDistributionFunction<T, N>
     ::GetNematicOrderParameter() const
   {
     // not yet implemented
     return 0;
   }
 
   /**
   * Calculate StdDev by MaxValue
   */
   template < typename T, unsigned int N >
   T itk::OrientationDistributionFunction<T, N>
     ::GetStdDevByMaxValue() const
   {
     double mean = 0;
     double std = 0;
     T max = NumericTraits<T>::NonpositiveMin();
     for( unsigned int i=0; i<InternalDimension; i++)
     {
       T val = (*this)[i];
       mean += val;
       max = (*this)[i] > max ? (*this)[i] : max;
     }
     mean /= InternalDimension;
     for( unsigned int i=0; i<InternalDimension; i++)
     {
       T val = (*this)[i];
       std += (val - mean) * (val - mean);
     }
     std /= InternalDimension-1;
 
     if(max == 0)
     {
       return 0;
     }
     else
     {
       return (sqrt(std)/max);
     }
   }
 
   template < typename T, unsigned int N >
   T itk::OrientationDistributionFunction<T, N>
     ::GetPrincipleCurvature(double alphaMinDegree, double alphaMaxDegree, int invert) const
   {
     // following loop only performed once
     // (computing indices of each angular range)
     m_MutexAngularRange.Lock();
     if(m_AngularRangeIdxs == NULL)
     {
       m_AngularRangeIdxs = new std::vector< std::vector<int>* >();
       for(unsigned int i=0; i<N; i++)
       {
         vnl_vector_fixed<double,3> pDir = GetDirection(i);
         std::vector<int> *idxs = new std::vector<int>();
         for(unsigned int j=0; j<N; j++)
         {
           vnl_vector_fixed<double,3> cDir = GetDirection(j);
           double angle = ( 180 / ODF_PI ) * acos( dot_product(pDir, cDir) );
           if( (angle < alphaMaxDegree) && (angle > alphaMinDegree) )
           {
             idxs->push_back(j);
           }
         }
         m_AngularRangeIdxs->push_back(idxs);
       }
     }
     m_MutexAngularRange.Unlock();
 
     // find the maximum (or minimum) direction (remember index and value)
     T mode;
     int pIdx = -1;
     if(invert == 0)
     {
       pIdx = GetPrincipleDiffusionDirection();
       mode = (*this)[pIdx];
     }
     else
     {
       mode = NumericTraits<T>::max();
       for( unsigned int i=0; i<N; i++)
       {
         if((*this)[i] < mode )
         {
           mode = (*this)[i];
           pIdx = i;
         }
       }
     }
     //////////////////////
     //////// compute median of mode and its neighbors to become more stable to noise
     //////// compared to simply using the mode
     //////////////////////
 
     //////// values of mode and its neighbors
     //////std::vector<int> nbs = GetNeighbors(pIdx);
     //////std::vector<T> modeAndNeighborVals;
     //////modeAndNeighborVals.push_back(mode);
     //////int numNeighbors = nbs.size();
     //////for(int i=0; i<numNeighbors; i++)
     //////{
     //////  modeAndNeighborVals.push_back((*this)[nbs[i]]);
     //////}
 
     //////// sort by value
     //////std::sort( modeAndNeighborVals.begin(), modeAndNeighborVals.end() );
 
     //////// median of mode and neighbors
     //////mode = modeAndNeighborVals[floor(0.5*(double)(numNeighbors+1)+0.5)];
 
     ////////////////
     // computing a quantile of the angular range
     ////////////////
 
     // define quantile
     double quantile = 0.00;
 
     // collect all values in angular range of mode
     std::vector<T> odfValuesInAngularRange;
     int numInRange = m_AngularRangeIdxs->at(pIdx)->size();
     for(int i=0; i<numInRange; i++)
     {
       odfValuesInAngularRange.push_back((*this)[(*m_AngularRangeIdxs->at(pIdx))[i] ]);
     }
 
     // sort them by value
     std::sort( odfValuesInAngularRange.begin(), odfValuesInAngularRange.end() );
 
     // median of angular range
     T median = odfValuesInAngularRange[floor(quantile*(double)numInRange+0.5)];
 
     // compute and return final value
     if(mode > median)
     {
       return mode/median - 1.0;
     }
     else
     {
       return median/mode - 1.0;
     }
   }
 
   /**
   * Calculate Normalized Entropy
   */
   template < typename T, unsigned int N >
   T itk::OrientationDistributionFunction<T, N>
     ::GetNormalizedEntropy() const
   {
     double mean = 0;
     for( unsigned int i=0; i<InternalDimension; i++)
     {
       T val = (*this)[i];
       if( val != 0 )
       {
         val = log(val);
       }
       else
       {
         val = log(0.0000001);
       }
       mean += val;
     }
     double _n = (double) InternalDimension;
     mean /= _n;
     return (T) (-_n / log(_n) * mean);
   }
 
   /**
   * Pre-multiply the Tensor by a Matrix
   */
   template<class T, unsigned int NOdfDirections>
   OrientationDistributionFunction<T, NOdfDirections>
     OrientationDistributionFunction<T, NOdfDirections>
     ::PreMultiply( const MatrixType & m ) const
   {
     Self result;
     typedef typename NumericTraits<T>::AccumulateType  AccumulateType;
     for(unsigned int r=0; r<NOdfDirections; r++)
     {
       for(unsigned int c=0; c<NOdfDirections; c++)
       {
         AccumulateType sum = NumericTraits<AccumulateType>::ZeroValue();
         for(unsigned int t=0; t<NOdfDirections; t++)
         {
           sum += m(r,t) * (*this)(t,c);
         }
         result(r,c) = static_cast<T>( sum );
       }
     }
     return result;
   }
 
   /**
   * Post-multiply the Tensor by a Matrix
   */
   template<class T, unsigned int NOdfDirections>
   OrientationDistributionFunction<T, NOdfDirections>
     OrientationDistributionFunction<T, NOdfDirections>
     ::PostMultiply( const MatrixType & m ) const
   {
     Self result;
     typedef typename NumericTraits<T>::AccumulateType  AccumulateType;
     for(unsigned int r=0; r<NOdfDirections; r++)
     {
       for(unsigned int c=0; c<NOdfDirections; c++)
       {
         AccumulateType sum = NumericTraits<AccumulateType>::ZeroValue();
         for(unsigned int t=0; t<NOdfDirections; t++)
         {
           sum += (*this)(r,t) * m(t,c);
         }
         result(r,c) = static_cast<T>( sum );
       }
     }
     return result;
   }
 
   /**
   * Print content to an ostream
   */
   template<class T, unsigned int NOdfDirections>
   std::ostream &
     operator<<(std::ostream& os,const OrientationDistributionFunction<T, NOdfDirections> & c )
   {
     for(unsigned int i=0; i<c.GetNumberOfComponents(); i++)
     {
       os <<  static_cast<typename NumericTraits<T>::PrintType>(c[i]) << "  ";
     }
     return os;
   }
 
   /**
   * Read content from an istream
   */
   template<class T, unsigned int NOdfDirections>
   std::istream &
     operator>>(std::istream& is, OrientationDistributionFunction<T, NOdfDirections> & dt )
   {
     for(unsigned int i=0; i < dt.GetNumberOfComponents(); i++)
     {
       is >> dt[i];
     }
     return is;
   }
 
 
 
 } // end namespace itk
 
 #endif
diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidRegistrationMethodHelper.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidRegistrationMethodHelper.h
index 69b96e4475..771a8f6d28 100644
--- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidRegistrationMethodHelper.h
+++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidRegistrationMethodHelper.h
@@ -1,122 +1,122 @@
 /*===================================================================
 
 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 MITKPYRAMIDREGISTRATIONMETHODHELPER_H
 #define MITKPYRAMIDREGISTRATIONMETHODHELPER_H
 
 #include <DiffusionCoreExports.h>
 
 #include <itkCommand.h>
 
 #include <itkRegularStepGradientDescentOptimizer.h>
 #include <itkMattesMutualInformationImageToImageMetric.h>
 #include <itkNormalizedCorrelationImageToImageMetric.h>
 
 #include <itkAffineTransform.h>
 #include <itkEuler3DTransform.h>
 
 #include <itkMultiResolutionImageRegistrationMethod.h>
 #include <itkImageMomentsCalculator.h>
 
 #include "mitkImageAccessByItk.h"
 
 /**
  * @brief Provides same functionality as \a AccessTwoImagesFixedDimensionByItk for a subset of types
  *
  *  For now, the subset defined is only short and float.
  */
 #define AccessTwoImagesFixedDimensionTypeSubsetByItk(mitkImage1, mitkImage2, itkImageTypeFunction, dimension) \
 {                                                                                                   \
   const mitk::PixelType& pixelType1 = mitkImage1->GetPixelType();                                   \
   const mitk::PixelType& pixelType2 = mitkImage2->GetPixelType();                                   \
   const mitk::Image* constImage1 = mitkImage1;                                                      \
   const mitk::Image* constImage2 = mitkImage2;                                                      \
   mitk::Image* nonConstImage1 = const_cast<mitk::Image*>(constImage1);                              \
   mitk::Image* nonConstImage2 = const_cast<mitk::Image*>(constImage2);                              \
   nonConstImage1->Update();                                                                         \
   nonConstImage2->Update();                                                                         \
   _checkSpecificDimension(mitkImage1, (dimension));                                                 \
   _checkSpecificDimension(mitkImage2, (dimension));                                                 \
   _accessTwoImagesByItkForEach(itkImageTypeFunction, ((short, dimension))((unsigned short, dimension))((float, dimension))((double, dimension)), ((short, dimension))((unsigned short, dimension))((float, dimension))((double, dimension)) ) \
   {                                                                                                 \
     std::string msg("Pixel type ");                                                                 \
     msg.append(pixelType1.GetComponentTypeAsString() );                                             \
     msg.append(" or pixel type ");                                                                  \
     msg.append(pixelType2.GetComponentTypeAsString() );                                             \
     msg.append(" is not in " MITK_PP_STRINGIZE(MITK_ACCESSBYITK_TYPES_DIMN_SEQ(dimension)));        \
     throw mitk::AccessByItkException(msg);                                                          \
   }                                                                                                 \
 }
 
 
 /**
  * @brief The PyramidOptControlCommand class stears the step lenght of the
  * gradient descent optimizer in multi-scale registration methods
  */
 template <typename RegistrationType >
 class PyramidOptControlCommand : public itk::Command
 {
 public:
 
   typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
 
   itkNewMacro( PyramidOptControlCommand )
 
   void Execute(itk::Object *caller, const itk::EventObject & /*event*/)
   {
     RegistrationType* registration = dynamic_cast< RegistrationType* >( caller );
 
     MITK_DEBUG << "\t - Pyramid level " << registration->GetCurrentLevel();
     if( registration->GetCurrentLevel() == 0 )
       return;
 
     OptimizerType* optimizer = dynamic_cast< OptimizerType* >(registration->GetOptimizer());
 
     MITK_INFO << optimizer->GetStopConditionDescription()  << "\n"
                << optimizer->GetValue() << " : " << optimizer->GetCurrentPosition();
 
     optimizer->SetMaximumStepLength( optimizer->GetMaximumStepLength() * 0.25f );
     optimizer->SetMinimumStepLength( optimizer->GetMinimumStepLength() * 0.1f );
    // optimizer->SetNumberOfIterations( optimizer->GetNumberOfIterations() * 1.5f );
   }
 
   void Execute(const itk::Object * /*object*/, const itk::EventObject & /*event*/){}
 };
 
 
 template <typename OptimizerType>
 class OptimizerIterationCommand : public itk::Command
 {
 public:
   itkNewMacro( OptimizerIterationCommand )
 
-  void Execute(itk::Object *caller, const itk::EventObject & event)
+  void Execute(itk::Object *caller, const itk::EventObject & /*event*/)
   {
     OptimizerType* optimizer = dynamic_cast< OptimizerType* >( caller );
 
     unsigned int currentIter = optimizer->GetCurrentIteration();
     MITK_INFO << "[" << currentIter << "] : " << optimizer->GetValue() << " : " << optimizer->GetCurrentPosition();
 
   }
 
-  void Execute(const itk::Object * object, const itk::EventObject & event)
+  void Execute(const itk::Object * /*object*/, const itk::EventObject & /*event*/)
   {
 
   }
 };
 
 
 #endif // MITKPYRAMIDREGISTRATIONMETHODHELPER_H
diff --git a/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.txx b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.txx
index 5decb35c7b..3140f00cf3 100644
--- a/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.txx
+++ b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.txx
@@ -1,424 +1,424 @@
 /*===================================================================
 
 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 "itkImageRegionIterator.h"
 #include "itkImageRegionConstIterator.h"
 #include "mitkImageCast.h"
 
 template<typename TPixelType>
 mitk::DiffusionImage<TPixelType>::DiffusionImage()
   : m_VectorImage(0), m_Directions(0), m_OriginalDirections(0), m_B_Value(-1.0), m_VectorImageAdaptor(0)
 {
   MeasurementFrameType mf;
   for(int i=0; i<3; i++)
     for(int j=0; j<3; j++)
       mf[i][j] = 0;
   for(int i=0; i<3; i++)
     mf[i][i] = 1;
   m_MeasurementFrame = mf;
 }
 
 template<typename TPixelType>
 mitk::DiffusionImage<TPixelType>::~DiffusionImage()
 {
   // Remove Observer for m_Directions
   RemoveDirectionsContainerObserver();
 }
 
 template<typename TPixelType>
 void mitk::DiffusionImage<TPixelType>
 ::InitializeFromVectorImage()
 {
   if(!m_VectorImage || !m_Directions || m_B_Value==-1.0)
   {
     MITK_INFO << "DiffusionImage could not be initialized. Set all members first!" << std::endl;
     return;
   }
 
   // find bzero index
   int firstZeroIndex = -1;
   for(GradientDirectionContainerType::ConstIterator it = m_Directions->Begin();
       it != m_Directions->End(); ++it)
   {
     firstZeroIndex++;
     GradientDirectionType g = it.Value();
     if(g[0] == 0 && g[1] == 0 && g[2] == 0 )
       break;
   }
 
   typedef itk::Image<TPixelType,3> ImgType;
   typename ImgType::Pointer img = ImgType::New();
   img->SetSpacing( m_VectorImage->GetSpacing() );   // Set the image spacing
   img->SetOrigin( m_VectorImage->GetOrigin() );     // Set the image origin
   img->SetDirection( m_VectorImage->GetDirection() );  // Set the image direction
   img->SetLargestPossibleRegion( m_VectorImage->GetLargestPossibleRegion());
   img->SetBufferedRegion( m_VectorImage->GetLargestPossibleRegion() );
   img->Allocate();
 
   int vecLength = m_VectorImage->GetVectorLength();
   InitializeByItk( img.GetPointer(), 1, vecLength );
 
   itk::ImageRegionIterator<ImgType> itw (img, img->GetLargestPossibleRegion() );
   itw.GoToBegin();
 
   itk::ImageRegionConstIterator<ImageType> itr (m_VectorImage, m_VectorImage->GetLargestPossibleRegion() );
   itr.GoToBegin();
 
   while(!itr.IsAtEnd())
   {
     itw.Set(itr.Get().GetElement(firstZeroIndex));
     ++itr;
     ++itw;
   }
 
   // init
   SetImportVolume(img->GetBufferPointer());
 
   m_DisplayIndex = firstZeroIndex;
   MITK_INFO << "Diffusion-Image successfully initialized.";
 }
 
 template<typename TPixelType>
 void mitk::DiffusionImage<TPixelType>
 ::SetDisplayIndexForRendering(int displayIndex)
 {
   int index = displayIndex;
   int vecLength = m_VectorImage->GetVectorLength();
   index = index > vecLength-1 ? vecLength-1 : index;
   if( m_DisplayIndex != index )
   {
     typedef itk::Image<TPixelType,3> ImgType;
     typename ImgType::Pointer img = ImgType::New();
     CastToItkImage<ImgType>(this, img);
 
     itk::ImageRegionIterator<ImgType> itw (img, img->GetLargestPossibleRegion() );
     itw.GoToBegin();
 
     itk::ImageRegionConstIterator<ImageType> itr (m_VectorImage, m_VectorImage->GetLargestPossibleRegion() );
     itr.GoToBegin();
 
     while(!itr.IsAtEnd())
     {
       itw.Set(itr.Get().GetElement(index));
       ++itr;
       ++itw;
     }
   }
 
   m_DisplayIndex = index;
 }
 
 template<typename TPixelType>
 bool mitk::DiffusionImage<TPixelType>::AreAlike(GradientDirectionType g1,
                                                 GradientDirectionType g2,
                                                 double precision)
 {
   GradientDirectionType diff = g1 - g2;
   GradientDirectionType diff2 = g1 + g2;
   return diff.two_norm() < precision || diff2.two_norm() < precision;
 }
 
 template<typename TPixelType>
 void mitk::DiffusionImage<TPixelType>::CorrectDKFZBrokenGradientScheme(double precision)
 {
   GradientDirectionContainerType::Pointer directionSet = CalcAveragedDirectionSet(precision, m_Directions);
   if(directionSet->size() < 7)
   {
     MITK_INFO << "Too few directions, assuming and correcting DKFZ-bogus sequence details.";
 
     double v [7][3] =
     {{ 0,         0,         0        },
      {-0.707057,  0,         0.707057 },
      { 0.707057,  0,         0.707057 },
      { 0,         0.707057,  0.707057 },
      { 0,         0.707057, -0.707057 },
      {-0.707057,  0.707057,  0        },
      { 0.707057,  0.707057,  0        } };
 
     int i=0;
 
     for(GradientDirectionContainerType::Iterator it = m_OriginalDirections->Begin();
         it != m_OriginalDirections->End(); ++it)
     {
       it.Value().set(v[i++%7]);
     }
     ApplyMeasurementFrame();
   }
 }
 
 template<typename TPixelType>
 mitk::DiffusionImage<TPixelType>::GradientDirectionContainerType::Pointer
 mitk::DiffusionImage<TPixelType>::CalcAveragedDirectionSet(double precision, GradientDirectionContainerType::Pointer directions)
 {
   // save old and construct new direction container
   GradientDirectionContainerType::Pointer newDirections = GradientDirectionContainerType::New();
 
   // fill new direction container
   for(GradientDirectionContainerType::ConstIterator gdcitOld = directions->Begin();
       gdcitOld != directions->End(); ++gdcitOld)
   {
     // already exists?
     bool found = false;
     for(GradientDirectionContainerType::ConstIterator gdcitNew = newDirections->Begin();
         gdcitNew != newDirections->End(); ++gdcitNew)
     {
       if(AreAlike(gdcitNew.Value(), gdcitOld.Value(), precision))
       {
         found = true;
         break;
       }
     }
 
     // if not found, add it to new container
     if(!found)
     {
       newDirections->push_back(gdcitOld.Value());
     }
   }
 
   return newDirections;
 }
 
 template<typename TPixelType>
 void mitk::DiffusionImage<TPixelType>::AverageRedundantGradients(double precision)
 {
 
   GradientDirectionContainerType::Pointer newDirs =
       CalcAveragedDirectionSet(precision, m_Directions);
 
   GradientDirectionContainerType::Pointer newOriginalDirs =
       CalcAveragedDirectionSet(precision, m_OriginalDirections);
 
   // if sizes equal, we do not need to do anything in this function
   if(m_Directions->size() == newDirs->size() || m_OriginalDirections->size() == newOriginalDirs->size())
     return;
 
   GradientDirectionContainerType::Pointer oldDirections = m_OriginalDirections;
   m_Directions = newDirs;
   m_OriginalDirections = newOriginalDirs;
 
   // new image
   typename ImageType::Pointer oldImage = m_VectorImage;
   m_VectorImage = ImageType::New();
   m_VectorImage->SetSpacing( oldImage->GetSpacing() );   // Set the image spacing
   m_VectorImage->SetOrigin( oldImage->GetOrigin() );     // Set the image origin
   m_VectorImage->SetDirection( oldImage->GetDirection() );  // Set the image direction
   m_VectorImage->SetLargestPossibleRegion( oldImage->GetLargestPossibleRegion() );
   m_VectorImage->SetVectorLength( m_Directions->size() );
   m_VectorImage->SetBufferedRegion( oldImage->GetLargestPossibleRegion() );
   m_VectorImage->Allocate();
 
   // average image data that corresponds to identical directions
   itk::ImageRegionIterator< ImageType > newIt(m_VectorImage, m_VectorImage->GetLargestPossibleRegion());
   newIt.GoToBegin();
   itk::ImageRegionIterator< ImageType > oldIt(oldImage, oldImage->GetLargestPossibleRegion());
   oldIt.GoToBegin();
 
   // initial new value of voxel
   typename ImageType::PixelType newVec;
   newVec.SetSize(m_Directions->size());
   newVec.AllocateElements(m_Directions->size());
 
   std::vector<std::vector<int> > dirIndices;
   for(GradientDirectionContainerType::ConstIterator gdcitNew = m_Directions->Begin();
       gdcitNew != m_Directions->End(); ++gdcitNew)
   {
     dirIndices.push_back(std::vector<int>(0));
     for(GradientDirectionContainerType::ConstIterator gdcitOld = oldDirections->Begin();
         gdcitOld != oldDirections->End(); ++gdcitOld)
     {
       if(AreAlike(gdcitNew.Value(), gdcitOld.Value(), precision))
       {
         //MITK_INFO << gdcitNew.Value() << "  " << gdcitOld.Value();
         dirIndices[gdcitNew.Index()].push_back(gdcitOld.Index());
       }
     }
   }
 
   //int ind1 = -1;
   while(!newIt.IsAtEnd())
   {
 
     // progress
     //typename ImageType::IndexType ind = newIt.GetIndex();
     //ind1 = ind.m_Index[2];
 
     // init new vector with zeros
     newVec.Fill(0.0);
 
     // the old voxel value with duplicates
     typename ImageType::PixelType oldVec = oldIt.Get();
 
     for(unsigned int i=0; i<dirIndices.size(); i++)
     {
       // do the averaging
       const unsigned int numavg = dirIndices[i].size();
       unsigned int sum = 0;
-      for(int j=0; j<numavg; j++)
+      for(unsigned int j=0; j<numavg; j++)
       {
         //MITK_INFO << newVec[i] << " << " << oldVec[dirIndices[i].at(j)];
         sum += oldVec[dirIndices[i].at(j)];
       }
       if(numavg == 0)
       {
         MITK_ERROR << "mitkDiffusionImage: Error on averaging. Possibly due to corrupted data";
         return;
       }
       newVec[i] = sum / numavg;
     }
 
     newIt.Set(newVec);
 
     ++newIt;
     ++oldIt;
   }
   ApplyMeasurementFrame();
   std::cout << std::endl;
 }
 
 template<typename TPixelType>
 void mitk::DiffusionImage<TPixelType>::ApplyMeasurementFrame()
 {
   RemoveDirectionsContainerObserver();
 
   m_Directions = GradientDirectionContainerType::New();
   int c = 0;
   for(GradientDirectionContainerType::ConstIterator gdcit = m_OriginalDirections->Begin();
       gdcit != m_OriginalDirections->End(); ++gdcit)
   {
     vnl_vector<double> vec = gdcit.Value();
     vec = vec.pre_multiply(m_MeasurementFrame);
     m_Directions->InsertElement(c, vec);
     c++;
   }
 
   UpdateBValueMap();
   AddDirectionsContainerObserver();
 }
 
 // returns number of gradients
 template<typename TPixelType>
 int mitk::DiffusionImage<TPixelType>::GetNumDirections()
 {
   int gradients = m_OriginalDirections->Size();
   for (int i=0; i<m_OriginalDirections->Size(); i++)
     if (GetB_Value(i)<=0)
     {
       gradients--;
     }
   return gradients;
 }
 
 // returns number of not diffusion weighted images
 template<typename TPixelType>
 int mitk::DiffusionImage<TPixelType>::GetNumB0()
 {
   int b0 = 0;
   for (int i=0; i<m_OriginalDirections->Size(); i++)
     if (GetB_Value(i)<=0)
     {
       b0++;
     }
   return b0;
 }
 
 // returns a list of indices belonging to the not diffusion weighted images
 template<typename TPixelType>
 typename mitk::DiffusionImage<TPixelType>::IndicesVector mitk::DiffusionImage<TPixelType>::GetB0Indices()
 {
   IndicesVector indices;
   for (int i=0; i<m_OriginalDirections->Size(); i++)
     if (GetB_Value(i)<=0)
     {
       indices.push_back(i);
     }
   return indices;
 }
 
 template<typename TPixelType>
 bool mitk::DiffusionImage<TPixelType>::IsMultiBval()
 {
   int gradients = m_OriginalDirections->Size();
 
   for (int i=0; i<gradients; i++)
     if (GetB_Value(i)>0 && std::fabs(m_B_Value-GetB_Value(i))>50)
       return true;
   return false;
 }
 
 template<typename TPixelType>
 void mitk::DiffusionImage<TPixelType>::UpdateBValueMap()
 {
   m_B_ValueMap.clear();
 
   GradientDirectionContainerType::ConstIterator gdcit;
   for( gdcit = this->m_Directions->Begin(); gdcit != this->m_Directions->End(); ++gdcit)
     m_B_ValueMap[GetB_Value(gdcit.Index())].push_back(gdcit.Index());
 }
 
 template<typename TPixelType>
 float mitk::DiffusionImage<TPixelType>::GetB_Value(unsigned int i)
 {
   if(i > m_Directions->Size()-1)
     return -1;
 
   if(m_Directions->ElementAt(i).one_norm() <= 0.0)
   {
     return 0;
   }
   else
   {
     double twonorm = m_Directions->ElementAt(i).two_norm();
     double bval = m_B_Value*twonorm*twonorm;
 
     if (bval<0)
         bval = ceil(bval - 0.5);
     else
         bval = floor(bval + 0.5);
 
     return bval;
   }
 }
 
 template<typename TPixelType>
 void mitk::DiffusionImage<TPixelType>::SetDirections(const std::vector<itk::Vector<double,3> > directions)
 {
   m_OriginalDirections = GradientDirectionContainerType::New();
   for(unsigned int i=0; i<directions.size(); i++)
   {
     m_OriginalDirections->InsertElement( i, directions[i].GetVnlVector() );
   }
   this->ApplyMeasurementFrame();
 }
 
 template<typename TPixelType>
 void mitk::DiffusionImage<TPixelType>::AddDirectionsContainerObserver()
 {
   // Add Modified Observer to invoke an UpdateBValueList by modifieng the DirectionsContainer (m_Directions)
   typedef DiffusionImage< TPixelType > Self;
   typedef itk::SimpleMemberCommand< Self >  DCCommand ;
   typename DCCommand::Pointer command = DCCommand::New();
   command->SetCallbackFunction(this, &Self::UpdateBValueMap);
 }
 
 template<typename TPixelType>
 void mitk::DiffusionImage<TPixelType>::RemoveDirectionsContainerObserver()
 {
   if(m_Directions){
     m_Directions->RemoveAllObservers();
   }
 }
 
diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkParticle.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkParticle.h
index 03cb4c09db..5d1f64928c 100644
--- a/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkParticle.h
+++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkParticle.h
@@ -1,94 +1,98 @@
 /*===================================================================
 
 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 _PARTICLE
 #define _PARTICLE
 
 #include <FiberTrackingExports.h>
 #include <vnl/vnl_vector_fixed.h>
 
 
 namespace mitk
 {
 
 /**
 * \brief A particle is the basic element of the Gibbs fiber tractography method.   */
 
 class FiberTracking_EXPORT Particle
 {
 public:
 
     Particle()
     {
         label = 0;
         pID = -1;
         mID = -1;
     }
 
     ~Particle()
     {
     }
 
     int gridindex;          // index in the grid where it is living
     int ID;                 // particle ID
     int pID;                // successor ID
     int mID;                // predecessor ID
     unsigned char label;    // label used in the fiber building process
 
     vnl_vector_fixed<float, 3>& GetPos()
     {
       return pos;
     }
 
     vnl_vector_fixed<float, 3>& GetDir()
     {
       return dir;
     }
 
 private:
-#pragma warning(push)
-#pragma warning(disable: 4251)
+#ifdef _MSC_VER
+ #pragma warning(push)
+ #pragma warning(disable: 4251)
+#endif
     // this pragma ignores the following warning:
     // warning C4251: 'mitk::Particle::pos' : class   'ATL::CStringT'   needs to have dll-interface to be used   by clients of class 'Particle'
     vnl_vector_fixed<float, 3> pos; // particle position (world coordinates. corner based voxels. not accounted for image rotation.
     vnl_vector_fixed<float, 3> dir; // normalized direction vector
-#pragma warning(pop)
+#ifdef _MSC_VER
+ #pragma warning(pop)
+#endif
 
 };
 
 class FiberTracking_EXPORT EndPoint
 {
 public:
     EndPoint()
     {}
 
     EndPoint(Particle *p,int ep)
     {
         this->p = p;
         this->ep = ep;
     }
     Particle *p;
     int ep;
 
     inline bool operator==(EndPoint P)
     {
         return (P.p == p) && (P.ep == ep);
     }
 };
 
 }
 
 #endif