diff --git a/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp b/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp
index 59a4cc84ea..4142e043ca 100644
--- a/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp
+++ b/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp
@@ -1,753 +1,753 @@
 /*=========================================================================
 
 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 <stdio.h>
 
 #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),
     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;
 
           for(int i=0; i<NrOdfDirections; i++)
           {
             sum += odf[i];
           }
           if(sum>0)
             odf /= sum;
 
           return odf;
           break;
         }
       case QBAR_B_ZERO_B_VALUE:
         {
           for(int i=0; i<NrOdfDirections; i++)
           {
             odf[i] = ((TOdfPixelType)log((TOdfPixelType)b0)-odf[i])/m_BValue;
           }
           return odf;
           break;
         }
       case QBAR_B_ZERO:
         {
           odf *= 1.0/b0;
           return odf;
           break;
         }
       case QBAR_NONE:
         {
           return odf;
           break;
         }
       case QBAR_ADC_ONLY:
         {
           for(int i=0; i<NrOdfDirections; i++)
           {
             odf[i] = ((TOdfPixelType)log((TOdfPixelType)b0)-odf[i])/m_BValue;
           }
           return odf;
           break;
         }
       case QBAR_RAW_SIGNAL:
         {
           return odf;
           break;
         }
       case QBAR_SOLID_ANGLE:
         {
           for(int i=0; i<NrOdfDirections; i++)
             odf[i] *= QBALL_ANAL_RECON_PI*4/NrOdfDirections;
 
           break;
         }
       case QBAR_NONNEG_SOLID_ANGLE:
         {
           break;
         }
       }
 
       return odf;
     }
 
 
     template<
       class TReferenceImagePixelType,
       class TGradientImagePixelType,
       class TOdfPixelType,
         int NOrderL,
         int NrOdfDirections>
         vnl_vector<TOdfPixelType>
         itk::AnalyticalDiffusionQballReconstructionImageFilter
         <TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType,
         NOrderL, NrOdfDirections>
         ::PreNormalize( vnl_vector<TOdfPixelType> vec,
         typename NumericTraits<ReferencePixelType>::AccumulateType b0 )
       {
         switch( m_NormalizationMethod )
         {
         case QBAR_STANDARD:
           {
             return vec;
             break;
           }
         case QBAR_B_ZERO_B_VALUE:
           {
             int n = vec.size();
             for(int i=0; i<n; i++)
             {
               if (vec[i]<=0)
                 vec[i] = 0.001;
 
               vec[i] = log(vec[i]);
             }
             return vec;
             break;
           }
         case QBAR_B_ZERO:
           {
             return vec;
             break;
           }
         case QBAR_NONE:
           {
             return vec;
             break;
           }
         case QBAR_ADC_ONLY:
           {
             int n = vec.size();
             for(int i=0; i<n; i++)
             {
               if (vec[i]<=0)
                 vec[i] = 0.001;
 
               vec[i] = log(vec[i]);
             }
             return vec;
             break;
           }
         case QBAR_RAW_SIGNAL:
           {
             return vec;
             break;
           }
         case QBAR_SOLID_ANGLE:
         case QBAR_NONNEG_SOLID_ANGLE:
           {
             int n = vec.size();
             double b0f = (double)b0;
             for(int i=0; i<n; i++)
             {
               vec[i] = vec[i]/b0f;
 
               if (vec[i]<0)
                 vec[i] = m_Delta1;
               else if (vec[i]<m_Delta1)
                 vec[i] = m_Delta1/2 + vec[i]*vec[i]/(2*m_Delta1);
               else if (vec[i]>=1)
                 vec[i] = 1-m_Delta2/2;
               else if (vec[i]>=1-m_Delta2)
                 vec[i] = 1-m_Delta2/2-(1-vec[i])*(1-vec[i])/(2*m_Delta2);
 
               vec[i] = log(-log(vec[i]));
             }
             return vec;
             break;
           }
         }
 
         return vec;
       }
 
       template< class T, class TG, class TO, int L, int NODF>
       void AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
         ::BeforeThreadedGenerateData()
       {
         // If we have more than 2 inputs, then each input, except the first is a
         // gradient image. The number of gradient images must match the number of
         // gradient directions.
         //const unsigned int numberOfInputs = this->GetNumberOfInputs();
 
         // There need to be at least 6 gradient directions to be able to compute the
         // tensor basis
         if( m_NumberOfGradientDirections < 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();
 
         typename GradientImagesType::Pointer img = static_cast< GradientImagesType * >(
           this->ProcessObject::GetInput(0) );
 
         m_BZeroImage = BZeroImageType::New();
         m_BZeroImage->SetSpacing( img->GetSpacing() );   // Set the image spacing
         m_BZeroImage->SetOrigin( img->GetOrigin() );     // Set the image origin
         m_BZeroImage->SetDirection( img->GetDirection() );  // Set the image direction
         m_BZeroImage->SetLargestPossibleRegion( img->GetLargestPossibleRegion());
         m_BZeroImage->SetBufferedRegion( img->GetLargestPossibleRegion() );
         m_BZeroImage->Allocate();
 
         m_ODFSumImage = BZeroImageType::New();
         m_ODFSumImage->SetSpacing( img->GetSpacing() );   // Set the image spacing
         m_ODFSumImage->SetOrigin( img->GetOrigin() );     // Set the image origin
         m_ODFSumImage->SetDirection( img->GetDirection() );  // Set the image direction
         m_ODFSumImage->SetLargestPossibleRegion( img->GetLargestPossibleRegion());
         m_ODFSumImage->SetBufferedRegion( img->GetLargestPossibleRegion() );
         m_ODFSumImage->Allocate();
 
         m_CoefficientImage = CoefficientImageType::New();
         m_CoefficientImage->SetSpacing( img->GetSpacing() );   // Set the image spacing
         m_CoefficientImage->SetOrigin( img->GetOrigin() );     // Set the image origin
         m_CoefficientImage->SetDirection( img->GetDirection() );  // Set the image direction
         m_CoefficientImage->SetLargestPossibleRegion( img->GetLargestPossibleRegion());
         m_CoefficientImage->SetBufferedRegion( img->GetLargestPossibleRegion() );
         m_CoefficientImage->Allocate();
 
         if(m_NormalizationMethod == QBAR_SOLID_ANGLE ||
           m_NormalizationMethod == QBAR_NONNEG_SOLID_ANGLE)
         {
           m_Lambda = 0.0;
         }
 
       }
 
       template< class T, class TG, class TO, int L, int NODF>
       void AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
         ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
         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< FloatImageType > oit3(m_ODFSumImage, outputRegionForThread);
         oit3.GoToBegin();
 
         ImageRegionIterator< CoefficientImageType > oit4(m_CoefficientImage, outputRegionForThread);
         oit4.GoToBegin();
 
         typedef ImageRegionConstIterator< GradientImagesType > GradientIteratorType;
         typedef typename GradientImagesType::PixelType         GradientVectorType;
         typename GradientImagesType::Pointer gradientImagePointer = NULL;
 
         // Would have liked a dynamic_cast here, but seems SGI doesn't like it
         // The enum will ensure that an inappropriate cast is not done
         gradientImagePointer = static_cast< GradientImagesType * >(
           this->ProcessObject::GetInput(0) );
 
         GradientIteratorType git(gradientImagePointer, outputRegionForThread );
         git.GoToBegin();
 
         // Compute the indicies of the baseline images and gradient images
         std::vector<unsigned int> baselineind; // contains the indicies of
         // the baseline images
         std::vector<unsigned int> gradientind; // contains the indicies of
         // the gradient images
 
         for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin();
           gdcit != this->m_GradientDirectionContainer->End(); ++gdcit)
         {
           if(gdcit.Value().one_norm() <= 0.0)
           {
             baselineind.push_back(gdcit.Index());
           }
           else
           {
             gradientind.push_back(gdcit.Index());
           }
         }
 
         if( m_DirectionsDuplicated )
         {
           int gradIndSize = gradientind.size();
           for(int i=0; i<gradIndSize; i++)
             gradientind.push_back(gradientind[i]);
         }
 
         while( !git.IsAtEnd() )
         {
           GradientVectorType b = git.Get();
 
           typename NumericTraits<ReferencePixelType>::AccumulateType b0 = NumericTraits<ReferencePixelType>::Zero;
 
           // Average the baseline image pixels
           for(unsigned int i = 0; i < baselineind.size(); ++i)
           {
             b0 += b[baselineind[i]];
           }
           b0 /= this->m_NumberOfBaselineImages;
 
           OdfPixelType odf(0.0);
           typename CoefficientImageType::PixelType coeffPixel(0.0);
           vnl_vector<TO> B(m_NumberOfGradientDirections);
 
           if( (b0 != 0) && (b0 >= m_Threshold) )
           {
             for( unsigned int i = 0; i< m_NumberOfGradientDirections; i++ )
             {
               B[i] = static_cast<TO>(b[gradientind[i]]);
             }
 
             B = PreNormalize(B, b0);
             if(m_NormalizationMethod == QBAR_SOLID_ANGLE)
             {
               vnl_vector<TO> coeffs(m_NumberCoefficients);
               coeffs = ( (*m_CoeffReconstructionMatrix) * B );
               coeffs[0] += 1.0/(2.0*sqrt(QBALL_ANAL_RECON_PI));
               odf = ( (*m_SphericalHarmonicBasisMatrix) * coeffs ).data_block();
               coeffPixel = coeffs.data_block();
             }
             else if(m_NormalizationMethod == QBAR_NONNEG_SOLID_ANGLE)
             {
               /** this would be the place to implement a non-negative
               * solver for quadratic programming problem:
               * min .5*|| Bc-s ||^2 subject to -CLPc <= 4*pi*ones
               * (refer to MICCAI 2009 Goh et al. "Estimating ODFs with PDF constraints")
               * .5*|| Bc-s ||^2 == .5*c'B'Bc - x'B's + .5*s's
               */
 
               itkExceptionMacro( << "Nonnegative Solid Angle not yet implemented");
 
             }
             else
             {
               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++)
             sum += (float) odf[k];
           oit3.Set( sum-1 );
           oit4.Set(coeffPixel);
           ++oit;  // odf image iterator
           ++oit3; // odf sum image iterator
           ++oit2; // b0 image iterator
           ++oit4; // coefficient image iterator
           ++git;  // Gradient  image iterator
         }
 
         std::cout << "One Thread finished reconstruction" << std::endl;
       }
 
       template< class T, class TG, class TO, int L, int NODF>
       void AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
         ::tofile2(vnl_matrix<double> *pA, std::string fname)
       {
         vnl_matrix<double> A = (*pA);
         ofstream myfile;
         std::locale C("C");
         std::locale originalLocale = myfile.getloc();
         myfile.imbue(C);
 
         myfile.open (fname.c_str());
         myfile << "A1=[";
         for(int i=0; i<A.rows(); i++)
         {
           for(int j=0; j<A.columns(); j++)
           {
             myfile << A(i,j) << " ";
             if(j==A.columns()-1 && i!=A.rows()-1)
               myfile << ";";
           }
         }
         myfile << "];";
         myfile.close();
 
         myfile.imbue( originalLocale );
       }
 
       template< class T, class TG, class TO, int L, int NODF>
       void AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
         ::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];
