diff --git a/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp b/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp
index 5b6e1922d8..a61ecaf99d 100644
--- a/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp
+++ b/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp
@@ -1,882 +1,749 @@
 /*=========================================================================
 
 Program:   Medical Imaging & Interaction Toolkit
 Language:  C++
 Date:      $Date: 2009-07-14 19:11:20 +0200 (Tue, 14 Jul 2009) $
 Version:   $Revision: 18127 $
 
 Copyright (c) German Cancer Research Center, Division of Medical and
 Biological Informatics. All rights reserved.
 See MITKCopyright.txt or http://www.mitk.org/copyright.html for details.
 
 This software is distributed WITHOUT ANY WARRANTY; without even
 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
 PURPOSE.  See the above copyright notices for more information.
 
 =========================================================================*/
 #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 "itkPointShell.h"
 
+#include "mitkSphericalHarmonicsFunctions.h"
+
 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)
   {
     // 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] = odf[i] < 0 ? 0 : odf[i];
             odf[i] *= QBALL_ANAL_RECON_PI*4/NrOdfDirections;
           }
           TOdfPixelType sum = 0;
           for(int i=0; i<NrOdfDirections; i++)
           {
             sum += odf[i];
           }
           if(sum>0)
             odf /= sum;
 
           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)
               {
                 MITK_ERROR << "AnalyticalDiffusionQballReconstructionImageFilter: negative dwi image value";
                 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)
               {
                 MITK_ERROR << "AnalyticalDiffusionQballReconstructionImageFilter: negative dwi image value";
                 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++)
             {
               if (vec[i]<=0)
               {
                 MITK_ERROR << "AnalyticalDiffusionQballReconstructionImageFilter: negative dwi image value";
                 vec[i] = 0.001;
               }
               if (b0f==0)
                 b0f = 0.01;
               if(vec[i] >= b0f)
                 vec[i] = b0f - 0.001;
               vec[i] = log(-log(vec[i]/b0f));
             }
             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 < 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 );
         }
 
         this->ComputeReconstructionMatrix();
 
         m_BZeroImage = BZeroImageType::New();
         typename GradientImagesType::Pointer img = static_cast< GradientImagesType * >(
           this->ProcessObject::GetInput(0) );
         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();
 
         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,
         int )
       {
         typename OutputImageType::Pointer outputImage =
           static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0));
 
         ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread);
         oit.GoToBegin();
 
         ImageRegionIterator< BZeroImageType > oit2(m_BZeroImage, outputRegionForThread);
         oit2.GoToBegin();
 
         ImageRegionIterator< BlaImage > oit3(m_ODFSumImage, outputRegionForThread);
         oit3.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);
           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();
             }
             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
             {
               odf = ( (*m_ReconstructionMatrix) * B ).data_block();
             }
             odf = Normalize(odf, b0);
 
           }
 
           oit.Set( odf );
           ++oit;
           oit2.Set( b0 );
           float sum = 0;
           for (int k=0; k<odf.Size(); k++)
             sum += (float) odf[k];
           oit3.Set( sum-1 );
           ++oit3;
           ++oit2;
           ++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(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;
-        }
-
-        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
-        }
-
-        //std::cout << "legendre(" << l << "," << m << "," << cos(theta) << ") = " << pml << std::endl;
-
-        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( 0 < m && m <= k )
-        {
-          return sqrt(2.0) * spherical_harmonic(m,k,theta,phi,true);
-        }
-
-        return 0;
-      }
-
-
-
       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 < 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();
             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);
+              mitk::ShericalHarmonicsFunctions::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);
+                mitk::ShericalHarmonicsFunctions::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);
+              (*B)(i,j) = mitk::ShericalHarmonicsFunctions::Yj(m,k,th,phi);
               //std::cout << "B(" << i << "," << j << ") = Yj(" << m << "," << k << "," << th << "," << phi << ") = " << (*B)(i,j) << std::endl;
             }
           }
         }
 
         //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]);
