diff --git a/Modules/DiffusionImaging/Connectomics/CMakeLists.txt b/Modules/DiffusionImaging/Connectomics/CMakeLists.txt
index bfc3988cba..c9be1479dc 100644
--- a/Modules/DiffusionImaging/Connectomics/CMakeLists.txt
+++ b/Modules/DiffusionImaging/Connectomics/CMakeLists.txt
@@ -1,9 +1,15 @@
+# Boost throws a lot of unused variable warnings on linux
+
+if(UNIX)
+  mitkFunctionCheckCAndCXXCompilerFlags("-Wno-error=unused-variable" CMAKE_C_FLAGS CMAKE_CXX_FLAGS)
+endif()
+
 MITK_CREATE_MODULE( Connectomics
   SUBPROJECTS MITK-DTI
   INCLUDE_DIRS Algorithms IODataStructures Rendering ${CMAKE_CURRENT_BINARY_DIR}
   DEPENDS DiffusionCore FiberTracking QmitkExt
   PACKAGE_DEPENDS VTK|vtkInfovisLayout
   WARNINGS_AS_ERRORS
 )
 
 add_subdirectory(Testing)
diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp
index f91701735b..f59a948b47 100644
--- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp
+++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp
@@ -1,817 +1,817 @@
 /*===================================================================
 
 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 __itkAnalyticalDiffusionQballReconstructionImageFilter_cpp
 #define __itkAnalyticalDiffusionQballReconstructionImageFilter_cpp
 
 #include <itkAnalyticalDiffusionQballReconstructionImageFilter.h>
 #include <itkImageRegionConstIterator.h>
 #include <itkImageRegionConstIteratorWithIndex.h>
 #include <itkImageRegionIterator.h>
 #include <itkArray.h>
 #include <vnl/vnl_vector.h>
 
 #include <boost/version.hpp>
 #include <stdio.h>
 #include <locale>
 
 #define _USE_MATH_DEFINES
 #include <math.h>
 
 #include <boost/math/special_functions.hpp>
 
 #include "itkPointShell.h"
 
 using namespace boost::math;
 
 namespace itk {
 
 #define QBALL_ANAL_RECON_PI       M_PI
 
 template< class T, class TG, class TO, int L, int NODF>
 AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
 ::AnalyticalDiffusionQballReconstructionImageFilter() :
     m_GradientDirectionContainer(NULL),
     m_NumberOfGradientDirections(0),
     m_NumberOfBaselineImages(1),
     m_Threshold(NumericTraits< ReferencePixelType >::NonpositiveMin()),
     m_BValue(1.0),
     m_Lambda(0.0),
     m_DirectionsDuplicated(false),
     m_Delta1(0.001),
     m_Delta2(0.001),
     m_UseMrtrixBasis(false)
 {
     // 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 );
 }
 
 
 template<
         class TReferenceImagePixelType,
         class TGradientImagePixelType,
         class TOdfPixelType,
         int NOrderL,
         int NrOdfDirections>
 typename itk::AnalyticalDiffusionQballReconstructionImageFilter<
 TReferenceImagePixelType,TGradientImagePixelType,TOdfPixelType,
 NOrderL,NrOdfDirections>::OdfPixelType
 itk::AnalyticalDiffusionQballReconstructionImageFilter
 <TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType,
 NOrderL, NrOdfDirections>
 ::Normalize( OdfPixelType odf,
              typename NumericTraits<ReferencePixelType>::AccumulateType b0 )
 {
     switch( m_NormalizationMethod )
     {
     case QBAR_STANDARD:
     {
         TOdfPixelType sum = 0;
 
         for(int i=0; i<NrOdfDirections; i++)
         {
             sum += odf[i];
         }
         if(sum>0)
             odf /= sum;
 
         return odf;
         break;
     }
     case QBAR_B_ZERO_B_VALUE:
     {
         for(int i=0; i<NrOdfDirections; i++)
         {
             odf[i] = ((TOdfPixelType)log((TOdfPixelType)b0)-odf[i])/m_BValue;
         }
         return odf;
         break;
     }
     case QBAR_B_ZERO:
     {
         odf *= 1.0/b0;
         return odf;
         break;
     }
     case QBAR_NONE:
     {
         return odf;
         break;
     }
     case QBAR_ADC_ONLY:
     {
         for(int i=0; i<NrOdfDirections; i++)
         {
             odf[i] = ((TOdfPixelType)log((TOdfPixelType)b0)-odf[i])/m_BValue;
         }
         return odf;
         break;
     }
     case QBAR_RAW_SIGNAL:
     {
         return odf;
         break;
     }
     case QBAR_SOLID_ANGLE:
     {
         for(int i=0; i<NrOdfDirections; i++)
             odf[i] *= QBALL_ANAL_RECON_PI*4/NrOdfDirections;
 
         break;
     }
     case QBAR_NONNEG_SOLID_ANGLE:
     {
         break;
     }
     }
 
     return odf;
 }
 
 
 template<
         class TReferenceImagePixelType,
         class TGradientImagePixelType,
         class TOdfPixelType,
         int NOrderL,
         int NrOdfDirections>
 vnl_vector<TOdfPixelType>
 itk::AnalyticalDiffusionQballReconstructionImageFilter
 <TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType,
 NOrderL, NrOdfDirections>
 ::PreNormalize( vnl_vector<TOdfPixelType> vec,
                 typename NumericTraits<ReferencePixelType>::AccumulateType b0 )
 {
     switch( m_NormalizationMethod )
     {
     case QBAR_STANDARD:
     {
         return vec;
         break;
     }
     case QBAR_B_ZERO_B_VALUE:
     {
         int n = vec.size();
         for(int i=0; i<n; i++)
         {
             if (vec[i]<=0)
                 vec[i] = 0.001;
 
             vec[i] = log(vec[i]);
         }
         return vec;
         break;
     }
     case QBAR_B_ZERO:
     {
         return vec;
         break;
     }
     case QBAR_NONE:
     {
         return vec;
         break;
     }
     case QBAR_ADC_ONLY:
     {
         int n = vec.size();
         for(int i=0; i<n; i++)
         {
             if (vec[i]<=0)
                 vec[i] = 0.001;
 
             vec[i] = log(vec[i]);
         }
         return vec;
         break;
     }
     case QBAR_RAW_SIGNAL:
     {
         return vec;
         break;
     }
     case QBAR_SOLID_ANGLE:
     case QBAR_NONNEG_SOLID_ANGLE:
     {
         int n = vec.size();
         double b0f = (double)b0;
         for(int i=0; i<n; i++)
         {
             vec[i] = vec[i]/b0f;
 
             if (vec[i]<0)
                 vec[i] = m_Delta1;
             else if (vec[i]<m_Delta1)
                 vec[i] = m_Delta1/2 + vec[i]*vec[i]/(2*m_Delta1);
             else if (vec[i]>=1)
                 vec[i] = 1-m_Delta2/2;
             else if (vec[i]>=1-m_Delta2)
                 vec[i] = 1-m_Delta2/2-(1-vec[i])*(1-vec[i])/(2*m_Delta2);
 
             vec[i] = log(-log(vec[i]));
         }
         return vec;
         break;
     }
     }
 
     return vec;
 }
 
 template< class T, class TG, class TO, int L, int NODF>
 void AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
 ::BeforeThreadedGenerateData()
 {
     // If we have more than 2 inputs, then each input, except the first is a
     // gradient image. The number of gradient images must match the number of
     // gradient directions.
     //const unsigned int numberOfInputs = this->GetNumberOfInputs();
 
     // There need to be at least 6 gradient directions to be able to compute the
     // tensor basis
     if( m_NumberOfGradientDirections < (L*L + L + 2)/2 + L )
     {
         itkExceptionMacro( << "Not enough gradient directions supplied (" << m_NumberOfGradientDirections << "). At least " << (L*L + L + 2)/2 + L << " needed for SH-order " << L);
     }
 
     // Input must be an itk::VectorImage.
     std::string gradientImageClassName(
                 this->ProcessObject::GetInput(0)->GetNameOfClass());
     if ( strcmp(gradientImageClassName.c_str(),"VectorImage") != 0 )
     {
         itkExceptionMacro( <<
                            "There is only one Gradient image. I expect that to be a VectorImage. "
                            << "But its of type: " << gradientImageClassName );
     }
 
     this->ComputeReconstructionMatrix();
 
     typename GradientImagesType::Pointer img = static_cast< GradientImagesType * >(
                 this->ProcessObject::GetInput(0) );
 
     m_BZeroImage = BZeroImageType::New();
     m_BZeroImage->SetSpacing( img->GetSpacing() );   // Set the image spacing
     m_BZeroImage->SetOrigin( img->GetOrigin() );     // Set the image origin
     m_BZeroImage->SetDirection( img->GetDirection() );  // Set the image direction
     m_BZeroImage->SetLargestPossibleRegion( img->GetLargestPossibleRegion());
     m_BZeroImage->SetBufferedRegion( img->GetLargestPossibleRegion() );
     m_BZeroImage->Allocate();
 
     m_ODFSumImage = BZeroImageType::New();
     m_ODFSumImage->SetSpacing( img->GetSpacing() );   // Set the image spacing
     m_ODFSumImage->SetOrigin( img->GetOrigin() );     // Set the image origin
     m_ODFSumImage->SetDirection( img->GetDirection() );  // Set the image direction
     m_ODFSumImage->SetLargestPossibleRegion( img->GetLargestPossibleRegion());
     m_ODFSumImage->SetBufferedRegion( img->GetLargestPossibleRegion() );
     m_ODFSumImage->Allocate();
 
     m_CoefficientImage = CoefficientImageType::New();
     m_CoefficientImage->SetSpacing( img->GetSpacing() );   // Set the image spacing
     m_CoefficientImage->SetOrigin( img->GetOrigin() );     // Set the image origin
     m_CoefficientImage->SetDirection( img->GetDirection() );  // Set the image direction
     m_CoefficientImage->SetLargestPossibleRegion( img->GetLargestPossibleRegion());
     m_CoefficientImage->SetBufferedRegion( img->GetLargestPossibleRegion() );
     m_CoefficientImage->Allocate();
 
     if(m_NormalizationMethod == QBAR_SOLID_ANGLE || m_NormalizationMethod == QBAR_NONNEG_SOLID_ANGLE)
         m_Lambda = 0.0;
 }
 
 template< class T, class TG, class TO, int L, int NODF>
 void AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
 ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
                        ThreadIdType )
 {
     typename OutputImageType::Pointer outputImage =
             static_cast< OutputImageType * >(this->ProcessObject::GetPrimaryOutput());
 
     ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread);
     oit.GoToBegin();
 
     ImageRegionIterator< BZeroImageType > oit2(m_BZeroImage, outputRegionForThread);
     oit2.GoToBegin();
 
     ImageRegionIterator< FloatImageType > oit3(m_ODFSumImage, outputRegionForThread);
     oit3.GoToBegin();
 
     ImageRegionIterator< CoefficientImageType > oit4(m_CoefficientImage, outputRegionForThread);
     oit4.GoToBegin();
 
     typedef ImageRegionConstIterator< GradientImagesType > GradientIteratorType;
     typedef typename GradientImagesType::PixelType         GradientVectorType;
     typename GradientImagesType::Pointer gradientImagePointer = NULL;
 
     // Would have liked a dynamic_cast here, but seems SGI doesn't like it
     // The enum will ensure that an inappropriate cast is not done
     gradientImagePointer = static_cast< GradientImagesType * >(
                 this->ProcessObject::GetInput(0) );
 
     GradientIteratorType git(gradientImagePointer, outputRegionForThread );
     git.GoToBegin();
 
     // Compute the indicies of the baseline images and gradient images
     std::vector<unsigned int> baselineind; // contains the indicies of
     // the baseline images
     std::vector<unsigned int> gradientind; // contains the indicies of
     // the gradient images
 
     for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin();
         gdcit != this->m_GradientDirectionContainer->End(); ++gdcit)
     {
         if(gdcit.Value().one_norm() <= 0.0)
             baselineind.push_back(gdcit.Index());
         else
             gradientind.push_back(gdcit.Index());
     }
 
     if( m_DirectionsDuplicated )
     {
         int gradIndSize = gradientind.size();
         for(int i=0; i<gradIndSize; i++)
             gradientind.push_back(gradientind[i]);
     }
 
     while( !git.IsAtEnd() )
     {
         GradientVectorType b = git.Get();
 
         typename NumericTraits<ReferencePixelType>::AccumulateType b0 = NumericTraits<ReferencePixelType>::Zero;
 
         // Average the baseline image pixels
         for(unsigned int i = 0; i < baselineind.size(); ++i)
         {
             b0 += b[baselineind[i]];
         }
         b0 /= this->m_NumberOfBaselineImages;
 
         OdfPixelType odf(0.0);
         typename CoefficientImageType::PixelType coeffPixel(0.0);
         vnl_vector<TO> B(m_NumberOfGradientDirections);
 
         if( (b0 != 0) && (b0 >= m_Threshold) )
         {
             for( unsigned int i = 0; i< m_NumberOfGradientDirections; i++ )
             {
                 B[i] = static_cast<TO>(b[gradientind[i]]);
             }
 
             B = PreNormalize(B, b0);
             if(m_NormalizationMethod == QBAR_SOLID_ANGLE)
             {
                 vnl_vector<TO> coeffs(m_NumberCoefficients);
                 coeffs = ( (*m_CoeffReconstructionMatrix) * B );
                 coeffs[0] += 1.0/(2.0*sqrt(QBALL_ANAL_RECON_PI));
                 odf = ( (*m_SphericalHarmonicBasisMatrix) * coeffs ).data_block();
                 coeffPixel = coeffs.data_block();
             }
             else if(m_NormalizationMethod == QBAR_NONNEG_SOLID_ANGLE)
             {
                 /** this would be the place to implement a non-negative
               * solver for quadratic programming problem:
               * min .5*|| Bc-s ||^2 subject to -CLPc <= 4*pi*ones
               * (refer to MICCAI 2009 Goh et al. "Estimating ODFs with PDF constraints")
               * .5*|| Bc-s ||^2 == .5*c'B'Bc - x'B's + .5*s's
               */
 
                 itkExceptionMacro( << "Nonnegative Solid Angle not yet implemented");
             }
             else
             {
                 vnl_vector<TO> coeffs(m_NumberCoefficients);
                 coeffs = ( (*m_CoeffReconstructionMatrix) * B );
                 coeffs[0] += 1.0/(2.0*sqrt(QBALL_ANAL_RECON_PI));
                 coeffPixel = coeffs.data_block();
                 odf = ( (*m_ReconstructionMatrix) * B ).data_block();
             }
             odf = Normalize(odf, b0);
 
         }
 
         oit.Set( odf );
         oit2.Set( b0 );
         float sum = 0;
         for (unsigned int k=0; k<odf.Size(); k++)
             sum += (float) odf[k];
         oit3.Set( sum-1 );
         oit4.Set(coeffPixel);
         ++oit;  // odf image iterator
         ++oit3; // odf sum image iterator
         ++oit2; // b0 image iterator
         ++oit4; // coefficient image iterator
         ++git;  // Gradient  image iterator
     }
 
     std::cout << "One Thread finished reconstruction" << std::endl;
 }
 
 template< class T, class TG, class TO, int L, int NODF>
 void AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
 ::tofile2(vnl_matrix<double> *pA, std::string fname)
 {
     vnl_matrix<double> A = (*pA);
     ofstream myfile;
     std::locale C("C");
     std::locale originalLocale = myfile.getloc();
     myfile.imbue(C);
 
     myfile.open (fname.c_str());
     myfile << "A1=[";
-    for(int i=0; i<A.rows(); i++)
+    for(unsigned int i=0; i<A.rows(); i++)
     {
-        for(int j=0; j<A.columns(); j++)
+        for(unsigned int j=0; j<A.columns(); j++)
         {
             myfile << A(i,j) << " ";
             if(j==A.columns()-1 && i!=A.rows()-1)
                 myfile << ";";
         }
     }
     myfile << "];";
     myfile.close();
 
     myfile.imbue( originalLocale );
 }
 
 template< class T, class TG, class TO, int L, int NODF>
 void AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
 ::Cart2Sph(double x, double y, double z, double *spherical)
 {
     double phi, theta, r;
     r = sqrt(x*x+y*y+z*z);
 
     if( r<mitk::eps )
     {
         theta = M_PI/2;
         phi = M_PI/2;
     }
     else
     {
         theta = acos(z/r);
         phi = atan2(y, x);
     }
     spherical[0] = phi;
     spherical[1] = theta;
     spherical[2] = r;
 }
 
 template< class T, class TG, class TO, int L, int NODF>
 double AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
 ::Yj(int m, int l, double theta, double phi, bool useMRtrixBasis)
 {
     if (!useMRtrixBasis)
     {
         if (m<0)
             return sqrt(2.0)*spherical_harmonic_r(l, -m, theta, phi);
         else if (m==0)
             return spherical_harmonic_r(l, m, theta, phi);
         else
             return pow(-1.0,m)*sqrt(2.0)*spherical_harmonic_i(l, m, theta, phi);
     }
     else
     {
         double plm = legendre_p<double>(l,abs(m),-cos(theta));
         double mag = sqrt((double)(2*l+1)/(4.0*M_PI)*factorial<double>(l-abs(m))/factorial<double>(l+abs(m)))*plm;
         if (m>0)
             return mag*cos(m*phi);
         else if (m==0)
             return mag;
         else
             return mag*sin(-m*phi);
     }
 
     return 0;
 }
 
 template< class T, class TG, class TO, int L, int NODF>
 double AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
 ::Legendre0(int l)
 {
     if( l%2 != 0 )
     {
         return 0;
     }
     else
     {
         double prod1 = 1.0;
         for(int i=1;i<l;i+=2) prod1 *= i;
         double prod2 = 1.0;
         for(int i=2;i<=l;i+=2) prod2 *= i;
         return pow(-1.0,l/2.0)*(prod1/prod2);
     }
 }
 
 template< class T, class TG, class TO, int L, int NODF>
 void AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
 ::ComputeReconstructionMatrix()
 {
 
     //for(int i=-6;i<7;i++)
     //  std::cout << boost::math::legendre_p<double>(6, i, 0.65657) << std::endl;
 
     if( m_NumberOfGradientDirections < (L*L + L + 2)/2 + L )
     {
         itkExceptionMacro( << "Not enough gradient directions supplied (" << m_NumberOfGradientDirections << "). At least " << (L*L + L + 2)/2 + L << " needed for SH-order " << L);
     }
 
     {
         // check for duplicate diffusion gradients
         bool warning = false;
         for(GradientDirectionContainerType::ConstIterator gdcit1 = this->m_GradientDirectionContainer->Begin();
             gdcit1 != this->m_GradientDirectionContainer->End(); ++gdcit1)
         {
             for(GradientDirectionContainerType::ConstIterator gdcit2 = this->m_GradientDirectionContainer->Begin();
                 gdcit2 != this->m_GradientDirectionContainer->End(); ++gdcit2)
             {
                 if(gdcit1.Value() == gdcit2.Value() && gdcit1.Index() != gdcit2.Index())
                 {
                     itkWarningMacro( << "Some of the Diffusion Gradients equal each other. Corresponding image data should be averaged before calling this filter." );
                     warning = true;
                     break;
                 }
             }
             if (warning) break;
         }
 
         // handle acquisition schemes where only half of the spherical
         // shell is sampled by the gradient directions. In this case,
         // each gradient direction is duplicated in negative direction.
         vnl_vector<double> centerMass(3);
         centerMass.fill(0.0);
         int count = 0;
         for(GradientDirectionContainerType::ConstIterator gdcit1 = this->m_GradientDirectionContainer->Begin();
             gdcit1 != this->m_GradientDirectionContainer->End(); ++gdcit1)
         {
             if(gdcit1.Value().one_norm() > 0.0)
             {
                 centerMass += gdcit1.Value();
                 count ++;
             }
         }
         centerMass /= count;
         if(centerMass.two_norm() > 0.1)
         {
             m_DirectionsDuplicated = true;
             m_NumberOfGradientDirections *= 2;
         }
     }
 
     vnl_matrix<double> *Q = new vnl_matrix<double>(3, m_NumberOfGradientDirections);
 
     {
         int i = 0;
         for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin();
             gdcit != this->m_GradientDirectionContainer->End(); ++gdcit)
         {
             if(gdcit.Value().one_norm() > 0.0)
             {
                 double x = gdcit.Value().get(0);
                 double y = gdcit.Value().get(1);
                 double z = gdcit.Value().get(2);
                 double cart[3];
                 Cart2Sph(x,y,z,cart);
                 (*Q)(0,i) = cart[0];
                 (*Q)(1,i) = cart[1];
                 (*Q)(2,i++) = cart[2];
             }
         }
         if(m_DirectionsDuplicated)
         {
             for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin();
                 gdcit != this->m_GradientDirectionContainer->End(); ++gdcit)
             {
                 if(gdcit.Value().one_norm() > 0.0)
                 {
                     double x = gdcit.Value().get(0);
                     double y = gdcit.Value().get(1);
                     double z = gdcit.Value().get(2);
                     double cart[3];
                     Cart2Sph(x,y,z,cart);
                     (*Q)(0,i) = cart[0];
                     (*Q)(1,i) = cart[1];
                     (*Q)(2,i++) = cart[2];
                 }
             }
         }
     }
 
     int l = L;
     m_NumberCoefficients = (int)(l*l + l + 2.0)/2.0 + l;
     vnl_matrix<double>* B = new vnl_matrix<double>(m_NumberOfGradientDirections,m_NumberCoefficients);
     vnl_matrix<double>* _L = new vnl_matrix<double>(m_NumberCoefficients,m_NumberCoefficients);
     _L->fill(0.0);
     vnl_matrix<double>* LL = new vnl_matrix<double>(m_NumberCoefficients,m_NumberCoefficients);
     LL->fill(0.0);
     vnl_matrix<double>* P = new vnl_matrix<double>(m_NumberCoefficients,m_NumberCoefficients);
     P->fill(0.0);
     vnl_matrix<double>* Inv = new vnl_matrix<double>(m_NumberCoefficients,m_NumberCoefficients);
     P->fill(0.0);
     vnl_vector<int>* lj = new vnl_vector<int>(m_NumberCoefficients);
     m_LP = new vnl_vector<double>(m_NumberCoefficients);
 
     for(unsigned int i=0; i<m_NumberOfGradientDirections; i++)
     {
         for(int k=0; k<=l; k+=2)
         {
             for(int m=-k; m<=k; m++)
             {
                 int j = (k*k + k + 2)/2 + m - 1;
                 (*_L)(j,j) = -k*(k+1);
                 (*m_LP)(j) = -k*(k+1);
                 (*lj)[j] = k;
                 double phi = (*Q)(0,i);
                 double th = (*Q)(1,i);
                 (*B)(i,j) = Yj(m,k,th,phi, m_UseMrtrixBasis);
             }
         }
     }
 
     //vnl_matrix<double> temp((*_L)*(*_L));
     LL->update(*_L);
     *LL *= *_L;
     //tofile2(LL,"LL");
 
     for(int i=0; i<m_NumberCoefficients; i++)
     {
         // here we leave out the 2*pi multiplication from Descoteaux
         // this is because we do not want the real integral value
         // but an average value over the equator.
 
         if(m_NormalizationMethod == QBAR_SOLID_ANGLE || m_NormalizationMethod == QBAR_NONNEG_SOLID_ANGLE)
         {
             (*P)(i,i) = 2.0*QBALL_ANAL_RECON_PI*Legendre0((*lj)[i]);
             (*m_LP)(i) *= (*P)(i,i);
         }
         else
         {
             (*P)(i,i) = Legendre0((*lj)[i]);
         }
     }
     m_B_t = new vnl_matrix<double>(B->transpose());
     //tofile2(&m_B_t,"m_B_t");
     vnl_matrix<double> B_t_B = (*m_B_t) * (*B);
     //tofile2(&B_t_B,"B_t_B");
     vnl_matrix<double> lambdaLL(m_NumberCoefficients,m_NumberCoefficients);
     lambdaLL.update((*LL));
     lambdaLL *= m_Lambda;
     //tofile2(&lambdaLL,"lLL");
 
     vnl_matrix<double> tmp( B_t_B + lambdaLL);
     vnl_matrix_inverse<double> *pseudoInverse
             = new vnl_matrix_inverse<double>( tmp );
 
     (*Inv) = pseudoInverse->pinverse();
     //tofile2(Inv,"Inv");
     vnl_matrix<double> temp((*Inv) * (*m_B_t));
     double fac1 = (1.0/(16.0*QBALL_ANAL_RECON_PI*QBALL_ANAL_RECON_PI));
     switch(m_NormalizationMethod)
     {
     case QBAR_ADC_ONLY:
     case QBAR_RAW_SIGNAL:
         break;
     case QBAR_STANDARD:
     case QBAR_B_ZERO_B_VALUE:
     case QBAR_B_ZERO:
     case QBAR_NONE:
         temp = (*P) * temp;
         break;
     case QBAR_SOLID_ANGLE:
         temp = fac1 * (*P) * (*_L) * temp;
         break;
     case QBAR_NONNEG_SOLID_ANGLE:
         break;
     }
 
     //tofile2(&temp,"A");
 
     m_CoeffReconstructionMatrix = new vnl_matrix<TO>(m_NumberCoefficients,m_NumberOfGradientDirections);
     for(int i=0; i<m_NumberCoefficients; i++)
     {
         for(unsigned int j=0; j<m_NumberOfGradientDirections; j++)
         {
             (*m_CoeffReconstructionMatrix)(i,j) = (float) temp(i,j);
         }
     }
 
     // this code goes to the image adapter coeffs->odfs later
 
     int NOdfDirections = NODF;
     vnl_matrix_fixed<double, 3, NODF>* U =
             itk::PointShell<NODF, vnl_matrix_fixed<double, 3, NODF> >::DistributePointShell();
 
     m_SphericalHarmonicBasisMatrix  = new vnl_matrix<TO>(NOdfDirections,m_NumberCoefficients);
     vnl_matrix<double>* sphericalHarmonicBasisMatrix2
             = new vnl_matrix<double>(NOdfDirections,m_NumberCoefficients);
     for(int i=0; i<NOdfDirections; i++)
     {
         double x = (*U)(0,i);
         double y = (*U)(1,i);
         double z = (*U)(2,i);
         double cart[3];
         Cart2Sph(x,y,z,cart);
         (*U)(0,i) = cart[0];
         (*U)(1,i) = cart[1];
         (*U)(2,i) = cart[2];
     }
 
     for(int i=0; i<NOdfDirections; i++)
     {
         for(int k=0; k<=l; k+=2)
         {
             for(int m=-k; m<=k; m++)
             {
                 int j = (k*k + k + 2)/2 + m - 1;
                 double phi = (*U)(0,i);
                 double th = (*U)(1,i);
                 (*m_SphericalHarmonicBasisMatrix)(i,j) = Yj(m,k,th,phi,m_UseMrtrixBasis);
                 (*sphericalHarmonicBasisMatrix2)(i,j) = (*m_SphericalHarmonicBasisMatrix)(i,j);
             }
         }
     }
 
     m_ReconstructionMatrix = new vnl_matrix<TO>(NOdfDirections,m_NumberOfGradientDirections);
     *m_ReconstructionMatrix = (*m_SphericalHarmonicBasisMatrix) * (*m_CoeffReconstructionMatrix);
 
 }
 
 template< class T, class TG, class TO, int L, int NODF>
 void AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
 ::SetGradientImage(const GradientDirectionContainerType *gradientDirection,
                    const GradientImagesType *gradientImage )
 {
     // Copy Gradient Direction Container
     this->m_GradientDirectionContainer = GradientDirectionContainerType::New();
     for(GradientDirectionContainerType::ConstIterator it = gradientDirection->Begin();
         it != gradientDirection->End(); it++)
     {
         this->m_GradientDirectionContainer->push_back(it.Value());
     }
 
 
     unsigned int numImages = gradientDirection->Size();
     this->m_NumberOfBaselineImages = 0;
     for(GradientDirectionContainerType::Iterator it = this->m_GradientDirectionContainer->Begin();
         it != this->m_GradientDirectionContainer->End(); it++)
     {
         if(it.Value().one_norm() <= 0.0)
         {
             this->m_NumberOfBaselineImages++;
         }
         else // Normalize non-zero gradient directions
         {
             it.Value() = it.Value() / it.Value().two_norm();
         }
     }
 
     this->m_NumberOfGradientDirections = numImages - this->m_NumberOfBaselineImages;
 
     // ensure that the gradient image we received has as many components as
     // the number of gradient directions
     if( gradientImage->GetVectorLength() != this->m_NumberOfBaselineImages + m_NumberOfGradientDirections )
     {
         itkExceptionMacro( << m_NumberOfGradientDirections << " gradients + " << this->m_NumberOfBaselineImages
                            << "baselines = " << m_NumberOfGradientDirections + this->m_NumberOfBaselineImages
                            << " directions specified but image has " << gradientImage->GetVectorLength()
                            << " components.");
     }
 
     this->ProcessObject::SetNthInput( 0, const_cast< GradientImagesType* >(gradientImage) );
 }
 
 template< class T, class TG, class TO, int L, int NODF>
 void AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
 ::PrintSelf(std::ostream& os, Indent indent) const
 {
     std::locale C("C");
     std::locale originalLocale = os.getloc();
     os.imbue(C);
 
     Superclass::PrintSelf(os,indent);
 
     os << indent << "OdfReconstructionMatrix: " << m_ReconstructionMatrix << std::endl;
     if ( m_GradientDirectionContainer )
         os << indent << "GradientDirectionContainer: " << m_GradientDirectionContainer << std::endl;
     else
         os << indent << "GradientDirectionContainer: (Gradient directions not set)" << std::endl;
 
     os << indent << "NumberOfGradientDirections: " <<  m_NumberOfGradientDirections << std::endl;
     os << indent << "NumberOfBaselineImages: " << m_NumberOfBaselineImages << std::endl;
     os << indent << "Threshold for reference B0 image: " << m_Threshold << std::endl;
     os << indent << "BValue: " << m_BValue << std::endl;
     os.imbue( originalLocale );
 }
 
 }
 
 #endif // __itkAnalyticalDiffusionQballReconstructionImageFilter_cpp
diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkOrientationDistributionFunction.txx b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkOrientationDistributionFunction.txx
index 6c3a17ed2c..269dc8e85a 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_t 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
diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.cpp
index fe4df6a770..c0a7253736 100644
--- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.cpp
@@ -1,832 +1,832 @@
 #include "itkTractsToVectorImageFilter.h"
 
 // VTK
 #include <vtkPolyLine.h>
 #include <vtkCellArray.h>
 #include <vtkCellData.h>
 
 // ITK
 #include <itkTimeProbe.h>
 #include <itkImageRegionIterator.h>
 
 // misc
 #define _USE_MATH_DEFINES
 #include <math.h>
 #include <boost/progress.hpp>
 
 
 namespace itk{
 
 static bool CompareVectorLengths(const vnl_vector_fixed< double, 3 >& v1, const vnl_vector_fixed< double, 3 >& v2)
 {
     return (v1.magnitude()>v2.magnitude());
 }
 
 template< class PixelType >
 TractsToVectorImageFilter< PixelType >::TractsToVectorImageFilter():
     m_AngularThreshold(0.7),
     m_Epsilon(0.999),
     m_MaskImage(NULL),
     m_NormalizeVectors(false),
     m_UseWorkingCopy(true),
     m_UseTrilinearInterpolation(false),
     m_MaxNumDirections(3),
     m_Thres(0.5),
     m_NumDirectionsImage(NULL)
 {
     this->SetNumberOfRequiredOutputs(1);
 }
 
 
 template< class PixelType >
 TractsToVectorImageFilter< PixelType >::~TractsToVectorImageFilter()
 {
 }
 
 
 template< class PixelType >
 vnl_vector_fixed<double, 3> TractsToVectorImageFilter< 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 >
 itk::Point<float, 3> TractsToVectorImageFilter< 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 >
 void TractsToVectorImageFilter< PixelType >::GenerateData()
 {
     mitk::Geometry3D::Pointer geometry = m_FiberBundle->GetGeometry();
 
     // calculate new image parameters
     itk::Vector<double> spacing;
     itk::Point<double> origin;
     itk::Matrix<double, 3, 3> direction;
     ImageRegion<3> imageRegion;
     if (!m_MaskImage.IsNull())
     {
         spacing = m_MaskImage->GetSpacing();
         imageRegion = m_MaskImage->GetLargestPossibleRegion();
         origin = m_MaskImage->GetOrigin();
         direction = m_MaskImage->GetDirection();
     }
     else
     {
         spacing = geometry->GetSpacing();
         origin = geometry->GetOrigin();
         mitk::Geometry3D::BoundsArrayType bounds = geometry->GetBounds();
         origin[0] += bounds.GetElement(0);
         origin[1] += bounds.GetElement(2);
         origin[2] += bounds.GetElement(4);
 
         for (int i=0; i<3; i++)
             for (int j=0; j<3; j++)
                 direction[j][i] = geometry->GetMatrixColumn(i)[j];
         imageRegion.SetSize(0, geometry->GetExtent(0));
         imageRegion.SetSize(1, geometry->GetExtent(1));
         imageRegion.SetSize(2, geometry->GetExtent(2));
 
 
         m_MaskImage = ItkUcharImgType::New();
         m_MaskImage->SetSpacing( spacing );
         m_MaskImage->SetOrigin( origin );
         m_MaskImage->SetDirection( direction );
         m_MaskImage->SetRegions( imageRegion );
         m_MaskImage->Allocate();
         m_MaskImage->FillBuffer(1);
     }
     OutputImageType::RegionType::SizeType outImageSize = imageRegion.GetSize();
     m_OutImageSpacing = m_MaskImage->GetSpacing();
     m_ClusteredDirectionsContainer = ContainerType::New();
 
     // initialize crossings image
     m_CrossingsImage = ItkUcharImgType::New();
     m_CrossingsImage->SetSpacing( spacing );
     m_CrossingsImage->SetOrigin( origin );
     m_CrossingsImage->SetDirection( direction );
     m_CrossingsImage->SetRegions( imageRegion );
     m_CrossingsImage->Allocate();
     m_CrossingsImage->FillBuffer(0);
 
     // initialize num directions image
     m_NumDirectionsImage = ItkUcharImgType::New();
     m_NumDirectionsImage->SetSpacing( spacing );
     m_NumDirectionsImage->SetOrigin( origin );
     m_NumDirectionsImage->SetDirection( direction );
     m_NumDirectionsImage->SetRegions( imageRegion );
     m_NumDirectionsImage->Allocate();
     m_NumDirectionsImage->FillBuffer(0);
 
     // resample fiber bundle
     float minSpacing = 1;
     if(m_OutImageSpacing[0]<m_OutImageSpacing[1] && m_OutImageSpacing[0]<m_OutImageSpacing[2])
         minSpacing = m_OutImageSpacing[0];
     else if (m_OutImageSpacing[1] < m_OutImageSpacing[2])
         minSpacing = m_OutImageSpacing[1];
     else
         minSpacing = m_OutImageSpacing[2];
 
     if (m_UseWorkingCopy)
         m_FiberBundle = m_FiberBundle->GetDeepCopy();
 
     // resample fiber bundle for sufficient voxel coverage
     m_FiberBundle->ResampleFibers(minSpacing/3);
 
     // iterate over all fibers
     vtkSmartPointer<vtkPolyData> fiberPolyData = m_FiberBundle->GetFiberPolyData();
     vtkSmartPointer<vtkCellArray> vLines = fiberPolyData->GetLines();
     vLines->InitTraversal();
     int numFibers = m_FiberBundle->GetNumFibers();
     itk::TimeProbe clock;
     m_DirectionsContainer = ContainerType::New();
 
     if (m_UseTrilinearInterpolation)
         MITK_INFO << "Generating directions from tractogram (trilinear interpolation)";
     else
         MITK_INFO << "Generating directions from tractogram";
 
     boost::progress_display disp(numFibers);
     for( int i=0; i<numFibers; i++ )
     {
         ++disp;
         clock.Start();
         vtkIdType   numPoints(0);
         vtkIdType*  points(NULL);
         vLines->GetNextCell ( numPoints, points );
         if (numPoints<2)
             continue;
 
         itk::Index<3> index; index.Fill(0);
         itk::ContinuousIndex<float, 3> contIndex;
         vnl_vector_fixed<double, 3> dir, wDir;
         itk::Point<float, 3> worldPos;
         vnl_vector<double> v;
         for( int j=0; j<numPoints-1; j++)
         {
             double* temp = fiberPolyData->GetPoint(points[j]);
             worldPos = GetItkPoint(temp);
             v = GetVnlVector(temp);
 
             dir = GetVnlVector(fiberPolyData->GetPoint(points[j+1]))-v;
             dir.normalize();
 
             m_MaskImage->TransformPhysicalPointToIndex(worldPos, index);
             m_MaskImage->TransformPhysicalPointToContinuousIndex(worldPos, contIndex);
 
             if (m_MaskImage->GetPixel(index)==0)
                 continue;
 
             if (!m_UseTrilinearInterpolation)
             {
                 if (index[0] < 0 || (unsigned long)index[0] >= outImageSize[0])
                     continue;
                 if (index[1] < 0 || (unsigned long)index[1] >= outImageSize[1])
                     continue;
                 if (index[2] < 0 || (unsigned long)index[2] >= outImageSize[2])
                     continue;
 
                 unsigned int idx = index[0] + outImageSize[0]*(index[1] + outImageSize[1]*index[2]);
                 DirectionContainerType::Pointer dirCont = DirectionContainerType::New();
                 if (m_DirectionsContainer->IndexExists(idx))
                 {
                     dirCont = m_DirectionsContainer->GetElement(idx);
                     if (dirCont.IsNull())
                     {
                         dirCont = DirectionContainerType::New();
                         dirCont->InsertElement(0, dir);
                         m_DirectionsContainer->InsertElement(idx, dirCont);
                     }
                     else
                         dirCont->InsertElement(dirCont->Size(), dir);
                 }
                 else
                 {
                     dirCont->InsertElement(0, dir);
                     m_DirectionsContainer->InsertElement(idx, dirCont);
                 }
 
                 continue;
             }
 
             float frac_x = contIndex[0] - index[0];
             float frac_y = contIndex[1] - index[1];
             float frac_z = contIndex[2] - index[2];
 
             if (frac_x<0)
             {
                 index[0] -= 1;
                 frac_x += 1;
             }
             if (frac_y<0)
             {
                 index[1] -= 1;
                 frac_y += 1;
             }
             if (frac_z<0)
             {
                 index[2] -= 1;
                 frac_z += 1;
             }
 
             frac_x = 1-frac_x;
             frac_y = 1-frac_y;
             frac_z = 1-frac_z;
 
             // int coordinates inside image?
             if (index[0] < 0 || (unsigned long)index[0] >= outImageSize[0]-1)
                 continue;
             if (index[1] < 0 || (unsigned long)index[1] >= outImageSize[1]-1)
                 continue;
             if (index[2] < 0 || (unsigned long)index[2] >= outImageSize[2]-1)
                 continue;
 
             DirectionContainerType::Pointer dirCont;
             int idx;
             wDir = dir;
             float weight = (  frac_x)*(  frac_y)*(  frac_z);
             if (weight>m_Thres)
             {
                 wDir *= weight;
                 idx = index[0]   + outImageSize[0]*(index[1]  + outImageSize[1]*index[2]  );
                 dirCont = DirectionContainerType::New();
                 if (m_DirectionsContainer->IndexExists(idx))
                 {
                     dirCont = m_DirectionsContainer->GetElement(idx);
                     if (dirCont.IsNull())
                     {
                         dirCont = DirectionContainerType::New();
                         dirCont->InsertElement(0, wDir);
                         m_DirectionsContainer->InsertElement(idx, dirCont);
                     }
                     else
                         dirCont->InsertElement(dirCont->Size(), wDir);
                 }
                 else
                 {
                     dirCont->InsertElement(0, wDir);
                     m_DirectionsContainer->InsertElement(idx, dirCont);
                 }
             }
 
             wDir = dir;
             weight = (  frac_x)*(1-frac_y)*(  frac_z);
             if (weight>m_Thres)
             {
                 wDir *= weight;
                 idx = index[0]   + outImageSize[0]*(index[1]+1+ outImageSize[1]*index[2]  );
                 dirCont = DirectionContainerType::New();
                 if (m_DirectionsContainer->IndexExists(idx))
                 {
                     dirCont = m_DirectionsContainer->GetElement(idx);
                     if (dirCont.IsNull())
                     {
                         dirCont = DirectionContainerType::New();
                         dirCont->InsertElement(0, wDir);
                         m_DirectionsContainer->InsertElement(idx, dirCont);
                     }
                     else
                         dirCont->InsertElement(dirCont->Size(), wDir);
                 }
                 else
                 {
                     dirCont->InsertElement(0, wDir);
                     m_DirectionsContainer->InsertElement(idx, dirCont);
                 }
             }
 
             wDir = dir;
             weight = (  frac_x)*(  frac_y)*(1-frac_z);
             if (weight>m_Thres)
             {
                 wDir *= weight;
                 idx = index[0]   + outImageSize[0]*(index[1]  + outImageSize[1]*index[2]+outImageSize[1]);
                 dirCont = DirectionContainerType::New();
                 if (m_DirectionsContainer->IndexExists(idx))
                 {
                     dirCont = m_DirectionsContainer->GetElement(idx);
                     if (dirCont.IsNull())
                     {
                         dirCont = DirectionContainerType::New();
                         dirCont->InsertElement(0, wDir);
                         m_DirectionsContainer->InsertElement(idx, dirCont);
                     }
                     else
                         dirCont->InsertElement(dirCont->Size(), wDir);
                 }
                 else
                 {
                     dirCont->InsertElement(0, wDir);
                     m_DirectionsContainer->InsertElement(idx, dirCont);
                 }
             }
 
             wDir = dir;
             weight = (  frac_x)*(1-frac_y)*(1-frac_z);
             if (weight>m_Thres)
             {
                 wDir *= weight;
                 idx = index[0]   + outImageSize[0]*(index[1]+1+ outImageSize[1]*index[2]+outImageSize[1]);
                 dirCont = DirectionContainerType::New();
                 if (m_DirectionsContainer->IndexExists(idx))
                 {
                     dirCont = m_DirectionsContainer->GetElement(idx);
                     if (dirCont.IsNull())
                     {
                         dirCont = DirectionContainerType::New();
                         dirCont->InsertElement(0, wDir);
                         m_DirectionsContainer->InsertElement(idx, dirCont);
                     }
                     else
                         dirCont->InsertElement(dirCont->Size(), wDir);
                 }
                 else
                 {
                     dirCont->InsertElement(0, wDir);
                     m_DirectionsContainer->InsertElement(idx, dirCont);
                 }
             }
 
             wDir = dir;
             weight = (1-frac_x)*(  frac_y)*(  frac_z);
             if (weight>m_Thres)
             {
                 wDir *= weight;
                 idx = index[0]+1 + outImageSize[0]*(index[1]  + outImageSize[1]*index[2]  );
                 dirCont = DirectionContainerType::New();
                 if (m_DirectionsContainer->IndexExists(idx))
                 {
                     dirCont = m_DirectionsContainer->GetElement(idx);
                     if (dirCont.IsNull())
                     {
                         dirCont = DirectionContainerType::New();
                         dirCont->InsertElement(0, wDir);
                         m_DirectionsContainer->InsertElement(idx, dirCont);
                     }
                     else
                         dirCont->InsertElement(dirCont->Size(), wDir);
                 }
                 else
                 {
                     dirCont->InsertElement(0, wDir);
                     m_DirectionsContainer->InsertElement(idx, dirCont);
                 }
             }
 
             wDir = dir;
             weight = (1-frac_x)*(  frac_y)*(1-frac_z);
             if (weight>m_Thres)
             {
                 wDir *= weight;
                 idx = index[0]+1 + outImageSize[0]*(index[1]  + outImageSize[1]*index[2]+outImageSize[1]);
                 dirCont = DirectionContainerType::New();
                 if (m_DirectionsContainer->IndexExists(idx))
                 {
                     dirCont = m_DirectionsContainer->GetElement(idx);
                     if (dirCont.IsNull())
                     {
                         dirCont = DirectionContainerType::New();
                         dirCont->InsertElement(0, wDir);
                         m_DirectionsContainer->InsertElement(idx, dirCont);
                     }
                     else
                         dirCont->InsertElement(dirCont->Size(), wDir);
                 }
                 else
                 {
                     dirCont->InsertElement(0, wDir);
                     m_DirectionsContainer->InsertElement(idx, dirCont);
                 }
             }
 
             wDir = dir;
             weight = (1-frac_x)*(1-frac_y)*(  frac_z);
             if (weight>m_Thres)
             {
                 wDir *= weight;
                 idx = index[0]+1 + outImageSize[0]*(index[1]+1+ outImageSize[1]*index[2]  );
                 dirCont = DirectionContainerType::New();
                 if (m_DirectionsContainer->IndexExists(idx))
                 {
                     dirCont = m_DirectionsContainer->GetElement(idx);
                     if (dirCont.IsNull())
                     {
                         dirCont = DirectionContainerType::New();
                         dirCont->InsertElement(0, wDir);
                         m_DirectionsContainer->InsertElement(idx, dirCont);
                     }
                     else
                         dirCont->InsertElement(dirCont->Size(), wDir);
                 }
                 else
                 {
                     dirCont->InsertElement(0, wDir);
                     m_DirectionsContainer->InsertElement(idx, dirCont);
                 }
             }
 
             wDir = dir;
             weight = (1-frac_x)*(1-frac_y)*(1-frac_z);
             if (weight>m_Thres)
             {
                 wDir *= weight;
                 idx = index[0]+1 + outImageSize[0]*(index[1]+1+ outImageSize[1]*index[2]+outImageSize[1]);
                 dirCont = DirectionContainerType::New();
                 if (m_DirectionsContainer->IndexExists(idx))
                 {
                     dirCont = m_DirectionsContainer->GetElement(idx);
                     if (dirCont.IsNull())
                     {
                         dirCont = DirectionContainerType::New();
                         dirCont->InsertElement(0, wDir);
                         m_DirectionsContainer->InsertElement(idx, dirCont);
                     }
                     else
                         dirCont->InsertElement(dirCont->Size(), wDir);
                 }
                 else
                 {
                     dirCont->InsertElement(0, wDir);
                     m_DirectionsContainer->InsertElement(idx, dirCont);
                 }
             }
         }
         clock.Stop();
     }
 
     vtkSmartPointer<vtkCellArray> m_VtkCellArray = vtkSmartPointer<vtkCellArray>::New();
     vtkSmartPointer<vtkPoints>    m_VtkPoints = vtkSmartPointer<vtkPoints>::New();
 
     itk::ImageRegionIterator<ItkUcharImgType> dirIt(m_NumDirectionsImage, m_NumDirectionsImage->GetLargestPossibleRegion());
     itk::ImageRegionIterator<ItkUcharImgType> crossIt(m_CrossingsImage, m_CrossingsImage->GetLargestPossibleRegion());
 
     m_DirectionImageContainer = DirectionImageContainerType::New();
     unsigned int maxNumDirections = 0;
 
     MITK_INFO << "Clustering directions";
     boost::progress_display disp2(outImageSize[0]*outImageSize[1]*outImageSize[2]);
     for(crossIt.GoToBegin(); !crossIt.IsAtEnd(); ++crossIt)
     {
         ++disp2;
         OutputImageType::IndexType index = crossIt.GetIndex();
         int idx = index[0]+(index[1]+index[2]*outImageSize[1])*outImageSize[0];
 
         if (!m_DirectionsContainer->IndexExists(idx))
         {
             ++dirIt;
             continue;
         }
         DirectionContainerType::Pointer dirCont = m_DirectionsContainer->GetElement(idx);
         if (dirCont.IsNull() || index[0] < 0 || (unsigned long)index[0] >= outImageSize[0] || index[1] < 0 || (unsigned long)index[1] >= outImageSize[1] || index[2] < 0 || (unsigned long)index[2] >= outImageSize[2])
         {
             ++dirIt;
             continue;
         }
 
         std::vector< DirectionType > directions;
 
         for (unsigned int i=0; i<dirCont->Size(); i++)
             if (dirCont->ElementAt(i).magnitude()>0.0001)
                 directions.push_back(dirCont->ElementAt(i));
 
         if (!directions.empty())
             directions = FastClustering(directions);
 
         std::sort( directions.begin(), directions.end(), CompareVectorLengths );
 
         if ( directions.size() > maxNumDirections )
         {
             for (unsigned int i=maxNumDirections; i<std::min<unsigned int>(directions.size(), m_MaxNumDirections); i++)
             {
                 ItkDirectionImageType::Pointer directionImage = ItkDirectionImageType::New();
                 directionImage->SetSpacing( spacing );
                 directionImage->SetOrigin( origin );
                 directionImage->SetDirection( direction );
                 directionImage->SetRegions( imageRegion );
                 directionImage->Allocate();
                 Vector< float, 3 > nullVec; nullVec.Fill(0.0);
                 directionImage->FillBuffer(nullVec);
                 m_DirectionImageContainer->InsertElement(i, directionImage);
             }
             maxNumDirections = std::min<unsigned int>(directions.size(), m_MaxNumDirections);
         }
 
         unsigned int numDir = directions.size();
         if (numDir>m_MaxNumDirections)
             numDir = m_MaxNumDirections;
 
         for (unsigned int i=0; i<numDir; i++)
         {
             vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
             itk::ContinuousIndex<float, 3> center;
             center[0] = index[0];
             center[1] = index[1];
             center[2] = index[2];
             itk::Point<float> worldCenter;
             m_MaskImage->TransformContinuousIndexToPhysicalPoint( center, worldCenter );
             DirectionType dir = directions.at(i);
 
             // set direction image pixel
             ItkDirectionImageType::Pointer directionImage = m_DirectionImageContainer->GetElement(i);
             Vector< float, 3 > pixel;
             pixel.SetElement(0, dir[0]);
             pixel.SetElement(1, dir[1]);
             pixel.SetElement(2, dir[2]);
             directionImage->SetPixel(index, pixel);
 
             // add direction to vector field (with spacing compensation)
             itk::Point<float> worldStart;
             worldStart[0] = worldCenter[0]-dir[0]/2*minSpacing;
             worldStart[1] = worldCenter[1]-dir[1]/2*minSpacing;
             worldStart[2] = worldCenter[2]-dir[2]/2*minSpacing;
             vtkIdType id = m_VtkPoints->InsertNextPoint(worldStart.GetDataPointer());
             container->GetPointIds()->InsertNextId(id);
             itk::Point<float> worldEnd;
             worldEnd[0] = worldCenter[0]+dir[0]/2*minSpacing;
             worldEnd[1] = worldCenter[1]+dir[1]/2*minSpacing;
             worldEnd[2] = worldCenter[2]+dir[2]/2*minSpacing;
             id = m_VtkPoints->InsertNextPoint(worldEnd.GetDataPointer());
             container->GetPointIds()->InsertNextId(id);
             m_VtkCellArray->InsertNextCell(container);
         }
         dirIt.Set(numDir);
         ++dirIt;
     }
 
     vtkSmartPointer<vtkPolyData> directionsPolyData = vtkSmartPointer<vtkPolyData>::New();
     directionsPolyData->SetPoints(m_VtkPoints);
     directionsPolyData->SetLines(m_VtkCellArray);
     m_OutputFiberBundle = mitk::FiberBundleX::New(directionsPolyData);
 }
 
 
 template< class PixelType >
 std::vector< vnl_vector_fixed< double, 3 > > TractsToVectorImageFilter< PixelType >::FastClustering(std::vector< vnl_vector_fixed< double, 3 > >& inDirs)
 {
     std::vector< vnl_vector_fixed< double, 3 > > outDirs;
     if (inDirs.empty())
         return outDirs;
     vnl_vector_fixed< double, 3 > oldMean, currentMean, workingMean;
 
     std::vector< vnl_vector_fixed< double, 3 > > normalizedDirs;
     std::vector< int > touched;
     for (unsigned int i=0; i<inDirs.size(); i++)
     {
         normalizedDirs.push_back(inDirs[i]);
         normalizedDirs.back().normalize();
     }
 
     // initialize
     float max = 0.0;
     touched.resize(inDirs.size(), 0);
     bool free = true;
     currentMean = inDirs[0];  // initialize first seed
     while (free)
     {
         oldMean.fill(0.0);
 
         // start mean-shift clustering
         float angle = 0.0;
         int counter = 0;
         while ((currentMean-oldMean).magnitude()>0.0001)
         {
             counter = 0;
             oldMean = currentMean;
             workingMean = oldMean;
             workingMean.normalize();
             currentMean.fill(0.0);
             for (unsigned int i=0; i<normalizedDirs.size(); i++)
             {
                 angle = dot_product(workingMean, normalizedDirs[i]);
                 if (angle>=m_AngularThreshold)
                 {
                     currentMean += inDirs[i];
                     touched[i] = 1;
                     counter++;
                 }
                 else if (-angle>=m_AngularThreshold)
                 {
                     currentMean -= inDirs[i];
                     touched[i] = 1;
                     counter++;
                 }
             }
         }
 
         // found stable mean
         if (counter>0)
         {
             currentMean /= counter;
             float mag = currentMean.magnitude();
 
             if (mag>0)
             {
                 if (mag>max)
                     max = mag;
 
                 outDirs.push_back(currentMean);
             }
         }
 
         // find next unused seed
         free = false;
         for (unsigned int i=0; i<touched.size(); i++)
             if (touched[i]==0)
             {
                 currentMean = inDirs[i];
                 free = true;
             }
     }
 
     if (m_NormalizeVectors)
         for (unsigned int i=0; i<outDirs.size(); i++)
             outDirs[i].normalize();
     else if (max>0)
         for (unsigned int i=0; i<outDirs.size(); i++)
             outDirs[i] /= max;
 
     if (inDirs.size()==outDirs.size())
         return outDirs;
     else
         return FastClustering(outDirs);
 }
 
 
 template< class PixelType >
 std::vector< vnl_vector_fixed< double, 3 > > TractsToVectorImageFilter< PixelType >::Clustering(std::vector< vnl_vector_fixed< double, 3 > >& inDirs)
 {
     std::vector< vnl_vector_fixed< double, 3 > > outDirs;
     if (inDirs.empty())
         return outDirs;
     vnl_vector_fixed< double, 3 > oldMean, currentMean, workingMean;
 
     std::vector< vnl_vector_fixed< double, 3 > > normalizedDirs;
     std::vector< int > touched;
-    for (int i=0; i<inDirs.size(); i++)
+    for (std::size_t i=0; i<inDirs.size(); i++)
     {
         normalizedDirs.push_back(inDirs[i]);
         normalizedDirs.back().normalize();
     }
 
     // initialize
     float max = 0.0;
     touched.resize(inDirs.size(), 0);
-    for (int j=0; j<inDirs.size(); j++)
+    for (std::size_t j=0; j<inDirs.size(); j++)
     {
         currentMean = inDirs[j];
         oldMean.fill(0.0);
 
         // start mean-shift clustering
         float angle = 0.0;
         int counter = 0;
         while ((currentMean-oldMean).magnitude()>0.0001)
         {
             counter = 0;
             oldMean = currentMean;
             workingMean = oldMean;
             workingMean.normalize();
             currentMean.fill(0.0);
-            for (int i=0; i<normalizedDirs.size(); i++)
+            for (std::size_t i=0; i<normalizedDirs.size(); i++)
             {
                 angle = dot_product(workingMean, normalizedDirs[i]);
                 if (angle>=m_AngularThreshold)
                 {
                     currentMean += inDirs[i];
                     counter++;
                 }
                 else if (-angle>=m_AngularThreshold)
                 {
                     currentMean -= inDirs[i];
                     counter++;
                 }
             }
         }
 
         // found stable mean
         if (counter>0)
         {
             bool add = true;
             vnl_vector_fixed< double, 3 > normMean = currentMean;
             normMean.normalize();
-            for (int i=0; i<outDirs.size(); i++)
+            for (std::size_t i=0; i<outDirs.size(); i++)
             {
                 vnl_vector_fixed< double, 3 > dir = outDirs[i];
                 dir.normalize();
                 if ((normMean-dir).magnitude()<=0.0001)
                 {
                     add = false;
                     break;
                 }
             }
 
             currentMean /= counter;
             if (add)
             {
                 float mag = currentMean.magnitude();
                 if (mag>0)
                 {
                     if (mag>max)
                         max = mag;
 
                     outDirs.push_back(currentMean);
                 }
             }
         }
     }
 
     if (m_NormalizeVectors)
-        for (int i=0; i<outDirs.size(); i++)
+        for (std::size_t i=0; i<outDirs.size(); i++)
             outDirs[i].normalize();
     else if (max>0)
-        for (int i=0; i<outDirs.size(); i++)
+        for (std::size_t i=0; i<outDirs.size(); i++)
             outDirs[i] /= max;
 
     if (inDirs.size()==outDirs.size())
         return outDirs;
     else
         return FastClustering(outDirs);
 }
 
 
 template< class PixelType >
 TractsToVectorImageFilter< PixelType >::DirectionContainerType::Pointer TractsToVectorImageFilter< PixelType >::MeanShiftClustering(DirectionContainerType::Pointer dirCont)
 {
     DirectionContainerType::Pointer container = DirectionContainerType::New();
 
     float max = 0;
     for (DirectionContainerType::ConstIterator it = dirCont->Begin(); it!=dirCont->End(); ++it)
     {
         vnl_vector_fixed<double, 3> mean = ClusterStep(dirCont, it.Value());
 
         if (mean.is_zero())
             continue;
         bool addMean = true;
 
         for (DirectionContainerType::ConstIterator it2 = container->Begin(); it2!=container->End(); ++it2)
         {
             vnl_vector_fixed<double, 3> dir = it2.Value();
             float angle = fabs(dot_product(mean, dir)/(mean.magnitude()*dir.magnitude()));
             if (angle>=m_Epsilon)
             {
                 addMean = false;
                 break;
             }
         }
 
         if (addMean)
         {
             if (m_NormalizeVectors)
                 mean.normalize();
             else if (mean.magnitude()>max)
                 max = mean.magnitude();
             container->InsertElement(container->Size(), mean);
         }
     }
 
     // max normalize voxel directions
     if (max>0 && !m_NormalizeVectors)
-        for (int i=0; i<container->Size(); i++)
+        for (std::size_t i=0; i<container->Size(); i++)
             container->ElementAt(i) /= max;
 
     if (container->Size()<dirCont->Size())
         return MeanShiftClustering(container);
     else
         return container;
 }
 
 
 template< class PixelType >
 vnl_vector_fixed<double, 3> TractsToVectorImageFilter< PixelType >::ClusterStep(DirectionContainerType::Pointer dirCont, vnl_vector_fixed<double, 3> currentMean)
 {
     vnl_vector_fixed<double, 3> newMean; newMean.fill(0);
 
     for (DirectionContainerType::ConstIterator it = dirCont->Begin(); it!=dirCont->End(); ++it)
     {
         vnl_vector_fixed<double, 3> dir = it.Value();
         float angle = dot_product(currentMean, dir)/(currentMean.magnitude()*dir.magnitude());
         if (angle>=m_AngularThreshold)
             newMean += dir;
         else if (-angle>=m_AngularThreshold)
             newMean -= dir;
     }
 
     if (fabs(dot_product(currentMean, newMean)/(currentMean.magnitude()*newMean.magnitude()))>=m_Epsilon || newMean.is_zero())
         return newMean;
     else
         return ClusterStep(dirCont, newMean);
 }
 }
 
 
 
