diff --git a/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.cpp b/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.cpp index 35d2f71350..a9129781e9 100644 --- a/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.cpp +++ b/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.cpp @@ -1,882 +1,214 @@ /*========================================================================= 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 #include #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include #if BOOST_VERSION / 100000 > 0 #if BOOST_VERSION / 100 % 1000 > 34 #include #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 - ::DiffusionMultiShellQballReconstructionImageFilter() : +template< class T, class TG, class TO, int L, int NODF> +DiffusionMultiShellQballReconstructionImageFilter +::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 - - ::Normalize( OdfPixelType odf, - typename NumericTraits::AccumulateType b0 ) - { - switch( m_NormalizationMethod ) - { - case QBAR_STANDARD: - { - TOdfPixelType sum = 0; - - for(int i=0; i0) - odf /= sum; - - return odf; - break; - } - case QBAR_B_ZERO_B_VALUE: - { - for(int i=0; i0) - odf /= sum; - - break; - } - case QBAR_NONNEG_SOLID_ANGLE: - { - break; - } - } - - return odf; - } - - - template< - class TReferenceImagePixelType, - class TGradientImagePixelType, - class TOdfPixelType, - int NOrderL, - int NrOdfDirections> - vnl_vector - itk::DiffusionMultiShellQballReconstructionImageFilter - - ::PreNormalize( vnl_vector vec, - typename NumericTraits::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= 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 - ::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 - ::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 baselineind; // contains the indicies of - // the baseline images - std::vector 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::AccumulateType b0 = NumericTraits::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 B(m_NumberOfGradientDirections); - - if( (b0 != 0) && (b0 >= m_Threshold) ) - { - for( unsigned int i = 0; i< m_NumberOfGradientDirections; i++ ) - { - B[i] = static_cast(b[gradientind[i]]); - } - - B = PreNormalize(B, b0); - if(m_NormalizationMethod == QBAR_SOLID_ANGLE) - { - vnl_vector 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 - void DiffusionMultiShellQballReconstructionImageFilter - ::tofile2(vnl_matrix *pA, std::string fname) - { - vnl_matrix 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 - double DiffusionMultiShellQballReconstructionImageFilter - ::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 - ::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 - ::legendre0(int l) - { - if( l%2 != 0 ) - { - return 0; - } - else - { - double prod1 = 1.0; - for(int i=1;i - double DiffusionMultiShellQballReconstructionImageFilter - ::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; - } +template< class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NOrderL, int NrOdfDirections> +typename itk::DiffusionMultiShellQballReconstructionImageFilter< TReferenceImagePixelType,TGradientImagePixelType,TOdfPixelType, NOrderL,NrOdfDirections>::OdfPixelType itk::DiffusionMultiShellQballReconstructionImageFilter +::Normalize( OdfPixelType odf, typename NumericTraits::AccumulateType b0 ) +{ - 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(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(l, m, cos(theta)); -#endif -#endif - } - //std::cout << "legendre(" << l << "," << m << "," << cos(theta) << ") = " << pml << std::endl; +template +vnl_vectoritk::DiffusionMultiShellQballReconstructionImageFilter +::PreNormalize( vnl_vector vec, typename NumericTraits::AccumulateType b0 ) +{ + // Threshold + double sigma = 0.001; - 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 + int n = vec.size(); + double b0f = (double)b0; + for(int i=0; i - double DiffusionMultiShellQballReconstructionImageFilter - ::Yj(int m, int k, double theta, double phi) - { - if( -k <= m && m < 0 ) + if(0 <= vec[i] && vec[i] < sigma ) { - return sqrt(2.0) * spherical_harmonic(m,k,theta,phi,false); + veci[i] = sigma / 2 + (vec[i] * vec[i]) / 2 * sigma; } - - 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 - ::ComputeReconstructionMatrix() - { - - //for(int i=-6;i<7;i++) - // std::cout << boost::math::legendre_p(6, i, 0.65657) << std::endl; - - if( m_NumberOfGradientDirections < 6 ) + if( 1 - sigma <= vec[i] && vec[i] < 1 ) { - itkExceptionMacro( << "Not enough gradient directions supplied. Need to supply at least 6" ); + vec[i] = 1-(sigma/2) - (( 1 - vec[i] * vec[i] ) / 2 * sigma ); } - + if( 1 <= vec[i] ) { - // 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 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; - } + vec[i] = 1 - (sigma / 2); } - vnl_matrix *Q - = new vnl_matrix(3, m_NumberOfGradientDirections); + if (b0f==0) + b0f = 0.01; + if(vec[i] >= b0f) + vec[i] = b0f - 0.001; - { - 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* B = new vnl_matrix(m_NumberOfGradientDirections,m_NumberCoefficients); - vnl_matrix* _L = new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients); - _L->fill(0.0); - vnl_matrix* LL = new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients); - LL->fill(0.0); - vnl_matrix* P = new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients); - P->fill(0.0); - vnl_matrix* Inv = new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients); - P->fill(0.0); - vnl_vector* lj = new vnl_vector(m_NumberCoefficients); - m_LP = new vnl_vector(m_NumberCoefficients); - - for(unsigned int i=0; i +void DiffusionMultiShellQballReconstructionImageFilter +::SetGradientImage( GradientDirectionContainerType *gradientDirection, const GradientImagesType *gradientImage ) +{ + this->m_GradientDirectionContainer = gradientDirection; + this->m_NumberOfBaselineImages = 0; - for(int i=0; i(B->transpose()); - //tofile2(&m_B_t,"m_B_t"); - vnl_matrix B_t_B = (*m_B_t) * (*B); - //tofile2(&B_t_B,"B_t_B"); - vnl_matrix lambdaLL(m_NumberCoefficients,m_NumberCoefficients); - lambdaLL.update((*LL)); - lambdaLL *= m_Lambda; - //tofile2(&lambdaLL,"lLL"); - - vnl_matrix tmp( B_t_B + lambdaLL); - vnl_matrix_inverse *pseudoInverse - = new vnl_matrix_inverse( tmp ); - - (*Inv) = pseudoInverse->pinverse(); - //tofile2(Inv,"Inv"); - vnl_matrix temp((*Inv) * (*m_B_t)); - double fac1 = (1.0/(16.0*QBALL_ANAL_RECON_PI*QBALL_ANAL_RECON_PI)); - switch(m_NormalizationMethod) + for(GradientDirectionContainerType::Iterator it = gradientDirection->Begin(); it != gradientDirection->End(); it++) + { + if( it.Value().one_norm() <= 0.0 ) { - 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; + this->m_NumberOfBaselineImages ++; + }else{ + it.Value() = it.Value() / it.Value().two_norm(); } + } - //tofile2(&temp,"A"); +} - m_CoeffReconstructionMatrix = new vnl_matrix(m_NumberCoefficients,m_NumberOfGradientDirections); - for(int i=0; i +void DiffusionMultiShellQballReconstructionImageFilter +::BeforeThreadedGenerateData() +{ + if( m_NumberOfGradientDirections < 6 ) + { + itkExceptionMacro( << "At least 6 gradient directions are required" ); + } +} - // this code goes to the image adapter coeffs->odfs later +template< class T, class TG, class TO, int L, int NODF> +void DiffusionMultiShellQballReconstructionImageFilter +::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int NumberOfThreads) +{ - int NOdfDirections = NODF; - vnl_matrix_fixed* U = - itk::PointShell >::DistributePointShell(); +} - m_SphericalHarmonicBasisMatrix = new vnl_matrix(NOdfDirections,m_NumberCoefficients); - vnl_matrix* sphericalHarmonicBasisMatrix2 - = new vnl_matrix(NOdfDirections,m_NumberCoefficients); - for(int i=0; i +void DiffusionMultiShellQballReconstructionImageFilter +::ComputeReconstructionMatrix() +{ - for(int i=0; i(NOdfDirections,m_NumberOfGradientDirections); - *m_ReconstructionMatrix = (*m_SphericalHarmonicBasisMatrix) * (*m_CoeffReconstructionMatrix); +} - } +template< class T, class TG, class TO, int L, int NODF> +void DiffusionMultiShellQballReconstructionImageFilter +::PrintSelf(std::ostream& os, Indent indent) const +{ + std::locale C("C"); + std::locale originalLocale = os.getloc(); + os.imbue(C); - template< class T, class TG, class TO, int L, int NODF> - void DiffusionMultiShellQballReconstructionImageFilter - ::SetGradientImage( GradientDirectionContainerType *gradientDirection, - const GradientImagesType *gradientImage ) - { - this->m_GradientDirectionContainer = gradientDirection; + Superclass::PrintSelf(os,indent); - 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(); - } - } + 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; - this->m_NumberOfGradientDirections = numImages - this->m_NumberOfBaselineImages; + os.imbue( originalLocale ); +} - // 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 ) +template< class T, class TG, class TO, int L, int NODF> +void DiffusionMultiShellQballReconstructionImageFilter +::tofile2(vnl_matrix *pA, std::string fname) +{ + vnl_matrix 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; im_NumberOfBaselineImages - << "baselines = " << m_NumberOfGradientDirections + this->m_NumberOfBaselineImages - << " directions specified but image has " << gradientImage->GetVectorLength() - << " components."); + myfile << A(i,j) << " "; + if(j==A.columns()-1 && i!=A.rows()-1) + myfile << ";"; } + } + myfile << "];"; + myfile.close(); - this->ProcessObject::SetNthInput( 0, - const_cast< GradientImagesType* >(gradientImage) ); - - } - - template< class T, class TG, class TO, int L, int NODF> - void DiffusionMultiShellQballReconstructionImageFilter - ::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; + myfile.imbue( originalLocale ); +} - os.imbue( originalLocale ); - } } #endif // __itkDiffusionMultiShellQballReconstructionImageFilter_cpp diff --git a/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h b/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h index 416a814de9..42fe6d40ce 100644 --- a/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h +++ b/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h @@ -1,300 +1,173 @@ /*========================================================================= 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 - * NRRD - * 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[1] - * Tuch DS, - * "Q-ball imaging", Magn Reson Med. 2004 Dec;52(6):1358-72. - * + Aganj_2010 */ -template< class TReferenceImagePixelType, - class TGradientImagePixelType, - class TOdfPixelType, - int NOrderL, - int NrOdfDirections> -class DiffusionMultiShellQballReconstructionImageFilter : - public ImageToImageFilter< Image< TReferenceImagePixelType, 3 >, - Image< Vector< TOdfPixelType, NrOdfDirections >, 3 > > +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 Pointer; - typedef SmartPointer 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 + typedef DiffusionMultiShellQballReconstructionImageFilter Self; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + typedef ImageToImageFilter< Image< TReferenceImagePixelType, 3>, Image< Vector< TOdfPixelType, NrOdfDirections >, 3 > > Superclass; + typedef TReferenceImagePixelType ReferencePixelType; + typedef TGradientImagePixelType GradientPixelType; + typedef Vector< TOdfPixelType, NrOdfDirections > OdfPixelType; + typedef typename Superclass::InputImageType ReferenceImageType; + typedef Image< OdfPixelType, 3 > OdfImageType; + typedef OdfImageType OutputImageType; + typedef TOdfPixelType BZeroPixelType; + typedef Image< BZeroPixelType, 3 > BZeroImageType; + typedef typename Superclass::OutputImageRegionType OutputImageRegionType; + /** Typedef defining one (of the many) gradient images. */ + typedef Image< GradientPixelType, 3 > GradientImageType; + /** An alternative typedef defining one (of the many) gradient images. + * It will be assumed that the vectorImage has the same dimension as the * Reference image and a vector length parameter of \c n (number of * gradient directions)*/ - typedef VectorImage< GradientPixelType, 3 > GradientImagesType; - - /** Holds the ODF reconstruction matrix */ - typedef vnl_matrix< TOdfPixelType >* - OdfReconstructionMatrixType; - - typedef vnl_matrix< double > CoefficientMatrixType; - - /** Holds each magnetic field gradient used to acquire one DWImage */ - typedef vnl_vector_fixed< double, 3 > GradientDirectionType; - - /** Container to hold gradient directions of the 'n' DW measurements */ - typedef VectorContainer< unsigned int, - GradientDirectionType > GradientDirectionContainerType; - - /** 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 + typedef VectorImage< GradientPixelType, 3 > GradientImagesType; + /** Holds the ODF reconstruction matrix */ + typedef vnl_matrix< TOdfPixelType >* OdfReconstructionMatrixType; + typedef vnl_matrix< double > CoefficientMatrixType; + /** Holds each magnetic field gradient used to acquire one DWImage */ + typedef vnl_vector_fixed< double, 3 > GradientDirectionType; + /** Container to hold gradient directions of the 'n' DW measurements */ + typedef VectorContainer< unsigned int, GradientDirectionType > GradientDirectionContainerType; + typedef Image BlaImage; + + // --------------------------------------------------------------------------------------------// + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Runtime information support. */ + itkTypeMacro(DiffusionMultiShellQballReconstructionImageFilter, ImageToImageFilter); + + /** set method to add gradient directions and its corresponding + * image. The image here is a VectorImage. The user is expected to pass the + * gradient directions in a container. The ith element of the container + * corresponds to the gradient direction of the ith component image the * VectorImage. For the baseline image, a vector of all zeros * should be set.*/ - void SetGradientImage( GradientDirectionContainerType *, - const GradientImagesType *image); - - /** 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 + 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 ); + if( idx >= m_NumberOfGradientDirections ) + { + itkExceptionMacro( << "Gradient direction " << idx << "does not exist" ); + } + return m_GradientDirectionContainer->ElementAt( idx+1 ); } - static void tofile2(vnl_matrix *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::AccumulateType b0 ); - OdfPixelType Normalize(OdfPixelType odf, typename NumericTraits::AccumulateType b0 ); + vnl_vector PreNormalize( vnl_vector vec, typename NumericTraits::AccumulateType b0 ); - vnl_vector PreNormalize( vnl_vector vec, typename NumericTraits::AccumulateType b0 ); - - /** Threshold on the reference image data. The output ODF will be a null + /** Threshold on the reference image data. The output ODF will be a null * pdf for pixels in the reference image that have a value less than this * threshold. */ - itkSetMacro( Threshold, ReferencePixelType ); - itkGetMacro( Threshold, ReferencePixelType ); + itkSetMacro( Threshold, ReferencePixelType ); + itkGetMacro( Threshold, ReferencePixelType ); - itkSetMacro( NormalizationMethod, Normalization); - itkGetMacro( NormalizationMethod, Normalization ); + itkSetMacro( NormalizationMethod, Normalization); + itkGetMacro( NormalizationMethod, Normalization ); - typedef Image BlaImage; - itkGetMacro( BZeroImage, typename BZeroImageType::Pointer); - itkGetMacro( ODFSumImage, typename BlaImage::Pointer); + itkGetMacro( BZeroImage, typename BZeroImageType::Pointer); + itkGetMacro( ODFSumImage, typename BlaImage::Pointer); + + itkSetMacro( BValue, TOdfPixelType); + + itkSetMacro( Lambda, double ); + itkGetMacro( Lambda, double ); + + itkGetConstReferenceMacro( BValue, TOdfPixelType); - 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)); - itkConceptMacro(TensorEqualityComparableCheck, - (Concept::EqualityComparable)); - itkConceptMacro(GradientConvertibleToDoubleCheck, - (Concept::Convertible)); - itkConceptMacro(DoubleConvertibleToTensorCheck, - (Concept::Convertible)); - itkConceptMacro(GradientReferenceAdditiveOperatorsCheck, - (Concept::AdditiveOperators)); - itkConceptMacro(ReferenceOStreamWritableCheck, - (Concept::OStreamWritable)); - itkConceptMacro(TensorOStreamWritableCheck, - (Concept::OStreamWritable)); - /** 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; + DiffusionMultiShellQballReconstructionImageFilter(); + ~DiffusionMultiShellQballReconstructionImageFilter() {}; + void PrintSelf(std::ostream& os, Indent indent) const; + + void ComputeReconstructionMatrix(); - OdfReconstructionMatrixType m_SphericalHarmonicBasisMatrix; - /** container to hold gradient directions */ - GradientDirectionContainerType::Pointer m_GradientDirectionContainer; + void BeforeThreadedGenerateData(); + void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, int NumberOfThreads ); + +private: + + OdfReconstructionMatrixType m_ReconstructionMatrix; + OdfReconstructionMatrixType m_CoeffReconstructionMatrix; + OdfReconstructionMatrixType m_SphericalHarmonicBasisMatrix; - /** Number of gradient measurements */ - unsigned int m_NumberOfGradientDirections; + /** container to hold gradient directions */ + GradientDirectionContainerType::Pointer m_GradientDirectionContainer; - /** Number of baseline images */ - unsigned int m_NumberOfBaselineImages; + /** Number of gradient measurements */ + unsigned int m_NumberOfGradientDirections; - /** Threshold on the reference image data */ - ReferencePixelType m_Threshold; + /** Number of baseline images */ + unsigned int m_NumberOfBaselineImages; - /** LeBihan's b-value for normalizing tensors */ - TOdfPixelType m_BValue; + /** Threshold on the reference image data */ + ReferencePixelType m_Threshold; - typename BZeroImageType::Pointer m_BZeroImage; + /** LeBihan's b-value for normalizing tensors */ + TOdfPixelType m_BValue; - double m_Lambda; + typename BZeroImageType::Pointer m_BZeroImage; - bool m_DirectionsDuplicated; + double m_Lambda; - Normalization m_NormalizationMethod; + bool m_DirectionsDuplicated; - int m_NumberCoefficients; + int m_NumberCoefficients; - vnl_matrix* m_B_t; - vnl_vector* m_LP; - BlaImage::Pointer m_ODFSumImage; + BlaImage::Pointer m_ODFSumImage; }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkDiffusionMultiShellQballReconstructionImageFilter.cpp" #endif #endif //__itkDiffusionMultiShellQballReconstructionImageFilter_h_