+            (*P)(i,i) = 2.0*QBALL_ANAL_RECON_PI*mitk::ShericalHarmonicsFunctions::legendre0((*lj)[i]);
             (*m_LP)(i) *= (*P)(i,i);
           }
           else
           {
-            (*P)(i,i) = legendre0((*lj)[i]);
+            (*P)(i,i) = mitk::ShericalHarmonicsFunctions::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);
+          mitk::ShericalHarmonicsFunctions::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_SphericalHarmonicBasisMatrix)(i,j) = mitk::ShericalHarmonicsFunctions::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);
 
       }
 
       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)
           {
             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/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.h b/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.h
index 1c89d2c9a2..145c19c8ce 100644
--- a/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.h
+++ b/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.h
@@ -1,300 +1,295 @@
 /*=========================================================================
 
 Program:   Medical Imaging & Interaction Toolkit
 Language:  C++
 Date:      $Date: 2009-07-14 19:11:20 +0200 (Tue, 14 Jul 2009) $
 Version:   $Revision: 18127 $
 
 Copyright (c) German Cancer Research Center, Division of Medical and
 Biological Informatics. All rights reserved.
 See MITKCopyright.txt or http://www.mitk.org/copyright.html for details.
 
 This software is distributed WITHOUT ANY WARRANTY; without even
 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
 PURPOSE.  See the above copyright notices for more information.
 
 =========================================================================*/
 
 #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:
 
   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;
   
    /** Method for creation through the object factory. */
   itkNewMacro(Self);  
 
   /** Runtime information support. */
   itkTypeMacro(AnalyticalDiffusionQballReconstructionImageFilter, 
                                                    ImageToImageFilter);
  
   typedef TReferenceImagePixelType                 ReferencePixelType;
 
   typedef TGradientImagePixelType                  GradientPixelType;
 
   typedef Vector< TOdfPixelType, NrOdfDirections >
                                                    OdfPixelType;
 
   /** Reference image data,  This image is aquired in the absence 
    * of a diffusion sensitizing field gradient */
   typedef typename Superclass::InputImageType      ReferenceImageType;
   
   typedef Image< OdfPixelType, 3 >                 OdfImageType;
 
   typedef OdfImageType                             OutputImageType;
 
   typedef TOdfPixelType                            BZeroPixelType;
 
   typedef Image< BZeroPixelType, 3 >      BZeroImageType;
 
   typedef typename Superclass::OutputImageRegionType
                                                    OutputImageRegionType;
 
   /** Typedef defining one (of the many) gradient images.  */
   typedef Image< GradientPixelType, 3 >            GradientImageType;
 
   /** 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;
 
   /** Holds the ODF reconstruction matrix */
   typedef vnl_matrix< TOdfPixelType >*
                                                    OdfReconstructionMatrixType;
   
   typedef vnl_matrix< double >                     CoefficientMatrixType;
 
   /** 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;
   
   /** 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);
   
   /** 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
     {
     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);
 
   OdfPixelType Normalize(OdfPixelType odf, 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
    * pdf for pixels in the reference image that have a value less than this
    * threshold. */
   itkSetMacro( Threshold, ReferencePixelType );
   itkGetMacro( Threshold, ReferencePixelType );
 
   itkSetMacro( NormalizationMethod, Normalization);
   itkGetMacro( NormalizationMethod, Normalization );
 
   typedef Image<float, 3> BlaImage;
   itkGetMacro( BZeroImage, typename BZeroImageType::Pointer);
   itkGetMacro( ODFSumImage, typename BlaImage::Pointer);
 
   itkSetMacro( BValue, TOdfPixelType);
 #ifdef GetBValue
 #undef GetBValue
 #endif
   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 */
 #endif
 
 protected:
   AnalyticalDiffusionQballReconstructionImageFilter();
   ~AnalyticalDiffusionQballReconstructionImageFilter() {};
   void PrintSelf(std::ostream& os, Indent indent) const;
 
   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;
   BlaImage::Pointer m_ODFSumImage;
 };
 
 }
 
 #ifndef ITK_MANUAL_INSTANTIATION
 #include "itkAnalyticalDiffusionQballReconstructionImageFilter.cpp"
 #endif
 
 #endif //__itkAnalyticalDiffusionQballReconstructionImageFilter_h_
 
