diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp
index 02a8ca2d37..261e98ffd5 100644
--- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp
+++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp
@@ -1,818 +1,824 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 #ifndef __itkAnalyticalDiffusionQballReconstructionImageFilter_cpp
 #define __itkAnalyticalDiffusionQballReconstructionImageFilter_cpp
 
 #include <itkAnalyticalDiffusionQballReconstructionImageFilter.h>
 #include <itkImageRegionConstIterator.h>
 #include <itkImageRegionConstIteratorWithIndex.h>
 #include <itkImageRegionIterator.h>
 #include <itkArray.h>
 #include <vnl/vnl_vector.h>
 
 #include <boost/version.hpp>
 #include <stdio.h>
 #include <locale>
 #include <fstream>
 
 #define _USE_MATH_DEFINES
 #include <math.h>
 
 #include <boost/math/special_functions.hpp>
 
 #include "itkPointShell.h"
 
 using namespace boost::math;
 
 namespace itk {
 
 #define QBALL_ANAL_RECON_PI       M_PI
 
 template< class T, class TG, class TO, int L, int NODF>
 AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
 ::AnalyticalDiffusionQballReconstructionImageFilter() :
     m_GradientDirectionContainer(NULL),
     m_NumberOfGradientDirections(0),
     m_NumberOfBaselineImages(1),
     m_Threshold(NumericTraits< ReferencePixelType >::NonpositiveMin()),
     m_BValue(1.0),
     m_Lambda(0.0),
     m_DirectionsDuplicated(false),
     m_Delta1(0.001),
     m_Delta2(0.001),
     m_UseMrtrixBasis(false)
 {
     // At least 1 inputs is necessary for a vector image.
     // For images added one at a time we need at least six
     this->SetNumberOfRequiredInputs( 1 );
 }
 
 
 template<
         class TReferenceImagePixelType,
         class TGradientImagePixelType,
         class TOdfPixelType,
         int NOrderL,
         int NrOdfDirections>
 typename itk::AnalyticalDiffusionQballReconstructionImageFilter<
 TReferenceImagePixelType,TGradientImagePixelType,TOdfPixelType,
 NOrderL,NrOdfDirections>::OdfPixelType
 itk::AnalyticalDiffusionQballReconstructionImageFilter
 <TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType,
 NOrderL, NrOdfDirections>
 ::Normalize( OdfPixelType odf,
              typename NumericTraits<ReferencePixelType>::AccumulateType b0 )
 {
     switch( m_NormalizationMethod )
     {
     case QBAR_STANDARD:
     {
         TOdfPixelType sum = 0;
 
         for(int i=0; i<NrOdfDirections; i++)
         {
             sum += odf[i];
         }
         if(sum>0)
             odf /= sum;
 
         return odf;
         break;
     }
     case QBAR_B_ZERO_B_VALUE:
     {
         for(int i=0; i<NrOdfDirections; i++)
         {
             odf[i] = ((TOdfPixelType)log((TOdfPixelType)b0)-odf[i])/m_BValue;
         }
         return odf;
         break;
     }
     case QBAR_B_ZERO:
     {
         odf *= 1.0/b0;
         return odf;
         break;
     }
     case QBAR_NONE:
     {
         return odf;
         break;
     }
     case QBAR_ADC_ONLY:
     {
         for(int i=0; i<NrOdfDirections; i++)
         {
             odf[i] = ((TOdfPixelType)log((TOdfPixelType)b0)-odf[i])/m_BValue;
         }
         return odf;
         break;
     }
     case QBAR_RAW_SIGNAL:
     {
         return odf;
         break;
     }
     case QBAR_SOLID_ANGLE:
     {
         for(int i=0; i<NrOdfDirections; i++)
             odf[i] *= QBALL_ANAL_RECON_PI*4/NrOdfDirections;
 
         break;
     }
     case QBAR_NONNEG_SOLID_ANGLE:
     {
         break;
     }
     }
 
     return odf;
 }
 
 
 template<
         class TReferenceImagePixelType,
         class TGradientImagePixelType,
         class TOdfPixelType,
         int NOrderL,
         int NrOdfDirections>
 vnl_vector<TOdfPixelType>
 itk::AnalyticalDiffusionQballReconstructionImageFilter
 <TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType,
 NOrderL, NrOdfDirections>
 ::PreNormalize( vnl_vector<TOdfPixelType> vec,
                 typename NumericTraits<ReferencePixelType>::AccumulateType b0 )
 {
     switch( m_NormalizationMethod )
     {
     case QBAR_STANDARD:
     {
+        int n = vec.size();
+        double b0f = (double)b0;
+        for(int i=0; i<n; i++)
+        {
+            vec[i] = vec[i]/b0f;
+        }
         return vec;
         break;
     }
     case QBAR_B_ZERO_B_VALUE:
     {
         int n = vec.size();
         for(int i=0; i<n; i++)
         {
             if (vec[i]<=0)
                 vec[i] = 0.001;
 
             vec[i] = log(vec[i]);
         }
         return vec;
         break;
     }
     case QBAR_B_ZERO:
     {
         return vec;
         break;
     }
     case QBAR_NONE:
     {
         return vec;
         break;
     }
     case QBAR_ADC_ONLY:
     {
         int n = vec.size();
         for(int i=0; i<n; i++)
         {
             if (vec[i]<=0)
                 vec[i] = 0.001;
 
             vec[i] = log(vec[i]);
         }
         return vec;
         break;
     }
     case QBAR_RAW_SIGNAL:
     {
         return vec;
         break;
     }
     case QBAR_SOLID_ANGLE:
     case QBAR_NONNEG_SOLID_ANGLE:
     {
         int n = vec.size();
         double b0f = (double)b0;
         for(int i=0; i<n; i++)
         {
             vec[i] = vec[i]/b0f;
 
             if (vec[i]<0)
                 vec[i] = m_Delta1;
             else if (vec[i]<m_Delta1)
                 vec[i] = m_Delta1/2 + vec[i]*vec[i]/(2*m_Delta1);
             else if (vec[i]>=1)
                 vec[i] = 1-m_Delta2/2;
             else if (vec[i]>=1-m_Delta2)
                 vec[i] = 1-m_Delta2/2-(1-vec[i])*(1-vec[i])/(2*m_Delta2);
 
             vec[i] = log(-log(vec[i]));
         }
         return vec;
         break;
     }
     }
 
     return vec;
 }
 
 template< class T, class TG, class TO, int L, int NODF>
 void AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
 ::BeforeThreadedGenerateData()
 {
     // If we have more than 2 inputs, then each input, except the first is a
     // gradient image. The number of gradient images must match the number of
     // gradient directions.
     //const unsigned int numberOfInputs = this->GetNumberOfInputs();
 
     // There need to be at least 6 gradient directions to be able to compute the
     // tensor basis
     if( m_NumberOfGradientDirections < (L*L + L + 2)/2 + L )
     {
         itkExceptionMacro( << "Not enough gradient directions supplied (" << m_NumberOfGradientDirections << "). At least " << (L*L + L + 2)/2 + L << " needed for SH-order " << L);
     }
 
     // Input must be an itk::VectorImage.
     std::string gradientImageClassName(
                 this->ProcessObject::GetInput(0)->GetNameOfClass());
     if ( strcmp(gradientImageClassName.c_str(),"VectorImage") != 0 )
     {
         itkExceptionMacro( <<
                            "There is only one Gradient image. I expect that to be a VectorImage. "
                            << "But its of type: " << gradientImageClassName );
     }
 
     this->ComputeReconstructionMatrix();
 
     typename GradientImagesType::Pointer img = static_cast< GradientImagesType * >(
                 this->ProcessObject::GetInput(0) );
 
     m_BZeroImage = BZeroImageType::New();
     m_BZeroImage->SetSpacing( img->GetSpacing() );   // Set the image spacing
     m_BZeroImage->SetOrigin( img->GetOrigin() );     // Set the image origin
     m_BZeroImage->SetDirection( img->GetDirection() );  // Set the image direction
     m_BZeroImage->SetLargestPossibleRegion( img->GetLargestPossibleRegion());
     m_BZeroImage->SetBufferedRegion( img->GetLargestPossibleRegion() );
     m_BZeroImage->Allocate();
 
     m_ODFSumImage = BZeroImageType::New();
     m_ODFSumImage->SetSpacing( img->GetSpacing() );   // Set the image spacing
     m_ODFSumImage->SetOrigin( img->GetOrigin() );     // Set the image origin
     m_ODFSumImage->SetDirection( img->GetDirection() );  // Set the image direction
     m_ODFSumImage->SetLargestPossibleRegion( img->GetLargestPossibleRegion());
     m_ODFSumImage->SetBufferedRegion( img->GetLargestPossibleRegion() );
     m_ODFSumImage->Allocate();
 
     m_CoefficientImage = CoefficientImageType::New();
     m_CoefficientImage->SetSpacing( img->GetSpacing() );   // Set the image spacing
     m_CoefficientImage->SetOrigin( img->GetOrigin() );     // Set the image origin
     m_CoefficientImage->SetDirection( img->GetDirection() );  // Set the image direction
     m_CoefficientImage->SetLargestPossibleRegion( img->GetLargestPossibleRegion());
     m_CoefficientImage->SetBufferedRegion( img->GetLargestPossibleRegion() );
     m_CoefficientImage->Allocate();
 
     if(m_NormalizationMethod == QBAR_SOLID_ANGLE || m_NormalizationMethod == QBAR_NONNEG_SOLID_ANGLE)
         m_Lambda = 0.0;
 }
 
 template< class T, class TG, class TO, int L, int NODF>
 void AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
 ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
                        ThreadIdType )
 {
     typename OutputImageType::Pointer outputImage =
             static_cast< OutputImageType * >(this->ProcessObject::GetPrimaryOutput());
 
     ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread);
     oit.GoToBegin();
 
     ImageRegionIterator< BZeroImageType > oit2(m_BZeroImage, outputRegionForThread);
     oit2.GoToBegin();
 
     ImageRegionIterator< FloatImageType > oit3(m_ODFSumImage, outputRegionForThread);
     oit3.GoToBegin();
 
     ImageRegionIterator< CoefficientImageType > oit4(m_CoefficientImage, outputRegionForThread);
     oit4.GoToBegin();
 
     typedef ImageRegionConstIterator< GradientImagesType > GradientIteratorType;
     typedef typename GradientImagesType::PixelType         GradientVectorType;
     typename GradientImagesType::Pointer gradientImagePointer = NULL;
 
     // Would have liked a dynamic_cast here, but seems SGI doesn't like it
     // The enum will ensure that an inappropriate cast is not done
     gradientImagePointer = static_cast< GradientImagesType * >(
                 this->ProcessObject::GetInput(0) );
 
     GradientIteratorType git(gradientImagePointer, outputRegionForThread );
     git.GoToBegin();
 
     // Compute the indicies of the baseline images and gradient images
     std::vector<unsigned int> baselineind; // contains the indicies of
     // the baseline images
     std::vector<unsigned int> gradientind; // contains the indicies of
     // the gradient images
 
     for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin();
         gdcit != this->m_GradientDirectionContainer->End(); ++gdcit)
     {
         if(gdcit.Value().one_norm() <= 0.0)
             baselineind.push_back(gdcit.Index());
         else
             gradientind.push_back(gdcit.Index());
     }
 
     if( m_DirectionsDuplicated )
     {
         int gradIndSize = gradientind.size();
         for(int i=0; i<gradIndSize; i++)
             gradientind.push_back(gradientind[i]);
     }
 
     while( !git.IsAtEnd() )
     {
         GradientVectorType b = git.Get();
 
         typename NumericTraits<ReferencePixelType>::AccumulateType b0 = NumericTraits<ReferencePixelType>::Zero;
 
         // Average the baseline image pixels
         for(unsigned int i = 0; i < baselineind.size(); ++i)
         {
             b0 += b[baselineind[i]];
         }
         b0 /= this->m_NumberOfBaselineImages;
 
         OdfPixelType odf(0.0);
         typename CoefficientImageType::PixelType coeffPixel(0.0);
         vnl_vector<TO> B(m_NumberOfGradientDirections);
 
         if( (b0 != 0) && (b0 >= m_Threshold) )
         {
             for( unsigned int i = 0; i< m_NumberOfGradientDirections; i++ )
             {
                 B[i] = static_cast<TO>(b[gradientind[i]]);
             }
 
             B = PreNormalize(B, b0);
             if(m_NormalizationMethod == QBAR_SOLID_ANGLE)
             {
                 vnl_vector<TO> coeffs(m_NumberCoefficients);
                 coeffs = ( (*m_CoeffReconstructionMatrix) * B );
                 coeffs[0] += 1.0/(2.0*sqrt(QBALL_ANAL_RECON_PI));
                 odf = ( (*m_SphericalHarmonicBasisMatrix) * coeffs ).data_block();
                 coeffPixel = coeffs.data_block();
             }
             else if(m_NormalizationMethod == QBAR_NONNEG_SOLID_ANGLE)
             {
                 /** this would be the place to implement a non-negative
               * solver for quadratic programming problem:
               * min .5*|| Bc-s ||^2 subject to -CLPc <= 4*pi*ones
               * (refer to MICCAI 2009 Goh et al. "Estimating ODFs with PDF constraints")
               * .5*|| Bc-s ||^2 == .5*c'B'Bc - x'B's + .5*s's
               */
 
                 itkExceptionMacro( << "Nonnegative Solid Angle not yet implemented");
             }
             else
             {
                 vnl_vector<TO> coeffs(m_NumberCoefficients);
                 coeffs = ( (*m_CoeffReconstructionMatrix) * B );
                 coeffs[0] += 1.0/(2.0*sqrt(QBALL_ANAL_RECON_PI));
                 coeffPixel = coeffs.data_block();
                 odf = ( (*m_ReconstructionMatrix) * B ).data_block();
             }
             odf = Normalize(odf, b0);
 
         }
 
         oit.Set( odf );
         oit2.Set( b0 );
         float sum = 0;
         for (unsigned int k=0; k<odf.Size(); k++)
             sum += (float) odf[k];
         oit3.Set( sum-1 );
         oit4.Set(coeffPixel);
         ++oit;  // odf image iterator
         ++oit3; // odf sum image iterator
         ++oit2; // b0 image iterator
         ++oit4; // coefficient image iterator
         ++git;  // Gradient  image iterator
     }
 
     std::cout << "One Thread finished reconstruction" << std::endl;
 }
 
 template< class T, class TG, class TO, int L, int NODF>
 void AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
 ::tofile2(vnl_matrix<double> *pA, std::string fname)
 {
     vnl_matrix<double> A = (*pA);
     std::ofstream myfile;
     std::locale C("C");
     std::locale originalLocale = myfile.getloc();
     myfile.imbue(C);
 
     myfile.open (fname.c_str());
     myfile << "A1=[";
     for(unsigned int i=0; i<A.rows(); i++)
     {
         for(unsigned int j=0; j<A.columns(); j++)
         {
             myfile << A(i,j) << " ";
             if(j==A.columns()-1 && i!=A.rows()-1)
                 myfile << ";";
         }
     }
     myfile << "];";
     myfile.close();
 
     myfile.imbue( originalLocale );
 }
 
 template< class T, class TG, class TO, int L, int NODF>
 void AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
 ::Cart2Sph(double x, double y, double z, double *spherical)
 {
     double phi, theta, r;
     r = sqrt(x*x+y*y+z*z);
 
     if( r<mitk::eps )
     {
         theta = M_PI/2;
         phi = M_PI/2;
     }
     else
     {
         theta = acos(z/r);
         phi = atan2(y, x);
     }
     spherical[0] = phi;
     spherical[1] = theta;
     spherical[2] = r;
 }
 
 template< class T, class TG, class TO, int L, int NODF>
 double AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
 ::Yj(int m, int l, double theta, double phi, bool useMRtrixBasis)
 {
     if (!useMRtrixBasis)
     {
         if (m<0)
             return sqrt(2.0)*spherical_harmonic_r(l, -m, theta, phi);
         else if (m==0)
             return spherical_harmonic_r(l, m, theta, phi);
         else
             return pow(-1.0,m)*sqrt(2.0)*spherical_harmonic_i(l, m, theta, phi);
     }
     else
     {
         double plm = legendre_p<double>(l,abs(m),-cos(theta));
         double mag = sqrt((double)(2*l+1)/(4.0*M_PI)*factorial<double>(l-abs(m))/factorial<double>(l+abs(m)))*plm;
         if (m>0)
             return mag*cos(m*phi);
         else if (m==0)
             return mag;
         else
             return mag*sin(-m*phi);
     }
 
     return 0;
 }
 
 template< class T, class TG, class TO, int L, int NODF>
 double AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
 ::Legendre0(int l)
 {
     if( l%2 != 0 )
     {
         return 0;
     }
     else
     {
         double prod1 = 1.0;
         for(int i=1;i<l;i+=2) prod1 *= i;
         double prod2 = 1.0;
         for(int i=2;i<=l;i+=2) prod2 *= i;
         return pow(-1.0,l/2.0)*(prod1/prod2);
     }
 }
 
 template< class T, class TG, class TO, int L, int NODF>
 void AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
 ::ComputeReconstructionMatrix()
 {
 
     //for(int i=-6;i<7;i++)
     //  std::cout << boost::math::legendre_p<double>(6, i, 0.65657) << std::endl;
 
     if( m_NumberOfGradientDirections < (L*L + L + 2)/2 + L )
     {
         itkExceptionMacro( << "Not enough gradient directions supplied (" << m_NumberOfGradientDirections << "). At least " << (L*L + L + 2)/2 + L << " needed for SH-order " << L);
     }
 
     {
         // check for duplicate diffusion gradients
         bool warning = false;
         for(GradientDirectionContainerType::ConstIterator gdcit1 = this->m_GradientDirectionContainer->Begin();
             gdcit1 != this->m_GradientDirectionContainer->End(); ++gdcit1)
         {
             for(GradientDirectionContainerType::ConstIterator gdcit2 = this->m_GradientDirectionContainer->Begin();
                 gdcit2 != this->m_GradientDirectionContainer->End(); ++gdcit2)
             {
                 if(gdcit1.Value() == gdcit2.Value() && gdcit1.Index() != gdcit2.Index())
                 {
                     itkWarningMacro( << "Some of the Diffusion Gradients equal each other. Corresponding image data should be averaged before calling this filter." );
                     warning = true;
                     break;
                 }
             }
             if (warning) break;
         }
 
         // handle acquisition schemes where only half of the spherical
         // shell is sampled by the gradient directions. In this case,
         // each gradient direction is duplicated in negative direction.
         vnl_vector<double> centerMass(3);
         centerMass.fill(0.0);
         int count = 0;
         for(GradientDirectionContainerType::ConstIterator gdcit1 = this->m_GradientDirectionContainer->Begin();
             gdcit1 != this->m_GradientDirectionContainer->End(); ++gdcit1)
         {
             if(gdcit1.Value().one_norm() > 0.0)
             {
                 centerMass += gdcit1.Value();
                 count ++;
             }
         }
         centerMass /= count;
         if(centerMass.two_norm() > 0.1)
         {
             m_DirectionsDuplicated = true;
             m_NumberOfGradientDirections *= 2;
         }
     }
 
     vnl_matrix<double> *Q = new vnl_matrix<double>(3, m_NumberOfGradientDirections);
 
     {
         int i = 0;
         for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin();
             gdcit != this->m_GradientDirectionContainer->End(); ++gdcit)
         {
             if(gdcit.Value().one_norm() > 0.0)
             {
                 double x = gdcit.Value().get(0);
                 double y = gdcit.Value().get(1);
                 double z = gdcit.Value().get(2);
                 double cart[3];
                 Cart2Sph(x,y,z,cart);
                 (*Q)(0,i) = cart[0];
                 (*Q)(1,i) = cart[1];
                 (*Q)(2,i++) = cart[2];
             }
         }
         if(m_DirectionsDuplicated)
         {
             for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin();
                 gdcit != this->m_GradientDirectionContainer->End(); ++gdcit)
             {
                 if(gdcit.Value().one_norm() > 0.0)
                 {
                     double x = gdcit.Value().get(0);
                     double y = gdcit.Value().get(1);
                     double z = gdcit.Value().get(2);
                     double cart[3];
                     Cart2Sph(x,y,z,cart);
                     (*Q)(0,i) = cart[0];
                     (*Q)(1,i) = cart[1];
                     (*Q)(2,i++) = cart[2];
                 }
             }
         }
     }
 
     int l = L;
     m_NumberCoefficients = (int)(l*l + l + 2.0)/2.0 + l;
     vnl_matrix<double>* B = new vnl_matrix<double>(m_NumberOfGradientDirections,m_NumberCoefficients);
     vnl_matrix<double>* _L = new vnl_matrix<double>(m_NumberCoefficients,m_NumberCoefficients);
     _L->fill(0.0);
     vnl_matrix<double>* LL = new vnl_matrix<double>(m_NumberCoefficients,m_NumberCoefficients);
     LL->fill(0.0);
     vnl_matrix<double>* P = new vnl_matrix<double>(m_NumberCoefficients,m_NumberCoefficients);
     P->fill(0.0);
     vnl_matrix<double>* Inv = new vnl_matrix<double>(m_NumberCoefficients,m_NumberCoefficients);
     P->fill(0.0);
     vnl_vector<int>* lj = new vnl_vector<int>(m_NumberCoefficients);
     m_LP = new vnl_vector<double>(m_NumberCoefficients);
 
     for(unsigned int i=0; i<m_NumberOfGradientDirections; i++)
     {
         for(int k=0; k<=l; k+=2)
         {
             for(int m=-k; m<=k; m++)
             {
                 int j = (k*k + k + 2)/2 + m - 1;
                 (*_L)(j,j) = -k*(k+1);
                 (*m_LP)(j) = -k*(k+1);
                 (*lj)[j] = k;
                 double phi = (*Q)(0,i);
                 double th = (*Q)(1,i);
                 (*B)(i,j) = Yj(m,k,th,phi, m_UseMrtrixBasis);
             }
         }
     }
 
     //vnl_matrix<double> temp((*_L)*(*_L));
     LL->update(*_L);
     *LL *= *_L;
     //tofile2(LL,"LL");
 
     for(int i=0; i<m_NumberCoefficients; i++)
     {
         // here we leave out the 2*pi multiplication from Descoteaux
         // this is because we do not want the real integral value
         // but an average value over the equator.
 
         if(m_NormalizationMethod == QBAR_SOLID_ANGLE || m_NormalizationMethod == QBAR_NONNEG_SOLID_ANGLE)
         {
             (*P)(i,i) = 2.0*QBALL_ANAL_RECON_PI*Legendre0((*lj)[i]);
             (*m_LP)(i) *= (*P)(i,i);
         }
         else
         {
             (*P)(i,i) = Legendre0((*lj)[i]);
         }
     }
     m_B_t = new vnl_matrix<double>(B->transpose());
     //tofile2(&m_B_t,"m_B_t");
     vnl_matrix<double> B_t_B = (*m_B_t) * (*B);
     //tofile2(&B_t_B,"B_t_B");
     vnl_matrix<double> lambdaLL(m_NumberCoefficients,m_NumberCoefficients);
     lambdaLL.update((*LL));
     lambdaLL *= m_Lambda;
     //tofile2(&lambdaLL,"lLL");
 
     vnl_matrix<double> tmp( B_t_B + lambdaLL);
     vnl_matrix_inverse<double> *pseudoInverse
             = new vnl_matrix_inverse<double>( tmp );
 
     (*Inv) = pseudoInverse->pinverse();
     //tofile2(Inv,"Inv");
     vnl_matrix<double> temp((*Inv) * (*m_B_t));
     double fac1 = (1.0/(16.0*QBALL_ANAL_RECON_PI*QBALL_ANAL_RECON_PI));
     switch(m_NormalizationMethod)
     {
     case QBAR_ADC_ONLY:
     case QBAR_RAW_SIGNAL:
         break;
     case QBAR_STANDARD:
     case QBAR_B_ZERO_B_VALUE:
     case QBAR_B_ZERO:
     case QBAR_NONE:
         temp = (*P) * temp;
         break;
     case QBAR_SOLID_ANGLE:
         temp = fac1 * (*P) * (*_L) * temp;
         break;
     case QBAR_NONNEG_SOLID_ANGLE:
         break;
     }
 
     //tofile2(&temp,"A");
 
     m_CoeffReconstructionMatrix = new vnl_matrix<TO>(m_NumberCoefficients,m_NumberOfGradientDirections);
     for(int i=0; i<m_NumberCoefficients; i++)
     {
         for(unsigned int j=0; j<m_NumberOfGradientDirections; j++)
         {
             (*m_CoeffReconstructionMatrix)(i,j) = (float) temp(i,j);
         }
     }
 
     // this code goes to the image adapter coeffs->odfs later
 
     int NOdfDirections = NODF;
     vnl_matrix_fixed<double, 3, NODF>* U =
             itk::PointShell<NODF, vnl_matrix_fixed<double, 3, NODF> >::DistributePointShell();
 
     m_SphericalHarmonicBasisMatrix  = new vnl_matrix<TO>(NOdfDirections,m_NumberCoefficients);
     vnl_matrix<double>* sphericalHarmonicBasisMatrix2
             = new vnl_matrix<double>(NOdfDirections,m_NumberCoefficients);
     for(int i=0; i<NOdfDirections; i++)
     {
         double x = (*U)(0,i);
         double y = (*U)(1,i);
         double z = (*U)(2,i);
         double cart[3];
         Cart2Sph(x,y,z,cart);
         (*U)(0,i) = cart[0];
         (*U)(1,i) = cart[1];
         (*U)(2,i) = cart[2];
     }
 
     for(int i=0; i<NOdfDirections; i++)
     {
         for(int k=0; k<=l; k+=2)
         {
             for(int m=-k; m<=k; m++)
             {
                 int j = (k*k + k + 2)/2 + m - 1;
                 double phi = (*U)(0,i);
                 double th = (*U)(1,i);
                 (*m_SphericalHarmonicBasisMatrix)(i,j) = Yj(m,k,th,phi,m_UseMrtrixBasis);
                 (*sphericalHarmonicBasisMatrix2)(i,j) = (*m_SphericalHarmonicBasisMatrix)(i,j);
             }
         }
     }
 
     m_ReconstructionMatrix = new vnl_matrix<TO>(NOdfDirections,m_NumberOfGradientDirections);
     *m_ReconstructionMatrix = (*m_SphericalHarmonicBasisMatrix) * (*m_CoeffReconstructionMatrix);
 
 }
 
 template< class T, class TG, class TO, int L, int NODF>
 void AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
 ::SetGradientImage(const GradientDirectionContainerType *gradientDirection,
                    const GradientImagesType *gradientImage )
 {
     // Copy Gradient Direction Container
     this->m_GradientDirectionContainer = GradientDirectionContainerType::New();
     for(GradientDirectionContainerType::ConstIterator it = gradientDirection->Begin();
         it != gradientDirection->End(); it++)
     {
         this->m_GradientDirectionContainer->push_back(it.Value());
     }
 
 
     unsigned int numImages = gradientDirection->Size();
     this->m_NumberOfBaselineImages = 0;
     for(GradientDirectionContainerType::Iterator it = this->m_GradientDirectionContainer->Begin();
         it != this->m_GradientDirectionContainer->End(); it++)
     {
         if(it.Value().one_norm() <= 0.0)
         {
             this->m_NumberOfBaselineImages++;
         }
         else // Normalize non-zero gradient directions
         {
             it.Value() = it.Value() / it.Value().two_norm();
         }
     }
 
     this->m_NumberOfGradientDirections = numImages - this->m_NumberOfBaselineImages;
 
     // ensure that the gradient image we received has as many components as
     // the number of gradient directions
     if( gradientImage->GetVectorLength() != this->m_NumberOfBaselineImages + m_NumberOfGradientDirections )
     {
         itkExceptionMacro( << m_NumberOfGradientDirections << " gradients + " << this->m_NumberOfBaselineImages
                            << "baselines = " << m_NumberOfGradientDirections + this->m_NumberOfBaselineImages
                            << " directions specified but image has " << gradientImage->GetVectorLength()
                            << " components.");
     }
 
     this->ProcessObject::SetNthInput( 0, const_cast< GradientImagesType* >(gradientImage) );
 }
 
 template< class T, class TG, class TO, int L, int NODF>
 void AnalyticalDiffusionQballReconstructionImageFilter<T,TG,TO,L,NODF>
 ::PrintSelf(std::ostream& os, Indent indent) const
 {
     std::locale C("C");
     std::locale originalLocale = os.getloc();
     os.imbue(C);
 
     Superclass::PrintSelf(os,indent);
 
     os << indent << "OdfReconstructionMatrix: " << m_ReconstructionMatrix << std::endl;
     if ( m_GradientDirectionContainer )
         os << indent << "GradientDirectionContainer: " << m_GradientDirectionContainer << std::endl;
     else
         os << indent << "GradientDirectionContainer: (Gradient directions not set)" << std::endl;
 
     os << indent << "NumberOfGradientDirections: " <<  m_NumberOfGradientDirections << std::endl;
     os << indent << "NumberOfBaselineImages: " << m_NumberOfBaselineImages << std::endl;
     os << indent << "Threshold for reference B0 image: " << m_Threshold << std::endl;
     os << indent << "BValue: " << m_BValue << std::endl;
     os.imbue( originalLocale );
 }
 
 }
 
 #endif // __itkAnalyticalDiffusionQballReconstructionImageFilter_cpp
diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleXReader.cpp b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleXReader.cpp
index 527341766a..139fc45445 100644
--- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleXReader.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleXReader.cpp
@@ -1,302 +1,193 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 
 #include "mitkFiberBundleXReader.h"
 #include <itkMetaDataObject.h>
 #include <vtkPolyData.h>
 #include <vtkDataReader.h>
 #include <vtkPolyDataReader.h>
 #include <vtkMatrix4x4.h>
 #include <vtkPolyLine.h>
 #include <vtkCellArray.h>
 #include <itksys/SystemTools.hxx>
 #include <tinyxml.h>
 #include <vtkCleanPolyData.h>
 #include <mitkTrackvis.h>
 
 namespace mitk
 {
 
 void FiberBundleXReader
 ::GenerateData()
 {
     if ( ( ! m_OutputCache ) )
     {
         Superclass::SetNumberOfRequiredOutputs(0);
         this->GenerateOutputInformation();
     }
 
     if (!m_OutputCache)
     {
         itkWarningMacro("Output cache is empty!");
     }
 
 
     Superclass::SetNumberOfRequiredOutputs(1);
     Superclass::SetNthOutput(0, m_OutputCache.GetPointer());
 }
 
 void FiberBundleXReader::GenerateOutputInformation()
 {
     try
     {
         const std::string& locale = "C";
         const std::string& currLocale = setlocale( LC_ALL, NULL );
         setlocale(LC_ALL, locale.c_str());
 
         std::string ext = itksys::SystemTools::GetFilenameLastExtension(m_FileName);
         ext = itksys::SystemTools::LowerCase(ext);
 
         if (ext==".trk")
         {
-            MITK_INFO << "Importing fiber bundle from TrackVis ...";
             m_OutputCache = OutputType::New();
-            TrackVis trk;
-            trk.open(m_FileName);
-            trk.read(m_OutputCache);
-            MITK_INFO << "... done";
+            TrackVisFiberReader reader;
+            reader.open(m_FileName);
+            reader.read(m_OutputCache);
+
+            mitk::Geometry3D::Pointer geometry = dynamic_cast<mitk::Geometry3D*>(m_OutputCache->GetGeometry());
+
+            mitk::Point3D origin;
+            origin[0]=reader.hdr.origin[0];
+            origin[1]=reader.hdr.origin[1];
+            origin[2]=reader.hdr.origin[2];
+            geometry->SetOrigin(origin);
+
+            mitk::Vector3D spacing;
+            spacing[0]=reader.hdr.voxel_size[0];
+            spacing[1]=reader.hdr.voxel_size[1];
+            spacing[2]=reader.hdr.voxel_size[2];
+            geometry->SetSpacing(spacing);
+
+            m_OutputCache->SetGeometry(geometry);
+
             return;
         }
 
         vtkSmartPointer<vtkDataReader> chooser=vtkSmartPointer<vtkDataReader>::New();
         chooser->SetFileName(m_FileName.c_str() );
         if( chooser->IsFilePolyData())
         {
-            MITK_INFO << "Reading vtk fiber bundle";
             vtkSmartPointer<vtkPolyDataReader> reader = vtkSmartPointer<vtkPolyDataReader>::New();
             reader->SetFileName( m_FileName.c_str() );
             reader->Update();
 
             if ( reader->GetOutput() != NULL )
             {
                 vtkSmartPointer<vtkPolyData> fiberPolyData = reader->GetOutput();
                 m_OutputCache = OutputType::New(fiberPolyData);
             }
         }
-        else // try to read deprecated fiber bundle file format
-        {
-            MITK_INFO << "Reading xml fiber bundle";
-
-            vtkSmartPointer<vtkPolyData> fiberPolyData = vtkSmartPointer<vtkPolyData>::New();
-            vtkSmartPointer<vtkCellArray> cellArray = vtkSmartPointer<vtkCellArray>::New();
-            vtkSmartPointer<vtkPoints>    points = vtkSmartPointer<vtkPoints>::New();
-            TiXmlDocument doc( m_FileName );
-            if(doc.LoadFile())
-            {
-                TiXmlHandle hDoc(&doc);
-                TiXmlElement* pElem;
-                TiXmlHandle hRoot(0);
-
-                pElem = hDoc.FirstChildElement().Element();
-
-                // save this for later
-                hRoot = TiXmlHandle(pElem);
-
-                pElem = hRoot.FirstChildElement("geometry").Element();
-
-                // read geometry
-                mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New();
-
-                // read origin
-                mitk::Point3D origin;
-                double temp = 0;
-                pElem->Attribute("origin_x", &temp);
-                origin[0] = temp;
-                pElem->Attribute("origin_y", &temp);
-                origin[1] = temp;
-                pElem->Attribute("origin_z", &temp);
-                origin[2] = temp;
-                geometry->SetOrigin(origin);
-
-                // read spacing
-          ScalarType spacing[3];
-                pElem->Attribute("spacing_x", &temp);
-                spacing[0] = temp;
-                pElem->Attribute("spacing_y", &temp);
-                spacing[1] = temp;
-                pElem->Attribute("spacing_z", &temp);
-                spacing[2] = temp;
-                geometry->SetSpacing(spacing);
-
-                // read transform
-                vtkMatrix4x4* m = vtkMatrix4x4::New();
-                pElem->Attribute("xx", &temp);
-                m->SetElement(0,0,temp);
-                pElem->Attribute("xy", &temp);
-                m->SetElement(1,0,temp);
-                pElem->Attribute("xz", &temp);
-                m->SetElement(2,0,temp);
-                pElem->Attribute("yx", &temp);
-                m->SetElement(0,1,temp);
-                pElem->Attribute("yy", &temp);
-                m->SetElement(1,1,temp);
-                pElem->Attribute("yz", &temp);
-                m->SetElement(2,1,temp);
-                pElem->Attribute("zx", &temp);
-                m->SetElement(0,2,temp);
-                pElem->Attribute("zy", &temp);
-                m->SetElement(1,2,temp);
-                pElem->Attribute("zz", &temp);
-                m->SetElement(2,2,temp);
-
-                m->SetElement(0,3,origin[0]);
-                m->SetElement(1,3,origin[1]);
-                m->SetElement(2,3,origin[2]);
-                m->SetElement(3,3,1);
-                geometry->SetIndexToWorldTransformByVtkMatrix(m);
-
-                // read bounds
-                float bounds[] = {0, 0, 0, 0, 0, 0};
-                pElem->Attribute("size_x", &temp);
-                bounds[1] = temp;
-                pElem->Attribute("size_y", &temp);
-                bounds[3] = temp;
-                pElem->Attribute("size_z", &temp);
-                bounds[5] = temp;
-                geometry->SetFloatBounds(bounds);
-                geometry->SetImageGeometry(true);
-
-                pElem = hRoot.FirstChildElement("fiber_bundle").FirstChild().Element();
-                for( ; pElem ; pElem=pElem->NextSiblingElement())
-                {
-                    TiXmlElement* pElem2 = pElem->FirstChildElement();
-
-                    vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
-
-                    for( ; pElem2; pElem2=pElem2->NextSiblingElement())
-                    {
-              Point3D point;
-                        pElem2->Attribute("pos_x", &temp);
-                        point[0] = temp;
-                        pElem2->Attribute("pos_y", &temp);
-                        point[1] = temp;
-                        pElem2->Attribute("pos_z", &temp);
-                        point[2] = temp;
-
-                        geometry->IndexToWorld(point, point);
-                        vtkIdType id = points->InsertNextPoint(point.GetDataPointer());
-                        container->GetPointIds()->InsertNextId(id);
-
-                    }
-                    cellArray->InsertNextCell(container);
-                }
-                fiberPolyData->SetPoints(points);
-                fiberPolyData->SetLines(cellArray);
-
-                vtkSmartPointer<vtkCleanPolyData> cleaner = vtkSmartPointer<vtkCleanPolyData>::New();
-                cleaner->SetInputData(fiberPolyData);
-                cleaner->Update();
-                fiberPolyData = cleaner->GetOutput();
-
-                m_OutputCache = OutputType::New(fiberPolyData);
-            }
-            else
-            {
-                MITK_INFO << "could not open xml file";
-                throw "could not open xml file";
-            }
-        }
         setlocale(LC_ALL, currLocale.c_str());
         MITK_INFO << "Fiber bundle read";
     }
     catch(...)
     {
         throw;
     }
 }
 
 void FiberBundleXReader::Update()
 {
     this->GenerateData();
 }
 
 const char* FiberBundleXReader
 ::GetFileName() const
 {
     return m_FileName.c_str();
 }
 
 
 void FiberBundleXReader
 ::SetFileName(const char* aFileName)
 {
     m_FileName = aFileName;
 }
 
 
 const char* FiberBundleXReader
 ::GetFilePrefix() const
 {
     return m_FilePrefix.c_str();
 }
 
 
 void FiberBundleXReader
 ::SetFilePrefix(const char* aFilePrefix)
 {
     m_FilePrefix = aFilePrefix;
 }
 
 
 const char* FiberBundleXReader
 ::GetFilePattern() const
 {
     return m_FilePattern.c_str();
 }
 
 
 void FiberBundleXReader
 ::SetFilePattern(const char* aFilePattern)
 {
     m_FilePattern = aFilePattern;
 }
 
 
 bool FiberBundleXReader
 ::CanReadFile(const std::string filename, const std::string /*filePrefix*/, const std::string /*filePattern*/)
 {
     // First check the extension
     if(  filename == "" )
     {
         return false;
     }
     std::string ext = itksys::SystemTools::GetFilenameLastExtension(filename);
     ext = itksys::SystemTools::LowerCase(ext);
 
     if (ext == ".fib" || ext == ".trk")
     {
         return true;
     }
 
     return false;
 }
 
 BaseDataSource::DataObjectPointer FiberBundleXReader::MakeOutput(const DataObjectIdentifierType &name)
 {
     itkDebugMacro("MakeOutput(" << name << ")");
     if( this->IsIndexedOutputName(name) )
     {
         return this->MakeOutput( this->MakeIndexFromOutputName(name) );
     }
     return static_cast<itk::DataObject*>(OutputType::New().GetPointer());
 }
 
 BaseDataSource::DataObjectPointer FiberBundleXReader::MakeOutput(DataObjectPointerArraySizeType /*idx*/)
 {
     return OutputType::New().GetPointer();
 }
 
 } //namespace MITK
diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleXReader.h b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleXReader.h
index dc864066d2..e05e54afe2 100644
--- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleXReader.h
+++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleXReader.h
@@ -1,82 +1,77 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 
 #ifndef __mitkFiberBundleXReader_h
 #define __mitkFiberBundleXReader_h
 
 #include <mitkCommon.h>
 #include <mitkFileReader.h>
 #include <mitkFiberBundleX.h>
 #include <vtkSmartPointer.h>
 #include <MitkFiberTrackingExports.h>
 
 namespace mitk
 {
 
   /** \brief
   */
 
   class MitkFiberTracking_EXPORT FiberBundleXReader : public FileReader, public BaseDataSource
   {
   public:
 
     /** Types for the standardized TractContainer **/
     /* direct linked includes of mitkFiberBundleX DataStructure */
 
     typedef mitk::FiberBundleX OutputType;
 
-    mitkClassMacro( FiberBundleXReader, BaseDataSource );
+    mitkClassMacro( FiberBundleXReader, BaseDataSource )
     itkFactorylessNewMacro(Self)
     itkCloneMacro(Self)
 
     const char* GetFileName() const;
     void SetFileName(const char* aFileName);
     const char* GetFilePrefix() const;
     void SetFilePrefix(const char* aFilePrefix);
     const char* GetFilePattern() const;
     void SetFilePattern(const char* aFilePattern);
 
     static bool CanReadFile(const std::string filename, const std::string filePrefix, const std::string filePattern);
 
-//    itkGetMacro(GroupFiberBundleX, FiberGroupType::Pointer);
-//    itkGetMacro(TractContainer, ContainerType::Pointer);
-
     virtual void Update();
 
     BaseDataSource::DataObjectPointer MakeOutput(const DataObjectIdentifierType &name);
-
     BaseDataSource::DataObjectPointer MakeOutput( DataObjectPointerArraySizeType idx);
 
   protected:
 
     /** Does the real work. */
     virtual void GenerateData();
     virtual void GenerateOutputInformation();
 
     OutputType::Pointer m_OutputCache;
 
     std::string m_FileName;
     std::string m_FilePrefix;
     std::string m_FilePattern;
 
-
   private:
     void operator=(const Self&); //purposely not implemented
   };
 
 } //namespace MITK
 
 #endif // __mitkFiberBundleXReader_h
diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleXWriter.cpp b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleXWriter.cpp
index c86fb17e9a..1a8d4ce53d 100644
--- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleXWriter.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleXWriter.cpp
@@ -1,135 +1,135 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 
 #include "mitkFiberBundleXWriter.h"
 #include <vtkSmartPointer.h>
 #include <vtkCleanPolyData.h>
 #include <itksys/SystemTools.hxx>
 #include <mitkTrackvis.h>
 #include <itkSize.h>
 
 mitk::FiberBundleXWriter::FiberBundleXWriter()
     : m_FileName(""), m_FilePrefix(""), m_FilePattern(""), m_Success(false)
 {
     this->SetNumberOfRequiredInputs( 1 );
 }
 
 
 mitk::FiberBundleXWriter::~FiberBundleXWriter()
 {}
 
 
 void mitk::FiberBundleXWriter::GenerateData()
 {
   try
   {
     const std::string& locale = "C";
     const std::string& currLocale = setlocale( LC_ALL, NULL );
     setlocale(LC_ALL, locale.c_str());
 
     m_Success = false;
     InputType* input = this->GetInput();
     if (input == NULL)
     {
         itkWarningMacro(<<"Sorry, input to FiberBundleXWriter is NULL!");
         return;
     }
     else if ( m_FileName == "" )
     {
         itkWarningMacro( << "Sorry, filename has not been set!" );
         return ;
     }
 
     std::string ext = itksys::SystemTools::GetFilenameLastExtension(m_FileName);
 
     if (ext==".fib" || ext==".vtk")
     {
         MITK_INFO << "Writing fiber bundle as binary VTK";
         vtkSmartPointer<vtkPolyDataWriter> writer = vtkSmartPointer<vtkPolyDataWriter>::New();
         writer->SetInputData(input->GetFiberPolyData());
         writer->SetFileName(m_FileName.c_str());
         writer->SetFileTypeToBinary();
         writer->Write();
     }
     else if (ext==".afib")
     {
         itksys::SystemTools::ReplaceString(m_FileName,".afib",".fib");
         MITK_INFO << "Writing fiber bundle as ascii VTK";
         vtkSmartPointer<vtkPolyDataWriter> writer = vtkSmartPointer<vtkPolyDataWriter>::New();
         writer->SetInputData(input->GetFiberPolyData());
         writer->SetFileName(m_FileName.c_str());
         writer->SetFileTypeToASCII();
         writer->Write();
     }
     else if (ext==".avtk")
     {
         itksys::SystemTools::ReplaceString(m_FileName,".avtk",".vtk");
         MITK_INFO << "Writing fiber bundle as ascii VTK";
         vtkSmartPointer<vtkPolyDataWriter> writer = vtkSmartPointer<vtkPolyDataWriter>::New();
         writer->SetInputData(input->GetFiberPolyData());
         writer->SetFileName(m_FileName.c_str());
         writer->SetFileTypeToASCII();
         writer->Write();
     }
     else if (ext==".trk")
     {
         MITK_INFO << "Writing fiber bundle as TRK";
-        TrackVis trk;
+        TrackVisFiberReader trk;
         trk.create(m_FileName, input);
         trk.writeHdr();
         trk.append(input);
     }
 
     setlocale(LC_ALL, currLocale.c_str());
     m_Success = true;
     MITK_INFO << "Fiber bundle written";
   }
   catch(...)
   {
     throw;
   }
 }
 
 
 void mitk::FiberBundleXWriter::SetInputFiberBundleX( InputType* diffVolumes )
 {
     this->ProcessObject::SetNthInput( 0, diffVolumes );
 }
 
 
 mitk::FiberBundleX* mitk::FiberBundleXWriter::GetInput()
 {
     if ( this->GetNumberOfInputs() < 1 )
     {
         return NULL;
     }
     else
     {
         return dynamic_cast<InputType*> ( this->ProcessObject::GetInput( 0 ) );
     }
 }
 
 
 std::vector<std::string> mitk::FiberBundleXWriter::GetPossibleFileExtensions()
 {
   std::vector<std::string> possibleFileExtensions;
   possibleFileExtensions.push_back(".fib");
   possibleFileExtensions.push_back(".afib");
   possibleFileExtensions.push_back(".vtk");
   possibleFileExtensions.push_back(".avtk");
   possibleFileExtensions.push_back(".trk");
   return possibleFileExtensions;
 }
diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkTrackvis.cpp b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkTrackvis.cpp
index 6b49d112dc..6100173580 100644
--- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkTrackvis.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkTrackvis.cpp
@@ -1,199 +1,198 @@
 #include <mitkTrackvis.h>
 
-TrackVis::TrackVis()  { filename = ""; fp = NULL; maxSteps = 20000; }
+TrackVisFiberReader::TrackVisFiberReader()  { filename = ""; fp = NULL; }
 
-TrackVis::~TrackVis() { if (fp) fclose( fp ); }
+TrackVisFiberReader::~TrackVisFiberReader() { if (fp) fclose( fp ); }
 
 
 // Create a TrackVis file and store standard metadata. The file is ready to append fibers.
 // ---------------------------------------------------------------------------------------
-short TrackVis::create(string filename , mitk::FiberBundleX *fib)
+short TrackVisFiberReader::create(string filename , mitk::FiberBundleX *fib)
 {
     // prepare the header
     for(int i=0; i<3 ;i++)
     {
         if (fib->GetReferenceImage().IsNotNull())
         {
             hdr.dim[i]            = fib->GetReferenceImage()->GetDimension(i);
             hdr.voxel_size[i]     = fib->GetReferenceImage()->GetGeometry()->GetSpacing().GetElement(i);
-            m_Origin[i]         = fib->GetReferenceImage()->GetGeometry()->GetOrigin().GetElement(i);
+            hdr.origin[i]         = fib->GetReferenceImage()->GetGeometry()->GetOrigin().GetElement(i);
         }
         else
         {
             hdr.dim[i]            = 1;
             hdr.voxel_size[i]     = 1;
-            m_Origin[i]         = 0;
+            hdr.origin[i]         = 0;
         }
-        hdr.origin[i] = 0;
     }
     hdr.n_scalars = 0;
     hdr.n_properties = 0;
     sprintf(hdr.voxel_order,"LPS");
     sprintf(hdr.pad2,"LPS");
     hdr.image_orientation_patient[0] = 1.0;
     hdr.image_orientation_patient[1] = 0.0;
     hdr.image_orientation_patient[2] = 0.0;
     hdr.image_orientation_patient[3] = 0.0;
     hdr.image_orientation_patient[4] = 1.0;
     hdr.image_orientation_patient[5] = 0.0;
     hdr.pad1[0] = 0;
     hdr.pad1[1] = 0;
     hdr.invert_x = 0;
     hdr.invert_y = 0;
     hdr.invert_z = 0;
     hdr.swap_xy = 0;
     hdr.swap_yz = 0;
     hdr.swap_zx = 0;
     hdr.n_count = 0;
     hdr.version = 1;
     hdr.hdr_size = 1000;
 
     // write the header to the file
     fp = fopen(filename.c_str(),"w+b");
     if (fp == NULL)
     {
         printf("[ERROR] Unable to create file '%s'\n",filename.c_str());
         return 0;
     }
     sprintf(hdr.id_string,"TRACK");
     if (fwrite((char*)&hdr, 1, 1000, fp) != 1000)
         MITK_ERROR << "TrackVis::create : Error occurding during writing fiber.";
 
     this->filename = filename;
 
     return 1;
 }
 
 
 // Open an existing TrackVis file and read metadata information.
 // The file pointer is positiond at the beginning of fibers data
 // -------------------------------------------------------------
