diff --git a/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp b/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp index 470297ca81..bc8dd5e7ce 100644 --- a/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp +++ b/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp @@ -1,879 +1,879 @@ /*========================================================================= -Program: Insight Segmentation & Registration Toolkit -Module: $RCSfile: itkDiffusionTensor3DReconstructionImageFilter.txx,v $ +Program: Medical Imaging & Interaction Toolkit Language: C++ -Date: $Date: 2006-07-19 15:11:41 $ -Version: $Revision: 1.11 $ +Date: $Date: 2009-07-14 19:11:20 +0200 (Tue, 14 Jul 2009) $ +Version: $Revision: 18127 $ -Copyright (c) Insight Software Consortium. All rights reserved. -See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details. +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 #include #if BOOST_VERSION / 100000 > 0 #if BOOST_VERSION / 100 % 1000 > 34 #include #endif #endif #include "itkPointShell.h" namespace itk { #define QBALL_ANAL_RECON_PI 3.14159265358979323846 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) { // 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; i0) odf /= sum; break; } case QBAR_NONNEG_SOLID_ANGLE: { break; } } return odf; } template< class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NOrderL, int NrOdfDirections> vnl_vector itk::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= b0f) vec[i] = b0f - 0.001; vec[i] = log(-log(vec[i]/b0f)); } return vec; break; } } return vec; } template< class T, class TG, class TO, int L, int NODF> void 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(); m_BZeroImage = BZeroImageType::New(); typename GradientImagesType::Pointer img = static_cast< GradientImagesType * >( this->ProcessObject::GetInput(0) ); m_BZeroImage->SetSpacing( img->GetSpacing() ); // Set the image spacing m_BZeroImage->SetOrigin( img->GetOrigin() ); // Set the image origin m_BZeroImage->SetDirection( img->GetDirection() ); // Set the image direction m_BZeroImage->SetLargestPossibleRegion( img->GetLargestPossibleRegion()); m_BZeroImage->SetBufferedRegion( img->GetLargestPossibleRegion() ); m_BZeroImage->Allocate(); m_ODFSumImage = BZeroImageType::New(); m_ODFSumImage->SetSpacing( img->GetSpacing() ); // Set the image spacing m_ODFSumImage->SetOrigin( img->GetOrigin() ); // Set the image origin m_ODFSumImage->SetDirection( img->GetDirection() ); // Set the image direction m_ODFSumImage->SetLargestPossibleRegion( img->GetLargestPossibleRegion()); m_ODFSumImage->SetBufferedRegion( img->GetLargestPossibleRegion() ); m_ODFSumImage->Allocate(); if(m_NormalizationMethod == QBAR_SOLID_ANGLE || m_NormalizationMethod == QBAR_NONNEG_SOLID_ANGLE) { m_Lambda = 0.0; } } template< class T, class TG, class TO, int L, int NODF> void 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< BlaImage > oit3(m_ODFSumImage, outputRegionForThread); oit3.GoToBegin(); typedef ImageRegionConstIterator< GradientImagesType > GradientIteratorType; typedef typename GradientImagesType::PixelType GradientVectorType; typename GradientImagesType::Pointer gradientImagePointer = NULL; // Would have liked a dynamic_cast here, but seems SGI doesn't like it // The enum will ensure that an inappropriate cast is not done gradientImagePointer = static_cast< GradientImagesType * >( this->ProcessObject::GetInput(0) ); GradientIteratorType git(gradientImagePointer, outputRegionForThread ); git.GoToBegin(); // Compute the indicies of the baseline images and gradient images std::vector baselineind; // contains the indicies of // the baseline images std::vector gradientind; // contains the indicies of // the gradient images for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin(); gdcit != this->m_GradientDirectionContainer->End(); ++gdcit) { if(gdcit.Value().one_norm() <= 0.0) { baselineind.push_back(gdcit.Index()); } else { gradientind.push_back(gdcit.Index()); } } if( m_DirectionsDuplicated ) { int gradIndSize = gradientind.size(); for(int i=0; i::AccumulateType b0 = NumericTraits::Zero; // Average the baseline image pixels for(unsigned int i = 0; i < baselineind.size(); ++i) { b0 += b[baselineind[i]]; } b0 /= this->m_NumberOfBaselineImages; OdfPixelType odf(0.0); vnl_vector B(m_NumberOfGradientDirections); if( (b0 != 0) && (b0 >= m_Threshold) ) { for( unsigned int i = 0; i< m_NumberOfGradientDirections; i++ ) { B[i] = static_cast(b[gradientind[i]]); } B = PreNormalize(B, b0); if(m_NormalizationMethod == QBAR_SOLID_ANGLE) { vnl_vector coeffs(m_NumberCoefficients); coeffs = ( (*m_CoeffReconstructionMatrix) * B ); coeffs[0] += 1.0/(2.0*sqrt(QBALL_ANAL_RECON_PI)); odf = ( (*m_SphericalHarmonicBasisMatrix) * coeffs ).data_block(); } else if(m_NormalizationMethod == QBAR_NONNEG_SOLID_ANGLE) { /** this would be the place to implement a non-negative * solver for quadratic programming problem: * min .5*|| Bc-s ||^2 subject to -CLPc <= 4*pi*ones * (refer to MICCAI 2009 Goh et al. "Estimating ODFs with PDF constraints") * .5*|| Bc-s ||^2 == .5*c'B'Bc - x'B's + .5*s's */ itkExceptionMacro( << "Nonnegative Solid Angle not yet implemented"); } else { odf = ( (*m_ReconstructionMatrix) * B ).data_block(); } odf = Normalize(odf, b0); } oit.Set( odf ); ++oit; oit2.Set( b0 ); float sum = 0; for (int k=0; k void 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 double AnalyticalDiffusionQballReconstructionImageFilter ::factorial(int number) { if(number <= 1) return 1; double result = 1.0; for(int i=1; i<=number; i++) result *= i; return result; } template< class T, class TG, class TO, int L, int NODF> void AnalyticalDiffusionQballReconstructionImageFilter ::Cart2Sph(double x, double y, double z, double *cart) { double phi, th, rad; rad = sqrt(x*x+y*y+z*z); th = atan2(z,sqrt(x*x+y*y)); phi = atan2(y,x); th = -th+QBALL_ANAL_RECON_PI/2; phi = -phi+QBALL_ANAL_RECON_PI; cart[0] = phi; cart[1] = th; cart[2] = rad; } template< class T, class TG, class TO, int L, int NODF> double AnalyticalDiffusionQballReconstructionImageFilter ::legendre0(int l) { if( l%2 != 0 ) { return 0; } else { double prod1 = 1.0; for(int i=1;i double AnalyticalDiffusionQballReconstructionImageFilter ::spherical_harmonic(int m,int l,double theta,double phi, bool complexPart) { if( theta < 0 || theta > QBALL_ANAL_RECON_PI ) { std::cout << "bad range" << std::endl; return 0; } if( phi < 0 || phi > 2*QBALL_ANAL_RECON_PI ) { std::cout << "bad range" << std::endl; return 0; } double pml = 0; double fac1 = factorial(l+m); double fac2 = factorial(l-m); if( m<0 ) { #if BOOST_VERSION / 100000 > 0 #if BOOST_VERSION / 100 % 1000 > 34 pml = ::boost::math::legendre_p(l, -m, cos(theta)); #else std::cout << "ERROR: Boost 1.35 minimum required" << std::endl; #endif #else std::cout << "ERROR: Boost 1.35 minimum required" << std::endl; #endif double mypow = pow(-1.0,-m); double myfac = (fac1/fac2); pml *= mypow*myfac; } else { #if BOOST_VERSION / 100000 > 0 #if BOOST_VERSION / 100 % 1000 > 34 pml = ::boost::math::legendre_p(l, m, cos(theta)); #endif #endif } //std::cout << "legendre(" << l << "," << m << "," << cos(theta) << ") = " << pml << std::endl; double retval = sqrt(((2.0*(double)l+1.0)/(4.0*QBALL_ANAL_RECON_PI))*(fac2/fac1)) * pml; if( !complexPart ) { retval *= cos(m*phi); } else { retval *= sin(m*phi); } //std::cout << retval << std::endl; return retval; } template< class T, class TG, class TO, int L, int NODF> double AnalyticalDiffusionQballReconstructionImageFilter ::Yj(int m, int k, double theta, double phi) { if( -k <= m && m < 0 ) { return sqrt(2.0) * spherical_harmonic(m,k,theta,phi,false); } if( m == 0 ) return spherical_harmonic(0,k,theta,phi,false); if( 0 < m && m <= k ) { return sqrt(2.0) * spherical_harmonic(m,k,theta,phi,true); } return 0; } template< class T, class TG, class TO, int L, int NODF> void 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]; Cart2Sph(x,y,z,cart); (*Q)(0,i) = cart[0]; (*Q)(1,i) = cart[1]; (*Q)(2,i++) = cart[2]; } } if(m_DirectionsDuplicated) { for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin(); gdcit != this->m_GradientDirectionContainer->End(); ++gdcit) { if(gdcit.Value().one_norm() > 0.0) { double x = gdcit.Value().get(0); double y = gdcit.Value().get(1); double z = gdcit.Value().get(2); double cart[3]; Cart2Sph(x,y,z,cart); (*Q)(0,i) = cart[0]; (*Q)(1,i) = cart[1]; (*Q)(2,i++) = cart[2]; } } } } int l = L; m_NumberCoefficients = (int)(l*l + l + 2.0)/2.0 + l; vnl_matrix* B = new vnl_matrix(m_NumberOfGradientDirections,m_NumberCoefficients); vnl_matrix* _L = new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients); _L->fill(0.0); vnl_matrix* LL = new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients); LL->fill(0.0); vnl_matrix* P = new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients); P->fill(0.0); vnl_matrix* Inv = new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients); P->fill(0.0); vnl_vector* lj = new vnl_vector(m_NumberCoefficients); m_LP = new vnl_vector(m_NumberCoefficients); for(unsigned int i=0; i(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/itkAnalyticalDiffusionQballReconstructionImageFilter.h b/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.h index b183891ff5..1c89d2c9a2 100644 --- a/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.h +++ b/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.h @@ -1,300 +1,300 @@ /*========================================================================= - Program: Insight Segmentation & Registration Toolkit - Module: $RCSfile: itkDiffusionTensor3DReconstructionImageFilter.h,v $ - Language: C++ - Date: $Date: 2006-03-27 17:01:06 $ - Version: $Revision: 1.12 $ +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) Insight Software Consortium. All rights reserved. - See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details. +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. +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_h_ #define __itkAnalyticalDiffusionQballReconstructionImageFilter_h_ #include "itkImageToImageFilter.h" -//#include "vnl/vnl_matrix.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 "QuadProg.h" + namespace itk{ /** \class AnalyticalDiffusionQballReconstructionImageFilter * \brief This class takes as input one or more reference image (acquired in the * absence of diffusion sensitizing gradients) and 'n' diffusion * weighted images and their gradient directions and computes an image of * orientation distribution function coefficients in a spherical harmonic basis. * * \par Inputs and Usage * \par * When you have the 'n' gradient and one or more reference images in a single * multi-component image (VectorImage), you can specify the images as * \code * filter->SetGradientImage( directionsContainer, vectorImage ); * \endcode * Note that this method is used to specify both the reference and gradient images. * This is convenient when the DWI images are read in using the * NRRD * format. Like the Nrrd format, the reference images are those components of the * vectorImage whose gradient direction is (0,0,0). If more than one reference image * is present, they are averaged prior to the reconstruction. * * \par Outputs * The output image is an image of vectors that must be understood as ODFs: * \code * Image< Vector< TPixelType, OdfNrDirections >, 3 > * \endcode * * \par Parameters * \li Threshold - Threshold on the reference image data. The output ODF will * be a null pdf for pixels in the reference image that have a value less * than this. * \li BValue - See the documentation of SetBValue(). * \li At least 6 gradient images must be specified for the filter to be able * to run. If the input gradient directions g_i are majorly sampled on one half * of the sqhere, then each input image I_i will be duplicated and assign -g_i * in order to guarantee stability of the algorithm. * \li OdfDirections - directions of resulting orientation distribution function * \li EquatorNrSamplingPoints - number of sampling points on equator when * performing Funk Radeon Transform (FRT) * \li BasisFunctionCenters - the centers of the basis functions are used for * the sRBF (spherical radial basis functions interpolation). If not set, they * will be defaulted to equal m_EquatorNrSamplingPoints * * \par Template parameters * The class is templated over * \li the pixel type of the reference and gradient images * (expected to be scalar data types) * \li the internal representation of the ODF pixels (double, float etc). * \li the number of OdfDirections * \li the number of basis function centers for the sRBF * * \par References: * \li[1] * Tuch DS, * "Q-ball imaging", Magn Reson Med. 2004 Dec;52(6):1358-72. * */ template< class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NOrderL, int NrOdfDirections> class AnalyticalDiffusionQballReconstructionImageFilter : public ImageToImageFilter< Image< TReferenceImagePixelType, 3 >, Image< Vector< TOdfPixelType, NrOdfDirections >, 3 > > { public: enum Normalization { QBAR_STANDARD, QBAR_B_ZERO_B_VALUE, QBAR_B_ZERO, QBAR_NONE, QBAR_ADC_ONLY, QBAR_RAW_SIGNAL, QBAR_SOLID_ANGLE, QBAR_NONNEG_SOLID_ANGLE }; typedef AnalyticalDiffusionQballReconstructionImageFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ImageToImageFilter< Image< TReferenceImagePixelType, 3>, Image< Vector< TOdfPixelType, NrOdfDirections >, 3 > > Superclass; /** Method for creation through the object factory. */ itkNewMacro(Self); /** Runtime information support. */ itkTypeMacro(AnalyticalDiffusionQballReconstructionImageFilter, ImageToImageFilter); typedef TReferenceImagePixelType ReferencePixelType; typedef TGradientImagePixelType GradientPixelType; typedef Vector< TOdfPixelType, NrOdfDirections > OdfPixelType; /** Reference image data, This image is aquired in the absence * of a diffusion sensitizing field gradient */ typedef typename Superclass::InputImageType ReferenceImageType; typedef Image< OdfPixelType, 3 > OdfImageType; typedef OdfImageType OutputImageType; typedef TOdfPixelType BZeroPixelType; typedef Image< BZeroPixelType, 3 > BZeroImageType; typedef typename Superclass::OutputImageRegionType OutputImageRegionType; /** Typedef defining one (of the many) gradient images. */ typedef Image< GradientPixelType, 3 > GradientImageType; /** An alternative typedef defining one (of the many) gradient images. * It will be assumed that the vectorImage has the same dimension as the * Reference image and a vector length parameter of \c n (number of * gradient directions)*/ typedef VectorImage< GradientPixelType, 3 > GradientImagesType; /** Holds the ODF reconstruction matrix */ typedef vnl_matrix< TOdfPixelType >* OdfReconstructionMatrixType; typedef vnl_matrix< double > CoefficientMatrixType; /** Holds each magnetic field gradient used to acquire one DWImage */ typedef vnl_vector_fixed< double, 3 > GradientDirectionType; /** Container to hold gradient directions of the 'n' DW measurements */ typedef VectorContainer< unsigned int, GradientDirectionType > GradientDirectionContainerType; /** set method to add gradient directions and its corresponding * image. The image here is a VectorImage. The user is expected to pass the * gradient directions in a container. The ith element of the container * corresponds to the gradient direction of the ith component image the * VectorImage. For the baseline image, a vector of all zeros * should be set.*/ void SetGradientImage( GradientDirectionContainerType *, const GradientImagesType *image); /** Get reference image */ virtual ReferenceImageType * GetReferenceImage() { return ( static_cast< ReferenceImageType *>(this->ProcessObject::GetInput(0)) ); } /** Return the gradient direction. idx is 0 based */ virtual GradientDirectionType GetGradientDirection( unsigned int idx) const { if( idx >= m_NumberOfGradientDirections ) { itkExceptionMacro( << "Gradient direction " << idx << "does not exist" ); } return m_GradientDirectionContainer->ElementAt( idx+1 ); } static void tofile2(vnl_matrix *A, std::string fname); static double factorial(int number); static void Cart2Sph(double x, double y, double z, double* cart); static double legendre0(int l); static double spherical_harmonic(int m,int l,double theta,double phi, bool complexPart); static double Yj(int m, int k, double theta, double phi); OdfPixelType Normalize(OdfPixelType odf, typename NumericTraits::AccumulateType b0 ); 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 ); itkSetMacro( NormalizationMethod, Normalization); itkGetMacro( NormalizationMethod, Normalization ); typedef Image BlaImage; itkGetMacro( BZeroImage, typename BZeroImageType::Pointer); itkGetMacro( ODFSumImage, typename BlaImage::Pointer); itkSetMacro( BValue, TOdfPixelType); #ifdef GetBValue #undef GetBValue #endif itkGetConstReferenceMacro( BValue, TOdfPixelType); itkSetMacro( Lambda, double ); itkGetMacro( Lambda, double ); #ifdef ITK_USE_CONCEPT_CHECKING /** Begin concept checking */ itkConceptMacro(ReferenceEqualityComparableCheck, (Concept::EqualityComparable)); itkConceptMacro(TensorEqualityComparableCheck, (Concept::EqualityComparable)); itkConceptMacro(GradientConvertibleToDoubleCheck, (Concept::Convertible)); itkConceptMacro(DoubleConvertibleToTensorCheck, (Concept::Convertible)); itkConceptMacro(GradientReferenceAdditiveOperatorsCheck, (Concept::AdditiveOperators)); itkConceptMacro(ReferenceOStreamWritableCheck, (Concept::OStreamWritable)); itkConceptMacro(TensorOStreamWritableCheck, (Concept::OStreamWritable)); /** End concept checking */ #endif protected: AnalyticalDiffusionQballReconstructionImageFilter(); ~AnalyticalDiffusionQballReconstructionImageFilter() {}; void PrintSelf(std::ostream& os, Indent indent) const; void ComputeReconstructionMatrix(); void BeforeThreadedGenerateData(); void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, int); private: OdfReconstructionMatrixType m_ReconstructionMatrix; OdfReconstructionMatrixType m_CoeffReconstructionMatrix; 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; double m_Lambda; bool m_DirectionsDuplicated; Normalization m_NormalizationMethod; int m_NumberCoefficients; vnl_matrix* m_B_t; vnl_vector* m_LP; BlaImage::Pointer m_ODFSumImage; }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkAnalyticalDiffusionQballReconstructionImageFilter.cpp" #endif #endif //__itkAnalyticalDiffusionQballReconstructionImageFilter_h_ diff --git a/Modules/DiffusionImaging/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.cpp b/Modules/DiffusionImaging/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.cpp index 9dc4809b85..328c39701a 100644 --- a/Modules/DiffusionImaging/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.cpp +++ b/Modules/DiffusionImaging/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.cpp @@ -1,769 +1,776 @@ -/*======================================================================= -Copyright (c) Insight Software Consortium. All rights reserved. -See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details. +/*========================================================================= -This software is distributed WITHOUT ANY WARRANTY; without even -the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +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 __itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter_cpp #define __itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter_cpp #include "itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h" #include "itkImageRegionConstIterator.h" #include "itkImageRegionConstIteratorWithIndex.h" #include "itkImageRegionIterator.h" #include "vnl/vnl_matrix.h" #include "vnl/algo/vnl_symmetric_eigensystem.h" #include "itkRegularizedIVIMReconstructionFilter.h" #include #define IVIM_FOO -100000 namespace itk { template< class TIn, class TOut> DiffusionIntravoxelIncoherentMotionReconstructionImageFilter ::DiffusionIntravoxelIncoherentMotionReconstructionImageFilter() : m_GradientDirectionContainer(NULL), m_Method(IVIM_DSTAR_FIX), m_FitDStar(true), m_Verbose(true) { this->SetNumberOfRequiredInputs( 1 ); this->SetNumberOfRequiredOutputs( 3 ); typename OutputImageType::Pointer outputPtr1 = OutputImageType::New(); this->SetNthOutput(0, outputPtr1.GetPointer()); typename OutputImageType::Pointer outputPtr2 = OutputImageType::New(); this->SetNthOutput(1, outputPtr2.GetPointer()); typename OutputImageType::Pointer outputPtr3 = OutputImageType::New(); this->SetNthOutput(2, outputPtr3.GetPointer()); } template< class TIn, class TOut> void DiffusionIntravoxelIncoherentMotionReconstructionImageFilter ::BeforeThreadedGenerateData() { // 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 ); } typename InputImageType::Pointer img = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) ); // Compute the indicies of the baseline images and gradient images GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin(); while( gdcit != this->m_GradientDirectionContainer->End() ) { if(gdcit.Value().one_norm() <= 0.0) { m_Snap.baselineind.push_back(gdcit.Index()); } else { m_Snap.gradientind.push_back(gdcit.Index()); double twonorm = gdcit.Value().two_norm(); m_Snap.bvals.push_back( m_BValue*twonorm*twonorm ); } ++gdcit; } // check sind die grad und base gleichlang? alle grad gerade und base ungerade? dann iterierende aufnahme!! m_Snap.iterated_sequence = false; if(m_Snap.baselineind.size() == m_Snap.gradientind.size()) { int size = m_Snap.baselineind.size(); int sum_b = 0, sum_g = 0; for(int i=0; im_BThres) { m_Snap.high_indices.push_back(i); } } } m_Snap.Nhigh = m_Snap.high_indices.size(); m_Snap.high_bvalues.set_size(m_Snap.Nhigh); m_Snap.high_meas.set_size(m_Snap.Nhigh); for(int i=0; i MeasAndBvals DiffusionIntravoxelIncoherentMotionReconstructionImageFilter ::ApplyS0Threshold(vnl_vector &meas, vnl_vector &bvals) { std::vector newmeas; std::vector newbvals; int N = meas.size(); for(int i=0; i void DiffusionIntravoxelIncoherentMotionReconstructionImageFilter ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int ) { typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread); oit.GoToBegin(); typename OutputImageType::Pointer dImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(1)); ImageRegionIterator< OutputImageType > oit1(dImage, outputRegionForThread); oit1.GoToBegin(); typename OutputImageType::Pointer dstarImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(2)); ImageRegionIterator< OutputImageType > oit2(dstarImage, outputRegionForThread); oit2.GoToBegin(); typedef ImageRegionConstIterator< InputImageType > InputIteratorType; typedef typename InputImageType::PixelType InputVectorType; typename InputImageType::Pointer inputImagePointer = NULL; // Would have liked a dynamic_cast here, but seems SGI doesn't like it // The enum will DiffusionIntravoxelIncoherentMotionReconstructionImageFilterensure that an inappropriate cast is not done inputImagePointer = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) ); InputIteratorType iit(inputImagePointer, outputRegionForThread ); iit.GoToBegin(); // init internal vector image for regularized fit m_InternalVectorImage = VectorImageType::New(); m_InternalVectorImage->SetSpacing( inputImagePointer->GetSpacing() ); // Set the image spacing m_InternalVectorImage->SetOrigin( inputImagePointer->GetOrigin() ); // Set the image origin m_InternalVectorImage->SetDirection( inputImagePointer->GetDirection() ); // Set the image direction m_InternalVectorImage->SetRegions( inputImagePointer->GetLargestPossibleRegion() ); m_InitialFitImage = InitialFitImageType::New(); m_InitialFitImage->SetSpacing( inputImagePointer->GetSpacing() ); // Set the image spacing m_InitialFitImage->SetOrigin( inputImagePointer->GetOrigin() ); // Set the image origin m_InitialFitImage->SetDirection( inputImagePointer->GetDirection() ); // Set the image direction m_InitialFitImage->SetRegions( inputImagePointer->GetLargestPossibleRegion() ); if(m_Method == IVIM_REGULARIZED) { m_InternalVectorImage->SetVectorLength(m_Snap.Nhigh); m_InternalVectorImage->Allocate(); VectorImageType::PixelType varvec(m_Snap.Nhigh); for(int i=0; iFillBuffer(varvec); m_InitialFitImage->Allocate(); InitialFitImageType::PixelType vec; vec[0] = 0.5; vec[1] = 0.01; vec[2]=0.001; m_InitialFitImage->FillBuffer(vec); } typedef itk::ImageRegionIterator VectorIteratorType; VectorIteratorType vecit(m_InternalVectorImage, outputRegionForThread ); vecit.GoToBegin(); typedef itk::ImageRegionIterator InitIteratorType; InitIteratorType initit(m_InitialFitImage, outputRegionForThread ); initit.GoToBegin(); while( !iit.IsAtEnd() ) { InputVectorType measvec = iit.Get(); typename NumericTraits::AccumulateType b0 = NumericTraits::Zero; m_Snap.meas.set_size(m_Snap.N); m_Snap.allmeas.set_size(m_Snap.N); if(!m_Snap.iterated_sequence) { // Average the baseline image pixels for(unsigned int i = 0; i < m_Snap.baselineind.size(); ++i) { b0 += measvec[m_Snap.baselineind[i]]; } if(m_Snap.baselineind.size()) b0 /= m_Snap.baselineind.size(); // measurement vector for(int i = 0; i < m_Snap.N; ++i) { m_Snap.allmeas[i] = measvec[m_Snap.gradientind[i]] / (b0+.0001); if(measvec[m_Snap.gradientind[i]] > m_S0Thres) { m_Snap.meas[i] = measvec[m_Snap.gradientind[i]] / (b0+.0001); } else { m_Snap.meas[i] = IVIM_FOO; } } } else { // measurement vector for(int i = 0; i < m_Snap.N; ++i) { b0 = measvec[m_Snap.baselineind[i]]; m_Snap.allmeas[i] = measvec[m_Snap.gradientind[i]] / (b0+.0001); if(measvec[m_Snap.gradientind[i]] > m_S0Thres) { m_Snap.meas[i] = measvec[m_Snap.gradientind[i]] / (b0+.0001); } else { m_Snap.meas[i] = IVIM_FOO; } } } m_Snap.currentF = 0; m_Snap.currentD = 0; m_Snap.currentDStar = 0; switch(m_Method) { case IVIM_D_THEN_DSTAR: { for(int i=0; i x_donly(2); x_donly[0] = 0.001; x_donly[1] = 0.1; // f 0.1 Dstar 0.01 D 0.001 vnl_levenberg_marquardt lm_donly(f_donly); lm_donly.set_f_tolerance(0.0001); lm_donly.minimize(x_donly); m_Snap.currentD = x_donly[0]; m_Snap.currentF = x_donly[1]; if(m_FitDStar) { MeasAndBvals input2 = ApplyS0Threshold(m_Snap.meas, m_Snap.bvalues); m_Snap.bvals2 = input2.bvals; m_Snap.meas2 = input2.meas; if (input2.N < 2) break; IVIM_dstar_only f_dstar_only(input2.N,m_Snap.currentD,m_Snap.currentF); f_dstar_only.set_bvalues(input2.bvals); f_dstar_only.set_measurements(input2.meas); vnl_vector< double > x_dstar_only(1); vnl_vector< double > fx_dstar_only(input2.N); double opt = 1111111111111111.0; int opt_idx = -1; int num_its = 100; double min_val = .001; double max_val = .15; for(int i=0; i x_fixd(2); // x_fixd[0] = 0.1; // x_fixd[1] = 0.01; // // f 0.1 Dstar 0.01 D 0.001 // vnl_levenberg_marquardt lm_fixd(f_fixd); // lm_fixd.set_f_tolerance(0.0001); // lm_fixd.minimize(x_fixd); // m_Snap.currentF = x_fixd[0]; // m_Snap.currentDStar = x_fixd[1]; } break; } case IVIM_DSTAR_FIX: { MeasAndBvals input = ApplyS0Threshold(m_Snap.meas, m_Snap.bvalues); m_Snap.bvals1 = input.bvals; m_Snap.meas1 = input.meas; if (input.N < 2) break; IVIM_fixdstar f_fixdstar(input.N,m_DStar); f_fixdstar.set_bvalues(input.bvals); f_fixdstar.set_measurements(input.meas); vnl_vector< double > x(2); x[0] = 0.1; x[1] = 0.001; // f 0.1 Dstar 0.01 D 0.001 vnl_levenberg_marquardt lm(f_fixdstar); lm.set_f_tolerance(0.0001); lm.minimize(x); m_Snap.currentF = x[0]; m_Snap.currentD = x[1]; m_Snap.currentDStar = m_DStar; break; } case IVIM_FIT_ALL: { MeasAndBvals input = ApplyS0Threshold(m_Snap.meas, m_Snap.bvalues); m_Snap.bvals1 = input.bvals; m_Snap.meas1 = input.meas; if (input.N < 3) break; IVIM_3param f_3param(input.N); f_3param.set_bvalues(input.bvals); f_3param.set_measurements(input.meas); vnl_vector< double > x(3); x[0] = 0.1; x[1] = 0.001; x[2] = 0.01; // f 0.1 Dstar 0.01 D 0.001 vnl_levenberg_marquardt lm(f_3param); lm.set_f_tolerance(0.0001); lm.minimize(x); m_Snap.currentF = x[0]; m_Snap.currentD = x[1]; m_Snap.currentDStar = x[2]; break; } case IVIM_LINEAR_D_THEN_F: { // // neglect zero-measurements // bool zero = false; // for(int i=0; i X(input.N,2); for(int i=0; i XX = X.transpose() * X; vnl_symmetric_eigensystem eigs(XX); vnl_vector eig; if(eigs.get_eigenvalue(0) > eigs.get_eigenvalue(1)) eig = eigs.get_eigenvector(0); else eig = eigs.get_eigenvector(1); m_Snap.currentF = 1 - exp( meas_m - bval_m*(eig(1)/eig(0)) ); m_Snap.currentD = -eig(1)/eig(0); if(m_FitDStar) { MeasAndBvals input2 = ApplyS0Threshold(m_Snap.meas, m_Snap.bvalues); m_Snap.bvals2 = input2.bvals; m_Snap.meas2 = input2.meas; if (input2.N < 2) break; IVIM_dstar_only f_dstar_only(input2.N,m_Snap.currentD,m_Snap.currentF); f_dstar_only.set_bvalues(input2.bvals); f_dstar_only.set_measurements(input2.meas); vnl_vector< double > x_dstar_only(1); vnl_vector< double > fx_dstar_only(input2.N); double opt = 1111111111111111.0; int opt_idx = -1; int num_its = 100; double min_val = .001; double max_val = .15; for(int i=0; i " << DStar; // x_dstar_only[0] = 0.01; // // f 0.1 Dstar 0.01 D 0.001 // vnl_levenberg_marquardt lm_dstar_only(f_dstar_only); // lm_dstar_only.set_f_tolerance(0.0001); // lm_dstar_only.minimize(x_dstar_only); // DStar = x_dstar_only[0]; break; } case IVIM_REGULARIZED: { //m_Snap.high_meas, m_Snap.high_bvalues; for(int i=0; i x_donly(2); x_donly[0] = 0.001; x_donly[1] = 0.1; if(input.N >= 2) { IVIM_d_and_f f_donly(input.N); f_donly.set_bvalues(input.bvals); f_donly.set_measurements(input.meas); //MITK_INFO << "initial fit N=" << input.N << ", min-b = " << input.bvals[0] << ", max-b = " << input.bvals[input.N-1]; vnl_levenberg_marquardt lm_donly(f_donly); lm_donly.set_f_tolerance(0.0001); lm_donly.minimize(x_donly); } typename InitialFitImageType::PixelType initvec; initvec[0] = x_donly[1]; initvec[1] = x_donly[0]; initit.Set(initvec); //MITK_INFO << "Init vox " << initit.GetIndex() << " with " << initvec[0] << "; " << initvec[1]; ++initit; int N = m_Snap.high_meas.size(); typename VectorImageType::PixelType vec(N); for(int i=0; i void DiffusionIntravoxelIncoherentMotionReconstructionImageFilter ::AfterThreadedGenerateData() { if(m_Method == IVIM_REGULARIZED) { typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); ImageRegionIterator< OutputImageType > oit0(outputImage, outputImage->GetLargestPossibleRegion()); oit0.GoToBegin(); typename OutputImageType::Pointer dImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(1)); ImageRegionIterator< OutputImageType > oit1(dImage, dImage->GetLargestPossibleRegion()); oit1.GoToBegin(); typename OutputImageType::Pointer dstarImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(2)); ImageRegionIterator< OutputImageType > oit2(dstarImage, dstarImage->GetLargestPossibleRegion()); oit2.GoToBegin(); typedef itk::RegularizedIVIMReconstructionFilter RegFitType; RegFitType::Pointer filter = RegFitType::New(); filter->SetInput(m_InitialFitImage); filter->SetReferenceImage(m_InternalVectorImage); filter->SetBValues(m_Snap.high_bvalues); filter->SetNumberIterations(m_NumberIterations); filter->SetNumberOfThreads(1); filter->SetLambda(m_Lambda); filter->Update(); typename RegFitType::OutputImageType::Pointer outimg = filter->GetOutput(); ImageRegionConstIterator< RegFitType::OutputImageType > iit(outimg, outimg->GetLargestPossibleRegion()); iit.GoToBegin(); while( !iit.IsAtEnd() ) { double f = iit.Get()[0]; IVIM_CEIL( f, 0.0, 1.0 ); oit0.Set( myround(f * 100.0) ); oit1.Set( myround(iit.Get()[1] * 10000.0) ); oit2.Set( myround(iit.Get()[2] * 1000.0) ); if(!m_Verbose) { // report the middle voxel if( iit.GetIndex()[0] == m_CrossPosition[0] && iit.GetIndex()[1] == m_CrossPosition[1] && iit.GetIndex()[2] == m_CrossPosition[2] ) { m_Snap.currentF = f; m_Snap.currentD = iit.Get()[1]; m_Snap.currentDStar = iit.Get()[2]; m_Snap.allmeas = m_tmp_allmeas; MITK_INFO << "setting " << f << ";" << iit.Get()[1] << ";" << iit.Get()[2]; } } ++oit0; ++oit1; ++oit2; ++iit; } } } template< class TIn, class TOut> double DiffusionIntravoxelIncoherentMotionReconstructionImageFilter ::myround(double number) { return number < 0.0 ? ceil(number - 0.5) : floor(number + 0.5); } template< class TIn, class TOut> void DiffusionIntravoxelIncoherentMotionReconstructionImageFilter ::SetGradientDirections( GradientDirectionContainerType *gradientDirection ) { this->m_GradientDirectionContainer = gradientDirection; this->m_NumberOfGradientDirections = gradientDirection->Size(); } template< class TIn, class TOut> void DiffusionIntravoxelIncoherentMotionReconstructionImageFilter ::PrintSelf(std::ostream& os, Indent indent) const { Superclass::PrintSelf(os,indent); if ( m_GradientDirectionContainer ) { os << indent << "GradientDirectionContainer: " << m_GradientDirectionContainer << std::endl; } else { os << indent << "GradientDirectionContainer: (Gradient directions not set)" << std::endl; } } } #endif // __itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter_cpp diff --git a/Modules/DiffusionImaging/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h b/Modules/DiffusionImaging/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h index bf87f6d2f4..8fb0ca16f7 100644 --- a/Modules/DiffusionImaging/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h +++ b/Modules/DiffusionImaging/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h @@ -1,381 +1,381 @@ /*========================================================================= - Program: Insight Segmentation & Registration Toolkit - Module: $RCSfile: itkDiffusionTensor3DReconstructionImageFilter.h,v $ - Language: C++ - Date: $Date: 2006-03-27 17:01:06 $ - Version: $Revision: 1.12 $ - - Copyright (c) Insight Software Consortium. All rights reserved. - See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm 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. +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 __itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter_h #define __itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter_h #include "itkImageToImageFilter.h" //#include "vnl/vnl_matrix.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 "QuadProg.h" #include "itkVectorImage.h" #include "vnl/vnl_least_squares_function.h" #include "vnl/algo/vnl_levenberg_marquardt.h" #include "vnl/vnl_math.h" #define IVIM_CEIL(val,u,o) (val) = \ ( (val) < (u) ) ? ( (u) ) : ( ( (val)>(o) ) ? ( (o) ) : ( (val) ) ); namespace itk{ /** baseclass for IVIM fitting algorithms */ struct IVIM_base { void set_measurements(const vnl_vector& x) { measurements.set_size(x.size()); measurements.copy_in(x.data_block()); } void set_bvalues(const vnl_vector& x) { bvalues.set_size(x.size()); bvalues.copy_in(x.data_block()); } vnl_vector measurements; vnl_vector bvalues; int N; }; /** Fitt all three parameters */ struct IVIM_3param : public IVIM_base, vnl_least_squares_function { IVIM_3param(unsigned int number_of_measurements) : vnl_least_squares_function(3, number_of_measurements, no_gradient) { N = get_number_of_residuals(); } void f(const vnl_vector& x, vnl_vector& fx) { double ef = x[0]; double D = x[1]; double Dstar = x[2]; for(int s=0; s& x, vnl_vector& fx) { double ef = x[0]; double D = x[1]; for(int s=0; s& x, vnl_vector& fx) { double D = x[0]; double f = x[1]; for(int s=0; s& x, vnl_vector& fx) { double ef = x[0]; double Dstar = x[1]; for(int s=0; s& x, vnl_vector& fx) { double Dstar = x[0]; for(int s=0; s meas; vnl_vector bvals; int N; }; /** \class DiffusionIntravoxelIncoherentMotionReconstructionImageFilter */ template< class TInputPixelType, class TOutputPixelType> class DiffusionIntravoxelIncoherentMotionReconstructionImageFilter : public ImageToImageFilter< VectorImage< TInputPixelType, 3 >, Image< TOutputPixelType, 3 > > { public: struct IVIMSnapshot { double currentF; double currentFunceiled; double currentD; double currentDStar; bool iterated_sequence; // wether each measurement has its own b0-acqu. std::vector baselineind; // baseline image indicies int N; // total number of measurements std::vector gradientind; // gradient image indicies std::vector bvals; // b-values != 0 vnl_vector bvalues; // copy of bvalues != 0 vnl_vector meas; // all measurements, thresholded blanked out vnl_vector allmeas; // all measurements int Nhigh; // number of used measurements std::vector high_indices; // indices of used measurements vnl_vector high_bvalues; // bvals of used measurements vnl_vector high_meas; // used measurements vnl_vector meas1; vnl_vector bvals1; vnl_vector meas2; vnl_vector bvals2; }; enum IVIM_Method { IVIM_FIT_ALL, IVIM_DSTAR_FIX, IVIM_D_THEN_DSTAR, IVIM_LINEAR_D_THEN_F, IVIM_REGULARIZED }; typedef DiffusionIntravoxelIncoherentMotionReconstructionImageFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ImageToImageFilter< VectorImage< TInputPixelType, 3>, Image< TOutputPixelType,3 > > Superclass; /** Method for creation through the object factory. */ itkNewMacro(Self); /** Runtime information support. */ itkTypeMacro(DiffusionIntravoxelIncoherentMotionReconstructionImageFilter, ImageToImageFilter); typedef TOutputPixelType OutputPixelType; typedef TInputPixelType InputPixelType; /** Reference image data, This image is aquired in the absence * of a diffusion sensitizing field gradient */ typedef typename Superclass::InputImageType InputImageType; typedef Image< OutputPixelType, 3 > OutputImageType; typedef typename Superclass::OutputImageRegionType OutputImageRegionType; /** 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; // vector image typedefs for regularized fit typedef itk::VectorImage VectorImageType; typedef itk::Image, 3> InitialFitImageType; /** 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 SetGradientDirections( GradientDirectionContainerType * ); void SetBValue(double bval){m_BValue = bval;} void SetBThres(double bval){m_BThres = bval;} void SetS0Thres(double val){m_S0Thres = val;} void SetDStar(double dstar){m_DStar = dstar;} void SetFitDStar(bool fit){m_FitDStar = fit;} void SetVerbose(bool verbose){m_Verbose = verbose;} void SetNumberIterations(int num){m_NumberIterations = num;} void SetLambda(double lambda){m_Lambda = lambda;} void SetCrossPosition(typename InputImageType::IndexType crosspos){this->m_CrossPosition = crosspos;} void SetMethod(IVIM_Method method){m_Method = method;} IVIMSnapshot GetSnapshot(){return m_Snap;} /** 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 ); } protected: DiffusionIntravoxelIncoherentMotionReconstructionImageFilter(); ~DiffusionIntravoxelIncoherentMotionReconstructionImageFilter() {}; void PrintSelf(std::ostream& os, Indent indent) const; void BeforeThreadedGenerateData(); void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, int); void AfterThreadedGenerateData(); MeasAndBvals ApplyS0Threshold(vnl_vector &meas, vnl_vector &bvals); private: double myround(double number); /** container to hold gradient directions */ GradientDirectionContainerType::Pointer m_GradientDirectionContainer; /** Number of gradient measurements */ unsigned int m_NumberOfGradientDirections; double m_BValue; double m_BThres; double m_S0Thres; double m_DStar; IVIM_Method m_Method; typename OutputImageType::Pointer m_DMap; typename OutputImageType::Pointer m_DStarMap; bool m_FitDStar; IVIMSnapshot m_Snap; bool m_Verbose; typename VectorImageType::Pointer m_InternalVectorImage; typename InitialFitImageType::Pointer m_InitialFitImage; int m_NumberIterations; // for total variation vnl_vector m_tmp_allmeas; double m_Lambda; typename InputImageType::IndexType m_CrossPosition; }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.cpp" #endif #endif //__itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter_h diff --git a/Modules/DiffusionImaging/Reconstruction/itkDiffusionQballReconstructionImageFilter.txx b/Modules/DiffusionImaging/Reconstruction/itkDiffusionQballReconstructionImageFilter.txx index a8e6683fd2..cdc17ee208 100644 --- a/Modules/DiffusionImaging/Reconstruction/itkDiffusionQballReconstructionImageFilter.txx +++ b/Modules/DiffusionImaging/Reconstruction/itkDiffusionQballReconstructionImageFilter.txx @@ -1,840 +1,841 @@ /*========================================================================= -Program: Insight Segmentation & Registration Toolkit +Program: Medical Imaging & Interaction Toolkit Language: C++ -Date: $Date: 2006-07-19 15:11:41 $ -Version: $Revision: 1.11 $ +Date: $Date: 2009-07-14 19:11:20 +0200 (Tue, 14 Jul 2009) $ +Version: $Revision: 18127 $ -Copyright (c) Insight Software Consortium. All rights reserved. -See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details. +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 +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 __itkDiffusionQballReconstructionImageFilter_txx #define __itkDiffusionQballReconstructionImageFilter_txx #include "itkDiffusionQballReconstructionImageFilter.h" #include "itkImageRegionConstIterator.h" #include "itkImageRegionConstIteratorWithIndex.h" #include "itkImageRegionIterator.h" #include "itkArray.h" #include "vnl/vnl_vector.h" #include "itkPointShell.h" namespace itk { #define QBALL_RECON_PI 3.14159265358979323846 template< class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NrOdfDirections, int NrBasisFunctionCenters> DiffusionQballReconstructionImageFilter< TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType, NrOdfDirections, NrBasisFunctionCenters> ::DiffusionQballReconstructionImageFilter() : m_GradientDirectionContainer(NULL), m_NumberOfGradientDirections(0), m_NumberOfEquatorSamplingPoints(0), m_NumberOfBaselineImages(1), m_Threshold(NumericTraits< ReferencePixelType >::NonpositiveMin()), m_BValue(1.0), m_GradientImageTypeEnumeration(Else), m_DirectionsDuplicated(false) { // At least 1 inputs is necessary for a vector image. // For images added one at a time we need at least six this->SetNumberOfRequiredInputs( 1 ); } template< class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NrOdfDirections, int NrBasisFunctionCenters> void DiffusionQballReconstructionImageFilter< TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType, NrOdfDirections, NrBasisFunctionCenters> ::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" ); } // If there is only 1 gradient image, it must be an itk::VectorImage. Otherwise // we must have a container of (numberOfInputs-1) itk::Image. Check to make sure if ( numberOfInputs == 1 && m_GradientImageTypeEnumeration != GradientIsInASingleImage ) { 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 ); } } // Compute reconstruction matrix that is multiplied to the data-vector // each voxel in order to reconstruct the ODFs this->ComputeReconstructionMatrix(); // Allocate the b-zero image m_BZeroImage = BZeroImageType::New(); if( m_GradientImageTypeEnumeration == GradientIsInManyImages ) { typename ReferenceImageType::Pointer img = static_cast< ReferenceImageType * >(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() ); } // The gradients are specified in a single multi-component image else if( m_GradientImageTypeEnumeration == GradientIsInASingleImage ) { 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 TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NrOdfDirections, int NrBasisFunctionCenters> typename itk::DiffusionQballReconstructionImageFilter::OdfPixelType itk::DiffusionQballReconstructionImageFilter::Normalize( OdfPixelType odf, typename NumericTraits::AccumulateType b0 ) { switch( m_NormalizationMethod ) { // divide by sum to retreive a PDF case QBR_STANDARD: { odf.Normalize(); return odf; break; } // ADC style case QBR_B_ZERO_B_VALUE: { for(int i=0; i vnl_vector itk::DiffusionQballReconstructionImageFilter::PreNormalize( vnl_vector vec ) { switch( m_NormalizationMethod ) { // standard: no normalization before reconstruction case QBR_STANDARD: { return vec; break; } // log of signal case QBR_B_ZERO_B_VALUE: { int n = vec.size(); for(int i=0; i void DiffusionQballReconstructionImageFilter< TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType, NrOdfDirections, NrBasisFunctionCenters> ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int ) { // init output and b-zero iterators 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(); vnl_vector B(m_NumberOfGradientDirections); // Two cases here: // 1. Gradients specified in multiple images // 'n' iterators for each of the gradient images // 2. Gradients specified in a single multi-component image // one iterator for all gradient directions if( m_GradientImageTypeEnumeration == GradientIsInManyImages ) { // b-zero reference image iterator ImageRegionConstIterator< ReferenceImageType > it(static_cast< ReferenceImageType * >(this->ProcessObject::GetInput(0)), outputRegionForThread); it.GoToBegin(); // fill vector with gradient iterators typedef ImageRegionConstIterator< GradientImageType > GradientIteratorType; std::vector< GradientIteratorType * > gradientItContainer; for( unsigned int i = 1; i<= m_NumberOfGradientDirections; i++ ) { // adapt index in case we have a duplicated shell int index = i; if(m_DirectionsDuplicated) index = i % (m_NumberOfGradientDirections/2); // init and pushback current input image iterator typename GradientImageType::Pointer gradientImagePointer = NULL; // dynamic_cast would be nice, static because of SGI gradientImagePointer = static_cast< GradientImageType * >( this->ProcessObject::GetInput(index) ); GradientIteratorType *git = new GradientIteratorType( gradientImagePointer, outputRegionForThread ); git->GoToBegin(); gradientItContainer.push_back(git); } // Following loop does the actual reconstruction work in each voxel // (Tuch, Q-Ball Reconstruction [1]) while( !it.IsAtEnd() ) { // b-zero reference value ReferencePixelType b0 = it.Get(); // init ODF OdfPixelType odf(0.0); // threshold on reference value to suppress noisy regions if( (b0 != 0) && (b0 >= m_Threshold) ) { // fill array of diffusion measurements for( unsigned int i = 0; i< m_NumberOfGradientDirections; i++ ) { GradientPixelType b = gradientItContainer[i]->Get(); B[i] = static_cast(b); ++(*gradientItContainer[i]); } // pre-normalization according to m_NormalizationMethod B = PreNormalize(B); // actual reconstruction odf = ( (*m_ReconstructionMatrix) * B ).data_block(); // post-normalization according to m_NormalizationMethod odf.Normalize(); } else { // in case we fall below threshold, we just increment to next voxel for( unsigned int i = 0; i< m_NumberOfGradientDirections; i++ ) { ++(*gradientItContainer[i]); } } // set and increment output iterators oit.Set( odf ); ++oit; oit2.Set( b0 ); ++oit2; ++it; } // clean up for( unsigned int i = 0; i< gradientItContainer.size(); i++ ) { delete gradientItContainer[i]; } } // The gradients are specified in a single multi-component image else if( m_GradientImageTypeEnumeration == GradientIsInASingleImage ) { // init input iterator typedef ImageRegionConstIterator< GradientImagesType > GradientIteratorType; typedef typename GradientImagesType::PixelType GradientVectorType; typename GradientImagesType::Pointer gradientImagePointer = NULL; // dynamic_cast would be nice, static because of SGI gradientImagePointer = static_cast< GradientImagesType * >( this->ProcessObject::GetInput(0) ); GradientIteratorType git(gradientImagePointer, outputRegionForThread ); git.GoToBegin(); // set of indicies each for the baseline images and gradient images std::vector baselineind; // contains baseline indicies std::vector gradientind; // contains gradient indicies 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()); } } // in case we have a duplicated shell, we also duplicate or indices if( m_DirectionsDuplicated ) { int gradIndSize = gradientind.size(); for(int i=0; i::AccumulateType b0 = NumericTraits::Zero; for(unsigned int i = 0; i < baselineind.size(); ++i) { b0 += b[baselineind[i]]; } b0 /= this->m_NumberOfBaselineImages; // init resulting ODF OdfPixelType odf(0.0); // threshold on reference value to suppress noisy regions if( (b0 != 0) && (b0 >= m_Threshold) ) { for( unsigned int i = 0; i< m_NumberOfGradientDirections; i++ ) { B[i] = static_cast(b[gradientind[i]]); } // pre-normalization according to m_NormalizationMethod B = PreNormalize(B); // actual reconstruction odf = ( (*m_ReconstructionMatrix) * B ).data_block(); // post-normalization according to m_NormalizationMethod odf = Normalize(odf, b0); } // set and increment output iterators oit.Set( odf ); ++oit; oit2.Set( b0 ); ++oit2; ++git; // Gradient image iterator } } std::cout << "One Thread finished reconstruction" << std::endl; } template< class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NrOdfDirections, int NrBasisFunctionCenters> void DiffusionQballReconstructionImageFilter< TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType, NrOdfDirections, NrBasisFunctionCenters> ::ComputeReconstructionMatrix() { 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; } } // set default number of equator sampling points if needed if(!this->m_NumberOfEquatorSamplingPoints) this->m_NumberOfEquatorSamplingPoints = (int) ceil((double)sqrt(8*QBALL_RECON_PI*this->m_NumberOfGradientDirections)); vnl_matrix* Q = new vnl_matrix(3, m_NumberOfGradientDirections); { // Fill matrix Q with gradient directions int i = 0; for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin(); gdcit != this->m_GradientDirectionContainer->End(); ++gdcit) { if(gdcit.Value().one_norm() > 0.0) { (*Q)(0,i) = gdcit.Value().get(0); (*Q)(1,i) = gdcit.Value().get(1); (*Q)(2,i++) = gdcit.Value().get(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) { (*Q)(0,i) = -gdcit.Value().get(0); (*Q)(1,i) = -gdcit.Value().get(1); (*Q)(2,i++) = -gdcit.Value().get(2); } } } } vnl_matrix_fixed* U = itk::PointShell >::DistributePointShell(); vnl_matrix_fixed* V = itk::PointShell >::DistributePointShell(); // calculate sampling points on the equator perpendicular to z-axis vnl_matrix *C = new vnl_matrix(3, m_NumberOfEquatorSamplingPoints); for(unsigned int i=0; i::Zero; } // rotate the sampling points to each directions of interest vnl_matrix *S = new vnl_matrix(3,m_NumberOfEquatorSamplingPoints*NOdfDirections); { vnl_vector_fixed z(NumericTraits::Zero); z.put(2,NumericTraits::One); vnl_matrix_fixed eye(NumericTraits::Zero); eye.fill_diagonal(NumericTraits::One); for(int i=0; i ui = (*U).get_column(i); vnl_matrix *RC = new vnl_matrix(3,m_NumberOfEquatorSamplingPoints); if( (z(0)*ui(0)+z(1)*ui(1)+z(2)*ui(2)+1) != 0 ) { vnl_matrix_fixed R; R.set_column(0, (z+ui)*(z(0)+ui(0))); R.set_column(1, (z+ui)*(z(1)+ui(1))); R.set_column(2, (z+ui)*(z(2)+ui(2))); R /= (z(0)*ui(0)+z(1)*ui(1)+z(2)*ui(2)+1); R -= eye; (*RC) = R*(*C); } else { RC = C; } (*S).set_columns(i*m_NumberOfEquatorSamplingPoints, *RC); } } // determine interpolation kernel width first // use to determine diffusion measurement contribution to each of the kernels vnl_matrix *H_plus = new vnl_matrix(NBasisFunctionCenters,m_NumberOfGradientDirections); double maxSigma = QBALL_RECON_PI/6; double bestSigma = maxSigma; { double stepsize = 0.01; double start = 0.01; double minCondition = NumericTraits::max(); vnl_matrix *H = new vnl_matrix(m_NumberOfGradientDirections,NBasisFunctionCenters,(double)0); { int increasing = 0; for( double sigma=start; sigma *tmpH = new vnl_matrix(m_NumberOfGradientDirections,NBasisFunctionCenters); for(unsigned int r=0; r1.0) ? 1.0 : qtv); double x = acos(qtv); (*tmpH)(r,c) = (1.0/(sigma*sqrt(2.0*QBALL_RECON_PI))) *exp((-x*x)/(2*sigma*sigma)); } } vnl_svd *solver; if(m_NumberOfGradientDirections>NBasisFunctionCenters) { solver = new vnl_svd(*tmpH); } else { solver = new vnl_svd(tmpH->transpose()); } double condition = solver->sigma_max() / solver->sigma_min(); std::cout << sigma << ": " << condition << std::endl; if( condition < minCondition ) { minCondition = condition; bestSigma = sigma; H->update(*tmpH); } else { // optimum assumed to be hit after condition increased 3 times if (++increasing>3) break; } } } vnl_matrix_inverse *pseudoInverse = new vnl_matrix_inverse(*H); (*H_plus) = pseudoInverse->pinverse(); std::cout << "choosing sigma = " << bestSigma << std::endl; } // this is the contribution of each kernel to each sampling point on the // equator vnl_matrix *G = new vnl_matrix(m_NumberOfEquatorSamplingPoints*NOdfDirections,NBasisFunctionCenters); { for(unsigned int r=0; r1.0) ? 1.0 : stv); double x = acos(stv); (*G)(r,c) = (1.0/(bestSigma*sqrt(2.0*QBALL_RECON_PI))) *exp((-x*x)/(2*bestSigma*bestSigma)); } } } vnl_matrix *GH_plus = new vnl_matrix(m_NumberOfEquatorSamplingPoints*NOdfDirections,m_NumberOfGradientDirections); // simple matrix multiplication, manual cause of stack overflow using operator for (unsigned i = 0; i < m_NumberOfEquatorSamplingPoints*NOdfDirections; ++i) for (unsigned j = 0; j < m_NumberOfGradientDirections; ++j) { double accum = (*G)(i,0) * (*H_plus)(0,j); for (unsigned k = 1; k < NOdfDirections; ++k) accum += (*G)(i,k) * (*H_plus)(k,j); (*GH_plus)(i,j) = accum; } typename vnl_matrix::iterator it3; for( it3 = (*GH_plus).begin(); it3 != (*GH_plus).end(); it3++) { if(*it3<0.0) *it3 = 0; } // this is an addition to the original tuch algorithm for(unsigned int i=0; i r = GH_plus->get_row(i); r /= r.sum(); GH_plus->set_row(i,r); } m_ReconstructionMatrix = new vnl_matrix(NOdfDirections,m_NumberOfGradientDirections,0.0); for(int i=0; i r = m_ReconstructionMatrix->get_row(i); r /= r.sum(); m_ReconstructionMatrix->set_row(i,r); } std::cout << "Reconstruction Matrix computed." << std::endl; } template< class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NrOdfDirections, int NrBasisFunctionCenters> void DiffusionQballReconstructionImageFilter< TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType, NrOdfDirections, NrBasisFunctionCenters> ::AddGradientImage( const GradientDirectionType &gradientDirection, const GradientImageType *gradientImage ) { // Make sure crazy users did not call both AddGradientImage and // SetGradientImage if( m_GradientImageTypeEnumeration == GradientIsInASingleImage) { itkExceptionMacro( << "Cannot call both methods:" << "AddGradientImage and SetGradientImage. Please call only one of them."); } // If the container to hold the gradient directions hasn't been allocated // yet, allocate it. if( !this->m_GradientDirectionContainer ) { this->m_GradientDirectionContainer = GradientDirectionContainerType::New(); } this->m_NumberOfGradientDirections = m_GradientDirectionContainer->Size(); m_GradientDirectionContainer->InsertElement( this->m_NumberOfGradientDirections, gradientDirection / gradientDirection.two_norm() ); this->ProcessObject::SetNthInput( this->m_NumberOfGradientDirections, const_cast< GradientImageType* >(gradientImage) ); m_GradientImageTypeEnumeration = GradientIsInManyImages; } template< class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NrOdfDirections, int NrBasisFunctionCenters> void DiffusionQballReconstructionImageFilter< TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType, NrOdfDirections, NrBasisFunctionCenters> ::SetGradientImage( GradientDirectionContainerType *gradientDirection, const GradientImagesType *gradientImage ) { // Make sure crazy users did not call both AddGradientImage and // SetGradientImage if( m_GradientImageTypeEnumeration == GradientIsInManyImages ) { itkExceptionMacro( << "Cannot call both methods:" << "AddGradientImage and SetGradientImage. Please call only one of them."); } 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) ); m_GradientImageTypeEnumeration = GradientIsInASingleImage; } template< class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NrOdfDirections, int NrBasisFunctionCenters> void DiffusionQballReconstructionImageFilter< TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType, NrOdfDirections, NrBasisFunctionCenters> ::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; if ( this->m_GradientImageTypeEnumeration == GradientIsInManyImages ) { os << indent << "Gradient images haven been supplied " << std::endl; } else if ( this->m_GradientImageTypeEnumeration == GradientIsInManyImages ) { os << indent << "A multicomponent gradient image has been supplied" << std::endl; } os.imbue( originalLocale ); } } #endif // __itkDiffusionQballReconstructionImageFilter_txx diff --git a/Modules/DiffusionImaging/Reconstruction/itkRegularizedIVIMLocalVariationImageFilter.txx b/Modules/DiffusionImaging/Reconstruction/itkRegularizedIVIMLocalVariationImageFilter.txx index a9f2b39aa4..dbeb2df29a 100644 --- a/Modules/DiffusionImaging/Reconstruction/itkRegularizedIVIMLocalVariationImageFilter.txx +++ b/Modules/DiffusionImaging/Reconstruction/itkRegularizedIVIMLocalVariationImageFilter.txx @@ -1,192 +1,193 @@ /*========================================================================= -Program: Insight Segmentation & Registration Toolkit +Program: Medical Imaging & Interaction Toolkit Language: C++ -Date: $Date: 2006-01-11 19:43:31 $ -Version: $Revision: x $ +Date: $Date: 2009-07-14 19:11:20 +0200 (Tue, 14 Jul 2009) $ +Version: $Revision: 18127 $ -Copyright (c) Insight Software Consortium. All rights reserved. -See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details. +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 +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 _itkRegularizedIVIMLocalVariationImageFilter_txx #define _itkRegularizedIVIMLocalVariationImageFilter_txx #include "itkConstShapedNeighborhoodIterator.h" #include "itkNeighborhoodInnerProduct.h" #include "itkImageRegionIterator.h" #include "itkNeighborhoodAlgorithm.h" #include "itkZeroFluxNeumannBoundaryCondition.h" #include "itkOffset.h" #include "itkProgressReporter.h" #include "itkVectorImage.h" #include #include namespace itk { template RegularizedIVIMLocalVariationImageFilter ::RegularizedIVIMLocalVariationImageFilter() {} template void RegularizedIVIMLocalVariationImageFilter ::GenerateInputRequestedRegion() throw (InvalidRequestedRegionError) { // call the superclass' implementation of this method Superclass::GenerateInputRequestedRegion(); // get pointers to the input and output typename Superclass::InputImagePointer inputPtr = const_cast< TInputImage * >( this->GetInput() ); typename Superclass::OutputImagePointer outputPtr = this->GetOutput(); if ( !inputPtr || !outputPtr ) { return; } // get a copy of the input requested region (should equal the output // requested region) typename TInputImage::RegionType inputRequestedRegion; inputRequestedRegion = inputPtr->GetRequestedRegion(); // pad the input requested region by 1 inputRequestedRegion.PadByRadius( 1 ); // crop the input requested region at the input's largest possible region if ( inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()) ) { inputPtr->SetRequestedRegion( inputRequestedRegion ); return; } else { // Couldn't crop the region (requested region is outside the largest // possible region). Throw an exception. // store what we tried to request (prior to trying to crop) inputPtr->SetRequestedRegion( inputRequestedRegion ); // build an exception InvalidRequestedRegionError e(__FILE__, __LINE__); e.SetLocation(ITK_LOCATION); e.SetDescription("Requested region outside possible region."); e.SetDataObject(inputPtr); throw e; } } template< class TInputImage, class TOutputImage> void RegularizedIVIMLocalVariationImageFilter< TInputImage, TOutputImage> ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId) { // Allocate output typename OutputImageType::Pointer output = this->GetOutput(); typename InputImageType::ConstPointer input = this->GetInput(); itk::Size size; for( int i=0; i bC; typename NeighborhoodAlgorithm:: ImageBoundaryFacesCalculator::FaceListType faceList = bC(input, outputRegionForThread, size); // support progress methods/callbacks ProgressReporter progress( this, threadId, outputRegionForThread.GetNumberOfPixels()); ZeroFluxNeumannBoundaryCondition nbc; std::vector pixels; // Process each of the boundary faces. These are N-d regions which border // the edge of the buffer. for ( typename NeighborhoodAlgorithm:: ImageBoundaryFacesCalculator::FaceListType::iterator fit=faceList.begin(); fit != faceList.end(); ++fit) { // iterators over output and input ImageRegionIterator output_image_it(output, *fit); ImageRegionConstIterator input_image_it(input.GetPointer(), *fit); // neighborhood iterator for input image ConstShapedNeighborhoodIterator input_image_neighbors_it(size, input, *fit); typename ConstShapedNeighborhoodIterator:: OffsetType offset; input_image_neighbors_it.OverrideBoundaryCondition(&nbc); input_image_neighbors_it.ClearActiveList(); for(int i=0; i:: ConstIterator input_neighbors_it; for (input_neighbors_it = input_image_neighbors_it.Begin(); ! input_neighbors_it.IsAtEnd(); input_neighbors_it++) { typename TInputImage::PixelType diffVec = input_neighbors_it.Get()-input_image_it.Get(); locVariation += IVIMSquaredEuclideanMetric ::Calc(diffVec); } locVariation = sqrt(locVariation + 0.0001); output_image_it.Set(locVariation); // update iterators ++input_image_neighbors_it; ++output_image_it; ++input_image_it; // report progress progress.CompletedPixel(); } } } /** * Standard "PrintSelf" method */ template void RegularizedIVIMLocalVariationImageFilter ::PrintSelf( std::ostream& os, Indent indent) const { Superclass::PrintSelf( os, indent ); } } // end namespace itk #endif //_itkRegularizedIVIMLocalVariationImageFilter_txx diff --git a/Modules/DiffusionImaging/Reconstruction/itkRegularizedIVIMReconstructionFilter.txx b/Modules/DiffusionImaging/Reconstruction/itkRegularizedIVIMReconstructionFilter.txx index eedc7ebdb7..fd0473a0bf 100644 --- a/Modules/DiffusionImaging/Reconstruction/itkRegularizedIVIMReconstructionFilter.txx +++ b/Modules/DiffusionImaging/Reconstruction/itkRegularizedIVIMReconstructionFilter.txx @@ -1,108 +1,109 @@ /*========================================================================= -Program: Insight Segmentation & Registration Toolkit +Program: Medical Imaging & Interaction Toolkit Language: C++ -Date: $Date: 2006-01-11 19:43:31 $ -Version: $Revision: x $ +Date: $Date: 2009-07-14 19:11:20 +0200 (Tue, 14 Jul 2009) $ +Version: $Revision: 18127 $ -Copyright (c) Insight Software Consortium. All rights reserved. -See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details. +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 +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 _itkRegularizedIVIMReconstructionFilter_txx #define _itkRegularizedIVIMReconstructionFilter_txx #include "itkConstShapedNeighborhoodIterator.h" #include "itkNeighborhoodInnerProduct.h" #include "itkImageRegionIterator.h" #include "itkImageRegionConstIterator.h" #include "itkNeighborhoodAlgorithm.h" #include "itkZeroFluxNeumannBoundaryCondition.h" #include "itkOffset.h" #include "itkProgressReporter.h" #include #include namespace itk { template RegularizedIVIMReconstructionFilter ::RegularizedIVIMReconstructionFilter() { m_Lambda = 1.0; } template void RegularizedIVIMReconstructionFilter ::GenerateData() { // first we cast the input image to match output type typename CastType::Pointer infilter = CastType::New(); infilter->SetInput(this->GetInput()); infilter->Update(); typename OutputImageType::Pointer image = infilter->GetOutput(); typename SingleIterationFilterType::Pointer filter; for(int i=0; iSetInput( image.GetPointer() ); filter->SetOriginalImage( m_ReferenceImage ); filter->SetBValues(m_BValues); filter->SetLambda(m_Lambda); filter->SetNumberOfThreads(this->GetNumberOfThreads()); filter->UpdateLargestPossibleRegion(); image = filter->GetOutput(); std::cout << "Iteration " << i+1 << "/" << m_NumberIterations << std::endl; } typename OutputImageType::Pointer output = this->GetOutput(); output->SetOrigin( image->GetOrigin() ); // Set the image origin output->SetDirection( image->GetDirection() ); // Set the image direction output->SetSpacing(image->GetSpacing()); output->SetLargestPossibleRegion( image->GetLargestPossibleRegion() ); output->SetBufferedRegion( image->GetLargestPossibleRegion() ); output->Allocate(); itk::ImageRegionIterator oit( output, output->GetLargestPossibleRegion()); oit.GoToBegin(); itk::ImageRegionConstIterator iit( image, image->GetLargestPossibleRegion()); iit.GoToBegin(); while(!oit.IsAtEnd()) { oit.Set(iit.Get()); ++iit; ++oit; } } /** * Standard "PrintSelf" method */ template void RegularizedIVIMReconstructionFilter ::PrintSelf( std::ostream& os, Indent indent) const { Superclass::PrintSelf( os, indent ); } } // end namespace itk #endif diff --git a/Modules/DiffusionImaging/Reconstruction/itkRegularizedIVIMReconstructionSingleIteration.txx b/Modules/DiffusionImaging/Reconstruction/itkRegularizedIVIMReconstructionSingleIteration.txx index 5564ab77fd..fdae01c8da 100644 --- a/Modules/DiffusionImaging/Reconstruction/itkRegularizedIVIMReconstructionSingleIteration.txx +++ b/Modules/DiffusionImaging/Reconstruction/itkRegularizedIVIMReconstructionSingleIteration.txx @@ -1,310 +1,311 @@ /*========================================================================= -Program: Insight Segmentation & Registration Toolkit +Program: Medical Imaging & Interaction Toolkit Language: C++ -Date: $Date: 2006-01-11 19:43:31 $ -Version: $Revision: x $ +Date: $Date: 2009-07-14 19:11:20 +0200 (Tue, 14 Jul 2009) $ +Version: $Revision: 18127 $ -Copyright (c) Insight Software Consortium. All rights reserved. -See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details. +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 +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 _itkRegularizedIVIMReconstructionSingleIteration_txx #define _itkRegularizedIVIMReconstructionSingleIteration_txx // itk includes #include "itkConstShapedNeighborhoodIterator.h" #include "itkNeighborhoodInnerProduct.h" #include "itkImageRegionIterator.h" #include "itkNeighborhoodAlgorithm.h" #include "itkZeroFluxNeumannBoundaryCondition.h" #include "itkOffset.h" #include "itkProgressReporter.h" #include "itkRegularizedIVIMLocalVariationImageFilter.h" // other includes #include #include #define IVIM_FOO -100000 namespace itk { /** * constructor */ template RegularizedIVIMReconstructionSingleIteration ::RegularizedIVIMReconstructionSingleIteration() { m_Lambda = 1.0; m_LocalVariation = LocalVariationImageType::New(); } /** * generate requested region */ template void RegularizedIVIMReconstructionSingleIteration ::GenerateInputRequestedRegion() throw (InvalidRequestedRegionError) { // call the superclass' implementation of this method Superclass::GenerateInputRequestedRegion(); // get pointers to the input and output typename Superclass::InputImagePointer inputPtr = const_cast< InputImageType * >( this->GetInput() ); typename Superclass::OutputImagePointer outputPtr = this->GetOutput(); if ( !inputPtr || !outputPtr ) { return; } // get a copy of the input requested region (should equal the output // requested region) typename InputImageType::RegionType inputRequestedRegion; inputRequestedRegion = inputPtr->GetRequestedRegion(); // pad the input requested region by 1 inputRequestedRegion.PadByRadius( 1 ); // crop the input requested region at the input's largest possible region if ( inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()) ) { inputPtr->SetRequestedRegion( inputRequestedRegion ); return; } else { // Couldn't crop the region (requested region is outside the largest // possible region). Throw an exception. // store what we tried to request (prior to trying to crop) inputPtr->SetRequestedRegion( inputRequestedRegion ); // build an exception InvalidRequestedRegionError e(__FILE__, __LINE__); e.SetLocation(ITK_LOCATION); e.SetDescription("Requested region outside possible region."); e.SetDataObject(inputPtr); throw e; } } /** * generate output */ template void RegularizedIVIMReconstructionSingleIteration ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId) { typename OutputImageType::Pointer output = this->GetOutput(); typename InputImageType::ConstPointer input = this->GetInput(); // Find the data-set boundary "faces" itk::Size size; for( int i=0; i bC; typename NeighborhoodAlgorithm:: ImageBoundaryFacesCalculator::FaceListType faceList = bC(input, outputRegionForThread, size); NeighborhoodAlgorithm:: ImageBoundaryFacesCalculator lv_bC; typename NeighborhoodAlgorithm:: ImageBoundaryFacesCalculator::FaceListType lv_faceList = lv_bC(m_LocalVariation, outputRegionForThread, size); ZeroFluxNeumannBoundaryCondition nbc; ZeroFluxNeumannBoundaryCondition lv_nbc; std::vector ws; std::vector hs; typename NeighborhoodAlgorithm:: ImageBoundaryFacesCalculator::FaceListType::iterator lv_fit=lv_faceList.begin(); // Process each of the boundary faces. These are N-d regions which border // the edge of the buffer. for ( typename NeighborhoodAlgorithm:: ImageBoundaryFacesCalculator::FaceListType::iterator fit=faceList.begin(); fit != faceList.end(); ++fit) { // iterators over output, input, original and local variation image ImageRegionIterator output_image_it = ImageRegionIterator(output, *fit); ImageRegionConstIterator input_image_it = ImageRegionConstIterator(input, *fit); ImageRegionConstIterator orig_image_it = ImageRegionConstIterator(m_OriginalImage, *fit); ImageRegionConstIterator loc_var_image_it = ImageRegionConstIterator( m_LocalVariation, *fit); // neighborhood in input image ConstShapedNeighborhoodIterator input_image_neighbors_it(size, input, *fit); typename ConstShapedNeighborhoodIterator:: OffsetType offset; input_image_neighbors_it.OverrideBoundaryCondition(&nbc); input_image_neighbors_it.ClearActiveList(); for(int i=0; i loc_var_image_neighbors_it(size, m_LocalVariation, *lv_fit); loc_var_image_neighbors_it.OverrideBoundaryCondition(&lv_nbc); loc_var_image_neighbors_it.ClearActiveList(); for(int i=0; i:: ConstIterator loc_var_neighbors_it; for (loc_var_neighbors_it = loc_var_image_neighbors_it.Begin(); ! loc_var_neighbors_it.IsAtEnd(); loc_var_neighbors_it++) { // w_alphabeta(u) = // 1 / ||nabla_alpha(u)||_a + 1 / ||nabla_beta(u)||_a ws[count] = locvar_alpha_inv + (1.0/(double)loc_var_neighbors_it.Get()); wsum += ws[count++]; //MITK_INFO << "nb var: " << count << ": " << loc_var_neighbors_it.Get(); } //MITK_INFO << "wsum: " << wsum; // h_alphaalpha * u_alpha^zero typename RefImageType::PixelType orig = orig_image_it.Get(); // vnl_vector estim(orig.GetSize()); // vnl_matrix diff(orig.GetSize(),1); // vnl_matrix estimdash(orig.GetSize(),2); vnl_vector_fixed step; step[0] = 0; step[1]=0; for(int ind=0; ind:: ConstIterator input_neighbors_it; for (input_neighbors_it = input_image_neighbors_it.Begin(); ! input_neighbors_it.IsAtEnd(); input_neighbors_it++) { step[1] += (input_neighbors_it.Get()[1] - input_image_it.Get()[1]) * (ws[count++] / (m_Lambda+wsum)); } //MITK_INFO << "stepfinal: " << step[0] << "; " << step[1]; // set output result OutputPixelType out; out[0] = input_image_it.Get()[0] + .001*step[0]; out[1] = input_image_it.Get()[1] + .00001*step[1]; output_image_it.Set( out ); //MITK_INFO << "(" << input_image_it.Get()[0] << " ; " << input_image_it.Get()[1] << ") => (" << out[0] << " ; " << out[1] << ")"; // increment iterators ++output_image_it; ++input_image_it; ++orig_image_it; ++loc_var_image_it; ++input_image_neighbors_it; ++loc_var_image_neighbors_it; } ++lv_fit; } } /** * first calculate local variation in the image */ template void RegularizedIVIMReconstructionSingleIteration ::BeforeThreadedGenerateData() { typedef typename itk::RegularizedIVIMLocalVariationImageFilter FilterType; typename FilterType::Pointer filter = FilterType::New(); filter->SetInput(this->GetInput(0)); filter->SetNumberOfThreads(this->GetNumberOfThreads()); filter->Update(); this->m_LocalVariation = filter->GetOutput(); } /** * Standard "PrintSelf" method */ template void RegularizedIVIMReconstructionSingleIteration ::PrintSelf( std::ostream& os, Indent indent) const { Superclass::PrintSelf( os, indent ); } } // end namespace itk #endif diff --git a/Modules/DiffusionImaging/Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.cpp b/Modules/DiffusionImaging/Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.cpp index 261d6c0c78..72032212fa 100644 --- a/Modules/DiffusionImaging/Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.cpp +++ b/Modules/DiffusionImaging/Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.cpp @@ -1,228 +1,224 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Language: C++ Date: $Date: 2007-12-11 14:46:19 +0100 (Di, 11 Dez 2007) $ Version: $Revision: 13129 $ 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 #include #include #include #include #include "mitkPixelType.h" #include "itkImageRegionIterator.h" #include "itkImageFileReader.h" #include "mitkNrrdDiffusionImageWriter.h" template< class D, class T > mitk::TeemDiffusionTensor3DReconstructionImageFilter ::TeemDiffusionTensor3DReconstructionImageFilter(): m_EstimateErrorImage(false),m_Sigma(-19191919), m_EstimationMethod(TeemTensorEstimationMethodsLLS), m_NumIterations(1),m_ConfidenceThreshold(-19191919.0), m_ConfidenceFuzzyness(0.0),m_MinPlausibleValue(1.0) { } template< class D, class T > mitk::TeemDiffusionTensor3DReconstructionImageFilter ::~TeemDiffusionTensor3DReconstructionImageFilter() { } void file_replace(std::string filename, std::string what, std::string with) { ofstream myfile2; std::locale C("C"); std::locale originalLocale2 = myfile2.getloc(); myfile2.imbue(C); char filename2[512]; sprintf(filename2, "%s2",filename.c_str()); myfile2.open (filename2); std::string line; ifstream myfile (filename.c_str()); std::locale originalLocale = myfile.getloc(); myfile.imbue(C); if (myfile.is_open()) { while (! myfile.eof() ) { getline (myfile,line); itksys::SystemTools::ReplaceString(line,what.c_str(),with.c_str()); myfile2 << line << std::endl; } myfile.close(); } myfile2.close(); itksys::SystemTools::RemoveFile(filename.c_str()); rename(filename2,filename.c_str()); myfile.imbue( originalLocale ); myfile2.imbue( originalLocale2 ); } // do the work template< class D, class T > void mitk::TeemDiffusionTensor3DReconstructionImageFilter ::Update() { // save input image to nrrd file in temp-folder char filename[512]; srand((unsigned)time(0)); int random_integer = rand(); sprintf( filename, "dwi_%d.nhdr",random_integer); typedef mitk::NrrdDiffusionImageWriter WriterType; typename WriterType::Pointer nrrdWriter = WriterType::New(); nrrdWriter->SetInput( m_Input ); //nrrdWriter->SetDirections(m_Input->GetDirections()); //nrrdWriter->SetB_Value(m_Input->GetB_Value()); nrrdWriter->SetFileName(filename); try { nrrdWriter->Update(); } catch (itk::ExceptionObject e) { std::cout << e << std::endl; } file_replace(filename,"vector","list"); // build up correct command from input params char command[4096]; sprintf( command, "tend estim -i %s -B kvp -o tensors_%d.nhdr -knownB0 true", filename, random_integer); //m_DiffusionImages; if(m_EstimateErrorImage) { sprintf( command, "%s -ee error_image_%d.nhdr", command, random_integer); } if(m_Sigma != -19191919) { sprintf( command, "%s -sigma %f", command, m_Sigma); } switch(m_EstimationMethod) { case TeemTensorEstimationMethodsLLS: sprintf( command, "%s -est lls", command); break; case TeemTensorEstimationMethodsMLE: sprintf( command, "%s -est mle", command); break; case TeemTensorEstimationMethodsNLS: sprintf( command, "%s -est nls", command); break; case TeemTensorEstimationMethodsWLS: sprintf( command, "%s -est wls", command); break; } sprintf( command, "%s -wlsi %d", command, m_NumIterations); if(m_ConfidenceThreshold != -19191919.0) { sprintf( command, "%s -t %f", command, m_ConfidenceThreshold); } sprintf( command, "%s -soft %f", command, m_ConfidenceFuzzyness); sprintf( command, "%s -mv %f", command, m_MinPlausibleValue); // call tend estim command std::cout << "Calling <" << command << ">" << std::endl; int success = system(command); if(!success) { MITK_ERROR << "system command could not be called!"; } remove(filename); sprintf( filename, "dwi_%d.raw", random_integer); remove(filename); // change kind from tensor to vector sprintf( filename, "tensors_%d.nhdr", random_integer); file_replace(filename,"3D-masked-symmetric-matrix","vector"); - //file_replace(filename,"3D-symmetric-matrix","vector"); - //itksys::SystemTools::ReplaceString(line,"3D-masked-symmetric-matrix","vector"); - //itksys::SystemTools::ReplaceString(line,"vector","domain"); - // read result as mitk::Image and provide it in m_Output typedef itk::ImageFileReader FileReaderType; typename FileReaderType::Pointer reader = FileReaderType::New(); reader->SetFileName(filename); reader->Update(); typename VectorImageType::Pointer vecImage = reader->GetOutput(); remove(filename); sprintf( filename, "tensors_%d.raw", random_integer); remove(filename); typename ItkTensorImageType::Pointer itkTensorImage = ItkTensorImageType::New(); itkTensorImage->SetSpacing( vecImage->GetSpacing() ); // Set the image spacing itkTensorImage->SetOrigin( vecImage->GetOrigin() ); // Set the image origin itkTensorImage->SetDirection( vecImage->GetDirection() ); // Set the image direction itkTensorImage->SetLargestPossibleRegion( vecImage->GetLargestPossibleRegion() ); itkTensorImage->SetBufferedRegion( vecImage->GetLargestPossibleRegion() ); itkTensorImage->SetRequestedRegion( vecImage->GetLargestPossibleRegion() ); itkTensorImage->Allocate(); itk::ImageRegionIterator it(vecImage, vecImage->GetLargestPossibleRegion()); itk::ImageRegionIterator it2(itkTensorImage, itkTensorImage->GetLargestPossibleRegion()); it2 = it2.Begin(); //#pragma omp parallel private (it) { for(it=it.Begin();!it.IsAtEnd(); ++it, ++it2) { //#pragma omp single nowait { VectorType vec = it.Get(); TensorType tensor; for(int i=1;i<7;i++) tensor[i-1] = vec[i] * vec[0]; it2.Set( tensor ); } } // end for } // end ompparallel m_OutputItk = mitk::TensorImage::New(); m_OutputItk->InitializeByItk(itkTensorImage.GetPointer()); m_OutputItk->SetVolume( itkTensorImage->GetBufferPointer() ); // in case: read resulting error-image and provide it in m_ErrorImage if(m_EstimateErrorImage) { // open error image here } }