diff --git a/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.cpp b/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.cpp index d39c168ca8..6cccb352e4 100644 --- a/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.cpp +++ b/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.cpp @@ -1,887 +1,1060 @@ /*========================================================================= 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 #include #define _USE_MATH_DEFINES #include #include "mitkSphericalHarmonicsFunctions.h" #include "itkPointShell.h" #include namespace itk { #define QBALL_ANAL_RECON_PI M_PI 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_IsHemisphericalArrangementOfGradientDirections(false), - m_IsArithmeticProgession(false), - m_ReconstructionType(Standard1Shell) + m_GradientDirectionContainer(NULL), + m_NumberOfGradientDirections(0), + m_NumberOfBaselineImages(1), + m_Threshold(NumericTraits< ReferencePixelType >::NonpositiveMin()), + m_BValue(1.0), + m_Lambda(0.0), + m_IsHemisphericalArrangementOfGradientDirections(false), + m_IsArithmeticProgession(false), + m_ReconstructionType(Standard1Shell) { - // At least 1 inputs is necessary for a vector image. - // For images added one at a time we need at least six - this->SetNumberOfRequiredInputs( 1 ); + // 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 ) { - for(int i=0; i0) + odf /= sum; + + return odf; +} + +template +void itk::DiffusionMultiShellQballReconstructionImageFilter +::Threshold(vnl_vector & vec, float sigma) +{ + for(int i = 0; i < vec.size(); i++) + { + if(vec[i] < 0) { - odf[i] = odf[i] < 0 ? 0 : odf[i]; - odf[i] *= QBALL_ANAL_RECON_PI*4/NrOdfDirections; + vec[i] = sigma / 2 ; } - TOdfPixelType sum = 0; - for(int i=0; i0) - odf /= sum; - - return odf; + if( 1 - sigma <= vec[i] && vec[i] < 1 ) + { + vec[i] = 1-(sigma/2) - (( 1 - vec[i] * vec[i] ) / 2 * sigma ); + } + if( 1 <= vec[i] ) + { + vec[i] = 1 - (sigma / 2); + } + } } template void itk::DiffusionMultiShellQballReconstructionImageFilter -::Threshold(vnl_vector & vec, float sigma) +::Threshold(vnl_matrix & mat, float sigma) { - for(int i = 0; i < vec.size(); i++) - { - if(vec[i] < 0) - { - vec[i] = sigma / 2 ; - } - if(0 <= vec[i] && vec[i] < sigma ) - { - vec[i] = sigma / 2 + (vec[i] * vec[i]) / 2 * sigma; - } - if( 1 - sigma <= vec[i] && vec[i] < 1 ) - { - vec[i] = 1-(sigma/2) - (( 1 - vec[i] * vec[i] ) / 2 * sigma ); - } - if( 1 <= vec[i] ) - { - vec[i] = 1 - (sigma / 2); - } + for(int i = 0; i < mat.rows(); i++) + { + for( int j = 0; j < mat.cols(); j++ ){ + + if(mat(i,j) < 0) + { + mat(i,j) = sigma / 2 ; + } + if(0 <= mat(i,j) && mat(i,j) < sigma ) + { + mat(i,j) = sigma / 2 + (mat(i,j) * mat(i,j)) / 2 * sigma; + } + if( 1 - sigma <= mat(i,j) && mat(i,j) < 1 ) + { + mat(i,j) = 1-(sigma/2) - (( 1 - mat(i,j) * mat(i,j) ) / 2 * sigma ); + } + if( 1 <= mat(i,j) ) + { + mat(i,j) = 1 - (sigma / 2); + } } + } } template void itk::DiffusionMultiShellQballReconstructionImageFilter -::PreNormalize( vnl_vector & vec, typename NumericTraits::AccumulateType b0 ) +::Projection1( vnl_matrix & E, float delta ) { - double b0f = (double)b0; - for(int i = 0; i < vec.size(); i++) + double sF = sqrt(5.0); + + vnl_vector vOnes(E.rows()); + vOnes.fill(1.0); + + vnl_matrix T0(E.rows(), E.cols()); + vnl_matrix C(E.rows(), 7); + vnl_matrix A(E.rows(), 7); + vnl_matrix B(E.rows(), 7); + + vnl_vector s0(E.rows()); + vnl_vector a0(E.rows()); + vnl_vector b0(E.rows()); + vnl_vector ta(E.rows()); + vnl_vector tb(E.rows()); + vnl_vector e(E.rows()); + vnl_vector m(E.rows()); + vnl_vector a(E.rows()); + vnl_vector b(E.rows()); + + + // logarithmierung aller werte in E + for(int i = 0 ; i < E.rows(); i++) + { + for(int j = 0 ; j < E.cols(); j++) + { + T0(i,j) = -log(E(i,j)); + } + } + + // Summeiere Zeilenweise über alle Shells sum = E1+E2+E3 + for(int i = 0 ; i < E.rows(); i++) + { + s0[i] = T0(i,0) + T0(i,1) + T0(i,2); + } + + + for(int i = 0; i < E.rows(); i ++) + { + // Alle Signal-Werte auf der Ersten shell E(N,0) normiert auf s0 + a0 = E(i,0) / s0[i]; + // Alle Signal-Werte auf der Zweiten shell E(N,1) normiert auf s0 + b0 = E(i,1) / s0[i]; + } + + ta = 3*a0; + tb = 3*b0; + e = tb - (2*ta); + m = (2*tb) + ta; + + for(int i = 0; i < E.rows(); i++) + { + C(i,0) = tb[i] < 1+3*delta && 0.5+1.5*(sF+1)*delta < ta[i] && ta[i] < 1-3* (sF+2) *delta; + C(i,1) = e[i] <= -1 +3*(2*sF+5)* delta && ta[i] >= 1-3*(sF+2)*delta; + C(i,2) = m[i] > 3 -3*sF*delta && -1+3*(2*sF+5)*d < e[i] && e<-3*sF*d; + C(i,3) = m[i] >= 3-3*sF*d && e >= -3 *sF * delta; + C(i,4) = 2.5 + 1.5*(5+sF)*delta < m[i] && m < 3-3*sF*delta && e> -3*sF*delta; + C(i,5) = ta[i] <= 0.5+1.5 *delta(sF+1)*delta && m[i] <= 2.5 + 1.5 *(5+sF) * delta; + C(i,6) = !( C(i,0) || C(i,1) || C(i,2) || C(i,3) || C(i,4) || C(i,5) ); // ~ANY(C(i,[0-5] ),2) + } + + A.set_column(0, a0); + A.set_column(1, (1/3 - (sF+2) * delta ) * vOnes); + A.set_column(2, 0.2+0.8 * a0 - 0.4 * b0 -d/sF); + A.set_column(3, (0.2 + delta /sF)*vOnes); + A.set_column(4, 0.2 * a0 + 0.4 * b0 + 2*d/sF); + A.set_column(5, (1/6+0.5*(sF+1)*delat) * vOnes); + A.set_column(6, a0); + + + B.set_column(0, (1/3 +delta) * vOnes ); + B.set_column(1, (1/3 +delta) * vOnes ); + B.set_column(2, 0.4 * 0.4 * a0 + 0.2 * b0 - 2*d/sF ); + B.set_column(3, (0.4 - 3* delta / sF) * vOnes ); + B.set_column(4, 0.4 * a0 + 0.8 * b0 - delta /sF); + B.set_column(5, (1/3+d) * vOnes); + B.set_column(6, b0 ); + + for(int i = 0 ; i < E.rows(); i++) + { + double sumA = 0; + double sumB = 0; + for(int j = 0 ; j < 7; j++) { - if (b0f==0) - b0f = 0.01; - if(vec[i] >= b0f) - vec[i] = b0f - 0.001; - vec[i] /= b0f; + if(C(i,j)) + { + sumA += A(i,j); + sumB += B(i,j); + } } + a[i] = sumA; + b[i] = sumB; + } + + for(int i = 0; i < E.rows(); i++) + { + E(i,0) = exp(-(a[i]*s0[i])); + E(i,1) = exp(-(b[i]*s0[i])); + E(i,2) = exp(-((1-a[i]-b[i])*s0[i])); + } +} + +template +void itk::DiffusionMultiShellQballReconstructionImageFilter +::S_S0Normalization( vnl_vector & vec, typename NumericTraits::AccumulateType b0 ) +{ + + double b0f = (double)b0; + for(int i = 0; i < vec.size(); i++) + { + if (b0f==0) + b0f = 0.01; + //if(vec[i] >= b0f) + // vec[i] = b0f - 0.001; + vec[i] /= b0f; + } } +template +void itk::DiffusionMultiShellQballReconstructionImageFilter +::S_S0Normalization( vnl_matrix & mat, typename NumericTraits::AccumulateType b0 = 0 ) +{ + double b0f = (double)b0; + for(int i = 0; i < mat.rows(); i++) + { + for( int j = 0; j < mat.cols(); j++ ){ + if (b0f==0) + b0f = 0.01; + //if(mat(i,j) >= b0f) + // mat(i,j) = b0f - 0.001; + mat(i,j) /= b0f; + } + } +} + + + template< class T, class TG, class TO, int L, int NODF> void DiffusionMultiShellQballReconstructionImageFilter ::DoubleLogarithm(vnl_vector & vec) { - for(int i = 0; i < vec.size(); i++) - { - vec[i] = log(-log(vec[i])); - } + for(int i = 0; i < vec.size(); i++) + { + vec[i] = log(-log(vec[i])); + } } template< class T, class TG, class TO, int L, int NODF> void DiffusionMultiShellQballReconstructionImageFilter ::SetGradientImage( GradientDirectionContainerType *gradientDirection, const GradientImagesType *gradientImage ) { - this->m_GradientDirectionContainer = gradientDirection; - this->m_NumberOfBaselineImages = 0; - - for(GradientDirectionContainerType::Iterator it = gradientDirection->Begin(); it != gradientDirection->End(); it++) - { - if( it.Value().one_norm() <= 0.0 ) - this->m_NumberOfBaselineImages ++; - else - it.Value() = it.Value() / it.Value().two_norm(); - } - - this->m_NumberOfGradientDirections = gradientDirection->Size() - this->m_NumberOfBaselineImages; - - // ensure that the gradient image we received has as many components as - // the number of gradient directions - if( gradientImage->GetVectorLength() != this->m_NumberOfBaselineImages + m_NumberOfGradientDirections ) - { - itkExceptionMacro( << m_NumberOfGradientDirections << " gradients + " << this->m_NumberOfBaselineImages - << "baselines = " << m_NumberOfGradientDirections + this->m_NumberOfBaselineImages - << " directions specified but image has " << gradientImage->GetVectorLength() - << " components."); - } + this->m_GradientDirectionContainer = gradientDirection; + this->m_NumberOfBaselineImages = 0; - this->ProcessObject::SetNthInput( 0, const_cast< GradientImagesType* >(gradientImage) ); + for(GradientDirectionContainerType::Iterator it = gradientDirection->Begin(); it != gradientDirection->End(); it++) + { + if( it.Value().one_norm() <= 0.0 ) + this->m_NumberOfBaselineImages ++; + else + it.Value() = it.Value() / it.Value().two_norm(); + } + + this->m_NumberOfGradientDirections = gradientDirection->Size() - this->m_NumberOfBaselineImages; + + // ensure that the gradient image we received has as many components as + // the number of gradient directions + if( gradientImage->GetVectorLength() != this->m_NumberOfBaselineImages + m_NumberOfGradientDirections ) + { + itkExceptionMacro( << m_NumberOfGradientDirections << " gradients + " << this->m_NumberOfBaselineImages + << "baselines = " << m_NumberOfGradientDirections + this->m_NumberOfBaselineImages + << " directions specified but image has " << gradientImage->GetVectorLength() + << " components."); + } + + this->ProcessObject::SetNthInput( 0, const_cast< GradientImagesType* >(gradientImage) ); } template< class T, class TG, class TO, int L, int NODF> void DiffusionMultiShellQballReconstructionImageFilter ::BeforeThreadedGenerateData() { - if( m_NumberOfGradientDirections < (((L+1)*(L+2))/2) /* && m_NumberOfGradientDirections < 6 */ ) + if( m_NumberOfGradientDirections < (((L+1)*(L+2))/2) /* && m_NumberOfGradientDirections < 6 */ ) + { + itkExceptionMacro( << "At least " << ((L+1)*(L+2))/2 << " gradient directions are required" ); + } + // Input must be an itk::VectorImage. + std::string gradientImageClassName(this->ProcessObject::GetInput(0)->GetNameOfClass()); + + if ( strcmp(gradientImageClassName.c_str(),"VectorImage") != 0 ) + itkExceptionMacro( << "There is only one Gradient image. I expect that to be a VectorImage. But its of type: " << gradientImageClassName ); + + m_BZeroImage = BZeroImageType::New(); + typename GradientImagesType::Pointer img = static_cast< GradientImagesType * >( this->ProcessObject::GetInput(0) ); + m_BZeroImage->SetSpacing( img->GetSpacing() ); // Set the image spacing + m_BZeroImage->SetOrigin( img->GetOrigin() ); // Set the image origin + m_BZeroImage->SetDirection( img->GetDirection() ); // Set the image direction + m_BZeroImage->SetLargestPossibleRegion( img->GetLargestPossibleRegion()); + m_BZeroImage->SetBufferedRegion( img->GetLargestPossibleRegion() ); + m_BZeroImage->Allocate(); + + // if no GradientIndexMap is set take all Gradients to GradientIndicies-Vector and add this to GradientIndexMap[ALL] + // add All Bzero + if(m_GradientIndexMap.size() == 0){ + + // split for all gradient directions the image in baseline-indicies and gradient-indicies for a fast access + GradientDirectionContainerType::ConstIterator gdcit; + for( gdcit = this->m_GradientDirectionContainer->Begin(); gdcit != this->m_GradientDirectionContainer->End(); ++gdcit) { - itkExceptionMacro( << "At least " << ((L+1)*(L+2))/2 << " gradient directions are required" ); + if(gdcit.Value().one_norm() <= 0.0) // if vector length at current position is 0, it is a basline signal + { + m_GradientIndexMap[0].push_back(gdcit.Index()); + }else{ // it is a gradient signal of length != 0 + m_GradientIndexMap[1].push_back(gdcit.Index()); + } } - // Input must be an itk::VectorImage. - std::string gradientImageClassName(this->ProcessObject::GetInput(0)->GetNameOfClass()); - - if ( strcmp(gradientImageClassName.c_str(),"VectorImage") != 0 ) - itkExceptionMacro( << "There is only one Gradient image. I expect that to be a VectorImage. But its of type: " << gradientImageClassName ); - - m_BZeroImage = BZeroImageType::New(); - typename GradientImagesType::Pointer img = static_cast< GradientImagesType * >( this->ProcessObject::GetInput(0) ); - m_BZeroImage->SetSpacing( img->GetSpacing() ); // Set the image spacing - m_BZeroImage->SetOrigin( img->GetOrigin() ); // Set the image origin - m_BZeroImage->SetDirection( img->GetDirection() ); // Set the image direction - m_BZeroImage->SetLargestPossibleRegion( img->GetLargestPossibleRegion()); - m_BZeroImage->SetBufferedRegion( img->GetLargestPossibleRegion() ); - m_BZeroImage->Allocate(); - - // if no GradientIndexMap is set take all Gradients to GradientIndicies-Vector and add this to GradientIndexMap[ALL] - // add All Bzero - if(m_GradientIndexMap.size() == 0){ - - // split for all gradient directions the image in baseline-indicies and gradient-indicies for a fast access - GradientDirectionContainerType::ConstIterator gdcit; - for( gdcit = this->m_GradientDirectionContainer->Begin(); gdcit != this->m_GradientDirectionContainer->End(); ++gdcit) - { - if(gdcit.Value().one_norm() <= 0.0) // if vector length at current position is 0, it is a basline signal - { - m_GradientIndexMap[0].push_back(gdcit.Index()); - }else{ // it is a gradient signal of length != 0 - m_GradientIndexMap[1].push_back(gdcit.Index()); - } - } - m_ReconstructionType = Standard1Shell; - }else if(m_GradientIndexMap.size() == 4){ + m_ReconstructionType = Standard1Shell; + }else if(m_GradientIndexMap.size() == 4){ - GradientIndexMapIteraotr it = m_GradientIndexMap.begin(); - it++; - int b1 = (*it).first; - it++; - int b2 = (*it).first; - it++; - int b3 = (*it).first; + GradientIndexMapIteraotr it = m_GradientIndexMap.begin(); + it++; + int b1 = (*it).first; + it++; + int b2 = (*it).first; + it++; + int b3 = (*it).first; - if(b2 - b1 == b1 && b3 - b2 == b1 ) - { - m_ReconstructionType = Analytical3Shells; - }else - { - m_ReconstructionType = NumericalNShells; - } - }else if(m_GradientIndexMap.size() > 2) + if(b2 - b1 == b1 && b3 - b2 == b1 ) + { + m_ReconstructionType = Analytical3Shells; + }else { - m_ReconstructionType = NumericalNShells; + m_ReconstructionType = NumericalNShells; } + }else if(m_GradientIndexMap.size() > 2) + { + m_ReconstructionType = NumericalNShells; + } - switch(m_ReconstructionType) - { - case Analytical3Shells: - { - GradientIndexMapIteraotr it = m_GradientIndexMap.begin(); - it++; - this->ComputeReconstructionMatrix((*it).second); - break; - } - case NumericalNShells: - this->ComputeReconstructionMatrix(m_GradientIndexMap[1]); - break; - case Standard1Shell: - this->ComputeReconstructionMatrix(m_GradientIndexMap[1]); - break; - } + switch(m_ReconstructionType) + { + case Analytical3Shells: + { + GradientIndexMapIteraotr it = m_GradientIndexMap.begin(); + it++; + this->ComputeReconstructionMatrix((*it).second); + break; + } + case NumericalNShells: + this->ComputeReconstructionMatrix(m_GradientIndexMap[1]); + break; + case Standard1Shell: + this->ComputeReconstructionMatrix(m_GradientIndexMap[1]); + break; + } } template< class T, class TG, class TO, int L, int NODF> void DiffusionMultiShellQballReconstructionImageFilter ::StandardOneShellReconstruction(const OutputImageRegionType& outputRegionForThread) { - // Get output image pointer - typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); - - // ImageRegionIterator for the output image - ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread); - oit.GoToBegin(); - - // ImageRegionIterator for the BZero (output) image - ImageRegionIterator< BZeroImageType > oit2(m_BZeroImage, outputRegionForThread); - oit2.GoToBegin(); - - IndiciesVector BZeroIndicies = m_GradientIndexMap[0]; - IndiciesVector SignalIndicies = m_GradientIndexMap[1]; - - // if the gradient directiosn aragement is hemispherical, duplicate all gradient directions - // alone, interested in the value, the direction can be neglected - if(m_IsHemisphericalArrangementOfGradientDirections){ - int NumbersOfGradientIndicies = SignalIndicies.size(); - for (int i = 0 ; i < NumbersOfGradientIndicies; i++) - SignalIndicies.push_back(SignalIndicies[i]); - } + // Get output image pointer + typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); + // ImageRegionIterator for the output image + ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread); + oit.GoToBegin(); - // Get input gradient image pointer - typename GradientImagesType::Pointer gradientImagePointer = static_cast< GradientImagesType * >( this->ProcessObject::GetInput(0) ); + // ImageRegionIterator for the BZero (output) image + ImageRegionIterator< BZeroImageType > oit2(m_BZeroImage, outputRegionForThread); + oit2.GoToBegin(); - // Const ImageRegionIterator for input gradient image - typedef ImageRegionConstIterator< GradientImagesType > GradientIteratorType; - GradientIteratorType git(gradientImagePointer, outputRegionForThread ); - git.GoToBegin(); + IndiciesVector BZeroIndicies = m_GradientIndexMap[0]; + IndiciesVector SignalIndicies = m_GradientIndexMap[1]; - typedef typename GradientImagesType::PixelType GradientVectorType; + // if the gradient directiosn aragement is hemispherical, duplicate all gradient directions + // alone, interested in the value, the direction can be neglected + if(m_IsHemisphericalArrangementOfGradientDirections){ + int NumbersOfGradientIndicies = SignalIndicies.size(); + for (int i = 0 ; i < NumbersOfGradientIndicies; i++) + SignalIndicies.push_back(SignalIndicies[i]); + } - // iterate overall voxels of the gradient image region - while( ! git.IsAtEnd() ) - { - GradientVectorType b = git.Get(); - // compute the average bzero signal - typename NumericTraits::AccumulateType b0 = NumericTraits::Zero; - for(unsigned int i = 0; i < BZeroIndicies.size(); ++i) - { - b0 += b[BZeroIndicies[i]]; - } - b0 /= this->m_NumberOfBaselineImages; + // Get input gradient image pointer + typename GradientImagesType::Pointer gradientImagePointer = static_cast< GradientImagesType * >( this->ProcessObject::GetInput(0) ); - // ODF Vector - OdfPixelType odf(0.0); + // Const ImageRegionIterator for input gradient image + typedef ImageRegionConstIterator< GradientImagesType > GradientIteratorType; + GradientIteratorType git(gradientImagePointer, outputRegionForThread ); + git.GoToBegin(); - // Create the Signal Vector - vnl_vector SignalVector(m_NumberOfGradientDirections); - if( (b0 != 0) && (b0 >= m_Threshold) ) - { + typedef typename GradientImagesType::PixelType GradientVectorType; - for( unsigned int i = 0; i< SignalIndicies.size(); i++ ) - { - SignalVector[i] = static_cast(b[SignalIndicies[i]]); - } + // iterate overall voxels of the gradient image region + while( ! git.IsAtEnd() ) + { + GradientVectorType b = git.Get(); - // apply threashold an generate ln(-ln(E)) signal - // Replace SignalVector with PreNormalized SignalVector - PreNormalize(SignalVector, b0); - DoubleLogarithm(SignalVector); + // compute the average bzero signal + typename NumericTraits::AccumulateType b0 = NumericTraits::Zero; + for(unsigned int i = 0; i < BZeroIndicies.size(); ++i) + { + b0 += b[BZeroIndicies[i]]; + } + b0 /= this->m_NumberOfBaselineImages; - // ODF coeffs-vector - vnl_vector coeffs(m_NumberCoefficients); + // ODF Vector + OdfPixelType odf(0.0); - // approximate ODF coeffs - coeffs = ( (*m_CoeffReconstructionMatrix) * SignalVector ); - coeffs[0] += 1.0/(2.0*sqrt(QBALL_ANAL_RECON_PI)); + // Create the Signal Vector + vnl_vector SignalVector(m_NumberOfGradientDirections); + if( (b0 != 0) && (b0 >= m_Threshold) ) + { - odf = ( (*m_ODFSphericalHarmonicBasisMatrix) * coeffs ).data_block(); - } - // set ODF to ODF-Image - oit.Set( odf ); - ++oit; + for( unsigned int i = 0; i< SignalIndicies.size(); i++ ) + { + SignalVector[i] = static_cast(b[SignalIndicies[i]]); + } + + // apply threashold an generate ln(-ln(E)) signal + // Replace SignalVector with PreNormalized SignalVector + S_S0Normalization(SignalVector, b0); + DoubleLogarithm(SignalVector); - // set b0(average) to b0-Image - oit2.Set( b0 ); - ++oit2; + // ODF coeffs-vector + vnl_vector coeffs(m_NumberCoefficients); - ++git; + // approximate ODF coeffs + coeffs = ( (*m_CoeffReconstructionMatrix) * SignalVector ); + coeffs[0] += 1.0/(2.0*sqrt(QBALL_ANAL_RECON_PI)); + odf = ( (*m_ODFSphericalHarmonicBasisMatrix) * coeffs ).data_block(); } + // set ODF to ODF-Image + oit.Set( odf ); + ++oit; - MITK_INFO << "One Thread finished reconstruction"; + // set b0(average) to b0-Image + oit2.Set( b0 ); + ++oit2; + + ++git; + + } + + MITK_INFO << "One Thread finished reconstruction"; } template< class T, class TG, class TO, int L, int NODF> void DiffusionMultiShellQballReconstructionImageFilter ::AnalyticalThreeShellReconstruction(const OutputImageRegionType& outputRegionForThread) { - int wrongODF = 0; + int wrongODF = 0; - typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); + typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); - ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread); - oit.GoToBegin(); + ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread); + oit.GoToBegin(); - ImageRegionIterator< BZeroImageType > oit2(m_BZeroImage, outputRegionForThread); - oit2.GoToBegin(); + ImageRegionIterator< BZeroImageType > oit2(m_BZeroImage, outputRegionForThread); + oit2.GoToBegin(); - IndiciesVector BZeroIndicies = m_GradientIndexMap[0]; + IndiciesVector BZeroIndicies = m_GradientIndexMap[0]; - GradientIndexMapIteraotr it = m_GradientIndexMap.begin(); - it++; - IndiciesVector Shell1Indiecies = (*it).second; - it++; - IndiciesVector Shell2Indiecies = (*it).second; - it++; - IndiciesVector Shell3Indiecies = (*it).second; + GradientIndexMapIteraotr it = m_GradientIndexMap.begin(); + it++; + IndiciesVector Shell1Indiecies = (*it).second; + it++; + IndiciesVector Shell2Indiecies = (*it).second; + it++; + IndiciesVector Shell3Indiecies = (*it).second; - //assert(Shell1Indiecies.size() == Shell2Indiecies.size() && Shell1Indiecies.size() == Shell3Indiecies.size()); + //assert(Shell1Indiecies.size() == Shell2Indiecies.size() && Shell1Indiecies.size() == Shell3Indiecies.size()); - if(m_IsHemisphericalArrangementOfGradientDirections){ - int NumbersOfGradientIndicies = Shell1Indiecies.size(); - for (int i = 0 ; i < NumbersOfGradientIndicies; i++){ - Shell1Indiecies.push_back(Shell1Indiecies[i]); - Shell2Indiecies.push_back(Shell2Indiecies[i]); - Shell3Indiecies.push_back(Shell3Indiecies[i]); - } + if(m_IsHemisphericalArrangementOfGradientDirections){ + int NumbersOfGradientIndicies = Shell1Indiecies.size(); + for (int i = 0 ; i < NumbersOfGradientIndicies; i++){ + Shell1Indiecies.push_back(Shell1Indiecies[i]); + Shell2Indiecies.push_back(Shell2Indiecies[i]); + Shell3Indiecies.push_back(Shell3Indiecies[i]); } + } - // Get input gradient image pointer - typename GradientImagesType::Pointer gradientImagePointer = static_cast< GradientImagesType * >( this->ProcessObject::GetInput(0) ); + // Get input gradient image pointer + typename GradientImagesType::Pointer gradientImagePointer = static_cast< GradientImagesType * >( this->ProcessObject::GetInput(0) ); - // Const ImageRegionIterator for input gradient image - typedef ImageRegionConstIterator< GradientImagesType > GradientIteratorType; - GradientIteratorType git(gradientImagePointer, outputRegionForThread ); - git.GoToBegin(); + // Const ImageRegionIterator for input gradient image + typedef ImageRegionConstIterator< GradientImagesType > GradientIteratorType; + GradientIteratorType git(gradientImagePointer, outputRegionForThread ); + git.GoToBegin(); - typedef typename GradientImagesType::PixelType GradientVectorType; + typedef typename GradientImagesType::PixelType GradientVectorType; - //------------------------- Preperations Zone ------------------------------------ - vnl_matrix SHBasis(m_SHBasisMatrix->rows(),m_SHBasisMatrix->cols()); - SHBasis.fill(0.0); + //------------------------- Preperations Zone ------------------------------------ + vnl_matrix SHBasis(m_SHBasisMatrix->rows(),m_SHBasisMatrix->cols()); + SHBasis.fill(0.0); - for(int i=0; i ReconstructionMatrix(m_SignalReonstructionMatrix->rows(),m_SignalReonstructionMatrix->cols()); - ReconstructionMatrix.fill(0.0); + vnl_matrix ReconstructionMatrix(m_SignalReonstructionMatrix->rows(),m_SignalReonstructionMatrix->cols()); + ReconstructionMatrix.fill(0.0); - for(int i=0; i < ReconstructionMatrix.rows(); i++) - for(unsigned int j=0; j < ReconstructionMatrix.cols(); j++) - ReconstructionMatrix(i,j) = (TO)(*m_SignalReonstructionMatrix)(i,j); + for(int i=0; i < ReconstructionMatrix.rows(); i++) + for(unsigned int j=0; j < ReconstructionMatrix.cols(); j++) + ReconstructionMatrix(i,j) = (TO)(*m_SignalReonstructionMatrix)(i,j); + // N x 3 - vnl_vector< TO > Shell1OriginalSignal(Shell1Indiecies.size()); - vnl_vector< TO > Shell2OriginalSignal(Shell1Indiecies.size()); - vnl_vector< TO > Shell3OriginalSignal(Shell1Indiecies.size()); + // x E1 , E2 , E3 + // + // N : : : + // : : : + // + vnl_matrix< TO > E(Shell1Indiecies.size(), 3); + vnl_vector< TO > Shell1OriginalSignal(Shell1Indiecies.size()); //E1 + vnl_vector< TO > Shell2OriginalSignal(Shell1Indiecies.size()); //E2 + vnl_vector< TO > Shell3OriginalSignal(Shell1Indiecies.size()); //E3 - vnl_vector< TO > Shell2Coefficients(Shell1Indiecies.size()); - vnl_vector< TO > Shell3Coefficients(Shell1Indiecies.size()); - vnl_vector< TO > SHApproximatedSignal2(Shell1Indiecies.size()); - vnl_vector< TO > SHApproximatedSignal3(Shell1Indiecies.size()); + vnl_vector< TO > Shell2Coefficients(Shell1Indiecies.size()); + vnl_vector< TO > Shell3Coefficients(Shell1Indiecies.size()); + vnl_vector< TO > SHApproximatedSignal2(Shell1Indiecies.size()); + vnl_vector< TO > SHApproximatedSignal3(Shell1Indiecies.size()); - //------------------------- Preperations Zone END --------------------------------- + //------------------------- Preperations Zone END --------------------------------- - // iterate overall voxels of the gradient image region - while( ! git.IsAtEnd() ) - { - Shell1OriginalSignal.fill(0.0); - Shell2OriginalSignal.fill(0.0); - Shell3OriginalSignal.fill(0.0); - SHApproximatedSignal3.fill(0.0); - SHApproximatedSignal2.fill(0.0); + // iterate overall voxels of the gradient image region + while( ! git.IsAtEnd() ) + { + Shell1OriginalSignal.fill(0.0); + Shell2OriginalSignal.fill(0.0); + Shell3OriginalSignal.fill(0.0); - Shell2Coefficients.fill(0.0); - Shell3Coefficients.fill(0.0); + SHApproximatedSignal3.fill(0.0); + SHApproximatedSignal2.fill(0.0); - GradientVectorType b = git.Get(); + Shell2Coefficients.fill(0.0); + Shell3Coefficients.fill(0.0); - // compute the average bzero signal - typename NumericTraits::AccumulateType b0 = NumericTraits::Zero; - for(unsigned int i = 0; i < BZeroIndicies.size(); ++i) - { - b0 += b[BZeroIndicies[i]]; - } - b0 /= this->m_NumberOfBaselineImages; + GradientVectorType b = git.Get(); - // ODF Vector - OdfPixelType odf(0.0); + // compute the average bzero signal + typename NumericTraits::AccumulateType b0 = NumericTraits::Zero; + for(unsigned int i = 0; i < BZeroIndicies.size(); ++i) + { + b0 += b[BZeroIndicies[i]]; + } + b0 /= this->m_NumberOfBaselineImages; - // Create the Signal Vector - //vnl_vector SignalVector(m_NumberOfGradientDirections); - if( (b0 != 0) && (b0 >= m_Threshold) ) - { + // ODF Vector + OdfPixelType odf(0.0); - for(int i = 0 ; i < Shell2Indiecies.size(); i++) - { - Shell1OriginalSignal[i] = static_cast(b[Shell1Indiecies[i]]); - Shell2OriginalSignal[i] = static_cast(b[Shell2Indiecies[i]]); - Shell3OriginalSignal[i] = static_cast(b[Shell3Indiecies[i]]); - } + // Create the Signal Vector + //vnl_vector SignalVector(m_NumberOfGradientDirections); + if( (b0 != 0) && (b0 >= m_Threshold) ) + { + for(int i = 0 ; i < Shell2Indiecies.size(); i++) + { + Shell1OriginalSignal[i] = static_cast(b[Shell1Indiecies[i]]); + Shell2OriginalSignal[i] = static_cast(b[Shell2Indiecies[i]]); + Shell3OriginalSignal[i] = static_cast(b[Shell3Indiecies[i]]); + } - Shell2Coefficients = ReconstructionMatrix * Shell2OriginalSignal; - Shell3Coefficients = ReconstructionMatrix * Shell3OriginalSignal; - SHApproximatedSignal2 = SHBasis * Shell2Coefficients; - SHApproximatedSignal3 = SHBasis * Shell3Coefficients; + Shell2Coefficients = ReconstructionMatrix * Shell2OriginalSignal; + Shell3Coefficients = ReconstructionMatrix * Shell3OriginalSignal; + SHApproximatedSignal2 = SHBasis * Shell2Coefficients; + SHApproximatedSignal3 = SHBasis * Shell3Coefficients; - PreNormalize(Shell1OriginalSignal,b0); - PreNormalize(SHApproximatedSignal2,b0); - PreNormalize(SHApproximatedSignal3,b0); + E.set_column(0,Shell1OriginalSignal); + E.set_column(1,SHApproximatedSignal2); + E.set_column(2,SHApproximatedSignal3); - vnl_vector AlphaValues(Shell1Indiecies.size()); - vnl_vector BetaValues(Shell1Indiecies.size()); - vnl_vector LAValues(Shell1Indiecies.size()); + // Si/S0 + S_S0Normalization(E,b0); - bool notCompliedWithConditions = false; - for( unsigned int i = 0; i< Shell1Indiecies.size(); i++ ) - { - float E1 = Shell1OriginalSignal[i]; - float E2 = SHApproximatedSignal2[i]; - float E3 = SHApproximatedSignal3[i]; + //Implements Eq. [19] and Fig. 4. + Threshold(E, m_Threshold); - if(!(0 < E3) || (!(E3 < E2)) || (!(E2 < E1)) || (!(E1 < 1)) || (!(E1 * E1 < E2)) || (!(E2 * E2 < E1 * E3)) || (!(E3-E1*E2 < E2 - E1*E1 + E1*E3 -E2*E2))) - { - notCompliedWithConditions = true; - break; - } + //inqualities [31]. Taking the lograithm of th first tree inqualities + //convert the quadratic inqualities to linear ones. + Projection1(E, m_Threshold); - float A = (E3 -E1*E2) / (2*(E2-E1*E1)) ; - float B = sqrt( A * A - ((E1*E3 - E2 * E2) / (E2-E1*E1)) ); - //MITK_INFO << A << " " << B; + vnl_vector AlphaValues(Shell1Indiecies.size()); + vnl_vector BetaValues(Shell1Indiecies.size()); + vnl_vector LAValues(Shell1Indiecies.size()); - AlphaValues[i] = A + B; - BetaValues[i] = A - B; - LAValues[i] = 0.5 + ((E1 - A)/(2*B)); + bool notCompliedWithConditions = false; - if(A != A || B != B || LAValues[i] != LAValues[i]) - { - // MITK_INFO << A << " " << B << " " << LAValues[i]; - } + for( unsigned int i = 0; i< Shell1Indiecies.size(); i++ ) + { + float E1 = Shell1OriginalSignal[i]; + float E2 = SHApproximatedSignal2[i]; + float E3 = SHApproximatedSignal3[i]; - } + if(!(0 < E3) || (!(E3 < E2)) || (!(E2 < E1)) || (!(E1 < 1)) || (!(E1 * E1 < E2)) || (!(E2 * E2 < E1 * E3)) || (!(E3-E1*E2 < E2 - E1*E1 + E1*E3 -E2*E2))) + { + notCompliedWithConditions = true; + break; + } - if(!notCompliedWithConditions){ - Threshold(AlphaValues); - Threshold(BetaValues); + float A = (E3 -E1*E2) / (2*(E2-E1*E1)) ; + float B = sqrt( A * A - ((E1*E3 - E2 * E2) / (E2-E1*E1)) ); + //MITK_INFO << A << " " << B; - DoubleLogarithm(AlphaValues); - DoubleLogarithm(BetaValues); + AlphaValues[i] = A + B; + BetaValues[i] = A - B; + LAValues[i] = 0.5 + ((E1 - A)/(2*B)); - vnl_vector SignalVector(Shell1Indiecies.size()); - for( unsigned int i = 0; i< AlphaValues.size(); i++ ) - { - SignalVector = LAValues[i] * AlphaValues + LAValues[i] * BetaValues; - } + if(A != A || B != B || LAValues[i] != LAValues[i]) + { + // MITK_INFO << A << " " << B << " " << LAValues[i]; + } - vnl_vector coeffs(m_NumberCoefficients); + } - coeffs = ( (*m_CoeffReconstructionMatrix) * SignalVector ); - coeffs[0] += 1.0/(2.0*sqrt(QBALL_ANAL_RECON_PI)); + if(!notCompliedWithConditions){ - odf = ( (*m_ODFSphericalHarmonicBasisMatrix) * coeffs ).data_block(); - }else{ - odf.Fill(1/252); - wrongODF ++; - } + Threshold(AlphaValues); + Threshold(BetaValues); + DoubleLogarithm(AlphaValues); + DoubleLogarithm(BetaValues); + + vnl_vector SignalVector(Shell1Indiecies.size()); + for( unsigned int i = 0; i< AlphaValues.size(); i++ ) + { + SignalVector = LAValues[i] * AlphaValues + LAValues[i] * BetaValues; } - // set ODF to ODF-Image - oit.Set( odf ); - ++oit; - // MITK_INFO << odf; - // set b0(average) to b0-Image - oit2.Set( b0 ); - ++oit2; - ++git; + vnl_vector coeffs(m_NumberCoefficients); + + coeffs = ( (*m_CoeffReconstructionMatrix) * SignalVector ); + coeffs[0] += 1.0/(2.0*sqrt(QBALL_ANAL_RECON_PI)); + + odf = ( (*m_ODFSphericalHarmonicBasisMatrix) * coeffs ).data_block(); + }else{ + odf.Fill(1/252); + wrongODF ++; + } } + // set ODF to ODF-Image + oit.Set( odf ); + ++oit; + // MITK_INFO << odf; + // set b0(average) to b0-Image + oit2.Set( b0 ); + ++oit2; + ++git; - MITK_INFO << "THREAD FINISHED with " << wrongODF << " WrongODFs"; + } + + + MITK_INFO << "THREAD FINISHED with " << wrongODF << " WrongODFs"; } template< class T, class TG, class TO, int L, int NODF> vnl_vector DiffusionMultiShellQballReconstructionImageFilter ::AnalyticalThreeShellParameterEstimation(const IndiciesVector * shell1Indicies,const IndiciesVector * shell2Indicies ,const IndiciesVector * shell3Indicies, vnl_vector) { } template< class T, class TG, class TO, int L, int NODF> void DiffusionMultiShellQballReconstructionImageFilter ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int NumberOfThreads) { - switch(m_ReconstructionType) - { - case Standard1Shell: - StandardOneShellReconstruction(outputRegionForThread); - break; - case Analytical3Shells: - AnalyticalThreeShellReconstruction(outputRegionForThread); - break; - case NumericalNShells: - break; - } + switch(m_ReconstructionType) + { + case Standard1Shell: + StandardOneShellReconstruction(outputRegionForThread); + break; + case Analytical3Shells: + AnalyticalThreeShellReconstruction(outputRegionForThread); + break; + case NumericalNShells: + break; + } } template< class T, class TG, class TO, int L, int NODF> void DiffusionMultiShellQballReconstructionImageFilter:: ComputeSphericalHarmonicsBasis(vnl_matrix * QBallReference, vnl_matrix *SHBasisOutput, vnl_matrix* LaplaciaBaltramiOutput, vnl_vector* SHOrderAssociation, vnl_matrix* SHEigenvalues ) { - for(unsigned int i=0; i< (*SHBasisOutput).rows(); i++) + for(unsigned int i=0; i< (*SHBasisOutput).rows(); i++) + { + for(int k = 0; k <= L; k += 2) { - for(int k = 0; k <= L; k += 2) - { - for(int m =- k; m <= k; m++) - { - int j = ( k * k + k + 2 ) / 2 + m - 1; + for(int m =- k; m <= k; m++) + { + int j = ( k * k + k + 2 ) / 2 + m - 1; - // Compute SHBasisFunctions - double phi = (*QBallReference)(0,i); - double th = (*QBallReference)(1,i); - (*SHBasisOutput)(i,j) = mitk::ShericalHarmonicsFunctions::Yj(m,k,th,phi); + // Compute SHBasisFunctions + double phi = (*QBallReference)(0,i); + double th = (*QBallReference)(1,i); + (*SHBasisOutput)(i,j) = mitk::ShericalHarmonicsFunctions::Yj(m,k,th,phi); - // Laplacian Baltrami Order Association - if(LaplaciaBaltramiOutput) - (*LaplaciaBaltramiOutput)(j,j) = k*k*(k + 1)*(k+1); + // Laplacian Baltrami Order Association + if(LaplaciaBaltramiOutput) + (*LaplaciaBaltramiOutput)(j,j) = k*k*(k + 1)*(k+1); - // SHEigenvalues with order Accosiation kj - if(SHEigenvalues) - (*SHEigenvalues)(j,j) = -k* (k+1); + // SHEigenvalues with order Accosiation kj + if(SHEigenvalues) + (*SHEigenvalues)(j,j) = -k* (k+1); - // Order Association - if(SHOrderAssociation) - (*SHOrderAssociation)[j] = k; + // Order Association + if(SHOrderAssociation) + (*SHOrderAssociation)[j] = k; - } - } + } } + } } template< class T, class TG, class TO, int L, int NODF> void DiffusionMultiShellQballReconstructionImageFilter ::ComputeFunkRadonTransformationMatrix(vnl_vector* SHOrderAssociationReference, vnl_matrix* FRTMatrixOutput ) { - for(int i=0; i bool DiffusionMultiShellQballReconstructionImageFilter ::CheckHemisphericalArrangementOfGradientDirections() { - // handle acquisition schemes where only half of the spherical - // shell is sampled by the gradient directions. In this case, - // each gradient direction is duplicated in negative direction. - vnl_vector centerMass(3); - centerMass.fill(0.0); - int count = 0; - GradientDirectionContainerType::ConstIterator gdcit1; - for( gdcit1 = this->m_GradientDirectionContainer->Begin(); gdcit1 != this->m_GradientDirectionContainer->End(); ++gdcit1) + // 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; + GradientDirectionContainerType::ConstIterator gdcit1; + for( gdcit1 = this->m_GradientDirectionContainer->Begin(); gdcit1 != this->m_GradientDirectionContainer->End(); ++gdcit1) + { + if(gdcit1.Value().one_norm() > 0.0) { - if(gdcit1.Value().one_norm() > 0.0) - { - centerMass += gdcit1.Value(); - count ++; - } - } - centerMass /= count; - if(centerMass.two_norm() > 0.1) - { - return false; + centerMass += gdcit1.Value(); + count ++; } - return true; + } + centerMass /= count; + if(centerMass.two_norm() > 0.1) + { + return false; + } + return true; } template< class T, class TG, class TO, int L, int NOdfDirections> void DiffusionMultiShellQballReconstructionImageFilter ::ComputeReconstructionMatrix(std::vector< unsigned int > gradientIndiciesVector) { - typedef std::auto_ptr< vnl_matrix< double> > MatrixDoublePtr; - typedef std::auto_ptr< vnl_vector< int > > VectorIntPtr; - typedef std::auto_ptr< vnl_matrix_inverse< double > > InverseMatrixDoublePtr; + typedef std::auto_ptr< vnl_matrix< double> > MatrixDoublePtr; + typedef std::auto_ptr< vnl_vector< int > > VectorIntPtr; + typedef std::auto_ptr< vnl_matrix_inverse< double > > InverseMatrixDoublePtr; - int numberOfGradientDirections = gradientIndiciesVector.size(); + int numberOfGradientDirections = gradientIndiciesVector.size(); - if( numberOfGradientDirections < (((L+1)*(L+2))/2) || numberOfGradientDirections < 6 ) - { - itkExceptionMacro( << "At least (L+1)(L+2)/2 gradient directions are required" ); - } + if( numberOfGradientDirections < (((L+1)*(L+2))/2) || numberOfGradientDirections < 6 ) + { + itkExceptionMacro( << "At least (L+1)(L+2)/2 gradient directions are required" ); + } - CheckDuplicateDiffusionGradients(); + CheckDuplicateDiffusionGradients(); - // check if gradient directions are arrangement as a hemisphere(true) or sphere(false) - m_IsHemisphericalArrangementOfGradientDirections = CheckHemisphericalArrangementOfGradientDirections(); - if(m_IsHemisphericalArrangementOfGradientDirections) numberOfGradientDirections *= 2; + // check if gradient directions are arrangement as a hemisphere(true) or sphere(false) + m_IsHemisphericalArrangementOfGradientDirections = CheckHemisphericalArrangementOfGradientDirections(); + if(m_IsHemisphericalArrangementOfGradientDirections) numberOfGradientDirections *= 2; - MatrixDoublePtr Q(new vnl_matrix(3, numberOfGradientDirections)); - Q->fill(0.0); - // Cartesian to spherical coordinates - { - int j = 0; + MatrixDoublePtr Q(new vnl_matrix(3, numberOfGradientDirections)); + Q->fill(0.0); + // Cartesian to spherical coordinates + { + int j = 0; - for(int i = 0; i < gradientIndiciesVector.size(); i++) - { + for(int i = 0; i < gradientIndiciesVector.size(); i++) + { - double x = this->m_GradientDirectionContainer->ElementAt(gradientIndiciesVector[i]).get(0); - double y = this->m_GradientDirectionContainer->ElementAt(gradientIndiciesVector[i]).get(1); - double z = this->m_GradientDirectionContainer->ElementAt(gradientIndiciesVector[i]).get(2); - double cart[3]; - mitk::ShericalHarmonicsFunctions::Cart2Sph(x,y,z,cart); - (*Q)(0,j) = cart[0]; - (*Q)(1,j) = cart[1]; - (*Q)(2,j++) = cart[2]; + double x = this->m_GradientDirectionContainer->ElementAt(gradientIndiciesVector[i]).get(0); + double y = this->m_GradientDirectionContainer->ElementAt(gradientIndiciesVector[i]).get(1); + double z = this->m_GradientDirectionContainer->ElementAt(gradientIndiciesVector[i]).get(2); + double cart[3]; + mitk::ShericalHarmonicsFunctions::Cart2Sph(x,y,z,cart); + (*Q)(0,j) = cart[0]; + (*Q)(1,j) = cart[1]; + (*Q)(2,j++) = cart[2]; - } - if(m_IsHemisphericalArrangementOfGradientDirections) - { - for(int i = 0; i < gradientIndiciesVector.size(); i++) - { - double x = this->m_GradientDirectionContainer->ElementAt(gradientIndiciesVector[i]).get(0); - double y = this->m_GradientDirectionContainer->ElementAt(gradientIndiciesVector[i]).get(1); - double z = this->m_GradientDirectionContainer->ElementAt(gradientIndiciesVector[i]).get(2); - double cart[3]; - mitk::ShericalHarmonicsFunctions::Cart2Sph(x,y,z,cart); - (*Q)(0,j) = cart[0]; - (*Q)(1,j) = cart[1]; - (*Q)(2,j++) = cart[2]; - } - } } + if(m_IsHemisphericalArrangementOfGradientDirections) + { + for(int i = 0; i < gradientIndiciesVector.size(); i++) + { + double x = this->m_GradientDirectionContainer->ElementAt(gradientIndiciesVector[i]).get(0); + double y = this->m_GradientDirectionContainer->ElementAt(gradientIndiciesVector[i]).get(1); + double z = this->m_GradientDirectionContainer->ElementAt(gradientIndiciesVector[i]).get(2); + double cart[3]; + mitk::ShericalHarmonicsFunctions::Cart2Sph(x,y,z,cart); + (*Q)(0,j) = cart[0]; + (*Q)(1,j) = cart[1]; + (*Q)(2,j++) = cart[2]; + } + } + } - int LOrder = L; - m_NumberCoefficients = (int)(LOrder*LOrder + LOrder + 2.0)/2.0 + LOrder; + int LOrder = L; + m_NumberCoefficients = (int)(LOrder*LOrder + LOrder + 2.0)/2.0 + LOrder; - m_SHBasisMatrix = new vnl_matrix(numberOfGradientDirections,m_NumberCoefficients); - m_SHBasisMatrix->fill(0.0); - VectorIntPtr SHOrderAssociation(new vnl_vector(m_NumberCoefficients)); - SHOrderAssociation->fill(0.0); - MatrixDoublePtr LaplacianBaltrami(new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients)); - LaplacianBaltrami->fill(0.0); - MatrixDoublePtr FRTMatrix(new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients)); - FRTMatrix->fill(0.0); - MatrixDoublePtr SHEigenvalues(new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients)); - SHEigenvalues->fill(0.0); + m_SHBasisMatrix = new vnl_matrix(numberOfGradientDirections,m_NumberCoefficients); + m_SHBasisMatrix->fill(0.0); + VectorIntPtr SHOrderAssociation(new vnl_vector(m_NumberCoefficients)); + SHOrderAssociation->fill(0.0); + MatrixDoublePtr LaplacianBaltrami(new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients)); + LaplacianBaltrami->fill(0.0); + MatrixDoublePtr FRTMatrix(new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients)); + FRTMatrix->fill(0.0); + MatrixDoublePtr SHEigenvalues(new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients)); + SHEigenvalues->fill(0.0); - // SHBasis-Matrix + LaplacianBaltrami-Matrix + SHOrderAssociationVector - ComputeSphericalHarmonicsBasis(Q.get() ,m_SHBasisMatrix, LaplacianBaltrami.get(), SHOrderAssociation.get(), SHEigenvalues.get()); + // SHBasis-Matrix + LaplacianBaltrami-Matrix + SHOrderAssociationVector + ComputeSphericalHarmonicsBasis(Q.get() ,m_SHBasisMatrix, LaplacianBaltrami.get(), SHOrderAssociation.get(), SHEigenvalues.get()); - // Compute FunkRadon Transformation Matrix Associated to SHBasis Order lj - ComputeFunkRadonTransformationMatrix(SHOrderAssociation.get() ,FRTMatrix.get()); + // Compute FunkRadon Transformation Matrix Associated to SHBasis Order lj + ComputeFunkRadonTransformationMatrix(SHOrderAssociation.get() ,FRTMatrix.get()); - MatrixDoublePtr temp(new vnl_matrix(((m_SHBasisMatrix->transpose()) * (*m_SHBasisMatrix)) + (m_Lambda * (*LaplacianBaltrami)))); + MatrixDoublePtr temp(new vnl_matrix(((m_SHBasisMatrix->transpose()) * (*m_SHBasisMatrix)) + (m_Lambda * (*LaplacianBaltrami)))); - InverseMatrixDoublePtr pseudo_inv(new vnl_matrix_inverse((*temp))); - MatrixDoublePtr inverse(new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients)); - inverse->fill(0.0); - (*inverse) = pseudo_inv->inverse(); + InverseMatrixDoublePtr pseudo_inv(new vnl_matrix_inverse((*temp))); + MatrixDoublePtr inverse(new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients)); + inverse->fill(0.0); + (*inverse) = pseudo_inv->inverse(); - // ODF Factor ( missing 1/4PI ?? ) - double factor = (1.0/(16.0*QBALL_ANAL_RECON_PI*QBALL_ANAL_RECON_PI)); + // ODF Factor ( missing 1/4PI ?? ) + double factor = (1.0/(16.0*QBALL_ANAL_RECON_PI*QBALL_ANAL_RECON_PI)); - m_SignalReonstructionMatrix = new vnl_matrix(); - m_SignalReonstructionMatrix->fill(0.0); - (*m_SignalReonstructionMatrix) = (*inverse) * (m_SHBasisMatrix->transpose()); + m_SignalReonstructionMatrix = new vnl_matrix(); + m_SignalReonstructionMatrix->fill(0.0); + (*m_SignalReonstructionMatrix) = (*inverse) * (m_SHBasisMatrix->transpose()); - MatrixDoublePtr TransformationMatrix(new vnl_matrix( factor * ((*FRTMatrix) * ((*SHEigenvalues) * (*m_SignalReonstructionMatrix))) ) ); + MatrixDoublePtr TransformationMatrix(new vnl_matrix( factor * ((*FRTMatrix) * ((*SHEigenvalues) * (*m_SignalReonstructionMatrix))) ) ); - m_CoeffReconstructionMatrix = new vnl_matrix(m_NumberCoefficients,numberOfGradientDirections); + m_CoeffReconstructionMatrix = new vnl_matrix(m_NumberCoefficients,numberOfGradientDirections); - // Cast double to float - for(int i=0; iodfs later + // this code goes to the image adapter coeffs->odfs later - vnl_matrix_fixed* U = itk::PointShell >::DistributePointShell(); + vnl_matrix_fixed* U = itk::PointShell >::DistributePointShell(); - m_ODFSphericalHarmonicBasisMatrix = new vnl_matrix(NOdfDirections,m_NumberCoefficients); + m_ODFSphericalHarmonicBasisMatrix = new vnl_matrix(NOdfDirections,m_NumberCoefficients); - for(int i=0; i((*m_ODFSphericalHarmonicBasisMatrix) * (*m_CoeffReconstructionMatrix)); + m_ReconstructionMatrix = new vnl_matrix((*m_ODFSphericalHarmonicBasisMatrix) * (*m_CoeffReconstructionMatrix)); } template< class T, class TG, class TO, int L, int NODF> template< class VNLType > void DiffusionMultiShellQballReconstructionImageFilter ::printMatrix( VNLType * mat ) { - std::stringstream stream; + std::stringstream stream; - for(int i = 0 ; i < mat->rows(); i++) + for(int i = 0 ; i < mat->rows(); i++) + { + stream.str(""); + for(int j = 0; j < mat->cols(); j++) { - stream.str(""); - for(int j = 0; j < mat->cols(); j++) - { - stream << (*mat)(i,j) << " "; - } - + stream << (*mat)(i,j) << " "; } + } - MITK_INFO << stream.str(); + + MITK_INFO << stream.str(); } template< class T, class TG, class TO, int L, int NODF> bool DiffusionMultiShellQballReconstructionImageFilter ::CheckDuplicateDiffusionGradients() { - GradientDirectionContainerType::ConstIterator gdcit1; - GradientDirectionContainerType::ConstIterator gdcit2; + GradientDirectionContainerType::ConstIterator gdcit1; + GradientDirectionContainerType::ConstIterator gdcit2; - for(gdcit1 = this->m_GradientDirectionContainer->Begin(); gdcit1 != this->m_GradientDirectionContainer->End(); ++gdcit1) + for(gdcit1 = this->m_GradientDirectionContainer->Begin(); gdcit1 != this->m_GradientDirectionContainer->End(); ++gdcit1) + { + for(gdcit2 = this->m_GradientDirectionContainer->Begin(); gdcit2 != this->m_GradientDirectionContainer->End(); ++gdcit2) { - for(gdcit2 = this->m_GradientDirectionContainer->Begin(); gdcit2 != this->m_GradientDirectionContainer->End(); ++gdcit2) - { - if(gdcit1.Value() == gdcit2.Value() && gdcit1.Index() != gdcit2.Index()) - { - itkWarningMacro( << "Some of the Diffusion Gradients equal each other. Corresponding image data should be averaged before calling this filter." ); - return true; - } - } + if(gdcit1.Value() == gdcit2.Value() && gdcit1.Index() != gdcit2.Index()) + { + itkWarningMacro( << "Some of the Diffusion Gradients equal each other. Corresponding image data should be averaged before calling this filter." ); + return true; + } } - return false; + } + return false; } 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; - - os.imbue( originalLocale ); + std::locale C("C"); + std::locale originalLocale = os.getloc(); + os.imbue(C); + + Superclass::PrintSelf(os,indent); + + os << indent << "OdfReconstructionMatrix: " << m_ReconstructionMatrix << std::endl; + if ( m_GradientDirectionContainer ) + { + os << indent << "GradientDirectionContainer: " + << m_GradientDirectionContainer << std::endl; + } + else + { + os << indent << + "GradientDirectionContainer: (Gradient directions not set)" << std::endl; + } + os << indent << "NumberOfGradientDirections: " << + m_NumberOfGradientDirections << std::endl; + os << indent << "NumberOfBaselineImages: " << + m_NumberOfBaselineImages << std::endl; + os << indent << "Threshold for reference B0 image: " << m_Threshold << std::endl; + os << indent << "BValue: " << m_BValue << std::endl; + + os.imbue( originalLocale ); } } #endif // __itkDiffusionMultiShellQballReconstructionImageFilter_cpp diff --git a/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h b/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h index 6c6a207d36..b50477affc 100644 --- a/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h +++ b/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h @@ -1,230 +1,236 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Language: C++ Date: $Date: 2009-07-14 19:11:20 +0200 (Tue, 14 Jul 2009) $ Version: $Revision: 18127 $ Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #ifndef __itkDiffusionMultiShellQballReconstructionImageFilter_h_ #define __itkDiffusionMultiShellQballReconstructionImageFilter_h_ #include "itkImageToImageFilter.h" #include "vnl/vnl_vector_fixed.h" #include "vnl/vnl_matrix.h" #include "vnl/algo/vnl_svd.h" #include "itkVectorContainer.h" #include "itkVectorImage.h" #include namespace itk{ /** \class DiffusionMultiShellQballReconstructionImageFilter Aganj_2010 */ template< class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NOrderL, int NrOdfDirections> class DiffusionMultiShellQballReconstructionImageFilter : public ImageToImageFilter< Image< TReferenceImagePixelType, 3 >, Image< Vector< TOdfPixelType, NrOdfDirections >, 3 > > { public: typedef DiffusionMultiShellQballReconstructionImageFilter Self; typedef SmartPointer 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; typedef std::map > GradientIndexMap; typedef std::map >::iterator GradientIndexMapIteraotr; typedef std::vector IndiciesVector; // --------------------------------------------------------------------------------------------// /** Method for creation through the object factory. */ itkNewMacro(Self); /** Runtime information support. */ itkTypeMacro(DiffusionMultiShellQballReconstructionImageFilter, ImageToImageFilter); /** set method to add gradient directions and its corresponding * image. The image here is a VectorImage. The user is expected to pass the * gradient directions in a container. The ith element of the container * corresponds to the gradient direction of the ith component image the * VectorImage. For the baseline image, a vector of all zeros * should be set.*/ void SetGradientImage( GradientDirectionContainerType *, const GradientImagesType *image); void SetGradientIndexMap(GradientIndexMap gradientIdexMap) { m_GradientIndexMap = gradientIdexMap; /* std::stringstream s1, s2, s3; for(int i = 0; i < m_GradientIndexMap[1000].size() ; i++){ s1 << m_GradientIndexMap[1000][i] << " "; s2 << m_GradientIndexMap[2000][i] << " "; s3 << m_GradientIndexMap[3000][i] << " "; } MITK_INFO << "1 SHELL " << std::endl << s1.str(); MITK_INFO << "2 SHELL " << std::endl << s2.str(); MITK_INFO << "3 SHELL " << std::endl << s3.str(); */ } /** Get reference image */ virtual ReferenceImageType * GetReferenceImage() { return ( static_cast< ReferenceImageType *>(this->ProcessObject::GetInput(0)) ); } /** Return the gradient direction. idx is 0 based */ virtual GradientDirectionType GetGradientDirection( unsigned int idx) const { if( idx >= m_GradientDirectionContainer->Size() ) { itkExceptionMacro( << "Gradient direction " << idx << "does not exist" ); } return m_GradientDirectionContainer->ElementAt( idx+1 ); } OdfPixelType Normalize(OdfPixelType odf, typename NumericTraits::AccumulateType b0 ); - void PreNormalize( vnl_vector & vec, typename NumericTraits::AccumulateType b0 = 0 ); + void S_S0Normalization( vnl_vector & vec, typename NumericTraits::AccumulateType b0 = 0 ); + void S_S0Normalization( vnl_matrix & mat, typename NumericTraits::AccumulateType b0 = 0 ); + void DoubleLogarithm(vnl_vector & vec); + void Threshold(vnl_vector & vec, float sigma = 0.0001); + void Threshold(vnl_matrix & mat, float sigma = 0.0001); + + void Projection1( vnl_matrix & mat, float sigma = 0.0001); /** Threshold on the reference image data. The output ODF will be a null * pdf for pixels in the reference image that have a value less than this * threshold. */ itkSetMacro( Threshold, ReferencePixelType ); itkGetMacro( Threshold, ReferencePixelType ); itkGetMacro( BZeroImage, typename BZeroImageType::Pointer); //itkGetMacro( ODFSumImage, typename BlaImage::Pointer); itkSetMacro( BValue, TOdfPixelType); itkSetMacro( Lambda, double ); itkGetMacro( Lambda, double ); itkGetConstReferenceMacro( BValue, TOdfPixelType); protected: DiffusionMultiShellQballReconstructionImageFilter(); ~DiffusionMultiShellQballReconstructionImageFilter() {}; void PrintSelf(std::ostream& os, Indent indent) const; void ComputeReconstructionMatrix(std::vector< unsigned int > gradientIndiciesVector); bool CheckDuplicateDiffusionGradients(); void ComputeSphericalHarmonicsBasis(vnl_matrix* QBallReference, vnl_matrix* SHBasisOutput, vnl_matrix* LaplaciaBaltramiOutput, vnl_vector* SHOrderAssociation , vnl_matrix * SHEigenvalues); void ComputeFunkRadonTransformationMatrix(vnl_vector* SHOrderAssociationReference, vnl_matrix* FRTMatrixOutput ); bool CheckHemisphericalArrangementOfGradientDirections(); void BeforeThreadedGenerateData(); void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, int NumberOfThreads ); vnl_vector AnalyticalThreeShellParameterEstimation(const IndiciesVector * shell1, const IndiciesVector * shell2, const IndiciesVector * shell3, vnl_vector b); void StandardOneShellReconstruction(const OutputImageRegionType& outputRegionForThread); void AnalyticalThreeShellReconstruction(const OutputImageRegionType& outputRegionForThread); private: enum ReconstructionType { Analytical3Shells, NumericalNShells, Standard1Shell }; OdfReconstructionMatrixType m_ReconstructionMatrix; OdfReconstructionMatrixType m_CoeffReconstructionMatrix; OdfReconstructionMatrixType m_ODFSphericalHarmonicBasisMatrix; CoefficientMatrixType m_SignalReonstructionMatrix; CoefficientMatrixType m_SHBasisMatrix; /** container to hold gradient directions */ GradientDirectionContainerType::Pointer m_GradientDirectionContainer; /** Number of gradient measurements */ unsigned int m_NumberOfGradientDirections; /** Number of baseline images */ unsigned int m_NumberOfBaselineImages; /** Threshold on the reference image data */ ReferencePixelType m_Threshold; /** LeBihan's b-value for normalizing tensors */ TOdfPixelType m_BValue; typename BZeroImageType::Pointer m_BZeroImage; GradientIndexMap m_GradientIndexMap; double m_Lambda; bool m_IsHemisphericalArrangementOfGradientDirections; bool m_IsArithmeticProgession; int m_NumberCoefficients; ReconstructionType m_ReconstructionType; template< class VNLType > void printMatrix( VNLType * mat ); }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkDiffusionMultiShellQballReconstructionImageFilter.cpp" #endif #endif //__itkDiffusionMultiShellQballReconstructionImageFilter_h_