-short TrackVis::open( string filename )
+short TrackVisFiberReader::open( string filename )
 {
     fp = fopen(filename.c_str(),"r+b");
     if (fp == NULL)
     {
         printf("[ERROR] Unable to open file '%s'\n",filename.c_str());
         return 0;
     }
     this->filename = filename;
 
     return fread((char*)(&hdr), 1, 1000, fp);
 }
 
 
 
 // Append a fiber to the file
 // --------------------------
-short TrackVis::append(mitk::FiberBundleX *fib)
+short TrackVisFiberReader::append(mitk::FiberBundleX *fib)
 {
     vtkPolyData* poly = fib->GetFiberPolyData();
     for (int i=0; i<fib->GetNumFibers(); i++)
     {
         vtkCell* cell = poly->GetCell(i);
         int numPoints = cell->GetNumberOfPoints();
         vtkPoints* points = cell->GetPoints();
 
         unsigned int     numSaved, pos = 0;
         //float* tmp = new float[3*maxSteps];
         std::vector< float > tmp;
         tmp.reserve(3*numPoints);
 
-        if ( numPoints > maxSteps )
-        {
-            printf( "[ERROR] Trying to write a fiber too long!\n" );
-            return 0;
-        }
-
         numSaved = numPoints;
         for(unsigned int i=0; i<numSaved ;i++)
         {
             double* p = points->GetPoint(i);
 
-            tmp[pos++] = p[0] - m_Origin[0];
-            tmp[pos++] = p[1] - m_Origin[1];
-            tmp[pos++] = p[2] - m_Origin[2];
+            tmp[pos++] = p[0];
+            tmp[pos++] = p[1];
+            tmp[pos++] = p[2];
         }
 
         // write the coordinates to the file
         if ( fwrite((char*)&numSaved, 1, 4, fp) != 4 )
         {
             printf( "[ERROR] Problems saving the fiber!\n" );
             return 1;
         }
         if ( fwrite((char*)&(tmp.front()), 1, 4*pos, fp) != 4*pos )
         {
             printf( "[ERROR] Problems saving the fiber!\n" );
             return 1;
         }
     }
 
     return 0;
 }
 