diff --git a/Modules/DiffusionImaging/Reconstruction/mitkSphericalHarmonicsFunctions.cpp b/Modules/DiffusionImaging/Reconstruction/mitkSphericalHarmonicsFunctions.cpp
new file mode 100644
index 0000000000..6317215e60
--- /dev/null
+++ b/Modules/DiffusionImaging/Reconstruction/mitkSphericalHarmonicsFunctions.cpp
@@ -0,0 +1,145 @@
+/*=========================================================================
+
+Program:   Medical Imaging & Interaction Toolkit
+Language:  C++
+Date:      $Date: 2009-07-14 19:11:20 +0200 (Tue, 14 Jul 2009) $
+Version:   $Revision: 18127 $
+
+Copyright (c) German Cancer Research Center, Division of Medical and
+Biological Informatics. All rights reserved.
+See MITKCopyright.txt or http://www.mitk.org/copyright.html for details.
+
+This software is distributed WITHOUT ANY WARRANTY; without even
+the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+PURPOSE.  See the above copyright notices for more information.
+
+=========================================================================*/
+
+#include <boost/version.hpp>
+
+#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 "mitkSphericalHarmonicsFunctions.h"
+namespace mitk{
+
+#define SPHERICAL_HARMONICS_PI       M_PI
+
+double mitk::ShericalHarmonicsFunctions::factorial(int number) {
+    if(number <= 1) return 1;
+    double result = 1.0;
+    for(int i=1; i<=number; i++)
+        result *= i;
+    return result;
+}
+
+void mitk::ShericalHarmonicsFunctions::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+SPHERICAL_HARMONICS_PI/2;
+    phi = -phi+SPHERICAL_HARMONICS_PI;
+    cart[0] = phi;
+    cart[1] = th;
+    cart[2] = rad;
+}
+
+double mitk::ShericalHarmonicsFunctions::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);
+    }
+}
+
+double mitk::ShericalHarmonicsFunctions::spherical_harmonic(int m,int l,double theta,double phi, bool complexPart)
+{
+    if( theta < 0 || theta > SPHERICAL_HARMONICS_PI )
+    {
+        std::cout << "bad range" << std::endl;
+        return 0;
+    }
+
+    if( phi < 0 || phi > 2*SPHERICAL_HARMONICS_PI )
+    {
+        std::cout << "bad range" << std::endl;
+        return 0;
+    }
+
+    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
+    }
+
+    //std::cout << "legendre(" << l << "," << m << "," << cos(theta) << ") = " << pml << std::endl;
+
+    double retval = sqrt(((2.0*(double)l+1.0)/(4.0*SPHERICAL_HARMONICS_PI))*(fac2/fac1)) * pml;
+    if( !complexPart )
+    {
+        retval *= cos(m*phi);
+    }
+    else
+    {
+        retval *= sin(m*phi);
+    }
+    //std::cout << retval << std::endl;
+    return retval;
+}
+
+double mitk::ShericalHarmonicsFunctions::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( 0 < m && m <= k )
+    {
+        return sqrt(2.0) * spherical_harmonic(m,k,theta,phi,true);
+    }
+
+    return 0;
+}
+}
diff --git a/Modules/DiffusionImaging/Reconstruction/mitkSphericalHarmonicsFunctions.h b/Modules/DiffusionImaging/Reconstruction/mitkSphericalHarmonicsFunctions.h
new file mode 100644
index 0000000000..9abc121e44
--- /dev/null
+++ b/Modules/DiffusionImaging/Reconstruction/mitkSphericalHarmonicsFunctions.h
@@ -0,0 +1,38 @@
+/*=========================================================================
+
+Program:   Medical Imaging & Interaction Toolkit
+Language:  C++
+Date:      $Date: 2009-07-14 19:11:20 +0200 (Tue, 14 Jul 2009) $
+Version:   $Revision: 18127 $
+
+Copyright (c) German Cancer Research Center, Division of Medical and
+Biological Informatics. All rights reserved.
+See MITKCopyright.txt or http://www.mitk.org/copyright.html for details.
+
+This software is distributed WITHOUT ANY WARRANTY; without even
+the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+PURPOSE.  See the above copyright notices for more information.
+
+=========================================================================*/
+
+#ifndef __mitkSphericalHarmonicsHandler_h_
+#define __mitkSphericalHarmonicsHandler_h_
+
+#include "MitkDiffusionImagingExports.h"
+
+namespace mitk{
+
+class MitkDiffusionImaging_EXPORT ShericalHarmonicsFunctions{
+
+public :
+   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);
+
+};
+}
+
+#endif //__mitkSphericalHarmonicsHandler_h_
+
diff --git a/Modules/DiffusionImaging/files.cmake b/Modules/DiffusionImaging/files.cmake
index e22043c204..c7d1299abe 100644
--- a/Modules/DiffusionImaging/files.cmake
+++ b/Modules/DiffusionImaging/files.cmake
@@ -1,167 +1,172 @@
 SET(CPP_FILES
 
   # DicomImport
   DicomImport/mitkDicomDiffusionImageReader.cpp
   DicomImport/mitkGroupDiffusionHeadersFilter.cpp
   DicomImport/mitkDicomDiffusionImageHeaderReader.cpp
   DicomImport/mitkGEDicomDiffusionImageHeaderReader.cpp
   DicomImport/mitkPhilipsDicomDiffusionImageHeaderReader.cpp
   DicomImport/mitkSiemensDicomDiffusionImageHeaderReader.cpp
   DicomImport/mitkSiemensMosaicDicomDiffusionImageHeaderReader.cpp
 
   # DataStructures
   IODataStructures/mitkDiffusionImagingObjectFactory.cpp
 
   # DataStructures -> DWI
   IODataStructures/DiffusionWeightedImages/mitkDiffusionImageHeaderInformation.cpp
   IODataStructures/DiffusionWeightedImages/mitkDiffusionImageSource.cpp
   IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageReader.cpp
   IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageWriter.cpp
   IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageIOFactory.cpp
   IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageWriterFactory.cpp
   IODataStructures/DiffusionWeightedImages/mitkDiffusionImageSerializer.cpp
 
   # DataStructures -> QBall
   IODataStructures/QBallImages/mitkQBallImageSource.cpp
   IODataStructures/QBallImages/mitkNrrdQBallImageReader.cpp
   IODataStructures/QBallImages/mitkNrrdQBallImageWriter.cpp
   IODataStructures/QBallImages/mitkNrrdQBallImageIOFactory.cpp
   IODataStructures/QBallImages/mitkNrrdQBallImageWriterFactory.cpp
   IODataStructures/QBallImages/mitkQBallImage.cpp
   IODataStructures/QBallImages/mitkQBallImageSerializer.cpp
 
   # DataStructures -> Tensor
   IODataStructures/TensorImages/mitkTensorImageSource.cpp
   IODataStructures/TensorImages/mitkNrrdTensorImageReader.cpp
   IODataStructures/TensorImages/mitkNrrdTensorImageWriter.cpp
   IODataStructures/TensorImages/mitkNrrdTensorImageIOFactory.cpp
   IODataStructures/TensorImages/mitkNrrdTensorImageWriterFactory.cpp
   IODataStructures/TensorImages/mitkTensorImage.cpp
   IODataStructures/TensorImages/mitkTensorImageSerializer.cpp
 
   IODataStructures/mitkParticle.cpp
   IODataStructures/mitkParticleGrid.cpp
 
 # DataStructures -> FiberBundleX
   IODataStructures/FiberBundleX/mitkFiberBundleX.cpp
   IODataStructures/FiberBundleX/mitkFiberBundleXWriter.cpp
   IODataStructures/FiberBundleX/mitkFiberBundleXReader.cpp
   IODataStructures/FiberBundleX/mitkFiberBundleXIOFactory.cpp
   IODataStructures/FiberBundleX/mitkFiberBundleXWriterFactory.cpp
   IODataStructures/FiberBundleX/mitkFiberBundleXSerializer.cpp
   IODataStructures/FiberBundleX/mitkFiberBundleXThreadMonitor.cpp
 
   # DataStructures -> PlanarFigureComposite
   IODataStructures/PlanarFigureComposite/mitkPlanarFigureComposite.cpp
 
   # DataStructures -> Tbss
   IODataStructures/TbssImages/mitkTbssImageSource.cpp
   IODataStructures/TbssImages/mitkTbssRoiImageSource.cpp
   IODataStructures/TbssImages/mitkNrrdTbssImageReader.cpp
   IODataStructures/TbssImages/mitkNrrdTbssImageIOFactory.cpp
   IODataStructures/TbssImages/mitkNrrdTbssRoiImageReader.cpp
   IODataStructures/TbssImages/mitkNrrdTbssRoiImageIOFactory.cpp
   IODataStructures/TbssImages/mitkTbssImage.cpp
   IODataStructures/TbssImages/mitkTbssRoiImage.cpp
   IODataStructures/TbssImages/mitkNrrdTbssImageWriter.cpp
   IODataStructures/TbssImages/mitkNrrdTbssImageWriterFactory.cpp
   IODataStructures/TbssImages/mitkNrrdTbssRoiImageWriter.cpp
   IODataStructures/TbssImages/mitkNrrdTbssRoiImageWriterFactory.cpp
   IODataStructures/TbssImages/mitkTbssImporter.cpp
 
   # Rendering
   Rendering/vtkMaskedProgrammableGlyphFilter.cpp
   Rendering/mitkCompositeMapper.cpp
   Rendering/mitkVectorImageVtkGlyphMapper3D.cpp
   Rendering/vtkOdfSource.cxx
   Rendering/vtkThickPlane.cxx
   Rendering/mitkOdfNormalizationMethodProperty.cpp
   Rendering/mitkOdfScaleByProperty.cpp
   Rendering/mitkFiberBundleXMapper2D.cpp
   Rendering/mitkFiberBundleXMapper3D.cpp
   Rendering/mitkFiberBundleXThreadMonitorMapper3D.cpp
   Rendering/mitkTbssImageMapper.cpp
   Rendering/mitkPlanarCircleMapper3D.cpp
   Rendering/mitkPlanarPolygonMapper3D.cpp
 
 # Interactions
   Interactions/mitkFiberBundleInteractor.cpp
 
   # Algorithms
   Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.cpp
   Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.cpp
   Algorithms/mitkTractAnalyzer.cpp
 
  # Tractography
  Tractography/itkStochasticTractographyFilter.h
+
+  # Reconstruction
+  Reconstruction/mitkSphericalHarmonicsFunctions.cpp
+
 )
 
 SET(H_FILES
   # Rendering
   Rendering/mitkDiffusionImageMapper.h
   Rendering/mitkTbssImageMapper.h
   Rendering/mitkOdfVtkMapper2D.h
   Rendering/mitkFiberBundleXMapper3D.h
   Rendering/mitkFiberBundleXMapper2D.h
   Rendering/mitkFiberBundleXThreadMonitorMapper3D.h
   Rendering/mitkPlanarCircleMapper3D.h
   Rendering/mitkPlanarPolygonMapper3D.h
 
   # Reconstruction
   Reconstruction/itkDiffusionQballReconstructionImageFilter.h
   Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.h
   Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.h
+  Reconstruction/mitkSphericalHarmonicsFunctions.h
   Reconstruction/itkPointShell.h
   Reconstruction/itkOrientationDistributionFunction.h
   Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h
   Reconstruction/itkRegularizedIVIMLocalVariationImageFilter.h
   Reconstruction/itkRegularizedIVIMReconstructionFilter.h
   Reconstruction/itkRegularizedIVIMReconstructionSingleIteration.h
 
   # IO Datastructures
   IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h
   IODataStructures/TbssImages/mitkTbssImporter.h
 
   # DataStructures -> FiberBundleX
   IODataStructures/FiberBundleX/mitkFiberBundleX.h
   IODataStructures/FiberBundleX/mitkFiberBundleXWriter.h
   IODataStructures/FiberBundleX/mitkFiberBundleXReader.h
   IODataStructures/FiberBundleX/mitkFiberBundleXIOFactory.h
   IODataStructures/FiberBundleX/mitkFiberBundleXWriterFactory.h
   IODataStructures/FiberBundleX/mitkFiberBundleXSerializer.h
   IODataStructures/FiberBundleX/mitkFiberBundleXThreadMonitor.h
 
   # Tractography
   Tractography/itkGibbsTrackingFilter.h
   Tractography/itkStochasticTractographyFilter.h
 
   # Algorithms
   Algorithms/itkDiffusionQballGeneralizedFaImageFilter.h
   Algorithms/itkDiffusionQballPrepareVisualizationImageFilter.h
   Algorithms/itkTensorDerivedMeasurementsFilter.h
   Algorithms/itkBrainMaskExtractionImageFilter.h
   Algorithms/itkB0ImageExtractionImageFilter.h
   Algorithms/itkTensorImageToDiffusionImageFilter.h
   Algorithms/itkTensorToL2NormImageFilter.h
   Algorithms/itkTractDensityImageFilter.h
   Algorithms/itkTractsToFiberEndingsImageFilter.h
   Algorithms/itkTractsToRgbaImageFilter.h
   Algorithms/itkGaussianInterpolateImageFunction.h
   Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.h
   Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.h
   Algorithms/itkDiffusionTensorPrincipleDirectionImageFilter.h
   Algorithms/itkCartesianToPolarVectorImageFilter.h
   Algorithms/itkPolarToCartesianVectorImageFilter.h
   Algorithms/itkDistanceMapFilter.h
   Algorithms/itkProjectionFilter.h
   Algorithms/itkSkeletonizationFilter.h
   Algorithms/itkReduceDirectionGradientsFilter.h
 )
 
 SET( TOOL_FILES
 )
 
 IF(WIN32)
 ENDIF(WIN32)
 
 #MITK_MULTIPLEX_PICTYPE( Algorithms/mitkImageRegistrationMethod-TYPE.cpp )