diff --git a/Modules/DiffusionImaging/Quantification/Algorithms/itkSkeletonizationFilter.txx b/Modules/DiffusionImaging/Quantification/Algorithms/itkSkeletonizationFilter.txx
index ead0836376..eb43aff6d5 100644
--- a/Modules/DiffusionImaging/Quantification/Algorithms/itkSkeletonizationFilter.txx
+++ b/Modules/DiffusionImaging/Quantification/Algorithms/itkSkeletonizationFilter.txx
@@ -1,306 +1,306 @@
 /*===================================================================
 
 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 _itkSkeletonizationFilter_txx
 #define _itkSkeletonizationFilter_txx
 
 #include "itkSkeletonizationFilter.h"
 #include "mitkProgressBar.h"
 #include <limits>
 #include <time.h>
 
 namespace itk
 {
 
   template< class TInputImage, class TOutputImage >
   SkeletonizationFilter<TInputImage, TOutputImage>::SkeletonizationFilter()
   {
     m_DirectionImage = VectorImageType::New();
   }
 
   template< class TInputImage, class TOutputImage >
   SkeletonizationFilter<TInputImage, TOutputImage>::~SkeletonizationFilter()
   {
 
   }
 
 
 
 
   template< class TInputImage, class TOutputImage >
   void SkeletonizationFilter<TInputImage, TOutputImage>::GenerateData()
   {
     //----------------------------------------------------------------------//
     //  Progress bar                                                        //
     //----------------------------------------------------------------------//
     mitk::ProgressBar::GetInstance()->AddStepsToDo( 3 );
 
     std::cout << "Skeletonize" << std::endl;
 
 
     const InputImageType* faImage = this->GetInput();
     typename InputImageType::SizeType size = faImage->GetRequestedRegion().GetSize();
 
     //typename RealImageType::SizeType size = m_FaImage->GetRequestedRegion().GetSize();
 
     m_DirectionImage->SetRegions(faImage->GetRequestedRegion());
     m_DirectionImage->SetDirection(faImage->GetDirection());
     m_DirectionImage->SetSpacing(faImage->GetSpacing());
     m_DirectionImage->SetOrigin(faImage->GetOrigin());
     m_DirectionImage->Allocate();
     m_DirectionImage->FillBuffer(0.0);
 
 
     for(unsigned int z=1; z<size[2]-1; z++) for(unsigned int y=1; y<size[1]-1; y++) for(unsigned int x=1; x<size[0]-1; x++)
     {
       typename InputImageType::IndexType ix;
       ix[0]=x; ix[1]=y; ix[2]=z;
       float theval = faImage->GetPixel(ix);
 
       if(theval != 0)
       {
 
         /* Calculate point of gravity. We will consider each 3x3x3 neighbourhood as a unit cube. The center
          * point of each voxel will be a multiplicative of 1/6. The center of the unit cube is 3/6 = 1/2/
          */
 
         float cogX = 0.0; float cogY = 0.0; float cogZ = 0.0; float sum = 0.0; float l;
         int vecX = 0; int vecY = 0; int vecZ = 0;
 
         for(int dz=-1; dz<=1; dz++) for(int dy=-1; dy<=1; dy++) for(int dx=-1; dx<=1;dx++)
         {
           typename InputImageType::IndexType p;
           p[0] = x+dx; p[1] = y+dy; p[2] = z+dz;
           float mass = faImage->GetPixel(p);
 
           sum += mass;
           cogX += (float)dx*mass; cogY += (float)dy*mass; cogZ += (float)dz*mass;
         }
 
         cogX /= sum; cogY /= sum; cogZ /= sum;
         l = sqrt(cogX*cogX + cogY*cogY + cogZ*cogZ);
 
         if (l > 0.1) /* is CofG far enough away from centre voxel? */
         {
           vecX = std::max(std::min(round(cogX/l),1),-1);
           vecY = std::max(std::min(round(cogY/l),1),-1);
           vecZ = std::max(std::min(round(cogZ/l),1),-1);
         }
         else
           // Find direction of max curvature
         {
 
           float maxcost=0, centreval=2*theval;
           for(int zz=0; zz<=1; zz++) // note - starts at zero as we're only searching half the voxels
           {
             for(int yy=-1; yy<=1; yy++)
             {
               for(int xx=-1; xx<=1; xx++)
               {
                 if ( (zz==1) || (yy==1) || ((yy==0)&&(xx==1)) )
                 {
                   float weighting = pow( (float)(xx*xx+yy*yy+zz*zz) , (float)-0.7 ); // power is arbitrary: maybe test other functions here
 
 
                   typename InputImageType::IndexType i,j;
                   i[0] = x+xx; i[1] = y+yy; i[2] = z+zz;
                   j[0] = x-xx; j[1] = y-yy; j[2] = z-zz;
                   float cost = weighting * ( centreval
                     - (float)faImage->GetPixel(i)
                     - (float)faImage->GetPixel(j));
 
                   if (cost>maxcost)
                   {
                     maxcost=cost;
                     vecX=xx;
                     vecY=yy;
                     vecZ=zz;
                   }
                 }
               }
             }
           }
         }
 
         VectorType vec;
         vec[0] = vecX; vec[1] = vecY; vec[2]=vecZ;
         m_DirectionImage->SetPixel(ix, vec);
 
       }
     }
 
     mitk::ProgressBar::GetInstance()->Progress();
 
 
     // Smooth m_DirectionImage and store in directionSmoothed by finding the
     // mode in a 3*3 neighbourhoud
     VectorImageType::Pointer directionSmoothed = VectorImageType::New();
     directionSmoothed->SetRegions(faImage->GetRequestedRegion());
     directionSmoothed->SetDirection(faImage->GetDirection());
     directionSmoothed->SetSpacing(faImage->GetSpacing());
     directionSmoothed->SetOrigin(faImage->GetOrigin());
     directionSmoothed->Allocate();
 
     VectorImageType::PixelType p;
     p[0]=0; p[1]=0; p[2]=0;
     directionSmoothed->FillBuffer(p);
 
 
 
     for(unsigned int z=1; z<size[2]-1; z++) for(unsigned int y=1; y<size[1]-1; y++) for(unsigned int x=1; x<size[0]-1; x++)
     {
       VectorImageType::IndexType ix;
       ix[0]=x; ix[1]=y; ix[2]=z;
 
       // Find the vector that occured most
       int* localsum = new int[27];
       int localmax=0, xxx, yyy, zzz;
 
       for(int zz=0; zz<27; zz++) localsum[zz]=0;
 
       for(int zz=-1; zz<=1; zz++) for(int yy=-1; yy<=1; yy++) for(int xx=-1; xx<=1; xx++)
       {
         VectorImageType::IndexType i;
         i[0] = x+xx; i[1] = y+yy; i[2] = z+zz;
         VectorType v = m_DirectionImage->GetPixel(i);
         xxx = v[0];
         yyy = v[1];
         zzz = v[2];
 
         localsum[(1+zzz)*9+(1+yyy)*3+1+xxx]++;
         localsum[(1-zzz)*9+(1-yyy)*3+1-xxx]++;
       }
 
       for(int zz=-1; zz<=1; zz++) for(int yy=-1; yy<=1; yy++) for(int xx=-1; xx<=1; xx++)
       {
         if (localsum[(1+zz)*9+(1+yy)*3+1+xx]>localmax)
         {
           localmax=localsum[(1+zz)*9+(1+yy)*3+1+xx];
           VectorType v;
           v[0] = xx; v[1] = yy; v[2] = zz;
           directionSmoothed->SetPixel(ix, v);
         }
       }
 
       delete localsum;
 
     }
 
     m_DirectionImage = directionSmoothed;
 
     mitk::ProgressBar::GetInstance()->Progress();
 
     // Do non-max-suppression in the direction of perp and set as output of the filter
     typename OutputImageType::Pointer outputImg = OutputImageType::New();
     outputImg->SetRegions(faImage->GetRequestedRegion());
     outputImg->SetDirection(faImage->GetDirection());
     outputImg->SetSpacing(faImage->GetSpacing());
     outputImg->SetOrigin(faImage->GetOrigin());
     outputImg->Allocate();
     outputImg->FillBuffer(0.0);
 
     for(unsigned int z=1; z<size[2]-1; z++) for(unsigned int y=1; y<size[1]-1; y++) for(unsigned int x=1; x<size[0]-1; x++)
     {
 
       typename InputImageType::IndexType ix;
       ix[0]=x; ix[1]=y; ix[2]=z;
 
       float theval = faImage->GetPixel(ix);
       VectorType v = directionSmoothed->GetPixel(ix);
 
       typename VectorImageType::IndexType i;
 
       i[0] = x-v[0]; i[1] = y-v[1]; i[2] = z-v[2];
       float min = faImage->GetPixel(i);
 
       i[0] = x+v[0]; i[1] = y+v[1]; i[2] = z+v[2];
       float plus = faImage->GetPixel(i);
 
       i[0] = x-2*v[0]; i[1] = y-2*v[1]; i[2] = z-2*v[2];
       float minmin = faImage->GetPixel(i);
 
       i[0] = x+2*v[0]; i[1] = y+2*v[1]; i[2] = z+2*v[2];
       float plusplus = faImage->GetPixel(i);
 
       if( ((v[0]!=0) || (v[1]!=0) || (v[2]!=0)) &&
            theval >= plus && theval >  min && theval >= plusplus && theval >  minmin  )
       {
         outputImg->SetPixel(ix, theval);
       }
 
     }
 
 
 
     Superclass::SetNthOutput( 0, outputImg );
     mitk::ProgressBar::GetInstance()->Progress();
 
   }
 
 
 
   // Can provide a vector image to visualize the gradient image used in the search for local maxima.
   template< class TInputImage, class TOutputImage >
   itk::VectorImage<int, 3>::Pointer SkeletonizationFilter<TInputImage, TOutputImage>::GetGradientImage()
   {
     GradientImageType::Pointer gradImg = GradientImageType::New();
 
     if(m_DirectionImage.IsNotNull())
     {
       gradImg->SetSpacing(m_DirectionImage->GetSpacing());
       gradImg->SetOrigin(m_DirectionImage->GetOrigin());
       gradImg->SetDirection(m_DirectionImage->GetDirection());
       gradImg->SetRegions(m_DirectionImage->GetLargestPossibleRegion().GetSize());
       gradImg->SetVectorLength(3);
       gradImg->Allocate();
 
 
       VectorImageType::SizeType size = m_DirectionImage->GetLargestPossibleRegion().GetSize();
 
-      for(int i=0; i<size[0]; i++)
+      for(std::size_t i=0; i<size[0]; i++)
       {
-        for(int j=0; j<size[1]; j++)
+        for(std::size_t j=0; j<size[1]; j++)
         {
-          for(int k=0; k<size[2]; k++)
+          for(std::size_t k=0; k<size[2]; k++)
           {
             itk::Index<3> ix;
             ix[0] = i;
             ix[1] = j;
             ix[2] = k;
 
             VectorType vec = m_DirectionImage->GetPixel(ix);
 
             itk::VariableLengthVector<int> pixel;
             pixel.SetSize(3);
             pixel.SetElement(0, vec.GetElement(0));
             pixel.SetElement(1, vec.GetElement(1));
             pixel.SetElement(2, vec.GetElement(2));
 
             gradImg->SetPixel(ix, pixel);
 
           }
         }
       }
 
     }
 
 
 
 
     return gradImg;
   }
 
 }
 #endif // _itkSkeletonizationFilter_txx