-
-
 //// Read one fiber from the file
 //// ----------------------------
-short TrackVis::read( mitk::FiberBundleX* fib )
+short TrackVisFiberReader::read( mitk::FiberBundleX* fib )
 {
     int numPoints;
     vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
     vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
 
     while (fread((char*)&numPoints, 1, 4, fp)==4)
     {
-        if ( numPoints >= maxSteps || numPoints <= 0 )
+        if ( numPoints <= 0 )
         {
             printf( "[ERROR] Trying to read a fiber with %d points!\n", numPoints );
             return -1;
         }
         vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
 
         float tmp[3];
         for(int i=0; i<numPoints; i++)
         {
             if (fread((char*)tmp, 1, 12, fp) == 0)
                 MITK_ERROR << "TrackVis::read: Error during read.";
 
+//            tmp[0] *= hdr.voxel_size[0];
+//            tmp[1] *= hdr.voxel_size[1];
+//            tmp[2] *= hdr.voxel_size[2];
+
+//            tmp[0] += hdr.origin[0];
+//            tmp[1] += hdr.origin[1];
+//            tmp[2] += hdr.origin[2];
+
             vtkIdType id = vtkNewPoints->InsertNextPoint(tmp);
             container->GetPointIds()->InsertNextId(id);
         }
         vtkNewCells->InsertNextCell(container);
     }
 
     vtkSmartPointer<vtkPolyData> fiberPolyData = vtkSmartPointer<vtkPolyData>::New();
     fiberPolyData->SetPoints(vtkNewPoints);
     fiberPolyData->SetLines(vtkNewCells);
     fib->SetFiberPolyData(fiberPolyData);
 
     return numPoints;
 }
 
 
 
 // Update the field in the header to the new FIBER TOTAL.
 // ------------------------------------------------------