-              mitk::ShericalHarmonicsFunctions::Cart2Sph(x,y,z,cart);
+              mitk::mitk_sh_functions::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];
-                mitk::ShericalHarmonicsFunctions::Cart2Sph(x,y,z,cart);
+                mitk::mitk_sh_functions::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) = mitk::ShericalHarmonicsFunctions::Yj(m,k,th,phi);
+              (*B)(i,j) = mitk::mitk_sh_functions::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*mitk::ShericalHarmonicsFunctions::legendre0((*lj)[i]);
+            (*P)(i,i) = 2.0*QBALL_ANAL_RECON_PI*mitk::mitk_sh_functions::legendre0((*lj)[i]);
             (*m_LP)(i) *= (*P)(i,i);
           }
           else
           {
-            (*P)(i,i) = mitk::ShericalHarmonicsFunctions::legendre0((*lj)[i]);
+            (*P)(i,i) = mitk::mitk_sh_functions::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];
-          mitk::ShericalHarmonicsFunctions::Cart2Sph(x,y,z,cart);
+          mitk::mitk_sh_functions::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) = mitk::ShericalHarmonicsFunctions::Yj(m,k,th,phi);
+              (*m_SphericalHarmonicBasisMatrix)(i,j) = mitk::mitk_sh_functions::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/itkDiffusionMultiShellQballReconstructionImageFilter.cpp b/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.cpp
index d221ea4d11..56b4b4a22a 100644
--- a/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.cpp
+++ b/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.cpp
@@ -1,1181 +1,1202 @@
 /*=========================================================================
 
 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 __itkDiffusionMultiShellQballReconstructionImageFilter_cpp
 #define __itkDiffusionMultiShellQballReconstructionImageFilter_cpp
 
 #include <itkDiffusionMultiShellQballReconstructionImageFilter.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>
 #include <strstream>
 
 #define _USE_MATH_DEFINES
 #include <math.h>
 
 #include "mitkSphericalHarmonicsFunctions.h"
 #include "mitkVNLVectorFunctions.h"
 #include "itkPointShell.h"
 #include <memory>
 
 
 namespace itk {
 
 #define QBALL_ANAL_RECON_PI       M_PI
 
 template< class T, class TG, class TO, int L, int NODF>
 DiffusionMultiShellQballReconstructionImageFilter<T,TG,TO,L,NODF>
 ::DiffusionMultiShellQballReconstructionImageFilter() :
   m_GradientDirectionContainer(NULL),
   m_NumberOfGradientDirections(0),
   m_NumberOfBaselineImages(1),
   m_Threshold(NumericTraits< ReferencePixelType >::NonpositiveMin()),
   m_BValue(1.0),
   m_Lambda(0.0),
   m_IsHemisphericalArrangementOfGradientDirections(false),
   m_IsArithmeticProgession(false),
   m_ReconstructionType(Standard1Shell)
 {
   // 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::DiffusionMultiShellQballReconstructionImageFilter< TReferenceImagePixelType,TGradientImagePixelType,TOdfPixelType, NOrderL,NrOdfDirections>::OdfPixelType itk::DiffusionMultiShellQballReconstructionImageFilter<TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType, NOrderL, NrOdfDirections>
 ::Normalize( OdfPixelType odf, typename NumericTraits<ReferencePixelType>::AccumulateType b0 )
 {
   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;
 
   return odf;
 }
 
 template<class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NOrderL, int NrOdfDirections>
 void itk::DiffusionMultiShellQballReconstructionImageFilter<TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType,NOrderL, NrOdfDirections>
 ::Threshold(vnl_vector<double> & vec, double sigma)
 {
   for(int i = 0; i < vec.size(); i++)
   {
     if(vec[i] < 0)
     {
       vec[i] = sigma / 2 ;
     }
     if(0 <= vec[i] && vec[i] < sigma )
     {
       vec[i] = sigma / 2 + (vec[i] * vec[i]) / (2 * sigma);
     }
     if( 1 - sigma <= vec[i] && vec[i] < 1 )
     {
       vec[i] = 1-(sigma/2) - ((( 1 - vec[i]) * (1 - vec[i] )) / (2 * sigma) );
     }
     if( 1 <= vec[i] )
     {
       vec[i] = 1 - (sigma / 2);
     }
   }
 }
 
 template<class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NOrderL, int NrOdfDirections>
 void itk::DiffusionMultiShellQballReconstructionImageFilter<TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType,NOrderL, NrOdfDirections>
 ::Threshold(vnl_matrix<double> & mat, double sigma)
 {
   for(int i = 0; i < mat.rows(); i++)
   {
     for( int j = 0; j < mat.cols(); j++ ){
 
       if(mat(i,j) < 0)
       {
         mat(i,j) = sigma / 2 ;
       }
       if(0 <= mat(i,j) && mat(i,j) < sigma )
       {
         mat(i,j) = sigma / 2 + (mat(i,j) * mat(i,j)) / (2 * sigma);
       }
       if( 1 - sigma <= mat(i,j) && mat(i,j) < 1 )
       {
         mat(i,j) = 1-(sigma/2) - (((1 - mat(i,j)) *( 1 - mat(i,j)))  / (2 * sigma) );
       }
       if( 1 <= mat(i,j) )
       {
         mat(i,j) = 1 - (sigma / 2);
       }
     }
   }
 }
 
 template<class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NOrderL, int NrOdfDirections>
 void itk::DiffusionMultiShellQballReconstructionImageFilter<TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType,NOrderL, NrOdfDirections>
 ::Projection1( vnl_matrix<double> & E, double delta )
 {
 
   const double sF = sqrt(5.0);
 
   vnl_vector<double> vOnes(E.rows());
   vOnes.fill(1.0);
 
   vnl_matrix<double> T0(E.rows(), E.cols());
   T0.fill(0.0);
   vnl_matrix<unsigned char> C(E.rows(), 7);
   C.fill(0);
   vnl_matrix<double> A(E.rows(), 7);
   A.fill(0.0);
   vnl_matrix<double> B(E.rows(), 7);
   B.fill(0.0);
 
   vnl_vector<double> s0(E.rows());
   s0.fill(0.0);
   vnl_vector<double> a0(E.rows());
   a0.fill(0.0);
   vnl_vector<double> b0(E.rows());
   b0.fill(0.0);
   vnl_vector<double> ta(E.rows());
   ta.fill(0.0);
   vnl_vector<double> tb(E.rows());
   tb.fill(0.0);
   vnl_vector<double> e(E.rows());
   e.fill(0.0);
   vnl_vector<double> m(E.rows());
   m.fill(0.0);
   vnl_vector<double> a(E.rows());
   a.fill(0.0);
   vnl_vector<double> b(E.rows());
   b.fill(0.0);
 
   // logarithmierung aller werte in E
   for(int i = 0 ; i < E.rows(); i++)
   {
     for(int j = 0 ; j < E.cols(); j++)
     {
       T0(i,j) = -log(E(i,j));
     }
   }
 
 
   //T0 = -T0.apply(std::log);
 
   // Summeiere Zeilenweise über alle Shells sum = E1+E2+E3
   for(int i = 0 ; i < E.rows(); i++)
   {
     s0[i] = T0(i,0) + T0(i,1) + T0(i,2);
   }
 
 
   for(int i = 0; i < E.rows(); i ++)
   {
     // Alle Signal-Werte auf der Ersten shell E(N,0) normiert auf s0
     a0 = E(i,0) / s0[i];
     // Alle Signal-Werte auf der Zweiten shell E(N,1) normiert auf s0
     b0 = E(i,1) / s0[i];
   }
 
   ta = a0 * 3.0;
   tb = b0 * 3.0;
   e = tb - (ta * 2.0);
   m = (tb *  2.0 ) + ta;
 
   for(int i = 0; i < E.rows(); i++)
   {
     C(i,0) = tb[i] < 1+3*delta && 0.5+1.5*(sF+1)*delta < ta[i] && ta[i] < 1-3* (sF+2) *delta;
     C(i,1) = e[i] <= -1 +3*(2*sF+5)* delta && ta[i] >= 1-3*(sF+2)*delta;
     C(i,2) = m[i] > 3 -3*sF*delta && -1+3*(2*sF+5)*delta < e[i] && e[i]<-3*sF*delta;
     C(i,3) = m[i] >= 3-3*sF*delta && e[i] >= -3 *sF * delta;
     C(i,4) = 2.5 + 1.5*(5+sF)*delta < m[i] && m[i] < 3-3*sF*delta && e[i] > -3*sF*delta;
     C(i,5) = ta[i] <= 0.5+1.5 *(sF+1)*delta && m[i] <= 2.5 + 1.5 *(5+sF) * delta;
     C(i,6) = !( C(i,0) || C(i,1) || C(i,2) || C(i,3) || C(i,4) || C(i,5) ); // ~ANY(C(i,[0-5] ),2)
 
     A(i,0)=(bool)C(i,0) * a0(i);
     A(i,1)=(bool)C(i,1) * (1.0/3.0-(sF+2)*delta);
     A(i,2)=(bool)C(i,2) * (0.2+0.8*a0(i)-0.4*b0(i)-delta/sF);
     A(i,3)=(bool)C(i,3) * (0.2+delta/sF);
     A(i,4)=(bool)C(i,4) * (0.2*a0(i)+0.4*b0(i)+2*delta/sF);
     A(i,5)=(bool)C(i,5) * (1.0/6.0+0.5*(sF+1)*delta);
     A(i,6)=(bool)C(i,6) * a0(i);
 
     B(i,0)=(bool)C(i,0) * (1.0/3.0+delta);
     B(i,1)=(bool)C(i,1) * (1.0/3.0+delta);
     B(i,2)=(bool)C(i,2) * (0.4-0.4*a0(i)+0.2*b0(i)-2*delta/sF);
     B(i,3)=(bool)C(i,3) * (0.4-3*delta/sF);
     B(i,4)=(bool)C(i,4) * (0.4*a0(i)+0.8*b0(i)-delta/sF);
     B(i,5)=(bool)C(i,5) * (1.0/3.0+delta);
     B(i,6)=(bool)C(i,6) * b0(i);
   }
 
 
 
   for(int i = 0 ; i < E.rows(); i++)
   {
     double sumA = 0;
     double sumB = 0;
     for(int j = 0 ; j < 7; j++)
     {
       sumA += A(i,j);
       sumB += B(i,j);
     }
     a[i] = sumA;
     b[i] = sumB;
   }
 
   for(int i = 0; i < E.rows(); i++)
   {
     E(i,0) = exp(-(a[i]*s0[i]));
     E(i,1) = exp(-(b[i]*s0[i]));
     E(i,2) = exp(-((1-a[i]-b[i])*s0[i]));
   }
 
 }
 
 
 
 
 template<class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NOrderL, int NrOdfDirections>
 void itk::DiffusionMultiShellQballReconstructionImageFilter<TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType,NOrderL, NrOdfDirections>
 ::Projection2( vnl_vector<double> & A, vnl_vector<double> & a, vnl_vector<double> & b, double delta0)
 {
 
   const double s6 = sqrt(6);
   const double s15 = s6/2.0;
 
   vnl_vector<double> delta(a.size());
   delta.fill(delta0);
 
   vnl_matrix<double> AM(a.size(), 15);
   AM.fill(0.0);
   vnl_matrix<double> aM(a.size(), 15);
   aM.fill(0.0);
   vnl_matrix<double> bM(a.size(), 15);
   bM.fill(0.0);
 
   vnl_matrix<unsigned char> B(a.size(), 15);
   B.fill(0);
 
 
   AM.set_column(0, A);
   AM.set_column(1, A);
   AM.set_column(2, A);
   AM.set_column(3, delta);
   AM.set_column(4, (A+a-b - (delta*s6))/3.0);
   AM.set_column(5, delta);
   AM.set_column(6, delta);
   AM.set_column(7, delta);
   AM.set_column(8, A);
   AM.set_column(9, 0.2*(a*2+A-2*(s6+1)*delta));
   AM.set_column(10,0.2*(b*(-2)+A+2-2*(s6+1)*delta));
   AM.set_column(11, delta);
   AM.set_column(12, delta);
   AM.set_column(13, delta);
   AM.set_column(14, 0.5-(1+s15)*delta);
 
 
   aM.set_column(0, a);
   aM.set_column(1, a);
   aM.set_column(2, -delta + 1);
   aM.set_column(3, a);
   aM.set_column(4, (A*2+a*5+b+s6*delta)/6.0);
   aM.set_column(5, a);
   aM.set_column(6, -delta + 1);
   aM.set_column(7, 0.5*(a+b)+(1+s15)*delta);
   aM.set_column(8, -delta + 1);
   aM.set_column(9, 0.2*(a*4+A*2+(s6+1)*delta));
   aM.set_column(10, -delta + 1);
   aM.set_column(11, (s6+3)*delta);
   aM.set_column(12, -delta + 1);
   aM.set_column(13, -delta + 1);
   aM.set_column(14, -delta + 1);
 
   bM.set_column(0, b);
   bM.set_column(1, delta);
   bM.set_column(2, b);
   bM.set_column(3, b);
   bM.set_column(4, (A*(-2)+a+b*5-s6*delta)/6.0);
   bM.set_column(5, delta);
   bM.set_column(6, b);
   bM.set_column(7, 0.5*(a+b)-(1+s15)*delta);
   bM.set_column(8, delta);
   bM.set_column(9, delta);
   bM.set_column(10, 0.2*(b*4-A*2+1-(s6+1)*delta));
   bM.set_column(11, delta);
   bM.set_column(12, delta);
   bM.set_column(13, -delta*(s6+3) + 1);
   bM.set_column(14, delta);
 
   delta0 *= 0.99;
 
   for(int i = 0 ; i < a.size(); i ++)
   {
     for(int j = 0 ; j < 15; j ++)
     {
       B(i,j) = delta0 < AM(i,j) && 2 * (AM(i,j) + delta0 * s15) < aM(i,j) - bM(i,j) && bM(i,j) > delta0 && aM(i,j) < 1- delta0;
     }
   }
 
   vnl_matrix<double> R2(a.size(), 15);
   R2.fill(0.0);
   vnl_matrix<double> A_(a.size(), 15);
   A_.fill(0.0);
   vnl_matrix<double> a_(a.size(), 15);
   a_.fill(0.0);
   vnl_matrix<double> b_(a.size(), 15);
   b_.fill(0.0);
 
 
 
   vnl_matrix<double> OnesVecMat(1, 15);
   OnesVecMat.fill(1.0);
 
   vnl_matrix<double> AVecMat(a.size(), 1);
   AVecMat.set_column(0,A);
 
   vnl_matrix<double> aVecMat(a.size(), 1);
   aVecMat.set_column(0,a);
 
   vnl_matrix<double> bVecMat(a.size(), 1);
   bVecMat.set_column(0,b);
 
   A_ = AM - (AVecMat * OnesVecMat);
   a_ = aM - (aVecMat * OnesVecMat);
   b_ = bM - (bVecMat * OnesVecMat);
 
   for(int i = 0 ; i < a.size(); i++)
     for(int j = 0 ; j < 15; j++)
     {
       A_(i,j) *= A_(i,j);
       a_(i,j) *= a_(i,j);
       b_(i,j) *= b_(i,j);
     }
 
 
   R2 = A_ + a_ + b_;
 
   for(int i = 0 ; i < a.size(); i ++)
   {
     for(int j = 0 ; j < 15; j ++)
     {
       if(B(i,j) == 0) R2(i,j) = 1e20;
     }
   }
 
   std::vector<unsigned int> indicies(a.size());
 
   // suche den spalten-index der zu der kleinsten Zahl einer Zeile korrespondiert
   for(int i = 0 ; i < a.size(); i++)
   {
     unsigned int index = 0;
     double minvalue = 999;
     for(int j = 0 ; j < 15 ; j++)
     {
       if(R2(i,j) < minvalue){
         minvalue = R2(i,j);
         index = j;
       }
     }
     indicies[i] = index;
   }
 
   for(int i = 0 ; i < a.size(); i++)
   {
     A[i] = AM(i,indicies[i]);
     a[i] = aM(i,indicies[i]);
     b[i] = bM(i,indicies[i]);
   }
 }
 
 
 template<class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NOrderL, int NrOdfDirections>
 void itk::DiffusionMultiShellQballReconstructionImageFilter<TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType,NOrderL, NrOdfDirections>
 ::S_S0Normalization( vnl_vector<double> & vec, typename NumericTraits<ReferencePixelType>::AccumulateType b0 )
 {
 
   double b0f = (double)b0;
   for(int i = 0; i < vec.size(); i++)
   {
     if (b0f==0)
       b0f = 0.01;
     if(vec[i] >= b0f)
       vec[i] = b0f - 0.001;
     vec[i] /= b0f;
   }
 
 }
 
 
 template<class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NOrderL, int NrOdfDirections>
 void itk::DiffusionMultiShellQballReconstructionImageFilter<TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType,NOrderL, NrOdfDirections>
 ::S_S0Normalization( vnl_matrix<double> & mat, typename NumericTraits<ReferencePixelType>::AccumulateType b0 )
 {
   double b0f = (double)b0;
   for(int i = 0; i < mat.rows(); i++)
   {
     for( int j = 0; j < mat.cols(); j++ ){
       if (b0f==0)
         b0f = 0.01;
       if(mat(i,j) >= b0f)
         mat(i,j) = b0f - 0.001;
       mat(i,j) /= b0f;
     }
   }
 }
 
 
 
 template< class T, class TG, class TO, int L, int NODF>
 void DiffusionMultiShellQballReconstructionImageFilter<T,TG,TO,L,NODF>
 ::DoubleLogarithm(vnl_vector<double> & vec)
 {
   for(int i = 0; i < vec.size(); i++)
   {
     vec[i] = log(-log(vec[i]));
   }
 }
 
 template< class T, class TG, class TO, int L, int NODF>
 void DiffusionMultiShellQballReconstructionImageFilter<T,TG,TO,L,NODF>
 ::SetGradientImage( GradientDirectionContainerType *gradientDirection, const GradientImagesType *gradientImage )
 {
   this->m_GradientDirectionContainer = gradientDirection;
   this->m_NumberOfBaselineImages = 0;
 
   for(GradientDirectionContainerType::Iterator it = gradientDirection->Begin(); it != gradientDirection->End(); it++)
   {
     if( it.Value().one_norm() <= 0.0 )
       this->m_NumberOfBaselineImages ++;
     else
       it.Value() = it.Value() / it.Value().two_norm();
   }
 
   this->m_NumberOfGradientDirections = gradientDirection->Size() - 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 DiffusionMultiShellQballReconstructionImageFilter<T,TG,TO,L,NODF>
 ::BeforeThreadedGenerateData()
 {
   itk::TimeProbe clock;
   clock.Start();
 
   if( m_NumberOfGradientDirections < (((L+1)*(L+2))/2) /* && m_NumberOfGradientDirections < 6 */ )
   {
     itkExceptionMacro( << "At least " << ((L+1)*(L+2))/2 << " 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 );
 
   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();
 
   // if no GradientIndexMap is set take all Gradients to GradientIndicies-Vector and add this to GradientIndexMap[ALL]
   // add All Bzero
   if(m_GradientIndexMap.size() == 0){
 
     // split for all gradient directions the image in baseline-indicies and gradient-indicies for a fast access
     GradientDirectionContainerType::ConstIterator gdcit;
     for( gdcit = this->m_GradientDirectionContainer->Begin(); gdcit != this->m_GradientDirectionContainer->End(); ++gdcit)
     {
       if(gdcit.Value().one_norm() <= 0.0) // if vector length at current position is 0, it is a basline signal
       {
         m_GradientIndexMap[0].push_back(gdcit.Index());
       }else{ // it is a gradient signal of length != 0
         m_GradientIndexMap[1].push_back(gdcit.Index());
       }
     }
     m_ReconstructionType = Standard1Shell;
   }else if(m_GradientIndexMap.size() == 4){
 
     GradientIndexMapIteraotr it = m_GradientIndexMap.begin();
     it++;
     int b1 = (*it).first;
     it++;
     int b2 = (*it).first;
     it++;
     int b3 = (*it).first;
 
 
     if(b2 - b1 == b1 && b3 - b2 == b1 )
     {
       m_ReconstructionType = Analytical3Shells;
     }else
     {
       m_ReconstructionType = NumericalNShells;
     }
   }else if(m_GradientIndexMap.size() > 2)
   {
     m_ReconstructionType = NumericalNShells;
   }
 
 
 
   switch(m_ReconstructionType)
   {
   case Analytical3Shells:
   {
     GradientIndexMapIteraotr it = m_GradientIndexMap.begin();
     it++;
     this->ComputeReconstructionMatrix((*it).second);
     break;
   }
   case NumericalNShells:
     this->ComputeReconstructionMatrix(m_GradientIndexMap[1]);
     break;
   case Standard1Shell:
     this->ComputeReconstructionMatrix(m_GradientIndexMap[1]);
     break;
   }
 
   clock.Stop();
   MITK_INFO << "Before GenerateData : " << clock.GetTotal();
 
 }
 
 
 
 template< class T, class TG, class TO, int L, int NODF>
 void DiffusionMultiShellQballReconstructionImageFilter<T,TG,TO,L,NODF>
 ::StandardOneShellReconstruction(const OutputImageRegionType& outputRegionForThread)
 {
   // Get output image pointer
   typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0));
 
   // ImageRegionIterator for the output image
   ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread);
   oit.GoToBegin();
 
   // ImageRegionIterator for the BZero (output) image
   ImageRegionIterator< BZeroImageType > oit2(m_BZeroImage, outputRegionForThread);
   oit2.GoToBegin();
 
   IndiciesVector BZeroIndicies = m_GradientIndexMap[0];
   IndiciesVector SignalIndicies = m_GradientIndexMap[1];
 
   // if the gradient directiosn aragement is hemispherical, duplicate all gradient directions
   // alone, interested in the value, the direction can be neglected
   if(m_IsHemisphericalArrangementOfGradientDirections){
     int NumbersOfGradientIndicies = SignalIndicies.size();
     for (int i = 0 ; i < NumbersOfGradientIndicies; i++)
       SignalIndicies.push_back(SignalIndicies[i]);
   }
 
 
   // Get input gradient image pointer
   typename GradientImagesType::Pointer gradientImagePointer = static_cast< GradientImagesType * >( this->ProcessObject::GetInput(0) );
 
   //  Const ImageRegionIterator for input gradient image
   typedef ImageRegionConstIterator< GradientImagesType > GradientIteratorType;
   GradientIteratorType git(gradientImagePointer, outputRegionForThread );
   git.GoToBegin();
 
   typedef typename GradientImagesType::PixelType         GradientVectorType;
 
   // iterate overall voxels of the gradient image region
   while( ! git.IsAtEnd() )
   {
     GradientVectorType b = git.Get();
 
     // compute the average bzero signal
     typename NumericTraits<ReferencePixelType>::AccumulateType b0 = NumericTraits<ReferencePixelType>::Zero;
     for(unsigned int i = 0; i < BZeroIndicies.size(); ++i)
     {
       b0 += b[BZeroIndicies[i]];
     }
     b0 /= BZeroIndicies.size();
 
     // ODF Vector
     OdfPixelType odf(0.0);
 
     // Create the Signal Vector
     vnl_vector<double> SignalVector(m_NumberOfGradientDirections);
     if( (b0 != 0) && (b0 >= m_Threshold) )
     {
 
       for( unsigned int i = 0; i< SignalIndicies.size(); i++ )
       {
         SignalVector[i] = static_cast<double>(b[SignalIndicies[i]]);
       }
 
       // apply threashold an generate ln(-ln(E)) signal
       // Replace SignalVector with PreNormalized SignalVector
       S_S0Normalization(SignalVector, b0);
       DoubleLogarithm(SignalVector);
 
       // ODF coeffs-vector
       vnl_vector<double> coeffs(m_NumberCoefficients);
 
       // approximate ODF coeffs
       coeffs = ( (*m_CoeffReconstructionMatrix) * SignalVector );
       coeffs[0] += 1.0/(2.0*sqrt(QBALL_ANAL_RECON_PI));
 
       odf = mitk::mitk_vnl_function::element_cast<double, TO>(( (*m_ODFSphericalHarmonicBasisMatrix) * coeffs )).data_block();
     }
     // set ODF to ODF-Image
     oit.Set( odf );
     ++oit;
 
     // set b0(average) to b0-Image
     oit2.Set( b0 );
     ++oit2;
 
     ++git;
 
   }
 
   MITK_INFO << "One Thread finished reconstruction";
 }
 
 #include <itkTimeProbe.h>
