diff --git a/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp b/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp
index 5e1ba61ac8..86c834df94 100644
--- a/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp
+++ b/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp
@@ -1,885 +1,791 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
-Copyright (c) German Cancer Research Center, 
+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 
+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>
 
-#if BOOST_VERSION / 100000 > 0
-#if BOOST_VERSION / 100 % 1000 > 34
-#include <boost/math/special_functions/legendre.hpp>
-#endif
-#endif
+#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() :
+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)
-  {
+{
     // 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;
+}
+
+
+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++)
-          {
+        for(int i=0; i<NrOdfDirections; i++)
+        {
             sum += odf[i];
-          }
-          if(sum>0)
+        }
+        if(sum>0)
             odf /= sum;
 
-          return odf;
-          break;
-        }
-      case QBAR_B_ZERO_B_VALUE:
+        return odf;
+        break;
+    }
+    case QBAR_B_ZERO_B_VALUE:
+    {
+        for(int i=0; i<NrOdfDirections; i++)
         {
-          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:
+        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++)
         {
-          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++)
+        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;
+        break;
+    }
+    case QBAR_NONNEG_SOLID_ANGLE:
+    {
+        break;
     }
+    }
+
+    return odf;
+}
 
 
-    template<
-      class TReferenceImagePixelType,
-      class TGradientImagePixelType,
-      class TOdfPixelType,
+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 )
+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++)
         {
-        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)
+            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] = 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;
+            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)
+            if (vec[i]<0)
                 vec[i] = m_Delta1;
-              else if (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)
+            else if (vec[i]>=1)
                 vec[i] = 1-m_Delta2/2;
-              else if (vec[i]>=1-m_Delta2)
+            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;
-          }
+            vec[i] = log(-log(vec[i]));
         }
-
         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 < 6 )
-        {
-          itkExceptionMacro( << "At least 6 gradient directions are required" );
-        }
+        break;
+    }
+    }
 
-        // 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 );
-        }
+    return vec;
+}
 
-        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>
+::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 < 6 )
+    {
+        itkExceptionMacro( << "At least 6 gradient directions are required" );
+    }
 
-      }
+    // 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 );
+    }
 
-      template< class T, class TG, class TO, int L, int NODF>
-      void AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
-        ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
-        int )
-      {
-        typename OutputImageType::Pointer outputImage =
-          static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0));
+    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;
+}
 
-        ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread);
-        oit.GoToBegin();
+template< class T, class TG, class TO, int L, int NODF>
+void AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
+::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
+                       int )
+{
+    typename OutputImageType::Pointer outputImage =
+            static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0));
 
-        ImageRegionIterator< BZeroImageType > oit2(m_BZeroImage, outputRegionForThread);
-        oit2.GoToBegin();
+    ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread);
+    oit.GoToBegin();
 
-        ImageRegionIterator< FloatImageType > oit3(m_ODFSumImage, outputRegionForThread);
-        oit3.GoToBegin();
+    ImageRegionIterator< BZeroImageType > oit2(m_BZeroImage, outputRegionForThread);
+    oit2.GoToBegin();
 
-        ImageRegionIterator< CoefficientImageType > oit4(m_CoefficientImage, outputRegionForThread);
-        oit4.GoToBegin();
+    ImageRegionIterator< FloatImageType > oit3(m_ODFSumImage, outputRegionForThread);
+    oit3.GoToBegin();
 
-        typedef ImageRegionConstIterator< GradientImagesType > GradientIteratorType;
-        typedef typename GradientImagesType::PixelType         GradientVectorType;
-        typename GradientImagesType::Pointer gradientImagePointer = NULL;
+    ImageRegionIterator< CoefficientImageType > oit4(m_CoefficientImage, outputRegionForThread);
+    oit4.GoToBegin();
 
-        // 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) );
+    typedef ImageRegionConstIterator< GradientImagesType > GradientIteratorType;
+    typedef typename GradientImagesType::PixelType         GradientVectorType;
+    typename GradientImagesType::Pointer gradientImagePointer = NULL;
 
-        GradientIteratorType git(gradientImagePointer, outputRegionForThread );
-        git.GoToBegin();
+    // 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) );
 
-        // 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
+    GradientIteratorType git(gradientImagePointer, outputRegionForThread );
+    git.GoToBegin();
 
-        for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin();
-          gdcit != this->m_GradientDirectionContainer->End(); ++gdcit)
-        {
-          if(gdcit.Value().one_norm() <= 0.0)
-          {
+    // 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
-          {
+        else
             gradientind.push_back(gdcit.Index());
-          }
-        }
+    }
 
