diff --git a/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp b/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp index 4142e043ca..925c7b2538 100644 --- a/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp +++ b/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp @@ -1,753 +1,753 @@ /*========================================================================= 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 __itkAnalyticalDiffusionQballReconstructionImageFilter_cpp #define __itkAnalyticalDiffusionQballReconstructionImageFilter_cpp #include #include #include #include #include #include #include #include "itkPointShell.h" #include "mitkSphericalHarmonicsFunctions.h" namespace itk { #define QBALL_ANAL_RECON_PI M_PI template< class T, class TG, class TO, int L, int NODF> AnalyticalDiffusionQballReconstructionImageFilter ::AnalyticalDiffusionQballReconstructionImageFilter() : m_GradientDirectionContainer(NULL), m_NumberOfGradientDirections(0), m_NumberOfBaselineImages(1), m_Threshold(NumericTraits< ReferencePixelType >::NonpositiveMin()), m_BValue(1.0), m_Lambda(0.0), m_DirectionsDuplicated(false), m_Delta1(0.001), m_Delta2(0.001) { // At least 1 inputs is necessary for a vector image. // For images added one at a time we need at least six this->SetNumberOfRequiredInputs( 1 ); } template< class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NOrderL, int NrOdfDirections> typename itk::AnalyticalDiffusionQballReconstructionImageFilter< TReferenceImagePixelType,TGradientImagePixelType,TOdfPixelType, NOrderL,NrOdfDirections>::OdfPixelType itk::AnalyticalDiffusionQballReconstructionImageFilter ::Normalize( OdfPixelType odf, typename NumericTraits::AccumulateType b0 ) { switch( m_NormalizationMethod ) { case QBAR_STANDARD: { TOdfPixelType sum = 0; for(int i=0; i0) odf /= sum; return odf; break; } case QBAR_B_ZERO_B_VALUE: { for(int i=0; i vnl_vector itk::AnalyticalDiffusionQballReconstructionImageFilter ::PreNormalize( vnl_vector vec, typename NumericTraits::AccumulateType b0 ) { switch( m_NormalizationMethod ) { case QBAR_STANDARD: { return vec; break; } case QBAR_B_ZERO_B_VALUE: { int n = vec.size(); for(int i=0; i=1) vec[i] = 1-m_Delta2/2; else if (vec[i]>=1-m_Delta2) vec[i] = 1-m_Delta2/2-(1-vec[i])*(1-vec[i])/(2*m_Delta2); vec[i] = log(-log(vec[i])); } return vec; break; } } return vec; } template< class T, class TG, class TO, int L, int NODF> void AnalyticalDiffusionQballReconstructionImageFilter ::BeforeThreadedGenerateData() { // If we have more than 2 inputs, then each input, except the first is a // gradient image. The number of gradient images must match the number of // gradient directions. //const unsigned int numberOfInputs = this->GetNumberOfInputs(); // There need to be at least 6 gradient directions to be able to compute the // tensor basis if( m_NumberOfGradientDirections < 6 ) { itkExceptionMacro( << "At least 6 gradient directions are required" ); } // Input must be an itk::VectorImage. std::string gradientImageClassName( this->ProcessObject::GetInput(0)->GetNameOfClass()); if ( strcmp(gradientImageClassName.c_str(),"VectorImage") != 0 ) { itkExceptionMacro( << "There is only one Gradient image. I expect that to be a VectorImage. " << "But its of type: " << gradientImageClassName ); } this->ComputeReconstructionMatrix(); typename GradientImagesType::Pointer img = static_cast< GradientImagesType * >( this->ProcessObject::GetInput(0) ); m_BZeroImage = BZeroImageType::New(); m_BZeroImage->SetSpacing( img->GetSpacing() ); // Set the image spacing m_BZeroImage->SetOrigin( img->GetOrigin() ); // Set the image origin m_BZeroImage->SetDirection( img->GetDirection() ); // Set the image direction m_BZeroImage->SetLargestPossibleRegion( img->GetLargestPossibleRegion()); m_BZeroImage->SetBufferedRegion( img->GetLargestPossibleRegion() ); m_BZeroImage->Allocate(); m_ODFSumImage = BZeroImageType::New(); m_ODFSumImage->SetSpacing( img->GetSpacing() ); // Set the image spacing m_ODFSumImage->SetOrigin( img->GetOrigin() ); // Set the image origin m_ODFSumImage->SetDirection( img->GetDirection() ); // Set the image direction m_ODFSumImage->SetLargestPossibleRegion( img->GetLargestPossibleRegion()); m_ODFSumImage->SetBufferedRegion( img->GetLargestPossibleRegion() ); m_ODFSumImage->Allocate(); m_CoefficientImage = CoefficientImageType::New(); m_CoefficientImage->SetSpacing( img->GetSpacing() ); // Set the image spacing m_CoefficientImage->SetOrigin( img->GetOrigin() ); // Set the image origin m_CoefficientImage->SetDirection( img->GetDirection() ); // Set the image direction m_CoefficientImage->SetLargestPossibleRegion( img->GetLargestPossibleRegion()); m_CoefficientImage->SetBufferedRegion( img->GetLargestPossibleRegion() ); m_CoefficientImage->Allocate(); if(m_NormalizationMethod == QBAR_SOLID_ANGLE || m_NormalizationMethod == QBAR_NONNEG_SOLID_ANGLE) { m_Lambda = 0.0; } } template< class T, class TG, class TO, int L, int NODF> void AnalyticalDiffusionQballReconstructionImageFilter ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int ) { typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread); oit.GoToBegin(); ImageRegionIterator< BZeroImageType > oit2(m_BZeroImage, outputRegionForThread); oit2.GoToBegin(); ImageRegionIterator< FloatImageType > oit3(m_ODFSumImage, outputRegionForThread); oit3.GoToBegin(); ImageRegionIterator< CoefficientImageType > oit4(m_CoefficientImage, outputRegionForThread); oit4.GoToBegin(); typedef ImageRegionConstIterator< GradientImagesType > GradientIteratorType; typedef typename GradientImagesType::PixelType GradientVectorType; typename GradientImagesType::Pointer gradientImagePointer = NULL; // Would have liked a dynamic_cast here, but seems SGI doesn't like it // The enum will ensure that an inappropriate cast is not done gradientImagePointer = static_cast< GradientImagesType * >( this->ProcessObject::GetInput(0) ); GradientIteratorType git(gradientImagePointer, outputRegionForThread ); git.GoToBegin(); // Compute the indicies of the baseline images and gradient images std::vector baselineind; // contains the indicies of // the baseline images std::vector gradientind; // contains the indicies of // the gradient images for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin(); gdcit != this->m_GradientDirectionContainer->End(); ++gdcit) { if(gdcit.Value().one_norm() <= 0.0) { baselineind.push_back(gdcit.Index()); } else { gradientind.push_back(gdcit.Index()); } } if( m_DirectionsDuplicated ) { int gradIndSize = gradientind.size(); for(int i=0; i::AccumulateType b0 = NumericTraits::Zero; // Average the baseline image pixels for(unsigned int i = 0; i < baselineind.size(); ++i) { b0 += b[baselineind[i]]; } b0 /= this->m_NumberOfBaselineImages; OdfPixelType odf(0.0); typename CoefficientImageType::PixelType coeffPixel(0.0); vnl_vector B(m_NumberOfGradientDirections); if( (b0 != 0) && (b0 >= m_Threshold) ) { for( unsigned int i = 0; i< m_NumberOfGradientDirections; i++ ) { B[i] = static_cast(b[gradientind[i]]); } B = PreNormalize(B, b0); if(m_NormalizationMethod == QBAR_SOLID_ANGLE) { vnl_vector coeffs(m_NumberCoefficients); coeffs = ( (*m_CoeffReconstructionMatrix) * B ); coeffs[0] += 1.0/(2.0*sqrt(QBALL_ANAL_RECON_PI)); odf = ( (*m_SphericalHarmonicBasisMatrix) * coeffs ).data_block(); coeffPixel = coeffs.data_block(); } else if(m_NormalizationMethod == QBAR_NONNEG_SOLID_ANGLE) { /** this would be the place to implement a non-negative * solver for quadratic programming problem: * min .5*|| Bc-s ||^2 subject to -CLPc <= 4*pi*ones * (refer to MICCAI 2009 Goh et al. "Estimating ODFs with PDF constraints") * .5*|| Bc-s ||^2 == .5*c'B'Bc - x'B's + .5*s's */ itkExceptionMacro( << "Nonnegative Solid Angle not yet implemented"); } else { odf = ( (*m_ReconstructionMatrix) * B ).data_block(); } odf = Normalize(odf, b0); } oit.Set( odf ); oit2.Set( b0 ); float sum = 0; for (int k=0; k void AnalyticalDiffusionQballReconstructionImageFilter ::tofile2(vnl_matrix *pA, std::string fname) { vnl_matrix A = (*pA); ofstream myfile; std::locale C("C"); std::locale originalLocale = myfile.getloc(); myfile.imbue(C); myfile.open (fname.c_str()); myfile << "A1=["; for(int i=0; i void AnalyticalDiffusionQballReconstructionImageFilter ::ComputeReconstructionMatrix() { //for(int i=-6;i<7;i++) // std::cout << boost::math::legendre_p(6, i, 0.65657) << std::endl; if( m_NumberOfGradientDirections < 6 ) { itkExceptionMacro( << "Not enough gradient directions supplied. Need to supply at least 6" ); } { // check for duplicate diffusion gradients bool warning = false; for(GradientDirectionContainerType::ConstIterator gdcit1 = this->m_GradientDirectionContainer->Begin(); gdcit1 != this->m_GradientDirectionContainer->End(); ++gdcit1) { for(GradientDirectionContainerType::ConstIterator gdcit2 = this->m_GradientDirectionContainer->Begin(); gdcit2 != this->m_GradientDirectionContainer->End(); ++gdcit2) { if(gdcit1.Value() == gdcit2.Value() && gdcit1.Index() != gdcit2.Index()) { itkWarningMacro( << "Some of the Diffusion Gradients equal each other. Corresponding image data should be averaged before calling this filter." ); warning = true; break; } } if (warning) break; } // handle acquisition schemes where only half of the spherical // shell is sampled by the gradient directions. In this case, // each gradient direction is duplicated in negative direction. vnl_vector centerMass(3); centerMass.fill(0.0); int count = 0; for(GradientDirectionContainerType::ConstIterator gdcit1 = this->m_GradientDirectionContainer->Begin(); gdcit1 != this->m_GradientDirectionContainer->End(); ++gdcit1) { if(gdcit1.Value().one_norm() > 0.0) { centerMass += gdcit1.Value(); count ++; } } centerMass /= count; if(centerMass.two_norm() > 0.1) { m_DirectionsDuplicated = true; m_NumberOfGradientDirections *= 2; } } vnl_matrix *Q = new vnl_matrix(3, m_NumberOfGradientDirections); { int i = 0; for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin(); gdcit != this->m_GradientDirectionContainer->End(); ++gdcit) { if(gdcit.Value().one_norm() > 0.0) { double x = gdcit.Value().get(0); double y = gdcit.Value().get(1); double z = gdcit.Value().get(2); double cart[3]; - mitk::mitk_sh_functions::Cart2Sph(x,y,z,cart); + mitk::sh::Cart2Sph(x,y,z,cart); (*Q)(0,i) = cart[0]; (*Q)(1,i) = cart[1]; (*Q)(2,i++) = cart[2]; } } if(m_DirectionsDuplicated) { for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin(); gdcit != this->m_GradientDirectionContainer->End(); ++gdcit) { if(gdcit.Value().one_norm() > 0.0) { double x = gdcit.Value().get(0); double y = gdcit.Value().get(1); double z = gdcit.Value().get(2); double cart[3]; - mitk::mitk_sh_functions::Cart2Sph(x,y,z,cart); + mitk::sh::Cart2Sph(x,y,z,cart); (*Q)(0,i) = cart[0]; (*Q)(1,i) = cart[1]; (*Q)(2,i++) = cart[2]; } } } } int l = L; m_NumberCoefficients = (int)(l*l + l + 2.0)/2.0 + l; vnl_matrix* B = new vnl_matrix(m_NumberOfGradientDirections,m_NumberCoefficients); vnl_matrix* _L = new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients); _L->fill(0.0); vnl_matrix* LL = new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients); LL->fill(0.0); vnl_matrix* P = new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients); P->fill(0.0); vnl_matrix* Inv = new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients); P->fill(0.0); vnl_vector* lj = new vnl_vector(m_NumberCoefficients); m_LP = new vnl_vector(m_NumberCoefficients); for(unsigned int i=0; i(B->transpose()); //tofile2(&m_B_t,"m_B_t"); vnl_matrix B_t_B = (*m_B_t) * (*B); //tofile2(&B_t_B,"B_t_B"); vnl_matrix lambdaLL(m_NumberCoefficients,m_NumberCoefficients); lambdaLL.update((*LL)); lambdaLL *= m_Lambda; //tofile2(&lambdaLL,"lLL"); vnl_matrix tmp( B_t_B + lambdaLL); vnl_matrix_inverse *pseudoInverse = new vnl_matrix_inverse( tmp ); (*Inv) = pseudoInverse->pinverse(); //tofile2(Inv,"Inv"); vnl_matrix temp((*Inv) * (*m_B_t)); double fac1 = (1.0/(16.0*QBALL_ANAL_RECON_PI*QBALL_ANAL_RECON_PI)); switch(m_NormalizationMethod) { case QBAR_ADC_ONLY: case QBAR_RAW_SIGNAL: break; case QBAR_STANDARD: case QBAR_B_ZERO_B_VALUE: case QBAR_B_ZERO: case QBAR_NONE: temp = (*P) * temp; break; case QBAR_SOLID_ANGLE: temp = fac1 * (*P) * (*_L) * temp; break; case QBAR_NONNEG_SOLID_ANGLE: break; } //tofile2(&temp,"A"); m_CoeffReconstructionMatrix = new vnl_matrix(m_NumberCoefficients,m_NumberOfGradientDirections); for(int i=0; iodfs later int NOdfDirections = NODF; vnl_matrix_fixed* U = itk::PointShell >::DistributePointShell(); m_SphericalHarmonicBasisMatrix = new vnl_matrix(NOdfDirections,m_NumberCoefficients); vnl_matrix* sphericalHarmonicBasisMatrix2 = new vnl_matrix(NOdfDirections,m_NumberCoefficients); for(int i=0; i(NOdfDirections,m_NumberOfGradientDirections); *m_ReconstructionMatrix = (*m_SphericalHarmonicBasisMatrix) * (*m_CoeffReconstructionMatrix); } template< class T, class TG, class TO, int L, int NODF> void AnalyticalDiffusionQballReconstructionImageFilter ::SetGradientImage( GradientDirectionContainerType *gradientDirection, const GradientImagesType *gradientImage ) { this->m_GradientDirectionContainer = gradientDirection; unsigned int numImages = gradientDirection->Size(); this->m_NumberOfBaselineImages = 0; for(GradientDirectionContainerType::Iterator it = this->m_GradientDirectionContainer->Begin(); it != this->m_GradientDirectionContainer->End(); it++) { if(it.Value().one_norm() <= 0.0) { this->m_NumberOfBaselineImages++; } else // Normalize non-zero gradient directions { it.Value() = it.Value() / it.Value().two_norm(); } } this->m_NumberOfGradientDirections = numImages - this->m_NumberOfBaselineImages; // ensure that the gradient image we received has as many components as // the number of gradient directions if( gradientImage->GetVectorLength() != this->m_NumberOfBaselineImages + m_NumberOfGradientDirections ) { itkExceptionMacro( << m_NumberOfGradientDirections << " gradients + " << this->m_NumberOfBaselineImages << "baselines = " << m_NumberOfGradientDirections + this->m_NumberOfBaselineImages << " directions specified but image has " << gradientImage->GetVectorLength() << " components."); } this->ProcessObject::SetNthInput( 0, const_cast< GradientImagesType* >(gradientImage) ); } template< class T, class TG, class TO, int L, int NODF> void AnalyticalDiffusionQballReconstructionImageFilter ::PrintSelf(std::ostream& os, Indent indent) const { std::locale C("C"); std::locale originalLocale = os.getloc(); os.imbue(C); Superclass::PrintSelf(os,indent); os << indent << "OdfReconstructionMatrix: " << m_ReconstructionMatrix << std::endl; if ( m_GradientDirectionContainer ) { os << indent << "GradientDirectionContainer: " << m_GradientDirectionContainer << std::endl; } else { os << indent << "GradientDirectionContainer: (Gradient directions not set)" << std::endl; } os << indent << "NumberOfGradientDirections: " << m_NumberOfGradientDirections << std::endl; os << indent << "NumberOfBaselineImages: " << m_NumberOfBaselineImages << std::endl; os << indent << "Threshold for reference B0 image: " << m_Threshold << std::endl; os << indent << "BValue: " << m_BValue << std::endl; os.imbue( originalLocale ); } } #endif // __itkAnalyticalDiffusionQballReconstructionImageFilter_cpp diff --git a/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.cpp b/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.cpp index be46230b5d..8bac58cbcf 100644 --- a/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.cpp +++ b/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.cpp @@ -1,1137 +1,1135 @@ /*========================================================================= 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 "mitkVNLVectorFunctions.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(Mode_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 ); } template void DiffusionMultiShellQballReconstructionImageFilter ::Normalize( OdfPixelType & out) { for(int i=0; i double DiffusionMultiShellQballReconstructionImageFilter ::CalculateThreashold(const double value, const double delta) { return (value<0)*(0.5*delta) + (value>=0 && value=delta && value<1-delta)*value+(value>=1-delta && value<1)*(1-0.5*delta-0.5*((1-value)*(1-value))/delta) + (value>=1)*(1-0.5*delta); } template void DiffusionMultiShellQballReconstructionImageFilter ::Threshold(vnl_vector & vec, double delta) { if (delta==0){ //Clip attenuation values. If att<0 => att=0, if att>1 => att=1 for (int i=0; i=0 && vec[i]<=1)*vec[i]+(vec[i]>1); } else{ //Use function from Aganj et al, MRM, 2010 for (int i=0; i< vec.size(); i++) vec[i]=CalculateThreashold(vec[i], delta); } } template void DiffusionMultiShellQballReconstructionImageFilter ::Threshold(vnl_matrix & mat, double delta) { if (delta==0){ //Clip attenuation values. If att<0 => att=0, if att>1 => att=1 for (int i=0; i=0 && mat(i,j)<=1)*mat(i,j)+(mat(i,j)>1); } else{ //Use function from Aganj et al, MRM, 2010 for (int i=0; i void DiffusionMultiShellQballReconstructionImageFilter ::Projection1( vnl_matrix & E, double delta ) { const 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)); } } //T0 = -T0.apply(std::log); // 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 = a0 * 3.0; tb = b0 * 3.0; e = tb - (ta * 2.0); m = (tb * 2.0 ) + 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)*delta < e[i] && e[i]<-3*sF*delta; C(i,3) = m[i] >= 3-3*sF*delta && e[i] >= -3 *sF * delta; C(i,4) = 2.5 + 1.5*(5+sF)*delta < m[i] && m[i] < 3-3*sF*delta && e[i] > -3*sF*delta; C(i,5) = ta[i] <= 0.5+1.5 *(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(i,0)=(bool)C(i,0) * a0(i); A(i,1)=(bool)C(i,1) * (1.0/3.0-(sF+2)*delta); A(i,2)=(bool)C(i,2) * (0.2+0.8*a0(i)-0.4*b0(i)-delta/sF); A(i,3)=(bool)C(i,3) * (0.2+delta/sF); A(i,4)=(bool)C(i,4) * (0.2*a0(i)+0.4*b0(i)+2*delta/sF); A(i,5)=(bool)C(i,5) * (1.0/6.0+0.5*(sF+1)*delta); A(i,6)=(bool)C(i,6) * a0(i); B(i,0)=(bool)C(i,0) * (1.0/3.0+delta); B(i,1)=(bool)C(i,1) * (1.0/3.0+delta); B(i,2)=(bool)C(i,2) * (0.4-0.4*a0(i)+0.2*b0(i)-2*delta/sF); B(i,3)=(bool)C(i,3) * (0.4-3*delta/sF); B(i,4)=(bool)C(i,4) * (0.4*a0(i)+0.8*b0(i)-delta/sF); B(i,5)=(bool)C(i,5) * (1.0/3.0+delta); B(i,6)=(bool)C(i,6) * b0(i); } for(int i = 0 ; i < E.rows(); i++) { double sumA = 0; double sumB = 0; for(int j = 0 ; j < 7; 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 DiffusionMultiShellQballReconstructionImageFilter ::Projection2( vnl_vector & A, vnl_vector & a, vnl_vector & b, double delta0) { const double s6 = sqrt(6); const double s15 = s6/2.0; vnl_vector delta(a.size()); delta.fill(delta0); vnl_matrix AM(a.size(), 15); vnl_matrix aM(a.size(), 15); vnl_matrix bM(a.size(), 15); vnl_matrix B(a.size(), 15); AM.set_column(0, A); AM.set_column(1, A); AM.set_column(2, A); AM.set_column(3, delta); AM.set_column(4, (A+a-b - (delta*s6))/3.0); AM.set_column(5, delta); AM.set_column(6, delta); AM.set_column(7, delta); AM.set_column(8, A); AM.set_column(9, 0.2*(a*2+A-2*(s6+1)*delta)); AM.set_column(10,0.2*(b*(-2)+A+2-2*(s6+1)*delta)); AM.set_column(11, delta); AM.set_column(12, delta); AM.set_column(13, delta); AM.set_column(14, 0.5-(1+s15)*delta); aM.set_column(0, a); aM.set_column(1, a); aM.set_column(2, -delta + 1); aM.set_column(3, a); aM.set_column(4, (A*2+a*5+b+s6*delta)/6.0); aM.set_column(5, a); aM.set_column(6, -delta + 1); aM.set_column(7, 0.5*(a+b)+(1+s15)*delta); aM.set_column(8, -delta + 1); aM.set_column(9, 0.2*(a*4+A*2+(s6+1)*delta)); aM.set_column(10, -delta + 1); aM.set_column(11, (s6+3)*delta); aM.set_column(12, -delta + 1); aM.set_column(13, -delta + 1); aM.set_column(14, -delta + 1); bM.set_column(0, b); bM.set_column(1, delta); bM.set_column(2, b); bM.set_column(3, b); bM.set_column(4, (A*(-2)+a+b*5-s6*delta)/6.0); bM.set_column(5, delta); bM.set_column(6, b); bM.set_column(7, 0.5*(a+b)-(1+s15)*delta); bM.set_column(8, delta); bM.set_column(9, delta); bM.set_column(10, 0.2*(b*4-A*2+1-(s6+1)*delta)); bM.set_column(11, delta); bM.set_column(12, delta); bM.set_column(13, -delta*(s6+3) + 1); bM.set_column(14, delta); delta0 *= 0.99; for(int i = 0 ; i < a.size(); i ++) { for(int j = 0 ; j < 15; j ++) { B(i,j) = delta0 < AM(i,j) && 2 * (AM(i,j) + delta0 * s15) < aM(i,j) - bM(i,j) && bM(i,j) > delta0 && aM(i,j) < 1- delta0; } } vnl_matrix R2(a.size(), 15); vnl_matrix A_(a.size(), 15); vnl_matrix a_(a.size(), 15); vnl_matrix b_(a.size(), 15); vnl_matrix OnesVecMat(1, 15); OnesVecMat.fill(1.0); vnl_matrix AVecMat(a.size(), 1); AVecMat.set_column(0,A); vnl_matrix aVecMat(a.size(), 1); aVecMat.set_column(0,a); vnl_matrix bVecMat(a.size(), 1); bVecMat.set_column(0,b); A_ = AM - (AVecMat * OnesVecMat); a_ = aM - (aVecMat * OnesVecMat); b_ = bM - (bVecMat * OnesVecMat); for(int i = 0 ; i < a.size(); i++) for(int j = 0 ; j < 15; j++) { A_(i,j) *= A_(i,j); a_(i,j) *= a_(i,j); b_(i,j) *= b_(i,j); } R2 = A_ + a_ + b_; for(int i = 0 ; i < a.size(); i ++) { for(int j = 0 ; j < 15; j ++) { if(B(i,j) == 0) R2(i,j) = 1e20; } } std::vector indicies(a.size()); // suche den spalten-index der zu der kleinsten Zahl einer Zeile korrespondiert for(int i = 0 ; i < a.size(); i++) { unsigned int index = 0; double minvalue = 999; for(int j = 0 ; j < 15 ; j++) { if(R2(i,j) < minvalue){ minvalue = R2(i,j); index = j; } } indicies[i] = index; } for(int i = 0 ; i < a.size(); i++) { A[i] = AM(i,indicies[i]); a[i] = aM(i,indicies[i]); b[i] = bM(i,indicies[i]); } } template void 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 DiffusionMultiShellQballReconstructionImageFilter ::S_S0Normalization( vnl_matrix & mat, typename NumericTraits::AccumulateType b0 ) { 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])); } } template< class T, class TG, class TO, int L, int NODF> void DiffusionMultiShellQballReconstructionImageFilter ::SetGradientImage( GradientDirectionContainerType *gradientDirection , const GradientImagesType *gradientImage , float bvalue) { this->m_BValue = bvalue; this->m_GradientDirectionContainer = gradientDirection; this->m_NumberOfBaselineImages = 0; this->m_ReconstructionType = Mode_Standard1Shell; GradientDirectionContainerType::ConstIterator gdcit; for( gdcit = this->m_GradientDirectionContainer->Begin(); gdcit != this->m_GradientDirectionContainer->End(); ++gdcit) { double bValueKey = int(((m_BValue* gdcit.Value().two_norm() * gdcit.Value().two_norm())+7.5)/10)*10; m_GradientIndexMap[bValueKey].push_back(gdcit.Index()); } //if(listOfUserSelctedBValues.size() == 0){ // itkExceptionMacro(<< "DiffusionMultiShellQballReconstructionImageFilter.cpp : No list Of User Selcted B Values available"); //} if(m_GradientIndexMap.size() == 0){ itkExceptionMacro(<< "DiffusionMultiShellQballReconstructionImageFilter.cpp : no GradientIndexMapAvalible"); } //if(listOfUserSelctedBValues.size() != m_GradientIndexMap.size()){ // itkExceptionMacro(<< "DiffusionMultiShellQballReconstructionImageFilter.cpp : The number of user selected B Values != number of Image BValues"); //} if(m_GradientIndexMap.size() == 4){ GradientIndexMapIteraotr it = m_GradientIndexMap.begin(); it++; const int b1 = (*it).first; it++; const int b2 = (*it).first; it++; const int b3 = (*it).first; if(b2 - b1 == b1 && b3 - b2 == b1 ) { m_ReconstructionType = Mode_Analytical3Shells; } } if(m_GradientIndexMap.size() > 2 && m_ReconstructionType != Mode_Analytical3Shells) { m_ReconstructionType = Mode_NumericalNShells; } this->m_NumberOfBaselineImages = m_GradientIndexMap[0].size(); 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() { itk::TimeProbe clock; clock.Start(); 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(); this->ComputeReconstructionMatrix(); clock.Stop(); MITK_INFO << "Before GenerateData : " << clock.GetTotal(); } 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 ImageRegionConstIterator< BZeroImageType > bzeroImageIterator(m_BZeroImage, outputRegionForThread); bzeroImageIterator.GoToBegin(); 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 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(); typedef typename GradientImagesType::PixelType GradientVectorType; // iterate overall voxels of the gradient image region while( ! git.IsAtEnd() ) { GradientVectorType b = git.Get(); // ODF Vector OdfPixelType odf(0.0); // Create the Signal Vector vnl_vector SignalVector(m_NumberOfGradientDirections); if( (bzeroImageIterator.Get() != 0) && (bzeroImageIterator.Get() >= m_Threshold) ) { 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, bzeroImageIterator.Get()); DoubleLogarithm(SignalVector); // ODF coeffs-vector vnl_vector coeffs(m_NumberCoefficients); // approximate ODF coeffs coeffs = ( (*m_CoeffReconstructionMatrix) * SignalVector ); coeffs[0] = 1.0/(2.0*sqrt(QBALL_ANAL_RECON_PI)); odf = mitk::mitk_vnl_function::element_cast(( (*m_ODFSphericalHarmonicBasisMatrix) * coeffs )).data_block(); odf *= (QBALL_ANAL_RECON_PI*4/NODF); } // set ODF to ODF-Image oit.Set( odf ); ++oit; ++bzeroImageIterator; ++git; } MITK_INFO << "One Thread finished reconstruction"; } #include //#include //#include template< class T, class TG, class TO, int L, int NODF> void DiffusionMultiShellQballReconstructionImageFilter ::NumericalNShellReconstruction(const OutputImageRegionType& outputRegionForThread) { // vnl_levenberg_marquardt LMOptimizer = new vnl_levenberg_marquardt(); } template< class T, class TG, class TO, int L, int NODF> void DiffusionMultiShellQballReconstructionImageFilter ::GenerateAveragedBZeroImage(const OutputImageRegionType& outputRegionForThread) { typedef typename GradientImagesType::PixelType GradientVectorType; ImageRegionIterator< BZeroImageType > bzeroIterator(m_BZeroImage, outputRegionForThread); bzeroIterator.GoToBegin(); IndiciesVector BZeroIndicies = m_GradientIndexMap[0]; 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(); 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 /= BZeroIndicies.size(); bzeroIterator.Set(b0); ++bzeroIterator; ++git; } } template< class T, class TG, class TO, int L, int NODF> void DiffusionMultiShellQballReconstructionImageFilter ::AnalyticalThreeShellReconstruction(const OutputImageRegionType& outputRegionForThread) { typedef typename GradientImagesType::PixelType GradientVectorType; // Input Gradient Image and Output ODF Image typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); typename GradientImagesType::Pointer gradientImagePointer = static_cast< GradientImagesType * >( this->ProcessObject::GetInput(0) ); // Define Image iterators ImageRegionIterator< OutputImageType > odfOutputImageIterator(outputImage, outputRegionForThread); ImageRegionConstIterator< BZeroImageType > bzeroImageIterator(m_BZeroImage, outputRegionForThread); ImageRegionConstIterator< GradientImagesType > gradientInputImageIterator(gradientImagePointer, outputRegionForThread ); // All iterators seht to Begin of the specific OutputRegion odfOutputImageIterator.GoToBegin(); bzeroImageIterator.GoToBegin(); gradientInputImageIterator.GoToBegin(); // Get Shell Indicies for all non-BZero Gradients // it MUST be a arithmetic progression eg.: 1000, 2000, 3000 GradientIndexMapIteraotr it = m_GradientIndexMap.begin(); it++; // it = b-value = 1000 IndiciesVector Shell1Indiecies = (*it).second; it++; // it = b-value = 2000 IndiciesVector Shell2Indiecies = (*it).second; it++; // it = b-value = 3000 IndiciesVector Shell3Indiecies = (*it).second; // if input data is a hemispherical arragement, duplicate eache gradient for each shell 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]); } } // Nx3 Signal Vector with E(0) = Shell 1, E(1) = Shell 2, E(2) = Shell 3 vnl_matrix< double > * E = new vnl_matrix(Shell1Indiecies.size(), 3); vnl_vector * AlphaValues = new vnl_vector(Shell1Indiecies.size()); vnl_vector * BetaValues = new vnl_vector(Shell1Indiecies.size()); vnl_vector * LAValues = new vnl_vector(Shell1Indiecies.size()); vnl_vector * PValues = new vnl_vector(Shell1Indiecies.size()); OdfPixelType odf(0.0); // iterate overall voxels of the gradient image region while( ! gradientInputImageIterator.IsAtEnd() ) { if( (bzeroImageIterator.Get() != 0) && (bzeroImageIterator.Get() >= m_Threshold) ) { // Get the Signal-Value for each Shell at each direction (specified in the ShellIndicies Vector .. this direction corresponse to this shell...) GradientVectorType b = gradientInputImageIterator.Get(); for(int i = 0 ; i < Shell1Indiecies.size(); i++) { E->put(i,0, static_cast(b[Shell1Indiecies[i]])); E->put(i,1, static_cast(b[Shell2Indiecies[i]])); E->put(i,2, static_cast(b[Shell3Indiecies[i]])); } //Approximated-Signal by SH fit - using the specific shell directions and values // approximated Signal : S = SHBasis * Coeffs // with Coeffs : C = (B_T * B + lambda * L) ^ -1 * B_T * OS // OS := Original-Signal E->set_column(1, (*m_SHBasisMatrix) * ((*m_SignalReonstructionMatrix) * (E->get_column(1)))); E->set_column(2, (*m_SHBasisMatrix) * ((*m_SignalReonstructionMatrix) * (E->get_column(2)))); // Normalize the Signal: Si/S0 S_S0Normalization(*E,bzeroImageIterator.Get()); //Implements Eq. [19] and Fig. 4. Threshold(*E); //inqualities [31]. Taking the lograithm of th first tree inqualities //convert the quadratic inqualities to linear ones. Projection1(*E); double E1, E2, E3, P2,A,B2,B,P,alpha,beta,lambda, ER1, ER2; for( unsigned int i = 0; i< Shell1Indiecies.size(); i++ ) { E1 = E->get(i,0); E2 = E->get(i,1); E3 = E->get(i,2); P2 = E2-E1*E1; A = (E3 -E1*E2) / ( 2* P2); B2 = A * A -(E1 * E3 - E2 * E2) /P2; B = 0; if(B2 > 0) B = sqrt(B2); P = 0; if(P2 > 0) P = sqrt(P2); alpha = A + B; beta = A - B; lambda = 0.5 + 0.5 * std::sqrt(1 - std::pow((P * 2 ) / (alpha - beta), 2));; ER1 = std::fabs(lambda * (alpha - beta) + (beta - E1 )) + std::fabs(lambda * (std::pow(alpha, 2) - std::pow(beta, 2)) + (std::pow(beta, 2) - E2 )) + std::fabs(lambda * (std::pow(alpha, 3) - std::pow(beta, 3)) + (std::pow(beta, 3) - E3 )); ER2 = std::fabs((lambda-1) * (alpha - beta) + (beta - E1 )) + std::fabs((lambda-1) * (std::pow(alpha, 2) - std::pow(beta, 2)) + (std::pow(beta, 2) - E2 )) + std::fabs((lambda-1) * (std::pow(alpha, 3) - std::pow(beta, 3)) + (std::pow(beta, 3) - E3 )); PValues->put(i, P); AlphaValues->put(i, alpha); BetaValues->put(i, beta); LAValues->put(i,(lambda * (ER1 < ER2)) + ((1-lambda) * (ER1 >= ER2))); } Projection2(*PValues, *AlphaValues, *BetaValues); //Threshold(*AlphaValues); //Threshold(*BetaValues); DoubleLogarithm(*AlphaValues); DoubleLogarithm(*BetaValues); vnl_vector SignalVector(element_product((*LAValues) , (*AlphaValues)-(*BetaValues)) + (*BetaValues)); vnl_vector coeffs((*m_CoeffReconstructionMatrix) *SignalVector ); // the first coeff is a fix value coeffs[0] = 1.0/(2.0*sqrt(QBALL_ANAL_RECON_PI)); // Cast the Signal-Type from double to float for the ODF-Image odf = mitk::mitk_vnl_function::element_cast( (*m_ODFSphericalHarmonicBasisMatrix) * coeffs ).data_block(); odf *= (QBALL_ANAL_RECON_PI*4/NODF); //Normalize(odf); } // set ODF to ODF-Image odfOutputImageIterator.Set( odf ); ++odfOutputImageIterator; // iterate ++bzeroImageIterator; ++gradientInputImageIterator; } MITK_INFO << "THREAD FINISHED"; delete E; delete AlphaValues; delete BetaValues; delete PValues; delete LAValues; } 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) { itk::TimeProbe clock; GenerateAveragedBZeroImage(outputRegionForThread); clock.Start(); switch(m_ReconstructionType) { case Mode_Standard1Shell: StandardOneShellReconstruction(outputRegionForThread); break; case Mode_Analytical3Shells: AnalyticalThreeShellReconstruction(outputRegionForThread); break; case Mode_NumericalNShells: break; } clock.Stop(); MITK_INFO << "Reconstruction in : " << clock.GetTotal() << " TU"; } 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++) { 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::mitk_sh_functions::Yj(m,k,th,phi); + (*SHBasisOutput)(i,j) = mitk::sh::Yj(m,k,th,phi); // 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); // 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) { if(gdcit1.Value().one_norm() > 0.0) { centerMass += gdcit1.Value(); count ++; } } 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() { 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; std::map >::const_iterator it = (m_GradientIndexMap.begin()); it++; const std::vector gradientIndiciesVector= (*it).second; int numberOfGradientDirections = gradientIndiciesVector.size(); if( numberOfGradientDirections < (((L+1)*(L+2))/2) || numberOfGradientDirections < 6 ) { itkExceptionMacro( << "At least (L+1)(L+2)/2 gradient directions for each shell are required; current : " << numberOfGradientDirections ); } CheckDuplicateDiffusionGradients(); // 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; 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::mitk_sh_functions::Cart2Sph(x,y,z,cart); + mitk::sh::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::mitk_sh_functions::Cart2Sph(x,y,z,cart); + mitk::sh::Cart2Sph(x,y,z,cart); (*Q)(0,j) = cart[0]; (*Q)(1,j) = cart[1]; (*Q)(2,j++) = cart[2]; } } } const int LOrder = L; m_NumberCoefficients = (int)(LOrder*LOrder + LOrder + 2.0)/2.0 + LOrder; MITK_INFO << m_NumberCoefficients; 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()); // 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)))); 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)); m_SignalReonstructionMatrix = new vnl_matrix((*inverse) * (m_SHBasisMatrix->transpose())); m_CoeffReconstructionMatrix = new vnl_matrix(( factor * ((*FRTMatrix) * ((*SHEigenvalues) * (*m_SignalReonstructionMatrix))) )); // this code goes to the image adapter coeffs->odfs later vnl_matrix_fixed* U = itk::PointShell >::DistributePointShell(); m_ODFSphericalHarmonicBasisMatrix = new vnl_matrix(NOdfDirections,m_NumberCoefficients); m_ODFSphericalHarmonicBasisMatrix->fill(0.0); for(int i=0; i 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(); } template< class T, class TG, class TO, int L, int NODF> bool DiffusionMultiShellQballReconstructionImageFilter ::CheckDuplicateDiffusionGradients() { bool value = false; GradientIndexMapIteraotr mapIterator = m_GradientIndexMap.begin(); while(mapIterator != m_GradientIndexMap.end()) { std::vector::const_iterator it1 = mapIterator->second.begin(); std::vector::const_iterator it2 = mapIterator->second.begin(); for(; it1 != mapIterator->second.end(); ++it1) { for(; it2 != mapIterator->second.end(); ++it2) { if(m_GradientDirectionContainer->ElementAt(*it1) == m_GradientDirectionContainer->ElementAt(*it2) && it1 != it2) { itkWarningMacro( << "Some of the Diffusion Gradients equal each other. Corresponding image data should be averaged before calling this filter." ); value = true; } } } ++mapIterator; } return value; } 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 ); } } #endif // __itkDiffusionMultiShellQballReconstructionImageFilter_cpp diff --git a/Modules/DiffusionImaging/Reconstruction/mitkSphericalHarmonicsFunctions.cpp b/Modules/DiffusionImaging/Reconstruction/mitkSphericalHarmonicsFunctions.cpp index 07d4d24ca6..3f931d4374 100644 --- a/Modules/DiffusionImaging/Reconstruction/mitkSphericalHarmonicsFunctions.cpp +++ b/Modules/DiffusionImaging/Reconstruction/mitkSphericalHarmonicsFunctions.cpp @@ -1,91 +1,91 @@ /*========================================================================= 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. =========================================================================*/ #include #define _USE_MATH_DEFINES #include #if BOOST_VERSION / 100000 > 0 #if BOOST_VERSION / 100 % 1000 > 34 #include #include #endif #endif +#include "mitkVector.h" #include "mitkSphericalHarmonicsFunctions.h" namespace mitk{ +using namespace ::boost::math; -#define SPHERICAL_HARMONICS_PI M_PI - -double mitk::mitk_sh_functions::factorial(int number) { +double mitk::sh::factorial(int number) { if(number <= 1) return 1; double result = 1.0; for(int i=1; i<=number; i++) result *= i; return result; } -void mitk::mitk_sh_functions::Cart2Sph(double x, double y, double z, double *cart) +void mitk::sh::Cart2Sph(double x, double y, double z, double *cart) { double phi, th, rad; rad = sqrt(x*x+y*y+z*z); - th = atan2(z,sqrt(x*x+y*y)); - phi = atan2(y,x); - th = -th+SPHERICAL_HARMONICS_PI/2; - phi = -phi+SPHERICAL_HARMONICS_PI; + if( rad < mitk::eps ) + { + th = M_PI/2; + phi = M_PI/2; + } + else + { + th = acos(z/rad); + phi = atan2(y, x); + } cart[0] = phi; cart[1] = th; cart[2] = rad; } -double mitk::mitk_sh_functions::legendre0(int l) +double mitk::sh::legendre0(int l) { if( l%2 != 0 ) { return 0; } else { double prod1 = 1.0; for(int i=1;i(k,m,theta,phi);; - } - - if( m == 0 ) - return ::boost::math::spherical_harmonic_r(k,m,theta,phi); - - if( 0 < m && m <= k ) - { - return sqrt(2.0) * ::boost::math::spherical_harmonic_i(k,m,theta,phi); - } - - return 0; + if (m<0) + return sqrt(2.0)*spherical_harmonic_r(l, -m, theta, phi); + else if (m==0) + return spherical_harmonic_r(l, m, theta, phi); + else + return pow(-1.0,m)*sqrt(2.0)*spherical_harmonic_i(l, m, theta, phi); + + return 0; } } diff --git a/Modules/DiffusionImaging/Reconstruction/mitkSphericalHarmonicsFunctions.h b/Modules/DiffusionImaging/Reconstruction/mitkSphericalHarmonicsFunctions.h index f185b00ceb..d553212471 100644 --- a/Modules/DiffusionImaging/Reconstruction/mitkSphericalHarmonicsFunctions.h +++ b/Modules/DiffusionImaging/Reconstruction/mitkSphericalHarmonicsFunctions.h @@ -1,37 +1,37 @@ /*========================================================================= 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 __mitkSphericalHarmonicsHandler_h_ #define __mitkSphericalHarmonicsHandler_h_ #include "MitkDiffusionImagingExports.h" namespace mitk{ -namespace mitk_sh_functions{ +namespace sh{ double factorial(int number); void Cart2Sph(double x, double y, double z, double* cart); double legendre0(int l); double spherical_harmonic(int m,int l,double theta,double phi, bool complexPart); double Yj(int m, int k, double theta, double phi); } } #endif //__mitkSphericalHarmonicsHandler_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 94c480a324..c565f526d1 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkQBallReconstructionView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkQBallReconstructionView.cpp @@ -1,1005 +1,1005 @@ /*========================================================================= 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); + m_View->m_Controls->m_groupQShellSelection->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); + m_Controls->m_groupQShellSelection->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()); + QVBoxLayout * layout = dynamic_cast(m_Controls->m_groupQShellSelection->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(); typedef itk::DiffusionMultiShellQballReconstructionImageFilter::GradientIndexMap GradientIndexMap; GradientIndexMap::iterator it; if(QString("DiffusionImage").compare(node->GetData()->GetNameOfClass())==0) { mitk::DiffusionImage* vols = static_cast*>(node->GetData()); vols->GetBvalueList(&m_GradientIndexMap); QString str = QString((dynamic_cast(node->GetPropertyList()->GetProperty("name")))->GetValue()); QLabel *label = new QLabel("Image: " + str ); layout->addWidget(label); for(it = m_GradientIndexMap.begin(); it != m_GradientIndexMap.end(); it ++){ if(it->first == 0) continue; QString checkboxName ; 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); + m_Controls->m_groupQShellSelection->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() , vols->GetB_Value() ); - filter->SetBValue(vols->GetB_Value()); + //filter->SetBValue(vols->GetB_Value()); filter->SetThreshold( m_Controls->m_QBallReconstructionThreasholdEdit->text().toFloat() ); filter->SetLambda(lambda); //filter->SetNumberOfThreads(1); //filter->SetGradientIndexMap(m_GradientIndexMap); //for(std::map >::iterator it = m_GradientIndexMap->begin() ; it != m_GradientIndexMap->end(); it++) // MITK_INFO << (*it).first << " Count " << (*it).second.size(); 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); } }