+#include <vnl/algo/vnl_levenberg_marquardt.h>
+#include <vnl/algo/vnl_lsqr.h>
+
+class LMFunction : public vnl_least_squares_function
+{
+    void f(vnl_vector<double> const& x, vnl_vector<double>& fx)
+    {
+
+    }
+};
+
+
+template< class T, class TG, class TO, int L, int NODF>
+void DiffusionMultiShellQballReconstructionImageFilter<T,TG,TO,L,NODF>
+::NumericalNShellReconstruction(const OutputImageRegionType& outputRegionForThread)
+{
+
+
+  vnl_levenberg_marquardt LMOptimizer = new vnl_levenberg_marquardt();
+
+}
 
 template< class T, class TG, class TO, int L, int NODF>
 void DiffusionMultiShellQballReconstructionImageFilter<T,TG,TO,L,NODF>
 ::AnalyticalThreeShellReconstruction(const OutputImageRegionType& outputRegionForThread)
 {
 
 
   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();
 
   IndiciesVector BZeroIndicies = m_GradientIndexMap[0];
 
   //MAP-Pair(B-value, SignalShell1Indiecies)
   GradientIndexMapIteraotr it = m_GradientIndexMap.begin();
   // it = b-vlaue = 0
 
   it++;
   // it = b-value = 1000
   IndiciesVector Shell1Indiecies = (*it).second;
 
   it++;
   // it = b-value = 2000
   IndiciesVector Shell2Indiecies = (*it).second;
 
   it++;
   // it = b-value = 3000
   IndiciesVector Shell3Indiecies = (*it).second;
 
   if(m_IsHemisphericalArrangementOfGradientDirections){
     int NumbersOfGradientIndicies = Shell1Indiecies.size();
     for (int i = 0 ; i < NumbersOfGradientIndicies; i++){
       Shell1Indiecies.push_back(Shell1Indiecies[i]);
       Shell2Indiecies.push_back(Shell2Indiecies[i]);
       Shell3Indiecies.push_back(Shell3Indiecies[i]);
     }
   }
 
   // Get input gradient image pointer
   typename GradientImagesType::Pointer gradientImagePointer = static_cast< GradientImagesType * >( this->ProcessObject::GetInput(0) );
 
   //  Const ImageRegionIterator for input gradient image
   typedef ImageRegionConstIterator< GradientImagesType > GradientIteratorType;
   GradientIteratorType git(gradientImagePointer, outputRegionForThread );
   git.GoToBegin();
 
   typedef typename GradientImagesType::PixelType         GradientVectorType;
 
   //------------------------- Preperations Zone ------------------------------------
 
   // N x 3
 
   //  x E1 , E2 , E3
   //
   //  N :    :    :
   //    :    :    :
   //
   vnl_matrix< double > * E = new vnl_matrix<double>(Shell1Indiecies.size(), 3);
   vnl_vector<double> * AlphaValues = new vnl_vector<double>(Shell1Indiecies.size());
   vnl_vector<double> * BetaValues = new vnl_vector<double>(Shell1Indiecies.size());
   vnl_vector<double> * LAValues = new vnl_vector<double>(Shell1Indiecies.size());
   vnl_vector<double> * PValues = new vnl_vector<double>(Shell1Indiecies.size());
   OdfPixelType odf(0.0);
   //------------------------- Preperations Zone END ---------------------------------
 
   // iterate overall voxels of the gradient image region
   while( ! git.IsAtEnd() )
   {
     GradientVectorType b = git.Get();
 
     // compute the average bzero signal
     typename NumericTraits<ReferencePixelType>::AccumulateType b0 = NumericTraits<ReferencePixelType>::Zero;
     for(unsigned int i = 0; i < BZeroIndicies.size(); ++i)
     {
       b0 += b[BZeroIndicies[i]];
     }
     b0 /= BZeroIndicies.size();
 
     if( (b0 != 0) && (b0 >= m_Threshold) )
     {
 
       for(int i = 0 ; i < Shell1Indiecies.size(); i++)
       {
         E->put(i,0, static_cast<double>(b[Shell1Indiecies[i]]));
         E->put(i,1, static_cast<double>(b[Shell2Indiecies[i]]));
         E->put(i,2, static_cast<double>(b[Shell3Indiecies[i]]));
       }
 
       //Approximated-Signal by SH fit - using the specific shell directions and values
       E->set_column(1, (*m_SHBasisMatrix) * ((*m_SignalReonstructionMatrix) * (E->get_column(1))));
       E->set_column(2, (*m_SHBasisMatrix) * ((*m_SignalReonstructionMatrix) * (E->get_column(2))));
 
       // Si/S0
       S_S0Normalization(*E,b0);
 
       //Implements Eq. [19] and Fig. 4.
       Threshold(*E);
 
       //inqualities [31]. Taking the lograithm of th first tree inqualities
       //convert the quadratic inqualities to linear ones.
       Projection1(*E);
 
       double E1, E2, E3, P2,A,B2,B,P,alpha,beta,lambda, ER1, ER2;
 
       for( unsigned int i = 0; i< Shell1Indiecies.size(); i++ )
       {
 
         E1 = E->get(i,0);
         E2 = E->get(i,1);
         E3 = E->get(i,2);
 
         P2 = E2-E1*E1;
         A = (E3 -E1*E2) / ( 2* P2);
         B2 = A * A -(E1 * E3 - E2 * E2) /P2;
         B = 0;
         if(B2 > 0) B = sqrt(B2);
 
         P = 0;
         if(P2 > 0) P = sqrt(P2);
 
         alpha = A + B;
         beta = A - B;
 
 
         lambda = 0.5 + 0.5 * std::sqrt(1 - std::pow((P * 2 ) / (alpha - beta), 2));;
         ER1 = std::fabs(lambda * (alpha - beta) + (beta - E1 ))
             + std::fabs(lambda * (std::pow(alpha, 2) - std::pow(beta, 2)) + (std::pow(beta, 2) - E2 ))
             + std::fabs(lambda * (std::pow(alpha, 3) - std::pow(beta, 3)) + (std::pow(beta, 3) - E3 ));
         ER2 = std::fabs((lambda-1) * (alpha - beta) + (beta - E1 ))
             + std::fabs((lambda-1) * (std::pow(alpha, 2) - std::pow(beta, 2)) + (std::pow(beta, 2) - E2 ))
             + std::fabs((lambda-1) * (std::pow(alpha, 3) - std::pow(beta, 3)) + (std::pow(beta, 3) - E3 ));
 
 
         PValues->put(i, P);
         AlphaValues->put(i, alpha);
         BetaValues->put(i, beta);
         LAValues->put(i,(lambda * (ER1 < ER2)) + ((1-lambda) * (ER1 >= ER2)));
 
       }
 
       Projection2(*PValues, *AlphaValues, *BetaValues);
 
       Threshold(*AlphaValues);
       Threshold(*BetaValues);
 
       DoubleLogarithm(*AlphaValues);
       DoubleLogarithm(*BetaValues);
 
       //L * A + (1-L) * B = S
       //L * A + B - L * B = S
       //L*(A - B) + B = S
 
       vnl_vector<double> SignalVector(element_product((*LAValues) , (*AlphaValues)-(*BetaValues)) + (*BetaValues));
       vnl_vector<double> coeffs((*m_CoeffReconstructionMatrix) *SignalVector );
 
       // the first coeff is a fix value
       coeffs[0] += 1.0/(2.0*sqrt(QBALL_ANAL_RECON_PI));
 
       // Cast the Signal-Type from double to float for the ODF-Image
       odf = mitk::mitk_vnl_function::element_cast<double, TO>( (*m_ODFSphericalHarmonicBasisMatrix) * coeffs ).data_block();
       odf *= (QBALL_ANAL_RECON_PI*4/NODF);
       Normalize(odf,b0);
     }
 
     // set ODF to ODF-Image
     oit.Set( odf );
     ++oit;
 
     // set b0(average) to b0-Image
     oit2.Set( b0 );
     ++oit2;
 
     ++git;
   }
   MITK_INFO << "THREAD FINISHED";
   delete E;
   delete AlphaValues;
   delete BetaValues;
   delete PValues;
   delete LAValues;
 }
 
 
 template< class T, class TG, class TO, int L, int NODF>
 vnl_vector<TO> DiffusionMultiShellQballReconstructionImageFilter<T,TG,TO,L,NODF>
 ::AnalyticalThreeShellParameterEstimation(const IndiciesVector * shell1Indicies,const IndiciesVector * shell2Indicies ,const IndiciesVector * shell3Indicies, vnl_vector<TO>)
 {
 
 
 
 }
 
 template< class T, class TG, class TO, int L, int NODF>
 void DiffusionMultiShellQballReconstructionImageFilter<T,TG,TO,L,NODF>
 ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int NumberOfThreads)
 {
   itk::TimeProbe clock;
   clock.Start();
   switch(m_ReconstructionType)
   {
   case Standard1Shell:
     StandardOneShellReconstruction(outputRegionForThread);
     break;
   case Analytical3Shells:
     AnalyticalThreeShellReconstruction(outputRegionForThread);
     break;
   case NumericalNShells:
     break;
   }
   clock.Stop();
   MITK_INFO << "Reconstruction in : " << clock.GetTotal() << " TU";
 
 
 
 }
 
 template< class T, class TG, class TO, int L, int NODF>
 void DiffusionMultiShellQballReconstructionImageFilter<T, TG, TO, L, NODF>::
 ComputeSphericalHarmonicsBasis(vnl_matrix<double> * QBallReference, vnl_matrix<double> *SHBasisOutput, vnl_matrix<double>* LaplaciaBaltramiOutput, vnl_vector<int>* SHOrderAssociation, vnl_matrix<double>* SHEigenvalues )
 {
   for(unsigned int i=0; i< (*SHBasisOutput).rows(); 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;
 
         // Compute SHBasisFunctions
         double phi = (*QBallReference)(0,i);
         double th = (*QBallReference)(1,i);
-        (*SHBasisOutput)(i,j) = mitk::ShericalHarmonicsFunctions::Yj(m,k,th,phi);
+        (*SHBasisOutput)(i,j) = mitk::mitk_sh_functions::Yj(m,k,th,phi);
 
         // Laplacian Baltrami Order Association
         if(LaplaciaBaltramiOutput)
           (*LaplaciaBaltramiOutput)(j,j) = k*k*(k + 1)*(k+1);
 
         // SHEigenvalues with order Accosiation kj
         if(SHEigenvalues)
           (*SHEigenvalues)(j,j) = -k* (k+1);
 
         // Order Association
         if(SHOrderAssociation)
           (*SHOrderAssociation)[j] = k;
 
       }
     }
   }
 }
 
 template< class T, class TG, class TO, int L, int NODF>
 void DiffusionMultiShellQballReconstructionImageFilter<T,TG,TO,L,NODF>
 ::ComputeFunkRadonTransformationMatrix(vnl_vector<int>* SHOrderAssociationReference, vnl_matrix<double>* FRTMatrixOutput )
 {
   for(int i=0; i<m_NumberCoefficients; i++)
   {
-    (*FRTMatrixOutput)(i,i) = 2.0 * QBALL_ANAL_RECON_PI * mitk::ShericalHarmonicsFunctions::legendre0((*SHOrderAssociationReference)[i]);
+    (*FRTMatrixOutput)(i,i) = 2.0 * QBALL_ANAL_RECON_PI * mitk::mitk_sh_functions::legendre0((*SHOrderAssociationReference)[i]);
   }
 }
 
 template< class T, class TG, class TO, int L, int NOdfDirections>
 bool DiffusionMultiShellQballReconstructionImageFilter<T,TG,TO,L,NOdfDirections>
 ::CheckHemisphericalArrangementOfGradientDirections()
 {
   // 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;
   GradientDirectionContainerType::ConstIterator gdcit1;
   for( 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)
   {
     return false;
   }
   return true;
 }
 
 template< class T, class TG, class TO, int L, int NOdfDirections>
 void DiffusionMultiShellQballReconstructionImageFilter<T,TG,TO,L,NOdfDirections>
 ::ComputeReconstructionMatrix(std::vector< unsigned int > gradientIndiciesVector)
 {
 
   typedef std::auto_ptr< vnl_matrix< double> >  MatrixDoublePtr;
   typedef std::auto_ptr< vnl_vector< int > >    VectorIntPtr;
   typedef std::auto_ptr< vnl_matrix_inverse< double > >  InverseMatrixDoublePtr;
 
   int numberOfGradientDirections = gradientIndiciesVector.size();
 
 
   if( numberOfGradientDirections < (((L+1)*(L+2))/2) || numberOfGradientDirections < 6  )
   {
     itkExceptionMacro( << "At least (L+1)(L+2)/2 gradient directions are required" );
   }
 
   CheckDuplicateDiffusionGradients();
 
 
   // check if gradient directions are arrangement as a hemisphere(true) or sphere(false)
   m_IsHemisphericalArrangementOfGradientDirections = CheckHemisphericalArrangementOfGradientDirections();
   if(m_IsHemisphericalArrangementOfGradientDirections) numberOfGradientDirections *= 2;
 
 
   MatrixDoublePtr Q(new vnl_matrix<double>(3, numberOfGradientDirections));
   Q->fill(0.0);
   // Cartesian to spherical coordinates
   {
     int j = 0;
 
     for(int i = 0; i < gradientIndiciesVector.size(); i++)
     {
 
       double x = this->m_GradientDirectionContainer->ElementAt(gradientIndiciesVector[i]).get(0);
       double y = this->m_GradientDirectionContainer->ElementAt(gradientIndiciesVector[i]).get(1);
       double z = this->m_GradientDirectionContainer->ElementAt(gradientIndiciesVector[i]).get(2);
       double cart[3];
-      mitk::ShericalHarmonicsFunctions::Cart2Sph(x,y,z,cart);
+      mitk::mitk_sh_functions::Cart2Sph(x,y,z,cart);
       (*Q)(0,j) = cart[0];
       (*Q)(1,j) = cart[1];
       (*Q)(2,j++) = cart[2];
 
     }
     if(m_IsHemisphericalArrangementOfGradientDirections)
     {
       for(int i = 0; i < gradientIndiciesVector.size(); i++)
       {
         double x = this->m_GradientDirectionContainer->ElementAt(gradientIndiciesVector[i]).get(0);
         double y = this->m_GradientDirectionContainer->ElementAt(gradientIndiciesVector[i]).get(1);
         double z = this->m_GradientDirectionContainer->ElementAt(gradientIndiciesVector[i]).get(2);
         double cart[3];
-        mitk::ShericalHarmonicsFunctions::Cart2Sph(x,y,z,cart);
+        mitk::mitk_sh_functions::Cart2Sph(x,y,z,cart);
         (*Q)(0,j) = cart[0];
         (*Q)(1,j) = cart[1];
         (*Q)(2,j++) = cart[2];
       }
     }
   }
 
   int LOrder = L;
   m_NumberCoefficients = (int)(LOrder*LOrder + LOrder + 2.0)/2.0 + LOrder;
 
   m_SHBasisMatrix = new vnl_matrix<double>(numberOfGradientDirections,m_NumberCoefficients);
   m_SHBasisMatrix->fill(0.0);
   VectorIntPtr SHOrderAssociation(new vnl_vector<int>(m_NumberCoefficients));
   SHOrderAssociation->fill(0.0);
   MatrixDoublePtr LaplacianBaltrami(new vnl_matrix<double>(m_NumberCoefficients,m_NumberCoefficients));
   LaplacianBaltrami->fill(0.0);
   MatrixDoublePtr FRTMatrix(new vnl_matrix<double>(m_NumberCoefficients,m_NumberCoefficients));
   FRTMatrix->fill(0.0);
   MatrixDoublePtr SHEigenvalues(new vnl_matrix<double>(m_NumberCoefficients,m_NumberCoefficients));
   SHEigenvalues->fill(0.0);
 
 
 
   // SHBasis-Matrix + LaplacianBaltrami-Matrix + SHOrderAssociationVector
   ComputeSphericalHarmonicsBasis(Q.get() ,m_SHBasisMatrix, LaplacianBaltrami.get(), SHOrderAssociation.get(), SHEigenvalues.get());
 
   // Compute FunkRadon Transformation Matrix Associated to SHBasis Order lj
   ComputeFunkRadonTransformationMatrix(SHOrderAssociation.get() ,FRTMatrix.get());
 
   MatrixDoublePtr temp(new vnl_matrix<double>(((m_SHBasisMatrix->transpose()) * (*m_SHBasisMatrix)) + (m_Lambda  * (*LaplacianBaltrami))));
 
   InverseMatrixDoublePtr pseudo_inv(new vnl_matrix_inverse<double>((*temp)));
   MatrixDoublePtr inverse(new vnl_matrix<double>(m_NumberCoefficients,m_NumberCoefficients));
   inverse->fill(0.0);
   (*inverse) = pseudo_inv->inverse();
 
   // ODF Factor ( missing 1/4PI ?? )
   double factor = (1.0/(16.0*QBALL_ANAL_RECON_PI*QBALL_ANAL_RECON_PI));
 
   m_SignalReonstructionMatrix = new vnl_matrix<double>((*inverse) * (m_SHBasisMatrix->transpose()));
 
   m_CoeffReconstructionMatrix = new vnl_matrix<double>(( factor * ((*FRTMatrix) * ((*SHEigenvalues) * (*m_SignalReonstructionMatrix))) ));
 
   // this code goes to the image adapter coeffs->odfs later
 
   vnl_matrix_fixed<double, 3, NOdfDirections>* U = itk::PointShell<NOdfDirections, vnl_matrix_fixed<double, 3, NOdfDirections> >::DistributePointShell();
 
   m_ODFSphericalHarmonicBasisMatrix  = new vnl_matrix<double>(NOdfDirections,m_NumberCoefficients);
   m_ODFSphericalHarmonicBasisMatrix->fill(0.0);
 
   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];
