diff --git a/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.cpp b/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.cpp index d31582eeeb..4b8d2f2ca3 100644 --- a/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.cpp +++ b/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.cpp @@ -1,580 +1,583 @@ /*========================================================================= 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_GradientDirectionContainer(NULL), + m_GradientIndexMap(NULL), + m_NumberOfGradientDirections(0), + m_NumberOfBaselineImages(1), + m_Threshold(NumericTraits< ReferencePixelType >::NonpositiveMin()), + m_BValue(1.0), + m_Lambda(0.0), + m_IsHemisphericalArrangementOfGradientDirections(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 ); + // 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; + for(int i=0; i0) + odf /= sum; + + return odf; } template vnl_vectoritk::DiffusionMultiShellQballReconstructionImageFilter ::PreNormalize( vnl_vector vec, typename NumericTraits::AccumulateType b0 ) { - // Threshold - double sigma = 0.001; - double temp; - int n = vec.size(); - double b0f = (double)b0; - for(int i=0; i " << vec[i]; + vec[i] /= b0f; + //temp = vec[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); + } + // MITK_INFO << temp << " <-> " << vec[i]; - vec[i] = log(-log(vec[i])); - } - return vec; + vec[i] = log(-log(vec[i])); + } + return vec; } 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; + this->m_GradientDirectionContainer = gradientDirection; + this->m_NumberOfBaselineImages = 0; - for(GradientDirectionContainerType::Iterator it = gradientDirection->Begin(); it != gradientDirection->End(); it++) - { - if( it.Value().one_norm() <= 0.0 ) - this->m_NumberOfBaselineImages ++; - else - it.Value() = it.Value() / it.Value().two_norm(); - } - - this->m_NumberOfGradientDirections = gradientDirection->Size() - this->m_NumberOfBaselineImages; - - // ensure that the gradient image we received has as many components as - // the number of gradient directions - if( gradientImage->GetVectorLength() != this->m_NumberOfBaselineImages + m_NumberOfGradientDirections ) - { - itkExceptionMacro( << m_NumberOfGradientDirections << " gradients + " << this->m_NumberOfBaselineImages - << "baselines = " << m_NumberOfGradientDirections + this->m_NumberOfBaselineImages - << " directions specified but image has " << gradientImage->GetVectorLength() - << " components."); - } - - this->ProcessObject::SetNthInput( 0, const_cast< GradientImagesType* >(gradientImage) ); + 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 */ ) - { - 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 ); - - 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(); + 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 ); + + 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(); } template< class T, class TG, class TO, int L, int NODF> void DiffusionMultiShellQballReconstructionImageFilter ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int NumberOfThreads) { - // 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 output image pointer + typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); - // ImageRegionIterator for the BZero (output) image - ImageRegionIterator< BZeroImageType > oit2(m_BZeroImage, outputRegionForThread); - oit2.GoToBegin(); + // 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(); - // define basline-signal indicies vector - std::vector baselineIndicies; - // define gradient-signal indicies vector - std::vector gradientIndicies; + // if no GradientIndexMap is set take all Gradients to GradientIndicies-Vector and add this to GradientIndexMap[ALL] + // add All Bzero + if(m_GradientIndexMap == 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) - { + m_GradientIndexMap = new GradientIndexMap(); - if(gdcit.Value().one_norm() <= 0.0) // if vector length at current position is 0, it is a basline signal - { - baselineIndicies.push_back(gdcit.Index()); - }else{ // it is a gradient signal of length != 0 - gradientIndicies.push_back(gdcit.Index()); + // 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)[BZero].push_back(gdcit.Index()); + }else{ // it is a gradient signal of length != 0 + (*m_GradientIndexMap)[AllGradients].push_back(gdcit.Index()); + } + } } - } - - // 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 = gradientIndicies.size(); - for (int i = 0 ; i < NumbersOfGradientIndicies; i++) - gradientIndicies.push_back(gradientIndicies[i]); - } + // 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 = (*m_GradientIndexMap)[AllGradients].size(); + for (int i = 0 ; i < NumbersOfGradientIndicies; i++) + (*m_GradientIndexMap)[AllGradients].push_back((*m_GradientIndexMap)[AllGradients][i]); + } - // Get input gradient image pointer - typename GradientImagesType::Pointer gradientImagePointer = static_cast< GradientImagesType * >( this->ProcessObject::GetInput(0) ); - // Const ImageRegionIterator for input gradient image - typedef ImageRegionConstIterator< GradientImagesType > GradientIteratorType; - GradientIteratorType git(gradientImagePointer, outputRegionForThread ); - git.GoToBegin(); + // Get input gradient image pointer + typename GradientImagesType::Pointer gradientImagePointer = static_cast< GradientImagesType * >( this->ProcessObject::GetInput(0) ); - typedef typename GradientImagesType::PixelType GradientVectorType; + // Const ImageRegionIterator for input gradient image + typedef ImageRegionConstIterator< GradientImagesType > GradientIteratorType; + GradientIteratorType git(gradientImagePointer, outputRegionForThread ); + git.GoToBegin(); - // iterate overall voxels of the gradient image region - while( ! git.IsAtEnd() ) - { - GradientVectorType b = git.Get(); + typedef typename GradientImagesType::PixelType GradientVectorType; - // compute the average bzero signal - typename NumericTraits::AccumulateType b0 = NumericTraits::Zero; - for(unsigned int i = 0; i < baselineIndicies.size(); ++i) + // iterate overall voxels of the gradient image region + while( ! git.IsAtEnd() ) { - b0 += b[baselineIndicies[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 < (*m_GradientIndexMap)[BZero].size(); ++i) + { + b0 += b[(*m_GradientIndexMap)[BZero][i]]; + } + b0 /= this->m_NumberOfBaselineImages; - // Create the Signal Vector - vnl_vector SignalVector(m_NumberOfGradientDirections); - if( (b0 != 0) && (b0 >= m_Threshold) ) - { - for( unsigned int i = 0; i< m_NumberOfGradientDirections; i++ ) - { - SignalVector[i] = static_cast(b[gradientIndicies[i]]); - } + // ODF Vector + OdfPixelType odf(0.0); - // apply threashold an generate ln(-ln(E)) signal - // Replace SignalVector with PreNormalized SignalVector - SignalVector = PreNormalize(SignalVector, b0); + // Create the Signal Vector + vnl_vector SignalVector(m_NumberOfGradientDirections); + if( (b0 != 0) && (b0 >= m_Threshold) ) + { + IndiciesVector & gradientIndicies = (*m_GradientIndexMap)[AllGradients]; + for( unsigned int i = 0; i< m_NumberOfGradientDirections; i++ ) + { + unsigned int index = gradientIndicies[i]; + SignalVector[i] = static_cast(b[index]); + } - // ODF coeffs-vector - vnl_vector coeffs(m_NumberCoefficients); + // apply threashold an generate ln(-ln(E)) signal + // Replace SignalVector with PreNormalized SignalVector + SignalVector = PreNormalize(SignalVector, b0); - // approximate ODF coeffs - coeffs = ( (*m_CoeffReconstructionMatrix) * SignalVector ); - coeffs[0] += 1.0/(2.0*sqrt(QBALL_ANAL_RECON_PI)); + // ODF coeffs-vector + vnl_vector coeffs(m_NumberCoefficients); - odf = ( (*m_SphericalHarmonicBasisMatrix) * coeffs ).data_block(); - } - // set ODF to ODF-Image - oit.Set( odf ); - ++oit; + // approximate ODF coeffs + coeffs = ( (*m_CoeffReconstructionMatrix) * SignalVector ); + coeffs[0] += 1.0/(2.0*sqrt(QBALL_ANAL_RECON_PI)); - // set b0(average) to b0-Image - oit2.Set( b0 ); - ++oit2; + odf = ( (*m_SphericalHarmonicBasisMatrix) * coeffs ).data_block(); + } + // set ODF to ODF-Image + oit.Set( odf ); + ++oit; - ++git; + // set b0(average) to b0-Image + oit2.Set( b0 ); + ++oit2; - } + ++git; - MITK_INFO << "One Thread finished reconstruction"; + } + + MITK_INFO << "One Thread finished reconstruction"; } - 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 ) - { +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(int k = 0; k <= L; k += 2) - { - for(int m =- k; m <= k; m++) + for(int k = 0; k <= L; k += 2) { - 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 ) - { +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() - { +template< class T, class TG, class TO, int L, int NOdfDirections> +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) { - if(gdcit1.Value().one_norm() > 0.0) - { - centerMass += gdcit1.Value(); - count ++; - } + if(gdcit1.Value().one_norm() > 0.0) + { + centerMass += gdcit1.Value(); + count ++; + } } centerMass /= count; if(centerMass.two_norm() > 0.1) { - return false; + return false; } return true; - } +} - template< class T, class TG, class TO, int L, int NOdfDirections> - void DiffusionMultiShellQballReconstructionImageFilter - ::ComputeReconstructionMatrix() - { +template< class T, class TG, class TO, int L, int NOdfDirections> +void DiffusionMultiShellQballReconstructionImageFilter +::ComputeReconstructionMatrix() +{ 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; if( m_NumberOfGradientDirections < (((L+1)*(L+2))/2) || m_NumberOfGradientDirections < 6 ) { - itkExceptionMacro( << "At least (L+1)(L+2)/2 gradient directions are required" ); + itkExceptionMacro( << "At least (L+1)(L+2)/2 gradient directions are required" ); } CheckDuplicateDiffusionGradients(); // check if gradient directions are arrangement as a hemisphere(true) or sphere(false) m_IsHemisphericalArrangementOfGradientDirections = CheckHemisphericalArrangementOfGradientDirections(); if(m_IsHemisphericalArrangementOfGradientDirections) m_NumberOfGradientDirections *= 2; MatrixDoublePtr Q(new vnl_matrix(3, m_NumberOfGradientDirections)); Q->fill(0.0); // Cartesian to spherical coordinates { - int i = 0; - GradientDirectionContainerType::ConstIterator gdcit; - for( gdcit = this->m_GradientDirectionContainer->Begin(); gdcit != this->m_GradientDirectionContainer->End(); ++gdcit) - { - if(gdcit.Value().one_norm() > 0.0) + int i = 0; + GradientDirectionContainerType::ConstIterator gdcit; + for( gdcit = this->m_GradientDirectionContainer->Begin(); gdcit != this->m_GradientDirectionContainer->End(); ++gdcit) { - double x = gdcit.Value().get(0); - double y = gdcit.Value().get(1); - double z = gdcit.Value().get(2); - double cart[3]; - mitk::ShericalHarmonicsFunctions::Cart2Sph(x,y,z,cart); - (*Q)(0,i) = cart[0]; - (*Q)(1,i) = cart[1]; - (*Q)(2,i++) = cart[2]; + if(gdcit.Value().one_norm() > 0.0) + { + double x = gdcit.Value().get(0); + double y = gdcit.Value().get(1); + double z = gdcit.Value().get(2); + double cart[3]; + mitk::ShericalHarmonicsFunctions::Cart2Sph(x,y,z,cart); + (*Q)(0,i) = cart[0]; + (*Q)(1,i) = cart[1]; + (*Q)(2,i++) = cart[2]; + } } - } - if(m_IsHemisphericalArrangementOfGradientDirections) - { - for( gdcit = this->m_GradientDirectionContainer->Begin(); gdcit != this->m_GradientDirectionContainer->End(); ++gdcit) + if(m_IsHemisphericalArrangementOfGradientDirections) { - if(gdcit.Value().one_norm() > 0.0) - { - double x = gdcit.Value().get(0); - double y = gdcit.Value().get(1); - double z = gdcit.Value().get(2); - double cart[3]; - mitk::ShericalHarmonicsFunctions::Cart2Sph(x,y,z,cart); - (*Q)(0,i) = cart[0]; - (*Q)(1,i) = cart[1]; - (*Q)(2,i++) = cart[2]; - } + for( gdcit = this->m_GradientDirectionContainer->Begin(); gdcit != this->m_GradientDirectionContainer->End(); ++gdcit) + { + if(gdcit.Value().one_norm() > 0.0) + { + double x = gdcit.Value().get(0); + double y = gdcit.Value().get(1); + double z = gdcit.Value().get(2); + double cart[3]; + mitk::ShericalHarmonicsFunctions::Cart2Sph(x,y,z,cart); + (*Q)(0,i) = cart[0]; + (*Q)(1,i) = cart[1]; + (*Q)(2,i++) = cart[2]; + } + } } - } } int LOrder = L; m_NumberCoefficients = (int)(LOrder*LOrder + LOrder + 2.0)/2.0 + LOrder; MatrixDoublePtr SHBasis(new vnl_matrix(m_NumberOfGradientDirections,m_NumberCoefficients)); SHBasis->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() ,SHBasis.get(), LaplacianBaltrami.get(), SHOrderAssociation.get(), SHEigenvalues.get()); // Compute FunkRadon Transformation Matrix Associated to SHBasis Order lj ComputeFunkRadonTransformationMatrix(SHOrderAssociation.get() ,FRTMatrix.get()); MatrixDoublePtr temp(new vnl_matrix(((SHBasis->transpose()) * (*SHBasis)) + (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(); // ODF Factor ( missing 1/4PI ?? ) double factor = (1.0/(16.0*QBALL_ANAL_RECON_PI*QBALL_ANAL_RECON_PI)); MatrixDoublePtr TransformationMatrix(new vnl_matrix( factor * ((*FRTMatrix) * ((*SHEigenvalues) * ((*inverse) * (SHBasis->transpose()) ) ) ))); m_CoeffReconstructionMatrix = new vnl_matrix(m_NumberCoefficients,m_NumberOfGradientDirections); // Cast double to float for(int i=0; iodfs later vnl_matrix_fixed* U = itk::PointShell >::DistributePointShell(); m_SphericalHarmonicBasisMatrix = new vnl_matrix(NOdfDirections,m_NumberCoefficients); 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> - template< class VNLType > - void DiffusionMultiShellQballReconstructionImageFilter - ::printMatrix( VNLType * mat ) - { +template< class T, class TG, class TO, int L, int NODF> +template< class VNLType > +void DiffusionMultiShellQballReconstructionImageFilter +::printMatrix( VNLType * mat ) +{ std::stringstream stream; for(int i = 0 ; i < mat->rows(); i++) { - stream.str(""); - for(int j = 0; j < mat->cols(); j++) - { - stream << (*mat)(i,j) << " "; - } - MITK_INFO << stream.str(); + stream.str(""); + for(int j = 0; j < mat->cols(); j++) + { + stream << (*mat)(i,j) << " "; + } + MITK_INFO << stream.str(); } - } +} - template< class T, class TG, class TO, int L, int NODF> - bool DiffusionMultiShellQballReconstructionImageFilter - ::CheckDuplicateDiffusionGradients() - { +template< class T, class TG, class TO, int L, int NODF> +bool DiffusionMultiShellQballReconstructionImageFilter +::CheckDuplicateDiffusionGradients() +{ GradientDirectionContainerType::ConstIterator gdcit1; GradientDirectionContainerType::ConstIterator gdcit2; for(gdcit1 = this->m_GradientDirectionContainer->Begin(); gdcit1 != this->m_GradientDirectionContainer->End(); ++gdcit1) { - for(gdcit2 = this->m_GradientDirectionContainer->Begin(); gdcit2 != this->m_GradientDirectionContainer->End(); ++gdcit2) - { - if(gdcit1.Value() == gdcit2.Value() && gdcit1.Index() != gdcit2.Index()) + for(gdcit2 = this->m_GradientDirectionContainer->Begin(); gdcit2 != this->m_GradientDirectionContainer->End(); ++gdcit2) { - 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; - } +} - template< class T, class TG, class TO, int L, int NODF> - void DiffusionMultiShellQballReconstructionImageFilter - ::PrintSelf(std::ostream& os, Indent indent) const - { +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; + os << indent << "GradientDirectionContainer: " + << m_GradientDirectionContainer << std::endl; } else { - os << indent << - "GradientDirectionContainer: (Gradient directions not set)" << std::endl; + 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 a2529bd31c..06a65027f1 100644 --- a/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h +++ b/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h @@ -1,174 +1,194 @@ /*========================================================================= 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(const GradientIndexMap * gradientIdexMap) + { + m_GradientIndexMap = gradientIdexMap; + } /** 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 ); } OdfPixelType Normalize(OdfPixelType odf, 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 * 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(); 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 ); + + private: OdfReconstructionMatrixType m_ReconstructionMatrix; OdfReconstructionMatrixType m_CoeffReconstructionMatrix; OdfReconstructionMatrixType m_SphericalHarmonicBasisMatrix; /** container to hold gradient directions */ GradientDirectionContainerType::Pointer m_GradientDirectionContainer; /** Number of gradient measurements */ unsigned int m_NumberOfGradientDirections; /** Number of baseline images */ unsigned int m_NumberOfBaselineImages; /** Threshold on the reference image data */ ReferencePixelType m_Threshold; /** LeBihan's b-value for normalizing tensors */ TOdfPixelType m_BValue; typename BZeroImageType::Pointer m_BZeroImage; + GradientIndexMap * m_GradientIndexMap; + double m_Lambda; bool m_IsHemisphericalArrangementOfGradientDirections; int m_NumberCoefficients; template< class VNLType > void printMatrix( VNLType * mat ); + enum GradientDirections + { + BZero = 0, + AllGradients = -1 + }; + }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkDiffusionMultiShellQballReconstructionImageFilter.cpp" #endif #endif //__itkDiffusionMultiShellQballReconstructionImageFilter_h_ diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkQBallReconstructionView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkQBallReconstructionView.cpp index 5bb914af18..30b85538e8 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkQBallReconstructionView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkQBallReconstructionView.cpp @@ -1,993 +1,999 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Module: $RCSfile$ Language: C++ Date: $Date: 2009-05-28 17:19:30 +0200 (Do, 28 Mai 2009) $ Version: $Revision: 17495 $ 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. =========================================================================*/ //#define MBILOG_ENABLE_DEBUG #include "QmitkQBallReconstructionView.h" #include "mitkDiffusionImagingConfigure.h" // qt includes #include // itk includes #include "itkTimeProbe.h" // mitk includes #include "mitkProgressBar.h" #include "mitkStatusBar.h" #include "mitkNodePredicateDataType.h" #include "QmitkDataStorageComboBox.h" #include "QmitkStdMultiWidget.h" #include "itkDiffusionQballReconstructionImageFilter.h" #include "itkAnalyticalDiffusionQballReconstructionImageFilter.h" #include "itkDiffusionMultiShellQballReconstructionImageFilter.h" #include "itkVectorContainer.h" #include "mitkQBallImage.h" #include "mitkProperties.h" #include "mitkVtkResliceInterpolationProperty.h" #include "mitkLookupTable.h" #include "mitkLookupTableProperty.h" #include "mitkTransferFunction.h" #include "mitkTransferFunctionProperty.h" #include "mitkDataNodeObject.h" #include "mitkOdfNormalizationMethodProperty.h" #include "mitkOdfScaleByProperty.h" #include "berryIStructuredSelection.h" #include "berryIWorkbenchWindow.h" #include "berryISelectionService.h" #include + const std::string QmitkQBallReconstructionView::VIEW_ID = "org.mitk.views.qballreconstruction"; #define DI_INFO MITK_INFO("DiffusionImaging") typedef float TTensorPixelType; const int QmitkQBallReconstructionView::nrconvkernels = 252; using namespace berry; struct QbrSelListener : ISelectionListener { berryObjectMacro(QbrSelListener); QbrSelListener(QmitkQBallReconstructionView* view) { m_View = view; } void DoSelectionChanged(ISelection::ConstPointer selection) { // save current selection in member variable m_View->m_CurrentSelection = selection.Cast(); // do something with the selected items if(m_View->m_CurrentSelection) { bool foundDwiVolume = false; // iterate selection for (IStructuredSelection::iterator i = m_View->m_CurrentSelection->Begin(); i != m_View->m_CurrentSelection->End(); ++i) { // extract datatree node if (mitk::DataNodeObject::Pointer nodeObj = i->Cast()) { mitk::DataNode::Pointer node = nodeObj->GetDataNode(); // only look at interesting types if(QString("DiffusionImage").compare(node->GetData()->GetNameOfClass())==0) { foundDwiVolume = true; } } } m_View->m_Controls->m_ButtonStandard->setEnabled(foundDwiVolume); m_View->m_Controls->m_groupQShellSelecton->setEnabled(foundDwiVolume); } } void SelectionChanged(IWorkbenchPart::Pointer part, ISelection::ConstPointer selection) { // check, if selection comes from datamanager if (part) { QString partname(part->GetPartName().c_str()); if(partname.compare("Datamanager")==0) { // apply selection DoSelectionChanged(selection); } } } QmitkQBallReconstructionView* m_View; }; QmitkQBallReconstructionView::QmitkQBallReconstructionView() : QmitkFunctionality(), m_Controls(NULL), m_MultiWidget(NULL) { } QmitkQBallReconstructionView::QmitkQBallReconstructionView(const QmitkQBallReconstructionView& other) { Q_UNUSED(other); throw std::runtime_error("Copy constructor not implemented"); } /*//void QmitkQBallReconstructionView::OpactiyChanged(int value) //{ // if (m_CurrentSelection) // { // if (mitk::DataNodeObject::Pointer nodeObj = m_CurrentSelection->Begin()->Cast()) // { // mitk::DataNode::Pointer node = nodeObj->GetDataNode(); // if(QString("DiffusionImage").compare(node->GetData()->GetNameOfClass())==0) // { // node->SetIntProperty("DisplayChannel", value); // mitk::RenderingManager::GetInstance()->RequestUpdateAll(); // } // } // } //} // //void QmitkQBallReconstructionView::OpactiyActionChanged() //{ // if (m_CurrentSelection) // { // if (mitk::DataNodeObject::Pointer nodeObj = m_CurrentSelection->Begin()->Cast()) // { // mitk::DataNode::Pointer node = nodeObj->GetDataNode(); // if(QString("DiffusionImage").compare(node->GetData()->GetNameOfClass())==0) // { // int displayChannel = 0.0; // if(node->GetIntProperty("DisplayChannel", displayChannel)) // { // m_OpacitySlider->setValue(displayChannel); // } // } // } // } // // MITK_INFO << "changed"; //} */ QmitkQBallReconstructionView::~QmitkQBallReconstructionView() { this->GetSite()->GetWorkbenchWindow()->GetSelectionService()->RemovePostSelectionListener(/*"org.mitk.views.datamanager",*/ m_SelListener); } void QmitkQBallReconstructionView::CreateQtPartControl(QWidget *parent) { if (!m_Controls) { // create GUI widgets m_Controls = new Ui::QmitkQBallReconstructionViewControls; m_Controls->setupUi(parent); this->CreateConnections(); QStringList items; items << "2" << "4" << "6" << "8" << "10" << "12"; m_Controls->m_QBallReconstructionMaxLLevelComboBox->addItems(items); m_Controls->m_QBallReconstructionMaxLLevelComboBox->setCurrentIndex(1); MethodChoosen(m_Controls->m_QBallReconstructionMethodComboBox->currentIndex()); #ifndef DIFFUSION_IMAGING_EXTENDED m_Controls->m_QBallReconstructionMethodComboBox->removeItem(3); #endif AdvancedCheckboxClicked(); // define data type for combobox //m_Controls->m_ImageSelector->SetDataStorage( this->GetDefaultDataStorage() ); //m_Controls->m_ImageSelector->SetPredicate( mitk::NodePredicateDataType::New("DiffusionImage") ); } m_SelListener = berry::ISelectionListener::Pointer(new QbrSelListener(this)); this->GetSite()->GetWorkbenchWindow()->GetSelectionService()->AddPostSelectionListener(/*"org.mitk.views.datamanager",*/ m_SelListener); berry::ISelection::ConstPointer sel( this->GetSite()->GetWorkbenchWindow()->GetSelectionService()->GetSelection("org.mitk.views.datamanager")); m_CurrentSelection = sel.Cast(); m_SelListener.Cast()->DoSelectionChanged(sel); } void QmitkQBallReconstructionView::StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget) { m_MultiWidget = &stdMultiWidget; } void QmitkQBallReconstructionView::StdMultiWidgetNotAvailable() { m_MultiWidget = NULL; } void QmitkQBallReconstructionView::CreateConnections() { if ( m_Controls ) { connect( (QObject*)(m_Controls->m_ButtonStandard), SIGNAL(clicked()), this, SLOT(ReconstructStandard()) ); connect( (QObject*)(m_Controls->m_AdvancedCheckbox), SIGNAL(clicked()), this, SLOT(AdvancedCheckboxClicked()) ); connect( (QObject*)(m_Controls->m_QBallReconstructionMethodComboBox), SIGNAL(currentIndexChanged(int)), this, SLOT(MethodChoosen(int)) ); } } void QmitkQBallReconstructionView::Activated() { QmitkFunctionality::Activated(); berry::ISelection::ConstPointer sel( this->GetSite()->GetWorkbenchWindow()->GetSelectionService()->GetSelection("org.mitk.views.datamanager")); m_CurrentSelection = sel.Cast(); m_SelListener.Cast()->DoSelectionChanged(sel); } void QmitkQBallReconstructionView::Deactivated() { QmitkFunctionality::Deactivated(); } void QmitkQBallReconstructionView::ReconstructStandard() { int index = m_Controls->m_QBallReconstructionMethodComboBox->currentIndex(); #ifndef DIFFUSION_IMAGING_EXTENDED if(index>=3) { index = index + 1; } #endif switch(index) { case 0: { // Numerical Reconstruct(0,0); break; } case 1: { // Standard Reconstruct(1,0); break; } case 2: { // Solid Angle Reconstruct(1,6); break; } case 3: { // Constrained Solid Angle Reconstruct(1,7); break; } case 4: { // ADC Reconstruct(1,4); break; } case 5: { // Raw Signal Reconstruct(1,5); break; } case 6: { // Multi Shell Reconstruct(2,7); break; } } } void QmitkQBallReconstructionView::MethodChoosen(int method) { #ifndef DIFFUSION_IMAGING_EXTENDED if(method>=3) { method = method + 1; } #endif m_Controls->m_groupQShellSelecton->setVisible(false); switch(method) { case 0: m_Controls->m_Description->setText("Numerical recon. (Tuch2004)"); break; case 1: m_Controls->m_Description->setText("Spherical harmonics recon. (Descoteaux2007)"); break; case 2: m_Controls->m_Description->setText("SH recon. with solid angle consideration (Aganj2009)"); break; case 3: m_Controls->m_Description->setText("SH solid angle with non-neg. constraint (Goh2009)"); break; case 4: m_Controls->m_Description->setText("SH recon. of the plain ADC-profiles"); break; case 5: m_Controls->m_Description->setText("SH recon. of the raw diffusion signal"); break; case 6: m_Controls->m_Description->setText("SH solid angle for multiple Q-Shell acquisitions"); if (m_CurrentSelection) { QVBoxLayout * layout = dynamic_cast(m_Controls->m_groupQShellSelecton->layout()); for(int i =0; i < layout->count(); i++){ layout->removeItem(layout->itemAt(i)); } for (IStructuredSelection::iterator i = m_CurrentSelection->Begin(); i != m_CurrentSelection->End(); ++i) { if (mitk::DataNodeObject::Pointer nodeObj = i->Cast()) { mitk::DataNode::Pointer node = nodeObj->GetDataNode(); - std::vector * bValues = new std::vector(); - std::vector::iterator it; + + typedef itk::DiffusionMultiShellQballReconstructionImageFilter::GradientIndexMap GradientIndexMap; + + GradientIndexMap * bValueMap = new GradientIndexMap(); + GradientIndexMap::iterator it; if(QString("DiffusionImage").compare(node->GetData()->GetNameOfClass())==0) { mitk::DiffusionImage* vols = static_cast*>(node->GetData()); - vols->GetBvalueList(bValues); + vols->GetBvalueList(bValueMap); QString str = QString((dynamic_cast(node->GetPropertyList()->GetProperty("name")))->GetValue()); QLabel *label = new QLabel("Image: " + str ); layout->addWidget(label); - for(it = bValues->begin(); it != bValues->end(); it ++){ + for(it = bValueMap->begin(); it != bValueMap->end(); it ++){ + if(it->first == 0) continue; QString checkboxName ; - checkboxName.setNum(*it); + checkboxName.setNum(it->first); QCheckBox *bValueCheckBox = new QCheckBox( "B-Value " + checkboxName ); + bValueCheckBox->setChecked(true); layout->addWidget(bValueCheckBox); } layout->addStretch(1); } } } } m_Controls->m_groupQShellSelecton->setVisible(true); break; } } void QmitkQBallReconstructionView::AdvancedCheckboxClicked() { bool check = m_Controls-> m_AdvancedCheckbox->isChecked(); m_Controls->m_QBallReconstructionMaxLLevelTextLabel_2->setVisible(check); m_Controls->m_QBallReconstructionMaxLLevelComboBox->setVisible(check); m_Controls->m_QBallReconstructionLambdaTextLabel_2->setVisible(check); m_Controls->m_QBallReconstructionLambdaLineEdit->setVisible(check); m_Controls->m_QBallReconstructionThresholdLabel_2->setVisible(check); m_Controls->m_QBallReconstructionThreasholdEdit->setVisible(check); m_Controls->m_OutputB0Image->setVisible(check); m_Controls->label_2->setVisible(check); //m_Controls->textLabel1_2->setVisible(check); //m_Controls->m_QBallReconstructionLambdaStepLineEdit->setVisible(check); //m_Controls->textLabel1_3->setVisible(check); m_Controls->frame_2->setVisible(check); } void QmitkQBallReconstructionView::Reconstruct(int method, int normalization) { if (m_CurrentSelection) { mitk::DataStorage::SetOfObjects::Pointer set = mitk::DataStorage::SetOfObjects::New(); int at = 0; for (IStructuredSelection::iterator i = m_CurrentSelection->Begin(); i != m_CurrentSelection->End(); ++i) { if (mitk::DataNodeObject::Pointer nodeObj = i->Cast()) { mitk::DataNode::Pointer node = nodeObj->GetDataNode(); if(QString("DiffusionImage").compare(node->GetData()->GetNameOfClass())==0) { set->InsertElement(at++, node); } } } if(method == 0) { NumericalQBallReconstruction(set, normalization); } else { #if BOOST_VERSION / 100000 > 0 #if BOOST_VERSION / 100 % 1000 > 34 if(method == 1) { AnalyticalQBallReconstruction(set, normalization); } if(method == 2) { MultiShellQBallReconstruction(set, normalization); } #else std::cout << "ERROR: Boost 1.35 minimum required" << std::endl; QMessageBox::warning(NULL,"ERROR","Boost 1.35 minimum required"); #endif #else std::cout << "ERROR: Boost 1.35 minimum required" << std::endl; QMessageBox::warning(NULL,"ERROR","Boost 1.35 minimum required"); #endif } } } void QmitkQBallReconstructionView::NumericalQBallReconstruction (mitk::DataStorage::SetOfObjects::Pointer inImages, int normalization) { try { itk::TimeProbe clock; int nrFiles = inImages->size(); if (!nrFiles) return; QString status; mitk::ProgressBar::GetInstance()->AddStepsToDo(nrFiles); mitk::DataStorage::SetOfObjects::const_iterator itemiter( inImages->begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( inImages->end() ); std::vector nodes; while ( itemiter != itemiterend ) // for all items { mitk::DiffusionImage* vols = static_cast*>( (*itemiter)->GetData()); std::string nodename; (*itemiter)->GetStringProperty("name", nodename); ++itemiter; // QBALL RECONSTRUCTION clock.Start(); MBI_INFO << "QBall reconstruction "; mitk::StatusBar::GetInstance()->DisplayText(status.sprintf( "QBall reconstruction for %s", nodename.c_str()).toAscii()); typedef itk::DiffusionQballReconstructionImageFilter QballReconstructionImageFilterType; QballReconstructionImageFilterType::Pointer filter = QballReconstructionImageFilterType::New(); filter->SetGradientImage( vols->GetDirections(), vols->GetVectorImage() ); filter->SetBValue(vols->GetB_Value()); filter->SetThreshold( m_Controls->m_QBallReconstructionThreasholdEdit->text().toFloat() ); switch(normalization) { case 0: { filter->SetNormalizationMethod(QballReconstructionImageFilterType::QBR_STANDARD); break; } case 1: { filter->SetNormalizationMethod(QballReconstructionImageFilterType::QBR_B_ZERO_B_VALUE); break; } case 2: { filter->SetNormalizationMethod(QballReconstructionImageFilterType::QBR_B_ZERO); break; } case 3: { filter->SetNormalizationMethod(QballReconstructionImageFilterType::QBR_NONE); break; } default: { filter->SetNormalizationMethod(QballReconstructionImageFilterType::QBR_STANDARD); } } filter->Update(); clock.Stop(); MBI_DEBUG << "took " << clock.GetMeanTime() << "s." ; // ODFs TO DATATREE mitk::QBallImage::Pointer image = mitk::QBallImage::New(); image->InitializeByItk( filter->GetOutput() ); //image->SetImportVolume( filter->GetOutput()->GetBufferPointer(), 0, 0, mitk::Image::ImportMemoryManagementType::ManageMemory ); image->SetVolume( filter->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( image ); QString newname; newname = newname.append(nodename.c_str()); newname = newname.append("_QN%1").arg(normalization); SetDefaultNodeProperties(node, newname.toStdString()); nodes.push_back(node); // B-Zero TO DATATREE if(m_Controls->m_OutputB0Image->isChecked()) { mitk::Image::Pointer image4 = mitk::Image::New(); image4->InitializeByItk( filter->GetBZeroImage().GetPointer() ); image4->SetVolume( filter->GetBZeroImage()->GetBufferPointer() ); mitk::DataNode::Pointer node4=mitk::DataNode::New(); node4->SetData( image4 ); node4->SetProperty( "name", mitk::StringProperty::New( QString(nodename.c_str()).append("_b0").toStdString()) ); nodes.push_back(node4); } mitk::ProgressBar::GetInstance()->Progress(); } std::vector::iterator nodeIt; for(nodeIt = nodes.begin(); nodeIt != nodes.end(); ++nodeIt) GetDefaultDataStorage()->Add(*nodeIt); mitk::StatusBar::GetInstance()->DisplayText(status.sprintf("Finished Processing %d Files", nrFiles).toAscii()); m_MultiWidget->RequestUpdate(); } catch (itk::ExceptionObject &ex) { MBI_INFO << ex ; return ; } } void QmitkQBallReconstructionView::AnalyticalQBallReconstruction( mitk::DataStorage::SetOfObjects::Pointer inImages, int normalization) { try { itk::TimeProbe clock; int nrFiles = inImages->size(); if (!nrFiles) return; std::vector lambdas; float minLambda = m_Controls->m_QBallReconstructionLambdaLineEdit->text().toFloat(); lambdas.push_back(minLambda); int nLambdas = lambdas.size(); QString status; mitk::ProgressBar::GetInstance()->AddStepsToDo(nrFiles*nLambdas); mitk::DataStorage::SetOfObjects::const_iterator itemiter( inImages->begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( inImages->end() ); std::vector* nodes = new std::vector(); while ( itemiter != itemiterend ) // for all items { mitk::DiffusionImage* vols = static_cast*>( (*itemiter)->GetData()); std::string nodename; (*itemiter)->GetStringProperty("name",nodename); itemiter++; // QBALL RECONSTRUCTION clock.Start(); MBI_INFO << "QBall reconstruction "; mitk::StatusBar::GetInstance()->DisplayText(status.sprintf( "QBall reconstruction for %s", nodename.c_str()).toAscii()); for(int i=0; im_QBallReconstructionMaxLLevelComboBox->currentIndex()) { case 0: { TemplatedAnalyticalQBallReconstruction<2>(vols, currentLambda, nodename, nodes, normalization); break; } case 1: { TemplatedAnalyticalQBallReconstruction<4>(vols, currentLambda, nodename, nodes, normalization); break; } case 2: { TemplatedAnalyticalQBallReconstruction<6>(vols, currentLambda, nodename, nodes, normalization); break; } case 3: { TemplatedAnalyticalQBallReconstruction<8>(vols, currentLambda, nodename, nodes, normalization); break; } case 4: { TemplatedAnalyticalQBallReconstruction<10>(vols, currentLambda, nodename, nodes, normalization); break; } case 5: { TemplatedAnalyticalQBallReconstruction<12>(vols, currentLambda, nodename, nodes, normalization); break; } } clock.Stop(); MBI_DEBUG << "took " << clock.GetMeanTime() << "s." ; mitk::ProgressBar::GetInstance()->Progress(); } } std::vector::iterator nodeIt; for(nodeIt = nodes->begin(); nodeIt != nodes->end(); ++nodeIt) GetDefaultDataStorage()->Add(*nodeIt); m_MultiWidget->RequestUpdate(); mitk::StatusBar::GetInstance()->DisplayText(status.sprintf("Finished Processing %d Files", nrFiles).toAscii()); } catch (itk::ExceptionObject &ex) { MBI_INFO << ex ; return ; } } template void QmitkQBallReconstructionView::TemplatedAnalyticalQBallReconstruction( mitk::DiffusionImage* vols, float lambda, std::string nodename, std::vector* nodes, int normalization) { typedef itk::AnalyticalDiffusionQballReconstructionImageFilter FilterType; typename FilterType::Pointer filter = FilterType::New(); filter->SetGradientImage( vols->GetDirections(), vols->GetVectorImage() ); filter->SetBValue(vols->GetB_Value()); filter->SetThreshold( m_Controls->m_QBallReconstructionThreasholdEdit->text().toFloat() ); filter->SetLambda(lambda); switch(normalization) { case 0: { filter->SetNormalizationMethod(FilterType::QBAR_STANDARD); break; } case 1: { filter->SetNormalizationMethod(FilterType::QBAR_B_ZERO_B_VALUE); break; } case 2: { filter->SetNormalizationMethod(FilterType::QBAR_B_ZERO); break; } case 3: { filter->SetNormalizationMethod(FilterType::QBAR_NONE); break; } case 4: { filter->SetNormalizationMethod(FilterType::QBAR_ADC_ONLY); break; } case 5: { filter->SetNormalizationMethod(FilterType::QBAR_RAW_SIGNAL); break; } case 6: { filter->SetNormalizationMethod(FilterType::QBAR_SOLID_ANGLE); break; } case 7: { filter->SetNormalizationMethod(FilterType::QBAR_NONNEG_SOLID_ANGLE); break; } default: { filter->SetNormalizationMethod(FilterType::QBAR_STANDARD); } } filter->Update(); // ODFs TO DATATREE mitk::QBallImage::Pointer image = mitk::QBallImage::New(); image->InitializeByItk( filter->GetOutput() ); image->SetVolume( filter->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( image ); QString newname; newname = newname.append(nodename.c_str()); newname = newname.append("_QA%1").arg(normalization); SetDefaultNodeProperties(node, newname.toStdString()); nodes->push_back(node); // mitk::Image::Pointer image5 = mitk::Image::New(); // image5->InitializeByItk( filter->GetODFSumImage().GetPointer() ); // image5->SetVolume( filter->GetODFSumImage()->GetBufferPointer() ); // mitk::DataNode::Pointer node5=mitk::DataNode::New(); // node5->SetData( image5 ); // node5->SetProperty( "name", mitk::StringProperty::New( // QString(nodename.c_str()).append("_ODF").toStdString()) ); // nodes->push_back(node5); // B-Zero TO DATATREE if(m_Controls->m_OutputB0Image->isChecked()) { mitk::Image::Pointer image4 = mitk::Image::New(); image4->InitializeByItk( filter->GetBZeroImage().GetPointer() ); image4->SetVolume( filter->GetBZeroImage()->GetBufferPointer() ); mitk::DataNode::Pointer node4=mitk::DataNode::New(); node4->SetData( image4 ); node4->SetProperty( "name", mitk::StringProperty::New( QString(nodename.c_str()).append("_b0").toStdString()) ); nodes->push_back(node4); } } void QmitkQBallReconstructionView::SetDefaultNodeProperties(mitk::DataNode::Pointer node, std::string name) { node->SetProperty( "ShowMaxNumber", mitk::IntProperty::New( 500 ) ); node->SetProperty( "Scaling", mitk::FloatProperty::New( 1.0 ) ); node->SetProperty( "Normalization", mitk::OdfNormalizationMethodProperty::New()); node->SetProperty( "ScaleBy", mitk::OdfScaleByProperty::New()); node->SetProperty( "IndexParam1", mitk::FloatProperty::New(2)); node->SetProperty( "IndexParam2", mitk::FloatProperty::New(1)); node->SetProperty( "visible", mitk::BoolProperty::New( true ) ); node->SetProperty( "VisibleOdfs", mitk::BoolProperty::New( false ) ); node->SetProperty ("layer", mitk::IntProperty::New(100)); node->SetProperty( "DoRefresh", mitk::BoolProperty::New( true ) ); //node->SetProperty( "opacity", mitk::FloatProperty::New(1.0f) ); node->SetProperty( "name", mitk::StringProperty::New(name) ); } /*//node->SetProperty( "volumerendering", mitk::BoolProperty::New( false ) ); //node->SetProperty( "use color", mitk::BoolProperty::New( true ) ); //node->SetProperty( "texture interpolation", mitk::BoolProperty::New( true ) ); //node->SetProperty( "reslice interpolation", mitk::VtkResliceInterpolationProperty::New() ); //node->SetProperty( "layer", mitk::IntProperty::New(0)); //node->SetProperty( "in plane resample extent by geometry", mitk::BoolProperty::New( false ) ); //node->SetOpacity(1.0f); //node->SetColor(1.0,1.0,1.0); //node->SetVisibility(true); //node->SetProperty( "IsQBallVolume", mitk::BoolProperty::New( true ) ); //mitk::LevelWindowProperty::Pointer levWinProp = mitk::LevelWindowProperty::New(); //mitk::LevelWindow levelwindow; //// levelwindow.SetAuto( image ); //levWinProp->SetLevelWindow( levelwindow ); //node->GetPropertyList()->SetProperty( "levelwindow", levWinProp ); //// add a default rainbow lookup table for color mapping //if(!node->GetProperty("LookupTable")) //{ // mitk::LookupTable::Pointer mitkLut = mitk::LookupTable::New(); // vtkLookupTable* vtkLut = mitkLut->GetVtkLookupTable(); // vtkLut->SetHueRange(0.6667, 0.0); // vtkLut->SetTableRange(0.0, 20.0); // vtkLut->Build(); // mitk::LookupTableProperty::Pointer mitkLutProp = mitk::LookupTableProperty::New(); // mitkLutProp->SetLookupTable(mitkLut); // node->SetProperty( "LookupTable", mitkLutProp ); //} //if(!node->GetProperty("binary")) // node->SetProperty( "binary", mitk::BoolProperty::New( false ) ); //// add a default transfer function //mitk::TransferFunction::Pointer tf = mitk::TransferFunction::New(); //node->SetProperty ( "TransferFunction", mitk::TransferFunctionProperty::New ( tf.GetPointer() ) ); //// set foldername as string property //mitk::StringProperty::Pointer nameProp = mitk::StringProperty::New( name ); //node->SetProperty( "name", nameProp );*/ void QmitkQBallReconstructionView::MultiShellQBallReconstruction(mitk::DataStorage::SetOfObjects::Pointer inImages, int normalization) { try { itk::TimeProbe clock; int nrFiles = inImages->size(); if (!nrFiles) return; std::vector lambdas; float minLambda = m_Controls->m_QBallReconstructionLambdaLineEdit->text().toFloat(); lambdas.push_back(minLambda); int nLambdas = lambdas.size(); QString status; mitk::ProgressBar::GetInstance()->AddStepsToDo(nrFiles*nLambdas); mitk::DataStorage::SetOfObjects::const_iterator itemiter( inImages->begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( inImages->end() ); std::vector* nodes = new std::vector(); while ( itemiter != itemiterend ) // for all items { mitk::DiffusionImage* vols = static_cast*>( (*itemiter)->GetData()); std::string nodename; (*itemiter)->GetStringProperty("name",nodename); itemiter++; // QBALL RECONSTRUCTION clock.Start(); MBI_INFO << "QBall reconstruction "; mitk::StatusBar::GetInstance()->DisplayText(status.sprintf( "QBall reconstruction for %s", nodename.c_str()).toAscii()); for(int i=0; im_QBallReconstructionMaxLLevelComboBox->currentIndex()) { case 0: { TemplatedMultiShellQBallReconstruction<2>(vols, currentLambda, nodename, nodes, normalization); break; } case 1: { TemplatedMultiShellQBallReconstruction<4>(vols, currentLambda, nodename, nodes, normalization); break; } case 2: { TemplatedMultiShellQBallReconstruction<6>(vols, currentLambda, nodename, nodes, normalization); break; } case 3: { TemplatedMultiShellQBallReconstruction<8>(vols, currentLambda, nodename, nodes, normalization); break; } case 4: { TemplatedMultiShellQBallReconstruction<10>(vols, currentLambda, nodename, nodes, normalization); break; } case 5: { TemplatedMultiShellQBallReconstruction<12>(vols, currentLambda, nodename, nodes, normalization); break; } } clock.Stop(); MBI_DEBUG << "took " << clock.GetMeanTime() << "s." ; mitk::ProgressBar::GetInstance()->Progress(); } } std::vector::iterator nodeIt; for(nodeIt = nodes->begin(); nodeIt != nodes->end(); ++nodeIt) GetDefaultDataStorage()->Add(*nodeIt); m_MultiWidget->RequestUpdate(); mitk::StatusBar::GetInstance()->DisplayText(status.sprintf("Finished Processing %d Files", nrFiles).toAscii()); } catch (itk::ExceptionObject &ex) { MBI_INFO << ex ; return ; } } template void QmitkQBallReconstructionView::TemplatedMultiShellQBallReconstruction(mitk::DiffusionImage* vols, float lambda, std::string nodename, std::vector* nodes, int normalization) { typedef itk::DiffusionMultiShellQballReconstructionImageFilter FilterType; typename FilterType::Pointer filter = FilterType::New(); filter->SetGradientImage( vols->GetDirections(), vols->GetVectorImage() ); filter->SetBValue(vols->GetB_Value()); filter->SetThreshold( m_Controls->m_QBallReconstructionThreasholdEdit->text().toFloat() ); filter->SetLambda(lambda); filter->Update(); // ODFs TO DATATREE mitk::QBallImage::Pointer image = mitk::QBallImage::New(); image->InitializeByItk( filter->GetOutput() ); image->SetVolume( filter->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( image ); QString newname; newname = newname.append(nodename.c_str()); newname = newname.append("_QA%1").arg(normalization); SetDefaultNodeProperties(node, newname.toStdString()); nodes->push_back(node); // B-Zero TO DATATREE if(m_Controls->m_OutputB0Image->isChecked()) { mitk::Image::Pointer image4 = mitk::Image::New(); image4->InitializeByItk( filter->GetBZeroImage().GetPointer() ); image4->SetVolume( filter->GetBZeroImage()->GetBufferPointer() ); mitk::DataNode::Pointer node4=mitk::DataNode::New(); node4->SetData( image4 ); node4->SetProperty( "name", mitk::StringProperty::New( QString(nodename.c_str()).append("_b0").toStdString()) ); nodes->push_back(node4); } }