-        if( m_DirectionsDuplicated )
-        {
-          int gradIndSize = gradientind.size();
-          for(int i=0; i<gradIndSize; i++)
+    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();
+    while( !git.IsAtEnd() )
+    {
+        GradientVectorType b = git.Get();
 
-          typename NumericTraits<ReferencePixelType>::AccumulateType b0 = NumericTraits<ReferencePixelType>::Zero;
+        typename NumericTraits<ReferencePixelType>::AccumulateType b0 = NumericTraits<ReferencePixelType>::Zero;
 
-          // Average the baseline image pixels
-          for(unsigned int i = 0; i < baselineind.size(); ++i)
-          {
+        // Average the baseline image pixels
+        for(unsigned int i = 0; i < baselineind.size(); ++i)
+        {
             b0 += b[baselineind[i]];
-          }
-          b0 /= this->m_NumberOfBaselineImages;
+        }
+        b0 /= this->m_NumberOfBaselineImages;
 
-          OdfPixelType odf(0.0);
-          typename CoefficientImageType::PixelType coeffPixel(0.0);
-          vnl_vector<TO> B(m_NumberOfGradientDirections);
+        OdfPixelType odf(0.0);
+        typename CoefficientImageType::PixelType coeffPixel(0.0);
+        vnl_vector<TO> B(m_NumberOfGradientDirections);
 
-          if( (b0 != 0) && (b0 >= m_Threshold) )
-          {
+        if( (b0 != 0) && (b0 >= m_Threshold) )
+        {
             for( unsigned int i = 0; i< m_NumberOfGradientDirections; i++ )
             {
-              B[i] = static_cast<TO>(b[gradientind[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();
+                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
+                /** 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");
-
+                itkExceptionMacro( << "Nonnegative Solid Angle not yet implemented");
             }
             else
             {
-              odf = ( (*m_ReconstructionMatrix) * B ).data_block();
+                odf = ( (*m_ReconstructionMatrix) * B ).data_block();
             }
             odf = Normalize(odf, b0);
 
-          }
+        }
 
-          oit.Set( odf );
-          oit2.Set( b0 );
-          float sum = 0;
-          for (int k=0; k<odf.Size(); k++)
+        oit.Set( odf );
+        oit2.Set( b0 );
+        float sum = 0;
+        for (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
-        }
+        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++)
+    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(int j=0; j<A.columns(); j++)
         {
-          for(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>
-      double AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
-        ::factorial(int number) {
-          if(number <= 1) return 1;
-          double result = 1.0;
-          for(int i=1; i<=number; i++)
-            result *= i;
-          return result;
-      }
-
-      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 *cart)
-      {
-        double phi, th, rad;
-        rad = sqrt(x*x+y*y+z*z);
-        th = atan2(z,sqrt(x*x+y*y));
-        phi = atan2(y,x);
-        th = -th+QBALL_ANAL_RECON_PI/2;
-        phi = -phi+QBALL_ANAL_RECON_PI;
-        cart[0] = phi;
-        cart[1] = th;
-        cart[2] = rad;
-      }
-
-      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>
-      double AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
-        ::spherical_harmonic(int m,int l,double theta,double phi, bool complexPart)
-      {
-        if( theta < 0 || theta > QBALL_ANAL_RECON_PI )
-        {
-          std::cout << "bad range" << std::endl;
-          return 0;
-        }
-
-        if( phi < 0 || phi > 2*QBALL_ANAL_RECON_PI )
-        {
-          std::cout << "bad range" << std::endl;
-          return 0;
+                myfile << ";";
         }
+    }
+    myfile << "];";
+    myfile.close();
 
-        double pml = 0;
-        double fac1 = factorial(l+m);
-        double fac2 = factorial(l-m);
-
-        if( m<0 )
-        {
-#if BOOST_VERSION / 100000 > 0
-#if BOOST_VERSION / 100 % 1000 > 34
-          pml = ::boost::math::legendre_p<double>(l, -m, cos(theta));
-#else
-          std::cout << "ERROR: Boost 1.35 minimum required" << std::endl;
-#endif
-#else
-          std::cout << "ERROR: Boost 1.35 minimum required" << std::endl;
-#endif
-          double mypow = pow(-1.0,-m);
-          double myfac = (fac1/fac2);
-          pml *= mypow*myfac;
-        }
-        else
-        {
-#if BOOST_VERSION / 100000 > 0
-#if BOOST_VERSION / 100 % 1000 > 34
-          pml = ::boost::math::legendre_p<double>(l, m, cos(theta));
-#endif
-#endif
-        }
+    myfile.imbue( originalLocale );
+}
 
-        //std::cout << "legendre(" << l << "," << m << "," << cos(theta) << ") = " << pml << std::endl;
+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);
 
-        double retval = sqrt(((2.0*(double)l+1.0)/(4.0*QBALL_ANAL_RECON_PI))*(fac2/fac1)) * pml;
-        if( !complexPart )
-        {
-          retval *= cos(m*phi);
-        }
-        else
-        {
-          retval *= sin(m*phi);
-        }
-        //std::cout << retval << std::endl;
-        return retval;
-      }
-
-      template< class T, class TG, class TO, int L, int NODF>
-      double AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
-        ::Yj(int m, int k, double theta, double phi)
-      {
-        if( -k <= m && m < 0 )
-        {
-          return sqrt(2.0) * spherical_harmonic(m,k,theta,phi,false);
-        }
-
-        if( m == 0 )
-          return spherical_harmonic(0,k,theta,phi,false);
+    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;
+}
 
-        if( 0 < m && m <= k )
-        {
-          return sqrt(2.0) * spherical_harmonic(m,k,theta,phi,true);
-        }
+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)
+{
+    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);
+
+    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()
-      {
+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;
+    //for(int i=-6;i<7;i++)
+    //  std::cout << boost::math::legendre_p<double>(6, i, 0.65657) << std::endl;
 
-        if( m_NumberOfGradientDirections < 6 )
-        {
-          itkExceptionMacro( << "Not enough gradient directions supplied. Need to supply at least 6" );
-        }
+    if( m_NumberOfGradientDirections < 6 )
+    {
+        itkExceptionMacro( << "Not enough gradient directions supplied. Need to supply at least 6" );
+    }
 
-        {
-          // check for duplicate diffusion gradients
-          bool warning = false;
-          for(GradientDirectionContainerType::ConstIterator gdcit1 = this->m_GradientDirectionContainer->Begin();
+    {
+        // 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)
+                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(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();
+        }
+
+        // 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 += gdcit1.Value();
+                count ++;
             }
-          }
-          centerMass /= count;
-          if(centerMass.two_norm() > 0.1)
-          {
+        }
+        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);
+    vnl_matrix<double> *Q = new vnl_matrix<double>(3, m_NumberOfGradientDirections);
 
-        {
-          int i = 0;
-          for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin();
+    {
+        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];
-              }
             }
-          }
         }
+        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++)
+    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 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);
-              //std::cout << "B(" << i << "," << j << ") = Yj(" << m << "," << k << "," << th << "," << phi << ") = " << (*B)(i,j) << std::endl;
+                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);
             }
-          }
         }
+    }
 
-        //vnl_matrix<double> temp((*_L)*(*_L));
-        LL->update(*_L);
-        *LL *= *_L;
-        //tofile2(LL,"LL");
+    //vnl_matrix<double> temp((*_L)*(*_L));
+    LL->update(*_L);
+    *LL *= *_L;
+    //tofile2(LL,"LL");
 
-        for(int i=0; i<m_NumberCoefficients; i++)
+    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)
         {
-          // 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]);
+            (*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)
+        else
         {
-        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;
+            (*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");
+    //tofile2(&temp,"A");
 
-        m_CoeffReconstructionMatrix = new vnl_matrix<TO>(m_NumberCoefficients,m_NumberOfGradientDirections);
-        for(int i=0; i<m_NumberCoefficients; i++)
+    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++)
         {
-          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
+    // 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();
+    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];
-        }
+    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 i=0; i<NOdfDirections; i++)
+    {
+        for(int k=0; k<=l; k+=2)
         {
-          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);
-              (*sphericalHarmonicBasisMatrix2)(i,j) = (*m_SphericalHarmonicBasisMatrix)(i,j);
+                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);
+                (*sphericalHarmonicBasisMatrix2)(i,j) = (*m_SphericalHarmonicBasisMatrix)(i,j);
             }
-          }
         }
+    }
 
-        m_ReconstructionMatrix = new vnl_matrix<TO>(NOdfDirections,m_NumberOfGradientDirections);
-        *m_ReconstructionMatrix = (*m_SphericalHarmonicBasisMatrix) * (*m_CoeffReconstructionMatrix);
-
-      }
+    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( GradientDirectionContainerType *gradientDirection,
-        const GradientImagesType *gradientImage )
-      {
-        this->m_GradientDirectionContainer = gradientDirection;
+}
 
-        unsigned int numImages = gradientDirection->Size();
-        this->m_NumberOfBaselineImages = 0;
-        for(GradientDirectionContainerType::Iterator it = this->m_GradientDirectionContainer->Begin();
-          it != this->m_GradientDirectionContainer->End(); it++)
+template< class T, class TG, class TO, int L, int NODF>
+void AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
+::SetGradientImage( GradientDirectionContainerType *gradientDirection,
+                    const GradientImagesType *gradientImage )
+{
+    this->m_GradientDirectionContainer = gradientDirection;
+
+    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)
         {
-          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 )
+        else // Normalize non-zero gradient directions
         {
-          itkExceptionMacro( << m_NumberOfGradientDirections << " gradients + " << this->m_NumberOfBaselineImages
-            << "baselines = " << m_NumberOfGradientDirections + this->m_NumberOfBaselineImages
-            << " directions specified but image has " << gradientImage->GetVectorLength()
-            << " components.");
+            it.Value() = it.Value() / it.Value().two_norm();
         }
+    }
 
-        this->ProcessObject::SetNthInput( 0,
-          const_cast< GradientImagesType* >(gradientImage) );
-
-      }
+    this->m_NumberOfGradientDirections = numImages - this->m_NumberOfBaselineImages;
 
-      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);
+    // 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.");
+    }
 
-        Superclass::PrintSelf(os,indent);
+    this->ProcessObject::SetNthInput( 0, const_cast< GradientImagesType* >(gradientImage) );
+}
 
-        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 );
-      }
+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/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.h b/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.h
index c6ab682894..2a837f283a 100644
--- a/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.h
+++ b/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.h
@@ -1,304 +1,282 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
-Copyright (c) German Cancer Research Center, 
+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 
+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_h_
 #define __itkAnalyticalDiffusionQballReconstructionImageFilter_h_
 
 #include "itkImageToImageFilter.h"
 #include "vnl/vnl_vector_fixed.h"
 #include "vnl/vnl_matrix.h"
 #include "vnl/algo/vnl_svd.h"
 #include "itkVectorContainer.h"
 #include "itkVectorImage.h"
 
 
 namespace itk{
 /** \class AnalyticalDiffusionQballReconstructionImageFilter
  * \brief This class takes as input one or more reference image (acquired in the
  * absence of diffusion sensitizing gradients) and 'n' diffusion
  * weighted images and their gradient directions and computes an image of
  * orientation distribution function coefficients in a spherical harmonic basis.
  *
  * \par Inputs and Usage
  * \par
  * When you have the 'n' gradient and one or more reference images in a single
  * multi-component image (VectorImage), you can specify the images as
  * \code
  *       filter->SetGradientImage( directionsContainer, vectorImage );
  * \endcode
  * Note that this method is used to specify both the reference and gradient images.
  * This is convenient when the DWI images are read in using the
  * <a href="http://wiki.na-mic.org/Wiki/index.php/NAMIC_Wiki:DTI:Nrrd_format">NRRD</a>
  * format. Like the Nrrd format, the reference images are those components of the
  * vectorImage whose gradient direction is (0,0,0). If more than one reference image
  * is present, they are averaged prior to the reconstruction.
  *
  * \par Outputs
  * The output image is an image of vectors that must be understood as ODFs:
  * \code
  *       Image< Vector< TPixelType, OdfNrDirections >, 3 >
  * \endcode
  *
  * \par Parameters
  * \li Threshold -  Threshold on the reference image data. The output ODF will
  * be a null pdf for pixels in the reference image that have a value less
  * than this.
  * \li BValue - See the documentation of SetBValue().
  * \li At least 6 gradient images must be specified for the filter to be able
  * to run. If the input gradient directions g_i are majorly sampled on one half
  * of the sqhere, then each input image I_i will be duplicated and assign -g_i
  * in order to guarantee stability of the algorithm.
  * \li OdfDirections - directions of resulting orientation distribution function
  * \li EquatorNrSamplingPoints - number of sampling points on equator when
  * performing Funk Radeon Transform (FRT)
  * \li BasisFunctionCenters - the centers of the basis functions are used for
  * the sRBF (spherical radial basis functions interpolation). If not set, they
  * will be defaulted to equal m_EquatorNrSamplingPoints
  *
  * \par Template parameters
  * The class is templated over
  * \li the pixel type of the reference and gradient images
  * (expected to be scalar data types)
  * \li the internal representation of the ODF pixels (double, float etc).
  * \li the number of OdfDirections
  * \li the number of basis function centers for the sRBF
  *
  * \par References:
  * \li<a href="http://www3.interscience.wiley.com/cgi-bin/fulltext/109800293/PDFSTART">[1]</a>
  * <em>Tuch DS,
  * "Q-ball imaging", Magn Reson Med. 2004 Dec;52(6):1358-72.</em>
  *
  */
 
 template< class TReferenceImagePixelType,
           class TGradientImagePixelType,
           class TOdfPixelType,
           int NOrderL,
           int NrOdfDirections>
 class AnalyticalDiffusionQballReconstructionImageFilter :
-  public ImageToImageFilter< Image< TReferenceImagePixelType, 3 >,
-                             Image< Vector< TOdfPixelType, NrOdfDirections >, 3 > >
+        public ImageToImageFilter< Image< TReferenceImagePixelType, 3 >,
+        Image< Vector< TOdfPixelType, NrOdfDirections >, 3 > >
 {
 
 public:
 
-  enum Normalization {
-    QBAR_STANDARD,
-    QBAR_B_ZERO_B_VALUE,
-    QBAR_B_ZERO,
-    QBAR_NONE,
-    QBAR_ADC_ONLY,
-    QBAR_RAW_SIGNAL,
-    QBAR_SOLID_ANGLE,
-    QBAR_NONNEG_SOLID_ANGLE
-  };
+    enum Normalization {
+        QBAR_STANDARD,
+        QBAR_B_ZERO_B_VALUE,
+        QBAR_B_ZERO,
+        QBAR_NONE,
+        QBAR_ADC_ONLY,
+        QBAR_RAW_SIGNAL,
+        QBAR_SOLID_ANGLE,
+        QBAR_NONNEG_SOLID_ANGLE
+    };
 
-  typedef AnalyticalDiffusionQballReconstructionImageFilter Self;
-  typedef SmartPointer<Self>                      Pointer;
-  typedef SmartPointer<const Self>                ConstPointer;
-  typedef ImageToImageFilter< Image< TReferenceImagePixelType, 3>,
-          Image< Vector< TOdfPixelType, NrOdfDirections >, 3 > >
-                          Superclass;
+    typedef AnalyticalDiffusionQballReconstructionImageFilter Self;
+    typedef SmartPointer<Self>                      Pointer;
+    typedef SmartPointer<const Self>                ConstPointer;
+    typedef ImageToImageFilter< Image< TReferenceImagePixelType, 3>,
+    Image< Vector< TOdfPixelType, NrOdfDirections >, 3 > > Superclass;
 
-   /** Method for creation through the object factory. */
-  itkNewMacro(Self);
+    /** Method for creation through the object factory. */
+    itkNewMacro(Self)
 
-  /** Runtime information support. */
-  itkTypeMacro(AnalyticalDiffusionQballReconstructionImageFilter,
-                                                   ImageToImageFilter);
+    /** Runtime information support. */
+    itkTypeMacro(AnalyticalDiffusionQballReconstructionImageFilter, ImageToImageFilter)
 
-  typedef TReferenceImagePixelType                 ReferencePixelType;
+    typedef TReferenceImagePixelType                 ReferencePixelType;
 
-  typedef TGradientImagePixelType                  GradientPixelType;
+    typedef TGradientImagePixelType                  GradientPixelType;
 
-  typedef Vector< TOdfPixelType, NrOdfDirections > OdfPixelType;
+    typedef Vector< TOdfPixelType, NrOdfDirections > OdfPixelType;
 
-  typedef TOdfPixelType                            BZeroPixelType;
+    typedef TOdfPixelType                            BZeroPixelType;
 
-  /** Reference image data,  This image is aquired in the absence
+    /** Reference image data,  This image is aquired in the absence
    * of a diffusion sensitizing field gradient */
-  typedef typename Superclass::InputImageType       ReferenceImageType;
+    typedef typename Superclass::InputImageType       ReferenceImageType;
 
-  typedef Image< OdfPixelType, 3 >                  OdfImageType;
+    typedef Image< OdfPixelType, 3 >                  OdfImageType;
 
-  typedef OdfImageType                              OutputImageType;
+    typedef OdfImageType                              OutputImageType;
 
-  typedef Image< Vector< TOdfPixelType, (NOrderL*NOrderL + NOrderL + 2)/2 + NOrderL >, 3 > CoefficientImageType;
+    typedef Image< Vector< TOdfPixelType, (NOrderL*NOrderL + NOrderL + 2)/2 + NOrderL >, 3 > CoefficientImageType;
 
-  typedef Image< BZeroPixelType, 3 >                BZeroImageType;
+    typedef Image< BZeroPixelType, 3 >                BZeroImageType;
 
-  typedef typename Superclass::OutputImageRegionType
-                                                   OutputImageRegionType;
+    typedef typename Superclass::OutputImageRegionType
+    OutputImageRegionType;
 
-  /** Typedef defining one (of the many) gradient images.  */
-  typedef Image< GradientPixelType, 3 >            GradientImageType;
+    /** Typedef defining one (of the many) gradient images.  */
+    typedef Image< GradientPixelType, 3 >            GradientImageType;
 
-  /** An alternative typedef defining one (of the many) gradient images.
+    /** An alternative typedef defining one (of the many) gradient images.
    * It will be assumed that the vectorImage has the same dimension as the
    * Reference image and a vector length parameter of \c n (number of
    * gradient directions)*/
-  typedef VectorImage< GradientPixelType, 3 >      GradientImagesType;
+    typedef VectorImage< GradientPixelType, 3 >      GradientImagesType;
 
-  /** Holds the ODF reconstruction matrix */
-  typedef vnl_matrix< TOdfPixelType >*
-                                                   OdfReconstructionMatrixType;
+    /** Holds the ODF reconstruction matrix */
+    typedef vnl_matrix< TOdfPixelType >*
+    OdfReconstructionMatrixType;
 
-  typedef vnl_matrix< double >                     CoefficientMatrixType;
+    typedef vnl_matrix< double >                     CoefficientMatrixType;
 
-  /** Holds each magnetic field gradient used to acquire one DWImage */
-  typedef vnl_vector_fixed< double, 3 >            GradientDirectionType;
+    /** Holds each magnetic field gradient used to acquire one DWImage */
+    typedef vnl_vector_fixed< double, 3 >            GradientDirectionType;
 
-  /** Container to hold gradient directions of the 'n' DW measurements */
-  typedef VectorContainer< unsigned int,
-          GradientDirectionType >                  GradientDirectionContainerType;
+    /** Container to hold gradient directions of the 'n' DW measurements */
+    typedef VectorContainer< unsigned int,
+    GradientDirectionType >                  GradientDirectionContainerType;
 
-  /** set method to add gradient directions and its corresponding
+    /** set method to add gradient directions and its corresponding
    * image. The image here is a VectorImage. The user is expected to pass the
    * gradient directions in a container. The ith element of the container
    * corresponds to the gradient direction of the ith component image the
    * VectorImage.  For the baseline image, a vector of all zeros
    * should be set.*/
-  void SetGradientImage( GradientDirectionContainerType *,
-                                             const GradientImagesType *image);
+    void SetGradientImage( GradientDirectionContainerType *,
+                           const GradientImagesType *image);
 
-  /** Get reference image */
-  virtual ReferenceImageType * GetReferenceImage()
-  { return ( static_cast< ReferenceImageType *>(this->ProcessObject::GetInput(0)) ); }
+    /** Get reference image */
+    virtual ReferenceImageType * GetReferenceImage()
+    { return ( static_cast< ReferenceImageType *>(this->ProcessObject::GetInput(0)) ); }
 
-  /** Return the gradient direction. idx is 0 based */
-  virtual GradientDirectionType GetGradientDirection( unsigned int idx) const
+    /** Return the gradient direction. idx is 0 based */
+    virtual GradientDirectionType GetGradientDirection( unsigned int idx) const
     {
-    if( idx >= m_NumberOfGradientDirections )
-      {
-      itkExceptionMacro( << "Gradient direction " << idx << "does not exist" );
-      }
-    return m_GradientDirectionContainer->ElementAt( idx+1 );
+        if( idx >= m_NumberOfGradientDirections )
+            itkExceptionMacro( << "Gradient direction " << idx << "does not exist" );
+        return m_GradientDirectionContainer->ElementAt( idx+1 );
     }
 
-  static void tofile2(vnl_matrix<double> *A, std::string fname);
-  static double factorial(int number);
-  static void Cart2Sph(double x, double y, double z, double* cart);
-  static double legendre0(int l);
-  static double spherical_harmonic(int m,int l,double theta,double phi, bool complexPart);
-  static double Yj(int m, int k, double theta, double phi);
+    static void tofile2(vnl_matrix<double> *A, std::string fname);
+    static void Cart2Sph(double x, double y, double z, double* cart);
+    static double Yj(int m, int k, double theta, double phi);
+    double Legendre0(int l);
 
-  OdfPixelType Normalize(OdfPixelType odf, typename NumericTraits<ReferencePixelType>::AccumulateType b0 );
+    OdfPixelType Normalize(OdfPixelType odf, typename NumericTraits<ReferencePixelType>::AccumulateType b0 );
+    vnl_vector<TOdfPixelType> PreNormalize( vnl_vector<TOdfPixelType> vec, typename NumericTraits<ReferencePixelType>::AccumulateType b0  );
 
-  vnl_vector<TOdfPixelType> PreNormalize( vnl_vector<TOdfPixelType> vec, typename NumericTraits<ReferencePixelType>::AccumulateType b0  );
-
-  /** Threshold on the reference image data. The output ODF will be a null
+    /** Threshold on the reference image data. The output ODF will be a null
    * pdf for pixels in the reference image that have a value less than this
    * threshold. */
-  itkSetMacro( Threshold, ReferencePixelType );
-  itkGetMacro( Threshold, ReferencePixelType );
+    itkSetMacro( Threshold, ReferencePixelType )
+    itkGetMacro( Threshold, ReferencePixelType )
 
-  itkSetMacro( NormalizationMethod, Normalization);
-  itkGetMacro( NormalizationMethod, Normalization );
+    itkSetMacro( NormalizationMethod, Normalization)
+    itkGetMacro( NormalizationMethod, Normalization )
 
-  typedef Image<float, 3> FloatImageType;
-  itkGetMacro( BZeroImage, typename BZeroImageType::Pointer);
-  itkGetMacro( ODFSumImage, typename FloatImageType::Pointer);
-  itkGetMacro( CoefficientImage, typename CoefficientImageType::Pointer);
+    typedef Image<float, 3> FloatImageType;
+    itkGetMacro( BZeroImage, typename BZeroImageType::Pointer)
+    itkGetMacro( ODFSumImage, typename FloatImageType::Pointer)
+    itkGetMacro( CoefficientImage, typename CoefficientImageType::Pointer)
 
-  itkSetMacro( BValue, TOdfPixelType);
+    itkSetMacro( BValue, TOdfPixelType)
 #ifdef GetBValue
 #undef GetBValue
 #endif
-  itkGetConstReferenceMacro( BValue, TOdfPixelType);
-
-  itkSetMacro( Lambda, double );
-  itkGetMacro( Lambda, double );
+    itkGetConstReferenceMacro( BValue, TOdfPixelType)
+    itkSetMacro( Lambda, double )
+    itkGetMacro( Lambda, double )
 
 #ifdef ITK_USE_CONCEPT_CHECKING
-  /** Begin concept checking */
-  itkConceptMacro(ReferenceEqualityComparableCheck,
-    (Concept::EqualityComparable<ReferencePixelType>));
-  itkConceptMacro(TensorEqualityComparableCheck,
-    (Concept::EqualityComparable<OdfPixelType>));
-  itkConceptMacro(GradientConvertibleToDoubleCheck,
-    (Concept::Convertible<GradientPixelType, double>));
-  itkConceptMacro(DoubleConvertibleToTensorCheck,
-    (Concept::Convertible<double, OdfPixelType>));
-  itkConceptMacro(GradientReferenceAdditiveOperatorsCheck,
-    (Concept::AdditiveOperators<GradientPixelType, GradientPixelType,
-                                ReferencePixelType>));
-  itkConceptMacro(ReferenceOStreamWritableCheck,
-    (Concept::OStreamWritable<ReferencePixelType>));
-  itkConceptMacro(TensorOStreamWritableCheck,
-    (Concept::OStreamWritable<OdfPixelType>));
-  /** End concept checking */
+    /** Begin concept checking */
+    itkConceptMacro(ReferenceEqualityComparableCheck,
+                    (Concept::EqualityComparable<ReferencePixelType>));
+    itkConceptMacro(TensorEqualityComparableCheck,
+                    (Concept::EqualityComparable<OdfPixelType>));
+    itkConceptMacro(GradientConvertibleToDoubleCheck,
+                    (Concept::Convertible<GradientPixelType, double>));
+    itkConceptMacro(DoubleConvertibleToTensorCheck,
+                    (Concept::Convertible<double, OdfPixelType>));
+    itkConceptMacro(GradientReferenceAdditiveOperatorsCheck,
+                    (Concept::AdditiveOperators<GradientPixelType, GradientPixelType,
+                     ReferencePixelType>));
+    itkConceptMacro(ReferenceOStreamWritableCheck,
+                    (Concept::OStreamWritable<ReferencePixelType>));
+    itkConceptMacro(TensorOStreamWritableCheck,
+                    (Concept::OStreamWritable<OdfPixelType>));
+    /** End concept checking */
 #endif
 
 protected:
-  AnalyticalDiffusionQballReconstructionImageFilter();
-  ~AnalyticalDiffusionQballReconstructionImageFilter() {};
-  void PrintSelf(std::ostream& os, Indent indent) const;
+    AnalyticalDiffusionQballReconstructionImageFilter();
+    ~AnalyticalDiffusionQballReconstructionImageFilter() {};
+    void PrintSelf(std::ostream& os, Indent indent) const;
 
-  void ComputeReconstructionMatrix();
-
-  void BeforeThreadedGenerateData();
-  void ThreadedGenerateData( const
-      OutputImageRegionType &outputRegionForThread, int);
+    void ComputeReconstructionMatrix();
+    void BeforeThreadedGenerateData();
+    void ThreadedGenerateData( const
+                               OutputImageRegionType &outputRegionForThread, int);
 
 private:
 
-  OdfReconstructionMatrixType                       m_ReconstructionMatrix;
-
-  OdfReconstructionMatrixType                       m_CoeffReconstructionMatrix;
-
-  OdfReconstructionMatrixType                       m_SphericalHarmonicBasisMatrix;
-
-  /** container to hold gradient directions */
-  GradientDirectionContainerType::Pointer           m_GradientDirectionContainer;
-
-  /** Number of gradient measurements */
-  unsigned int                                      m_NumberOfGradientDirections;
-
-  /** Number of baseline images */
-  unsigned int                                      m_NumberOfBaselineImages;
-
-  /** Threshold on the reference image data */
-  ReferencePixelType                                m_Threshold;
-
-  /** LeBihan's b-value for normalizing tensors */
-  TOdfPixelType                                     m_BValue;
-
-  typename BZeroImageType::Pointer                  m_BZeroImage;
-
-  double                                            m_Lambda;
-
-  bool                                              m_DirectionsDuplicated;
-
-  Normalization                                     m_NormalizationMethod;
-
-  int                                               m_NumberCoefficients;
-
-  vnl_matrix<double>*                               m_B_t;
-  vnl_vector<double>*                               m_LP;
-  FloatImageType::Pointer m_ODFSumImage;
-  typename CoefficientImageType::Pointer            m_CoefficientImage;
-  TOdfPixelType                            m_Delta1;
-  TOdfPixelType                            m_Delta2;
+    OdfReconstructionMatrixType                       m_ReconstructionMatrix;
+    OdfReconstructionMatrixType                       m_CoeffReconstructionMatrix;
+    OdfReconstructionMatrixType                       m_SphericalHarmonicBasisMatrix;
+    /** container to hold gradient directions */
+    GradientDirectionContainerType::Pointer           m_GradientDirectionContainer;
+    /** Number of gradient measurements */
+    unsigned int                                      m_NumberOfGradientDirections;
+    /** Number of baseline images */
+    unsigned int                                      m_NumberOfBaselineImages;
+    /** Threshold on the reference image data */
+    ReferencePixelType                                m_Threshold;
+    /** LeBihan's b-value for normalizing tensors */
+    TOdfPixelType                                     m_BValue;
+    typename BZeroImageType::Pointer                  m_BZeroImage;
+    double                                            m_Lambda;
+    bool                                              m_DirectionsDuplicated;
+    Normalization                                     m_NormalizationMethod;
+    int                                               m_NumberCoefficients;
+    vnl_matrix<double>*                               m_B_t;
+    vnl_vector<double>*                               m_LP;
+    FloatImageType::Pointer                           m_ODFSumImage;
+    typename CoefficientImageType::Pointer            m_CoefficientImage;
+    TOdfPixelType                                     m_Delta1;
+    TOdfPixelType                                     m_Delta2;
 };
 
 }
 
 #ifndef ITK_MANUAL_INSTANTIATION
 #include "itkAnalyticalDiffusionQballReconstructionImageFilter.cpp"
 #endif
 
 #endif //__itkAnalyticalDiffusionQballReconstructionImageFilter_h_