-    mitk::ShericalHarmonicsFunctions::Cart2Sph(x,y,z,cart);
+    mitk::mitk_sh_functions::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<=LOrder; 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_ODFSphericalHarmonicBasisMatrix)(i,j) = mitk::ShericalHarmonicsFunctions::Yj(m,k,th,phi);
+        (*m_ODFSphericalHarmonicBasisMatrix)(i,j) = mitk::mitk_sh_functions::Yj(m,k,th,phi);
       }
     }
   }
 
   m_ReconstructionMatrix = new vnl_matrix<double>((*m_ODFSphericalHarmonicBasisMatrix) * (*m_CoeffReconstructionMatrix));
 
 }
 
 template< class T, class TG, class TO, int L, int NODF>
 template< class VNLType >
 void DiffusionMultiShellQballReconstructionImageFilter<T,TG,TO,L,NODF>
 ::printMatrix( VNLType * mat  )
 {
 
   std::stringstream stream;
 
   for(int i = 0 ; i < mat->rows(); i++)
   {
     stream.str("");
     for(int j = 0; j < mat->cols(); j++)
     {
       stream << (*mat)(i,j) << "  ";
     }
 
   }
 
 
   MITK_INFO << stream.str();
 
 }
 
 template< class T, class TG, class TO, int L, int NODF>
 bool DiffusionMultiShellQballReconstructionImageFilter<T,TG,TO,L,NODF>
 ::CheckDuplicateDiffusionGradients()
 {
 
   GradientDirectionContainerType::ConstIterator gdcit1;
   GradientDirectionContainerType::ConstIterator gdcit2;
 
   for(gdcit1 = this->m_GradientDirectionContainer->Begin(); gdcit1 != this->m_GradientDirectionContainer->End(); ++gdcit1)
   {
     for(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." );
         return true;
       }
     }
   }
   return false;
 }
 
 
 template< class T, class TG, class TO, int L, int NODF>
 void DiffusionMultiShellQballReconstructionImageFilter<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 // __itkDiffusionMultiShellQballReconstructionImageFilter_cpp
