diff --git a/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.cpp b/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.cpp
new file mode 100644
index 0000000000..35d2f71350
--- /dev/null
+++ b/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.cpp
@@ -0,0 +1,882 @@
+/*=========================================================================
+
+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>
+
+#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"
+
+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_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::DiffusionMultiShellQballReconstructionImageFilter<
+      TReferenceImagePixelType,TGradientImagePixelType,TOdfPixelType,
+      NOrderL,NrOdfDirections>::OdfPixelType
+      itk::DiffusionMultiShellQballReconstructionImageFilter
+      <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::DiffusionMultiShellQballReconstructionImageFilter
+        <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 << "DiffusionMultiShellQballReconstructionImageFilter: 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 << "DiffusionMultiShellQballReconstructionImageFilter: 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 << "DiffusionMultiShellQballReconstructionImageFilter: 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 DiffusionMultiShellQballReconstructionImageFilter<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 DiffusionMultiShellQballReconstructionImageFilter<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 DiffusionMultiShellQballReconstructionImageFilter<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 DiffusionMultiShellQballReconstructionImageFilter<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 DiffusionMultiShellQballReconstructionImageFilter<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 DiffusionMultiShellQballReconstructionImageFilter<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 DiffusionMultiShellQballReconstructionImageFilter<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 DiffusionMultiShellQballReconstructionImageFilter<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 DiffusionMultiShellQballReconstructionImageFilter<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);
+              (*Q)(0,i) = cart[0];
+              (*Q)(1,i) = cart[1];
+              (*Q)(2,i++) = cart[2];
+            }
+          }
+          if(m_DirectionsDuplicated)
+          {
+            for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin();
+              gdcit != this->m_GradientDirectionContainer->End(); ++gdcit)
+            {
+              if(gdcit.Value().one_norm() > 0.0)
+              {
+                double x = gdcit.Value().get(0);
+                double y = gdcit.Value().get(1);
+                double z = gdcit.Value().get(2);
+                double cart[3];
+                Cart2Sph(x,y,z,cart);
+                (*Q)(0,i) = cart[0];
+                (*Q)(1,i) = cart[1];
+                (*Q)(2,i++) = cart[2];
+              }
+            }
+          }
+        }
+
+        int l = L;
+        m_NumberCoefficients = (int)(l*l + l + 2.0)/2.0 + l;
+        vnl_matrix<double>* B = new vnl_matrix<double>(m_NumberOfGradientDirections,m_NumberCoefficients);
+        vnl_matrix<double>* _L = new vnl_matrix<double>(m_NumberCoefficients,m_NumberCoefficients);
+        _L->fill(0.0);
+        vnl_matrix<double>* LL = new vnl_matrix<double>(m_NumberCoefficients,m_NumberCoefficients);
+        LL->fill(0.0);
+        vnl_matrix<double>* P = new vnl_matrix<double>(m_NumberCoefficients,m_NumberCoefficients);
+        P->fill(0.0);
+        vnl_matrix<double>* Inv = new vnl_matrix<double>(m_NumberCoefficients,m_NumberCoefficients);
+        P->fill(0.0);
+        vnl_vector<int>* lj = new vnl_vector<int>(m_NumberCoefficients);
+        m_LP = new vnl_vector<double>(m_NumberCoefficients);
+
+        for(unsigned int i=0; i<m_NumberOfGradientDirections; i++)
+        {
+          for(int k=0; k<=l; k+=2)
+          {
+            for(int m=-k; m<=k; m++)
+            {
+              int j = (k*k + k + 2)/2 + m - 1;
+              (*_L)(j,j) = -k*(k+1);
+              (*m_LP)(j) = -k*(k+1);
+              (*lj)[j] = k;
+              double phi = (*Q)(0,i);
+              double th = (*Q)(1,i);
+              (*B)(i,j) = Yj(m,k,th,phi);
+              //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]);
+            (*m_LP)(i) *= (*P)(i,i);
+          }
+          else
+          {
+            (*P)(i,i) = legendre0((*lj)[i]);
+          }
+        }
+        m_B_t = new vnl_matrix<double>(B->transpose());
+        //tofile2(&m_B_t,"m_B_t");
+        vnl_matrix<double> B_t_B = (*m_B_t) * (*B);
+        //tofile2(&B_t_B,"B_t_B");
+        vnl_matrix<double> lambdaLL(m_NumberCoefficients,m_NumberCoefficients);
+        lambdaLL.update((*LL));
+        lambdaLL *= m_Lambda;
+        //tofile2(&lambdaLL,"lLL");
+
+        vnl_matrix<double> tmp( B_t_B + lambdaLL);
+        vnl_matrix_inverse<double> *pseudoInverse
+          = new vnl_matrix_inverse<double>( tmp );
+
+        (*Inv) = pseudoInverse->pinverse();
+        //tofile2(Inv,"Inv");
+        vnl_matrix<double> temp((*Inv) * (*m_B_t));
+        double fac1 = (1.0/(16.0*QBALL_ANAL_RECON_PI*QBALL_ANAL_RECON_PI));
+        switch(m_NormalizationMethod)
+        {
+        case QBAR_ADC_ONLY:
+        case QBAR_RAW_SIGNAL:
+          break;
+        case QBAR_STANDARD:
+        case QBAR_B_ZERO_B_VALUE:
+        case QBAR_B_ZERO:
+        case QBAR_NONE:
+          temp = (*P) * temp;
+          break;
+        case QBAR_SOLID_ANGLE:
+          temp = fac1 * (*P) * (*_L) * temp;
+          break;
+        case QBAR_NONNEG_SOLID_ANGLE:
+          break;
+        }
+
+        //tofile2(&temp,"A");
+
+        m_CoeffReconstructionMatrix = new vnl_matrix<TO>(m_NumberCoefficients,m_NumberOfGradientDirections);
+        for(int i=0; i<m_NumberCoefficients; i++)
+        {
+          for(unsigned int j=0; j<m_NumberOfGradientDirections; j++)
+          {
+            (*m_CoeffReconstructionMatrix)(i,j) = (float) temp(i,j);
+          }
+        }
+
+        // this code goes to the image adapter coeffs->odfs later
+
+        int NOdfDirections = NODF;
+        vnl_matrix_fixed<double, 3, NODF>* U =
+          itk::PointShell<NODF, vnl_matrix_fixed<double, 3, NODF> >::DistributePointShell();
+
+        m_SphericalHarmonicBasisMatrix  = new vnl_matrix<TO>(NOdfDirections,m_NumberCoefficients);
+        vnl_matrix<double>* sphericalHarmonicBasisMatrix2
+          = new vnl_matrix<double>(NOdfDirections,m_NumberCoefficients);
+        for(int i=0; i<NOdfDirections; i++)
+        {
+          double x = (*U)(0,i);
+          double y = (*U)(1,i);
+          double z = (*U)(2,i);
+          double cart[3];
+          Cart2Sph(x,y,z,cart);
+          (*U)(0,i) = cart[0];
+          (*U)(1,i) = cart[1];
+          (*U)(2,i) = cart[2];
+        }
+
+        for(int i=0; i<NOdfDirections; i++)
+        {
+          for(int k=0; k<=l; k+=2)
+          {
+            for(int m=-k; m<=k; m++)
+            {
+              int j = (k*k + k + 2)/2 + m - 1;
+              double phi = (*U)(0,i);
+              double th = (*U)(1,i);
+              (*m_SphericalHarmonicBasisMatrix)(i,j) = Yj(m,k,th,phi);
+              (*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 DiffusionMultiShellQballReconstructionImageFilter<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 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
new file mode 100644
index 0000000000..416a814de9
--- /dev/null
+++ b/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h
@@ -0,0 +1,300 @@
+/*=========================================================================
+
+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"
+
+
+namespace itk{
+/** \class DiffusionMultiShellQballReconstructionImageFilter
+ * \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 DiffusionMultiShellQballReconstructionImageFilter :
+  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 DiffusionMultiShellQballReconstructionImageFilter 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(DiffusionMultiShellQballReconstructionImageFilter,
+                                                   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:
+  DiffusionMultiShellQballReconstructionImageFilter();
+  ~DiffusionMultiShellQballReconstructionImageFilter() {};
+  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 "itkDiffusionMultiShellQballReconstructionImageFilter.cpp"
+#endif
+
+#endif //__itkDiffusionMultiShellQballReconstructionImageFilter_h_
+
diff --git a/Modules/DiffusionImaging/files.cmake b/Modules/DiffusionImaging/files.cmake
index c7d1299abe..594c8f32c4 100644
--- a/Modules/DiffusionImaging/files.cmake
+++ b/Modules/DiffusionImaging/files.cmake
@@ -1,172 +1,173 @@
 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/itkDiffusionMultiShellQballReconstructionImageFilter.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 )