-void TrackVis::updateTotal( int totFibers )
+void TrackVisFiberReader::updateTotal( int totFibers )
 {
     fseek(fp, 1000-12, SEEK_SET);
     if (fwrite((char*)&totFibers, 1, 4, fp) != 4)
         MITK_ERROR << "[ERROR] Problems saving the fiber!";
 }
 
 
-void TrackVis::writeHdr()
+void TrackVisFiberReader::writeHdr()
 {
     fseek(fp, 0, SEEK_SET);
     if (fwrite((char*)&hdr, 1, 1000, fp) != 1000)
         MITK_ERROR << "[ERROR] Problems saving the fiber!";
 }
 
 
 // Close the TrackVis file, but keep the metadata in the header.
 // -------------------------------------------------------------
-void TrackVis::close()
+void TrackVisFiberReader::close()
 {
     fclose(fp);
     fp = NULL;
 }
diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkTrackvis.h b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkTrackvis.h
index 56d0111c19..2bca7b9198 100644
--- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkTrackvis.h
+++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkTrackvis.h
@@ -1,83 +1,82 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 #ifndef _TRACKVIS
 #define _TRACKVIS
 
 #include <mitkCommon.h>
 #include <mitkFiberBundleX.h>
 #include <vtkSmartPointer.h>
 #include <vtkPolyData.h>
 #include <vtkCellArray.h>
 #include <vtkPoints.h>
 #include <vtkPolyLine.h>
 #include <itkSize.h>
 
 using namespace std;
 
 // Structure to hold metadata of a TrackVis file
 // ---------------------------------------------
 struct TrackVis_header
 {
     char                id_string[6];
     short int           dim[3];
     float               voxel_size[3];
     float               origin[3];
     short int           n_scalars;
     char                scalar_name[10][20];
     short int           n_properties;
     char                property_name[10][20];
     char                reserved[508];
     char                voxel_order[4];
     char                pad2[4];
     float               image_orientation_patient[6];
     char                pad1[2];
     unsigned char       invert_x;
     unsigned char       invert_y;
     unsigned char       invert_z;
     unsigned char       swap_xy;
     unsigned char       swap_yz;
     unsigned char       swap_zx;
     int                 n_count;
     int                 version;
     int                 hdr_size;
 };
 
 // Class to handle TrackVis files.
 // -------------------------------
-class MitkFiberTracking_EXPORT TrackVis
+class MitkFiberTracking_EXPORT TrackVisFiberReader
 {
 private:
     string              filename;
     FILE*               fp;
-    int                 maxSteps;  // [TODO] should be related to the variable defined for fiber-tracking
 
 public:
     TrackVis_header     hdr;
     float               m_Origin[3];
 
-    short     create(string filename, mitk::FiberBundleX* fib);
-    short     open( string filename );
-    short     read( mitk::FiberBundleX* fib );
-    short    append( mitk::FiberBundleX* fib );
+    short   create(string filename, mitk::FiberBundleX* fib);
+    short   open( string filename );
+    short   read( mitk::FiberBundleX* fib );
+    short   append( mitk::FiberBundleX* fib );
     void    writeHdr();
     void    updateTotal( int totFibers );
     void    close();
 
-    TrackVis();
-    ~TrackVis();
+    TrackVisFiberReader();
+    ~TrackVisFiberReader();
 };
 
 #endif
diff --git a/Modules/DiffusionImaging/MiniApps/CMakeLists.txt b/Modules/DiffusionImaging/MiniApps/CMakeLists.txt
index 6a26aa8477..5f79886390 100755
--- a/Modules/DiffusionImaging/MiniApps/CMakeLists.txt
+++ b/Modules/DiffusionImaging/MiniApps/CMakeLists.txt
@@ -1,65 +1,64 @@
 option(BUILD_DiffusionMiniApps "Build commandline tools for diffusion" OFF)
 
 if(BUILD_DiffusionMiniApps OR MITK_BUILD_ALL_APPS)
 
   # needed include directories
   include_directories(
     ${CMAKE_CURRENT_SOURCE_DIR}
     ${CMAKE_CURRENT_BINARY_DIR}
     )
 
   project( MitkDiffusionMiniApps )
 
   # fill in the standalone executables here
   set(DIFFUSIONMINIAPPS
     mitkDiffusionMiniApps
   )
 
   # set additional files here
   set(DIFFUSIONCORE_ADDITIONAL_FILES
     MiniAppManager.cpp
     FileFormatConverter.cpp
     TensorReconstruction.cpp
     QballReconstruction.cpp
     DiffusionIndices.cpp
     CopyGeometry.cpp
     GibbsTracking.cpp
     StreamlineTracking.cpp
     FiberProcessing.cpp
     LocalDirectionalFiberPlausibility.cpp
     #TractogramAngularError.cpp
     FiberDirectionExtraction.cpp
     PeakExtraction.cpp
     PeaksAngularError.cpp
     MultishellMethods.cpp
     Fiberfox.cpp
     ExportShImage.cpp
-    ImportShImage.cpp
     NetworkCreation.cpp
     NetworkStatistics.cpp
     DwiDenoising.cpp
   )
 
   # deprecated
 #  FOREACH(tool ${DIFFUSIONMINIAPPS})
 #    ADD_EXECUTABLE(
 #    ${tool}
 #    ${tool}.cpp
 #    ${DIFFUSIONCORE_ADDITIONAL_FILES}
 #    )
 
 #    TARGET_LINK_LIBRARIES(
 #      ${tool}
 #      ${ALL_LIBRARIES} )
 #  ENDFOREACH(tool)
 
   mitk_create_executable(DiffusionMiniApps
     DEPENDS MitkDiffusionCore MitkFiberTracking MitkConnectomics
     PACKAGE_DEPENDS ITK|ITKDiffusionTensorImage
   )
 
   if(EXECUTABLE_IS_ENABLED)
     MITK_INSTALL_TARGETS(EXECUTABLES ${EXECUTABLE_TARGET})
   endif()
 
 endif()
diff --git a/Modules/DiffusionImaging/MiniApps/files.cmake b/Modules/DiffusionImaging/MiniApps/files.cmake
index a2bcaa1eba..98428da34c 100644
--- a/Modules/DiffusionImaging/MiniApps/files.cmake
+++ b/Modules/DiffusionImaging/MiniApps/files.cmake
@@ -1,32 +1,31 @@
 set(CPP_FILES
     mitkDiffusionMiniApps.cpp
     MiniAppManager.cpp
 
     BatchedFolderRegistration.cpp
     DicomFolderDump.cpp
     FileFormatConverter.cpp
     TensorReconstruction.cpp
     TensorDerivedMapsExtraction.cpp
     QballReconstruction.cpp
     DiffusionIndices.cpp
     ExtractImageStatistics.cpp
     CopyGeometry.cpp
     GibbsTracking.cpp
     StreamlineTracking.cpp
     FiberProcessing.cpp
     LocalDirectionalFiberPlausibility.cpp
     #TractogramAngularError.cpp
     FiberDirectionExtraction.cpp
     ImageResampler.cpp
     PeakExtraction.cpp
     PeaksAngularError.cpp
     MultishellMethods.cpp
     Fiberfox.cpp
     ExportShImage.cpp
-    ImportShImage.cpp
     NetworkCreation.cpp
     NetworkStatistics.cpp
     DwiDenoising.cpp
     FiberExtraction.cpp
     FiberJoin.cpp
 )