diff --git a/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h b/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h
index ad3324735a..f2118fdcd5 100644
--- a/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h
+++ b/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h
@@ -1,235 +1,236 @@
 /*=========================================================================
 
 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 __itkDiffusionMultiShellQballReconstructionImageFilter_h_
 #define __itkDiffusionMultiShellQballReconstructionImageFilter_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"
 #include <iomanip>
 
 
 
 namespace itk{
 /** \class DiffusionMultiShellQballReconstructionImageFilter
     Aganj_2010
  */
 
 template< class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NOrderL, int NrOdfDirections>
 class DiffusionMultiShellQballReconstructionImageFilter : public ImageToImageFilter< Image< TReferenceImagePixelType, 3 >, Image< Vector< TOdfPixelType, NrOdfDirections >, 3 > >
 {
 
 public:
 
     typedef DiffusionMultiShellQballReconstructionImageFilter Self;
     typedef SmartPointer<Self>                      Pointer;
     typedef SmartPointer<const Self>                ConstPointer;
     typedef ImageToImageFilter< Image< TReferenceImagePixelType, 3>, Image< Vector< TOdfPixelType, NrOdfDirections >, 3 > > Superclass;
     typedef TReferenceImagePixelType                  ReferencePixelType;
     typedef TGradientImagePixelType                   GradientPixelType;
     typedef Vector< TOdfPixelType, NrOdfDirections >  OdfPixelType;
     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;
 
     typedef std::map<double, std::vector<unsigned int> > GradientIndexMap;
     typedef std::map<double, std::vector<unsigned int> >::iterator GradientIndexMapIteraotr;
     typedef std::vector<unsigned int> IndiciesVector;
 
     // --------------------------------------------------------------------------------------------//
 
     /** Method for creation through the object factory. */
     itkNewMacro(Self);
 
     /** Runtime information support. */
     itkTypeMacro(DiffusionMultiShellQballReconstructionImageFilter, ImageToImageFilter);
 
     /** 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 SetGradientIndexMap(GradientIndexMap gradientIdexMap)
     {
 
 
         m_GradientIndexMap = gradientIdexMap;
 
         /*
         std::stringstream s1, s2, s3;
         for(int i = 0; i < m_GradientIndexMap[1000].size() ; i++){
             s1 << m_GradientIndexMap[1000][i] << " ";
             s2 << m_GradientIndexMap[2000][i] << " ";
             s3 << m_GradientIndexMap[3000][i] << " ";
         }
 
         MITK_INFO << "1 SHELL " << std::endl << s1.str();
         MITK_INFO << "2 SHELL " << std::endl << s2.str();
         MITK_INFO << "3 SHELL " << std::endl << s3.str();
         */
 
 
     }
 
     /** 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_GradientDirectionContainer->Size() )
         {
             itkExceptionMacro( << "Gradient direction " << idx << "does not exist" );
         }
         return m_GradientDirectionContainer->ElementAt( idx+1 );
     }
 
     OdfPixelType Normalize(OdfPixelType odf, typename NumericTraits<ReferencePixelType>::AccumulateType b0 );
 
     void S_S0Normalization( vnl_vector<double> & vec, typename NumericTraits<ReferencePixelType>::AccumulateType b0  = 0 );
     void S_S0Normalization( vnl_matrix<double> & mat, typename NumericTraits<ReferencePixelType>::AccumulateType b0  = 0 );
 
     void DoubleLogarithm(vnl_vector<double> & vec);
 
     void Threshold(vnl_vector<double> & vec, double sigma = 0.0001);
     void Threshold(vnl_matrix<double> & mat, double sigma = 0.0001);
 
     void Projection1( vnl_matrix<double> & mat, double delta = 0.0001);
     void Projection2( vnl_vector<double> & A, vnl_vector<double> & alpha, vnl_vector<double> & beta, double delta = 0.0001);
 
     /** 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 );
 
     itkGetMacro( BZeroImage, typename BZeroImageType::Pointer);
     //itkGetMacro( ODFSumImage, typename BlaImage::Pointer);
 
     itkSetMacro( BValue, TOdfPixelType);
 
     itkSetMacro( Lambda, double );
     itkGetMacro( Lambda, double );
 
     itkGetConstReferenceMacro( BValue, TOdfPixelType);
 
 
 protected:
     DiffusionMultiShellQballReconstructionImageFilter();
     ~DiffusionMultiShellQballReconstructionImageFilter() {};
     void PrintSelf(std::ostream& os, Indent indent) const;
 
     void ComputeReconstructionMatrix(std::vector< unsigned int > gradientIndiciesVector);
     bool CheckDuplicateDiffusionGradients();
     void ComputeSphericalHarmonicsBasis(vnl_matrix<double>* QBallReference, vnl_matrix<double>* SHBasisOutput, vnl_matrix<double>* LaplaciaBaltramiOutput, vnl_vector<int>* SHOrderAssociation , vnl_matrix<double> * SHEigenvalues);
     void ComputeFunkRadonTransformationMatrix(vnl_vector<int>* SHOrderAssociationReference, vnl_matrix<double>* FRTMatrixOutput );
     bool CheckHemisphericalArrangementOfGradientDirections();
 
     void BeforeThreadedGenerateData();
     void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, int NumberOfThreads );
 
     vnl_vector<TOdfPixelType> AnalyticalThreeShellParameterEstimation(const IndiciesVector * shell1, const IndiciesVector * shell2, const IndiciesVector * shell3, vnl_vector<TOdfPixelType> b);
     void StandardOneShellReconstruction(const OutputImageRegionType& outputRegionForThread);
     void AnalyticalThreeShellReconstruction(const OutputImageRegionType& outputRegionForThread);
+    void NumericalNShellReconstruction(const OutputImageRegionType& outputRegionForThread);
 
 private:
 
     enum ReconstructionType
     {
         Analytical3Shells,
         NumericalNShells,
         Standard1Shell
     };
 
 
 
     CoefficientMatrixType m_ReconstructionMatrix;
     CoefficientMatrixType m_CoeffReconstructionMatrix;
     CoefficientMatrixType m_ODFSphericalHarmonicBasisMatrix;
     CoefficientMatrixType m_SignalReonstructionMatrix;
     CoefficientMatrixType m_SHBasisMatrix;
 
     /** 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;
 
     GradientIndexMap m_GradientIndexMap;
 
     double m_Lambda;
 
     bool m_IsHemisphericalArrangementOfGradientDirections;
 
     bool m_IsArithmeticProgession;
 
     int m_NumberCoefficients;
 
     ReconstructionType m_ReconstructionType;
 
 
 
     template< class VNLType >
     void printMatrix( VNLType * mat );
 
 
 
 };
 
 }
 
 #ifndef ITK_MANUAL_INSTANTIATION
 #include "itkDiffusionMultiShellQballReconstructionImageFilter.cpp"
 #endif
 
 #endif //__itkDiffusionMultiShellQballReconstructionImageFilter_h_
 
diff --git a/Modules/DiffusionImaging/Reconstruction/mitkSphericalHarmonicsFunctions.cpp b/Modules/DiffusionImaging/Reconstruction/mitkSphericalHarmonicsFunctions.cpp
index 6317215e60..07d4d24ca6 100644
--- a/Modules/DiffusionImaging/Reconstruction/mitkSphericalHarmonicsFunctions.cpp
+++ b/Modules/DiffusionImaging/Reconstruction/mitkSphericalHarmonicsFunctions.cpp
@@ -1,145 +1,91 @@
 /*=========================================================================
 
 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>
+#include <boost/math/special_functions/spherical_harmonic.hpp>
 #endif
 #endif
 
 #include "mitkSphericalHarmonicsFunctions.h"
 namespace mitk{
 
 #define SPHERICAL_HARMONICS_PI       M_PI
 
-double mitk::ShericalHarmonicsFunctions::factorial(int number) {
+double mitk::mitk_sh_functions::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)
+void mitk::mitk_sh_functions::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)
+double mitk::mitk_sh_functions::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)
+double mitk::mitk_sh_functions::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);
+        return sqrt(2.0) * ::boost::math::spherical_harmonic_r<double, double>(k,m,theta,phi);;
     }
 
     if( m == 0 )
-        return spherical_harmonic(0,k,theta,phi,false);
+        return  ::boost::math::spherical_harmonic_r<double, double>(k,m,theta,phi);
 
     if( 0 < m && m <= k )
     {
-        return sqrt(2.0) * spherical_harmonic(m,k,theta,phi,true);
+        return sqrt(2.0) * ::boost::math::spherical_harmonic_i<double, double>(k,m,theta,phi);
     }
 
     return 0;
 }
+
 }
diff --git a/Modules/DiffusionImaging/Reconstruction/mitkSphericalHarmonicsFunctions.h b/Modules/DiffusionImaging/Reconstruction/mitkSphericalHarmonicsFunctions.h
index 9abc121e44..f185b00ceb 100644
--- a/Modules/DiffusionImaging/Reconstruction/mitkSphericalHarmonicsFunctions.h
+++ b/Modules/DiffusionImaging/Reconstruction/mitkSphericalHarmonicsFunctions.h
@@ -1,38 +1,37 @@
 /*=========================================================================
 
 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{
+namespace mitk_sh_functions{
 
-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);
+   double factorial(int number);
+   void Cart2Sph(double x, double y, double z, double* cart);
+   double legendre0(int l);
+   double spherical_harmonic(int m,int l,double theta,double phi, bool complexPart);
+   double Yj(int m, int k, double theta, double phi);
 
-};
+}
 }
 
 #endif //__mitkSphericalHarmonicsHandler_h_