diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkFiniteDiffOdfMaximaExtractionFilter.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkFiniteDiffOdfMaximaExtractionFilter.cpp index ef720b67d9..579a772c7e 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkFiniteDiffOdfMaximaExtractionFilter.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkFiniteDiffOdfMaximaExtractionFilter.cpp @@ -1,480 +1,498 @@ #ifndef __itkFiniteDiffOdfMaximaExtractionFilter_cpp #define __itkFiniteDiffOdfMaximaExtractionFilter_cpp #include "itkFiniteDiffOdfMaximaExtractionFilter.h" #include #include #include #include #include #include #include #include #include #include #include #include #include using namespace boost::math; namespace itk { static bool CompareVectors(const vnl_vector_fixed< double, 3 >& v1, const vnl_vector_fixed< double, 3 >& v2) { return (v1.magnitude()>v2.magnitude()); } template< class PixelType, int ShOrder, int NrOdfDirections > FiniteDiffOdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections> ::FiniteDiffOdfMaximaExtractionFilter() - : m_MaxNumPeaks(2) + : m_Toolkit(FSL) + , m_MaxNumPeaks(2) , m_PeakThreshold(0.4) , m_ClusteringThreshold(0.9) , m_AngularThreshold(0.7) , m_NumCoeffs((ShOrder*ShOrder + ShOrder + 2)/2 + ShOrder) , m_NormalizationMethod(MAX_VEC_NORM) , m_AbsolutePeakThreshold(0) { this->SetNumberOfRequiredInputs(1); } template< class PixelType, int ShOrder, int NrOdfDirections > void FiniteDiffOdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections> ::FindCandidatePeaks(OdfType& odf, double thr, std::vector< DirectionType >& container) { double gfa = odf.GetGeneralizedFractionalAnisotropy(); //Find the peaks using a finite difference method bool flag = true; vnl_vector_fixed< bool, NrOdfDirections > used; used.fill(false); //Find the peaks for (int i=0; ithr && val*gfa>m_AbsolutePeakThreshold) // limit to one hemisphere ??? { flag = true; std::vector< int > neighbours = odf.GetNeighbors(i); for (int j=0; j std::vector< vnl_vector_fixed< double, 3 > > FiniteDiffOdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections>::MeanShiftClustering(std::vector< vnl_vector_fixed< double, 3 > >& inDirs) { std::vector< DirectionType > outDirs; if (inDirs.empty()) return inDirs; DirectionType oldMean, currentMean, workingMean; std::vector< int > touched; // initialize touched.resize(inDirs.size(), 0); bool free = true; currentMean = inDirs[0]; // initialize first seed while (free) { oldMean.fill(0.0); // start mean-shift clustering float angle = 0.0; int counter = 0; while ((currentMean-oldMean).magnitude()>0.0001) { counter = 0; oldMean = currentMean; workingMean = oldMean; workingMean.normalize(); currentMean.fill(0.0); for (int i=0; i=m_ClusteringThreshold) { currentMean += inDirs[i]; touched[i] = 1; counter++; } else if (-angle>=m_ClusteringThreshold) { currentMean -= inDirs[i]; touched[i] = 1; counter++; } } } // found stable mean if (counter>0) { float mag = currentMean.magnitude(); if (mag>0) { currentMean /= mag; outDirs.push_back(currentMean); } } // find next unused seed free = false; for (int i=0; i void FiniteDiffOdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections> ::BeforeThreadedGenerateData() { typename CoefficientImageType::Pointer ShCoeffImage = static_cast< CoefficientImageType* >( this->ProcessObject::GetInput(0) ); mitk::Vector3D spacing = ShCoeffImage->GetSpacing(); double minSpacing = spacing[0]; if (spacing[1]GetOrigin(); itk::Matrix direction = ShCoeffImage->GetDirection(); ImageRegion<3> imageRegion = ShCoeffImage->GetLargestPossibleRegion(); m_DirectionImageContainer = ItkDirectionImageContainer::New(); for (int i=0; i nullVec; nullVec.Fill(0.0); ItkDirectionImage::Pointer img = ItkDirectionImage::New(); img->SetSpacing( spacing ); img->SetOrigin( origin ); img->SetDirection( direction ); img->SetRegions( imageRegion ); img->Allocate(); img->FillBuffer(nullVec); m_DirectionImageContainer->InsertElement(m_DirectionImageContainer->Size(), img); } if (m_MaskImage.IsNull()) { m_MaskImage = ItkUcharImgType::New(); m_MaskImage->SetSpacing( spacing ); m_MaskImage->SetOrigin( origin ); m_MaskImage->SetDirection( direction ); m_MaskImage->SetRegions( imageRegion ); m_MaskImage->Allocate(); m_MaskImage->FillBuffer(1); } m_NumDirectionsImage = ItkUcharImgType::New(); m_NumDirectionsImage->SetSpacing( spacing ); m_NumDirectionsImage->SetOrigin( origin ); m_NumDirectionsImage->SetDirection( direction ); m_NumDirectionsImage->SetRegions( imageRegion ); m_NumDirectionsImage->Allocate(); m_NumDirectionsImage->FillBuffer(0); this->SetNumberOfOutputs(m_MaxNumPeaks); // calculate SH basis OdfType odf; vnl_matrix_fixed* directions = odf.GetDirections(); vnl_matrix< double > sphCoords; std::vector< DirectionType > dirs; for (int i=0; iget_column(i)); Cart2Sph(dirs, sphCoords); // convert candidate peaks to spherical angles m_ShBasis = CalcShBasis(sphCoords); // evaluate spherical harmonics at each peak MITK_INFO << "Starting finite differences maximum extraction"; MITK_INFO << "ODF sampling points: " << NrOdfDirections; MITK_INFO << "SH order: " << ShOrder; MITK_INFO << "Maximum peaks: " << m_MaxNumPeaks; MITK_INFO << "Relative threshold: " << m_PeakThreshold; MITK_INFO << "Absolute threshold: " << m_AbsolutePeakThreshold; MITK_INFO << "Clustering threshold: " << m_ClusteringThreshold; MITK_INFO << "Angular threshold: " << m_AngularThreshold; } template< class PixelType, int ShOrder, int NrOdfDirections > void FiniteDiffOdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections> ::AfterThreadedGenerateData() { MITK_INFO << "Generating vector field"; vtkSmartPointer m_VtkCellArray = vtkSmartPointer::New(); vtkSmartPointer m_VtkPoints = vtkSmartPointer::New(); typename CoefficientImageType::Pointer ShCoeffImage = static_cast< CoefficientImageType* >( this->ProcessObject::GetInput(0) ); ImageRegionConstIterator< CoefficientImageType > cit(ShCoeffImage, ShCoeffImage->GetLargestPossibleRegion() ); mitk::Vector3D spacing = ShCoeffImage->GetSpacing(); double minSpacing = spacing[0]; if (spacing[1]GetLargestPossibleRegion().GetSize()[0]*ShCoeffImage->GetLargestPossibleRegion().GetSize()[1]*ShCoeffImage->GetLargestPossibleRegion().GetSize()[2]; boost::progress_display disp(maxProgress); while( !cit.IsAtEnd() ) { ++disp; typename CoefficientImageType::IndexType index = cit.GetIndex(); if (m_MaskImage->GetPixel(index)==0) { ++cit; continue; } for (int i=0; iSize(); i++) { ItkDirectionImage::Pointer img = m_DirectionImageContainer->GetElement(i); itk::Vector< float, 3 > pixel = img->GetPixel(index); DirectionType dir; dir[0] = pixel[0]; dir[1] = pixel[1]; dir[2] = pixel[2]; vtkSmartPointer container = vtkSmartPointer::New(); itk::ContinuousIndex center; center[0] = index[0]; center[1] = index[1]; center[2] = index[2]; itk::Point worldCenter; ShCoeffImage->TransformContinuousIndexToPhysicalPoint( center, worldCenter ); itk::Point worldStart; worldStart[0] = worldCenter[0]-dir[0]/2 * minSpacing; worldStart[1] = worldCenter[1]-dir[1]/2 * minSpacing; worldStart[2] = worldCenter[2]-dir[2]/2 * minSpacing; vtkIdType id = m_VtkPoints->InsertNextPoint(worldStart.GetDataPointer()); container->GetPointIds()->InsertNextId(id); itk::Point worldEnd; worldEnd[0] = worldCenter[0]+dir[0]/2 * minSpacing; worldEnd[1] = worldCenter[1]+dir[1]/2 * minSpacing; worldEnd[2] = worldCenter[2]+dir[2]/2 * minSpacing; id = m_VtkPoints->InsertNextPoint(worldEnd.GetDataPointer()); container->GetPointIds()->InsertNextId(id); m_VtkCellArray->InsertNextCell(container); } ++cit; } vtkSmartPointer directionsPolyData = vtkSmartPointer::New(); directionsPolyData->SetPoints(m_VtkPoints); directionsPolyData->SetLines(m_VtkCellArray); m_OutputFiberBundle = mitk::FiberBundleX::New(directionsPolyData); for (int i=0; iSize(); i++) { ItkDirectionImage::Pointer img = m_DirectionImageContainer->GetElement(i); this->SetNthOutput(i, img); } } template< class PixelType, int ShOrder, int NrOdfDirections > void FiniteDiffOdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections> ::ThreadedGenerateData( const OutputImageRegionType& outputRegionForThread, int threadID ) { typename CoefficientImageType::Pointer ShCoeffImage = static_cast< CoefficientImageType* >( this->ProcessObject::GetInput(0) ); ImageRegionConstIterator< CoefficientImageType > cit(ShCoeffImage, outputRegionForThread ); OdfType odf; while( !cit.IsAtEnd() ) { typename CoefficientImageType::IndexType index = cit.GetIndex(); if (m_MaskImage->GetPixel(index)==0) { ++cit; continue; } CoefficientPixelType c = cit.Get(); // calculate ODF double max = 0; odf.Fill(0.0); for (int i=0; imax) max = odf[i]; } if (max<0.0001) { ++cit; continue; } std::vector< DirectionType > candidates, peaks, temp; peaks.clear(); max *= m_PeakThreshold; // relative threshold FindCandidatePeaks(odf, max, candidates); // find all local maxima candidates = MeanShiftClustering(candidates); // cluster maxima vnl_matrix< double > shBasis, sphCoords; Cart2Sph(candidates, sphCoords); // convert candidate peaks to spherical angles shBasis = CalcShBasis(sphCoords); // evaluate spherical harmonics at each peak max = 0.0; for (int i=0; imax) max = val; peaks.push_back(candidates[i]*val); } std::sort( peaks.begin(), peaks.end(), CompareVectors ); // sort peaks // kick out directions to close to a larger direction (too far away to cluster but too close to keep) int m = peaks.size(); if ( m>m_DirectionImageContainer->Size() ) m = m_DirectionImageContainer->Size(); for (int i=0; im_AngularThreshold && valm_DirectionImageContainer->Size() ) num = m_DirectionImageContainer->Size(); for (int i=0; i dir = peaks.at(i); ItkDirectionImage::Pointer img = m_DirectionImageContainer->GetElement(i); switch (m_NormalizationMethod) { case NO_NORM: break; case SINGLE_VEC_NORM: dir.normalize(); break; case MAX_VEC_NORM: dir /= max; break; } dir = ShCoeffImage->GetDirection()*dir; itk::Vector< float, 3 > pixel; pixel.SetElement(0, dir[0]); pixel.SetElement(1, dir[1]); pixel.SetElement(2, dir[2]); img->SetPixel(index, pixel); } m_NumDirectionsImage->SetPixel(index, num); ++cit; } MITK_INFO << "Thread " << threadID << " finished extraction"; } // convert cartesian to spherical coordinates template< class PixelType, int ShOrder, int NrOdfDirections > void FiniteDiffOdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections> ::Cart2Sph(const std::vector< DirectionType >& dir, vnl_matrix& sphCoords) { sphCoords.set_size(dir.size(), 2); for (int i=0; i vnl_matrix FiniteDiffOdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections> ::CalcShBasis(vnl_matrix& sphCoords) { int M = sphCoords.rows(); int j, m; double mag, plm; vnl_matrix shBasis; shBasis.set_size(M, m_NumCoeffs); for (int p=0; p(l,abs(m),cos(sphCoords(p,0))); mag = sqrt((double)(2*l+1)/(4.0*M_PI)*factorial(l-abs(m))/factorial(l+abs(m)))*plm; if (m<0) shBasis(p,j) = sqrt(2.0)*mag*cos(fabs((double)m)*sphCoords(p,1)); else if (m==0) shBasis(p,j) = mag; else shBasis(p,j) = pow(-1.0, m)*sqrt(2.0)*mag*sin(m*sphCoords(p,1)); + break; + case MRTRIX: + + plm = legendre_p(l,abs(m),-cos(sphCoords(p,0))); + mag = sqrt((double)(2*l+1)/(4.0*M_PI)*factorial(l-abs(m))/factorial(l+abs(m)))*plm; + if (m>0) + shBasis(p,j) = mag*cos(m*sphCoords(p,1)); + else if (m==0) + shBasis(p,j) = mag; + else + shBasis(p,j) = mag*sin(-m*sphCoords(p,1)); + break; + } + j++; } } return shBasis; } } #endif // __itkFiniteDiffOdfMaximaExtractionFilter_cpp diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkFiniteDiffOdfMaximaExtractionFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkFiniteDiffOdfMaximaExtractionFilter.h index 6b104af504..948bff45ad 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkFiniteDiffOdfMaximaExtractionFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkFiniteDiffOdfMaximaExtractionFilter.h @@ -1,134 +1,144 @@ /*========================================================================= 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. =========================================================================*/ #ifndef __itkFiniteDiffOdfMaximaExtractionFilter_h_ #define __itkFiniteDiffOdfMaximaExtractionFilter_h_ #include "itkImageToImageFilter.h" #include "vnl/vnl_vector_fixed.h" #include "vnl/vnl_matrix.h" #include "vnl/algo/vnl_svd.h" #include "itkVectorContainer.h" #include "itkVectorImage.h" #include #include namespace itk{ /** * \brief Extract ODF peaks by searching all local maxima on a densely sampled ODF und clustering these maxima to get the underlying fiber direction. * NrOdfDirections: number of sampling points on the ODF surface (about 20000 is a good value) */ template< class PixelType, int ShOrder, int NrOdfDirections > class FiniteDiffOdfMaximaExtractionFilter : public ImageToImageFilter< Image< Vector< PixelType, (ShOrder*ShOrder + ShOrder + 2)/2 + ShOrder >, 3 >, Image< Vector< PixelType, 3 >, 3 > > { public: + enum Toolkit { ///< SH coefficient convention (depends on toolkit) + FSL, + MRTRIX + }; + enum NormalizationMethods { NO_NORM, ///< no length normalization of the output peaks SINGLE_VEC_NORM, ///< normalize the single peaks to length 1 MAX_VEC_NORM ///< normalize all peaks according to their length in comparison to the largest peak (0-1) }; typedef FiniteDiffOdfMaximaExtractionFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ImageToImageFilter< Image< Vector< PixelType, (ShOrder*ShOrder + ShOrder + 2)/2 + ShOrder >, 3 >, Image< Vector< PixelType, 3 >, 3 > > Superclass; /** Method for creation through the object factory. */ itkNewMacro(Self) /** Runtime information support. */ itkTypeMacro(FiniteDiffOdfMaximaExtractionFilter, ImageToImageFilter) typedef typename Superclass::InputImageType CoefficientImageType; typedef typename CoefficientImageType::RegionType CoefficientImageRegionType; typedef typename CoefficientImageType::PixelType CoefficientPixelType; typedef typename Superclass::OutputImageType OutputImageType; typedef typename Superclass::OutputImageRegionType OutputImageRegionType; typedef OrientationDistributionFunction OdfType; typedef itk::Image ItkUcharImgType; typedef vnl_vector_fixed< double, 3 > DirectionType; typedef Image< Vector< float, 3 >, 3> ItkDirectionImage; typedef VectorContainer< unsigned int, ItkDirectionImage::Pointer > ItkDirectionImageContainer; // input itkSetMacro( MaxNumPeaks, unsigned int) ///< maximum number of peaks per voxel. if more peaks are detected, only the largest are kept. itkSetMacro( PeakThreshold, double) ///< threshold on the peak length relative to the largest peak inside the current voxel itkSetMacro( AbsolutePeakThreshold, double) ///< hard threshold on the peak length of all local maxima itkSetMacro( ClusteringThreshold, double) ///< directions closer together than the specified angular threshold will be clustered (in rad) itkSetMacro( AngularThreshold, double) ///< directions closer together than the specified threshold that remain after clustering are discarded (largest is kept) (in rad) itkSetMacro( MaskImage, ItkUcharImgType::Pointer) ///< only voxels inside the binary mask are processed itkSetMacro( NormalizationMethod, NormalizationMethods) ///< normalization method of ODF peaks // output itkGetMacro( OutputFiberBundle, mitk::FiberBundleX::Pointer) ///< vector field (peak sizes rescaled for visualization purposes) itkGetMacro( DirectionImageContainer, ItkDirectionImageContainer::Pointer) ///< container for output peaks itkGetMacro( NumDirectionsImage, ItkUcharImgType::Pointer) ///< number of peaks per voxel + itkSetMacro( Toolkit, Toolkit) ///< define SH coefficient convention (depends on toolkit) + itkGetMacro( Toolkit, Toolkit) ///< SH coefficient convention (depends on toolkit) + protected: FiniteDiffOdfMaximaExtractionFilter(); ~FiniteDiffOdfMaximaExtractionFilter(){} void BeforeThreadedGenerateData(); void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, int threadID ); void AfterThreadedGenerateData(); /** Extract all local maxima from the densely sampled ODF surface. Thresholding possible. **/ void FindCandidatePeaks(OdfType& odf, double odfMax, std::vector< DirectionType >& inDirs); /** Cluster input directions within a certain angular threshold **/ std::vector< DirectionType > MeanShiftClustering(std::vector< DirectionType >& inDirs); /** Convert cartesian to spherical coordinates **/ void Cart2Sph(const std::vector< DirectionType >& dir, vnl_matrix& sphCoords); /** Calculate spherical harmonic basis of the defined order **/ vnl_matrix CalcShBasis(vnl_matrix& sphCoords); private: NormalizationMethods m_NormalizationMethod; ///< normalization method of ODF peaks unsigned int m_MaxNumPeaks; ///< maximum number of peaks per voxel. if more peaks are detected, only the largest are kept. double m_PeakThreshold; ///< threshold on the peak length relative to the largest peak inside the current voxel double m_AbsolutePeakThreshold;///< hard threshold on the peak length of all local maxima vnl_matrix< double > m_ShBasis; ///< container for evaluated SH base functions double m_ClusteringThreshold; ///< directions closer together than the specified angular threshold will be clustered (in rad) double m_AngularThreshold; ///< directions closer together than the specified threshold that remain after clustering are discarded (largest is kept) (in rad) const int m_NumCoeffs; ///< number of spherical harmonics coefficients mitk::FiberBundleX::Pointer m_OutputFiberBundle; ///< vector field (peak sizes rescaled for visualization purposes) ItkDirectionImageContainer::Pointer m_DirectionImageContainer;///< container for output peaks ItkUcharImgType::Pointer m_NumDirectionsImage; ///< number of peaks per voxel ItkUcharImgType::Pointer m_MaskImage; ///< only voxels inside the binary mask are processed + + Toolkit m_Toolkit; }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkFiniteDiffOdfMaximaExtractionFilter.cpp" #endif #endif //__itkFiniteDiffOdfMaximaExtractionFilter_h_ diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkFslShCoefficientImageConverter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkFslShCoefficientImageConverter.h deleted file mode 100644 index 2d734ee75f..0000000000 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkFslShCoefficientImageConverter.h +++ /dev/null @@ -1,85 +0,0 @@ -/*========================================================================= - - 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. - -=========================================================================*/ -#ifndef __itkFslShCoefficientImageConverter_h_ -#define __itkFslShCoefficientImageConverter_h_ - -#include - -namespace itk{ -/** \class FslShCoefficientImageConverter - Converts FSL reconstructions of diffusionweighted images (4D images containing the sh coefficients) to qball images or 3D sh-coefficient images. -*/ - -template< class PixelType, int ShOrder > -class FslShCoefficientImageConverter : public ProcessObject -{ - -public: - - enum NormalizationMethods { - NO_NORM, - SINGLE_VEC_NORM, - SPACING_COMPENSATION - }; - - typedef FslShCoefficientImageConverter Self; - typedef SmartPointer Pointer; - typedef SmartPointer ConstPointer; - typedef ProcessObject Superclass; - typedef itk::Image< float, 4 > InputImageType; - typedef Image< Vector< PixelType, (ShOrder*ShOrder + ShOrder + 2)/2 + ShOrder >, 3 > CoefficientImageType; - typedef Image< Vector< PixelType, QBALL_ODFSIZE >, 3 > QballImageType; - - /** Method for creation through the object factory. */ - itkNewMacro(Self) - - /** Runtime information support. */ - itkTypeMacro(FslShCoefficientImageConverter, ProcessObject) - - // input - itkSetMacro( InputImage, InputImageType::Pointer) ///< sh coefficient image in FSL file format - - // output - itkGetMacro( CoefficientImage, typename CoefficientImageType::Pointer) ///< mitk style image containing the SH coefficients - itkGetMacro( QballImage, typename QballImageType::Pointer) ///< mitk Q-Ball image generated from the coefficients - - void GenerateData(); - -protected: - FslShCoefficientImageConverter(); - ~FslShCoefficientImageConverter(){} - - void CalcShBasis(); - vnl_matrix_fixed GetSphericalOdfDirections(); - - InputImageType::Pointer m_InputImage; - typename CoefficientImageType::Pointer m_CoefficientImage; ///< mitk style image containing the SH coefficients - typename QballImageType::Pointer m_QballImage; ///< mitk Q-Ball image generated from the coefficients - vnl_matrix m_ShBasis; - -private: - -}; - -} - -#ifndef ITK_MANUAL_INSTANTIATION -#include "itkFslShCoefficientImageConverter.cpp" -#endif - -#endif //__itkFslShCoefficientImageConverter_h_ - diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkMrtrixPeakImageConverter.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkMrtrixPeakImageConverter.cpp index 2c33adc7cb..4580f4f34c 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkMrtrixPeakImageConverter.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkMrtrixPeakImageConverter.cpp @@ -1,157 +1,168 @@ #ifndef __itkMrtrixPeakImageConverter_cpp #define __itkMrtrixPeakImageConverter_cpp #include #include #include #include "itkMrtrixPeakImageConverter.h" #include #include #include #include #include #include #include #include #include namespace itk { template< class PixelType > MrtrixPeakImageConverter< PixelType >::MrtrixPeakImageConverter(): m_NormalizationMethod(NO_NORM) { } template< class PixelType > void MrtrixPeakImageConverter< PixelType > ::GenerateData() { // output vector field vtkSmartPointer m_VtkCellArray = vtkSmartPointer::New(); vtkSmartPointer m_VtkPoints = vtkSmartPointer::New(); Vector spacing4 = m_InputImage->GetSpacing(); Point origin4 = m_InputImage->GetOrigin(); Matrix direction4 = m_InputImage->GetDirection(); ImageRegion<4> imageRegion4 = m_InputImage->GetLargestPossibleRegion(); Vector spacing3; Point origin3; Matrix direction3; ImageRegion<3> imageRegion3; spacing3[0] = spacing4[0]; spacing3[1] = spacing4[1]; spacing3[2] = spacing4[2]; origin3[0] = origin4[0]; origin3[1] = origin4[1]; origin3[2] = origin4[2]; for (int r=0; r<3; r++) for (int c=0; c<3; c++) direction3[r][c] = direction4[r][c]; imageRegion3.SetSize(0, imageRegion4.GetSize()[0]); imageRegion3.SetSize(1, imageRegion4.GetSize()[1]); imageRegion3.SetSize(2, imageRegion4.GetSize()[2]); double minSpacing = spacing3[0]; if (spacing3[1] InputIteratorType; int x = m_InputImage->GetLargestPossibleRegion().GetSize(0); int y = m_InputImage->GetLargestPossibleRegion().GetSize(1); int z = m_InputImage->GetLargestPossibleRegion().GetSize(2); int numDirs = m_InputImage->GetLargestPossibleRegion().GetSize(3)/3; + m_NumDirectionsImage = ItkUcharImgType::New(); + m_NumDirectionsImage->SetSpacing( spacing3 ); + m_NumDirectionsImage->SetOrigin( origin3 ); + m_NumDirectionsImage->SetDirection( direction3 ); + m_NumDirectionsImage->SetRegions( imageRegion3 ); + m_NumDirectionsImage->Allocate(); + m_NumDirectionsImage->FillBuffer(0); + for (int i=0; iSetSpacing( spacing3 ); directionImage->SetOrigin( origin3 ); directionImage->SetDirection( direction3 ); directionImage->SetRegions( imageRegion3 ); directionImage->Allocate(); Vector< PixelType, 3 > nullVec; nullVec.Fill(0.0); directionImage->FillBuffer(nullVec); for (int a=0; a dirVec; dirVec.set_size(4); for (int k=0; k<3; k++) { index.SetElement(3,k+i*3); dirVec[k] = m_InputImage->GetPixel(index); } dirVec[3] = 0; vtkSmartPointer container = vtkSmartPointer::New(); itk::ContinuousIndex center; center[0] = index[0]; center[1] = index[1]; center[2] = index[2]; center[3] = 0; itk::Point worldCenter; m_InputImage->TransformContinuousIndexToPhysicalPoint( center, worldCenter ); switch (m_NormalizationMethod) { case NO_NORM: break; case SINGLE_VEC_NORM: dirVec.normalize(); break; } dirVec.normalize(); dirVec = m_InputImage->GetDirection()*dirVec; itk::Point worldStart; worldStart[0] = worldCenter[0]-dirVec[0]/2 * minSpacing; worldStart[1] = worldCenter[1]-dirVec[1]/2 * minSpacing; worldStart[2] = worldCenter[2]-dirVec[2]/2 * minSpacing; vtkIdType id = m_VtkPoints->InsertNextPoint(worldStart.GetDataPointer()); container->GetPointIds()->InsertNextId(id); itk::Point worldEnd; worldEnd[0] = worldCenter[0]+dirVec[0]/2 * minSpacing; worldEnd[1] = worldCenter[1]+dirVec[1]/2 * minSpacing; worldEnd[2] = worldCenter[2]+dirVec[2]/2 * minSpacing; id = m_VtkPoints->InsertNextPoint(worldEnd.GetDataPointer()); container->GetPointIds()->InsertNextId(id); m_VtkCellArray->InsertNextCell(container); // generate direction image typename ItkDirectionImageType::IndexType index2; index2[0] = a; index2[1] = b; index2[2] = c; Vector< PixelType, 3 > pixel; pixel.SetElement(0, dirVec[0]); pixel.SetElement(1, dirVec[1]); pixel.SetElement(2, dirVec[2]); directionImage->SetPixel(index2, pixel); + + if (dirVec.magnitude()>0.0001) + m_NumDirectionsImage->SetPixel(index2, m_NumDirectionsImage->GetPixel(index2)+1); } m_DirectionImageContainer->InsertElement(m_DirectionImageContainer->Size(), directionImage); } vtkSmartPointer directionsPolyData = vtkSmartPointer::New(); directionsPolyData->SetPoints(m_VtkPoints); directionsPolyData->SetLines(m_VtkCellArray); m_OutputFiberBundle = mitk::FiberBundleX::New(directionsPolyData); } } #endif // __itkMrtrixPeakImageConverter_cpp diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkMrtrixPeakImageConverter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkMrtrixPeakImageConverter.h index c8f18b4aed..01aa33e2ad 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkMrtrixPeakImageConverter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkMrtrixPeakImageConverter.h @@ -1,87 +1,90 @@ /*========================================================================= 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. =========================================================================*/ #ifndef __itkMrtrixPeakImageConverter_h_ #define __itkMrtrixPeakImageConverter_h_ #include #include #include #include namespace itk{ /** \class MrtrixPeakImageConverter Converts a series of 4D images containing directional (3D vector) information into a vector field stored as mitkFiberBundleX. These 4D files are for example generated by the FSL qboot command. */ template< class PixelType > class MrtrixPeakImageConverter : public ProcessObject { public: enum NormalizationMethods { NO_NORM, ///< don't normalize peaks SINGLE_VEC_NORM ///< normalize peaks to length 1 }; typedef MrtrixPeakImageConverter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ProcessObject Superclass; /** Method for creation through the object factory. */ itkNewMacro(Self) /** Runtime information support. */ itkTypeMacro(MrtrixPeakImageConverter, ImageToImageFilter) typedef vnl_vector_fixed< double, 3 > DirectionType; typedef VectorContainer< int, DirectionType > DirectionContainerType; typedef VectorContainer< int, DirectionContainerType::Pointer > ContainerType; typedef Image< float, 4 > InputImageType; typedef Image< Vector< float, 3 >, 3> ItkDirectionImageType; typedef VectorContainer< int, ItkDirectionImageType::Pointer > DirectionImageContainerType; + typedef itk::Image ItkUcharImgType; itkSetMacro( NormalizationMethod, NormalizationMethods) ///< normalization method of ODF peaks itkSetMacro( InputImage, InputImageType::Pointer) ///< MRtrix direction image of type itk::Image< float, 4 > itkGetMacro( OutputFiberBundle, mitk::FiberBundleX::Pointer) ///< vector field (peak sizes rescaled for visualization purposes) itkGetMacro( DirectionImageContainer, DirectionImageContainerType::Pointer) ///< container for output peaks + itkGetMacro( NumDirectionsImage, ItkUcharImgType::Pointer) ///< number of peaks per voxel void GenerateData(); protected: MrtrixPeakImageConverter(); ~MrtrixPeakImageConverter(){} NormalizationMethods m_NormalizationMethod; ///< normalization method of ODF peaks mitk::FiberBundleX::Pointer m_OutputFiberBundle; ///< vector field (peak sizes rescaled for visualization purposes) InputImageType::Pointer m_InputImage; ///< MRtrix direction image of type itk::Image< float, 4 > DirectionImageContainerType::Pointer m_DirectionImageContainer; ///< container for output peaks + ItkUcharImgType::Pointer m_NumDirectionsImage; ///< number of peaks per voxel private: }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkMrtrixPeakImageConverter.cpp" #endif #endif //__itkMrtrixPeakImageConverter_h_ diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkFslShCoefficientImageConverter.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkShCoefficientImageImporter.cpp similarity index 74% rename from Modules/DiffusionImaging/DiffusionCore/Algorithms/itkFslShCoefficientImageConverter.cpp rename to Modules/DiffusionImaging/DiffusionCore/Algorithms/itkShCoefficientImageImporter.cpp index 6582911678..6a21735284 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkFslShCoefficientImageConverter.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkShCoefficientImageImporter.cpp @@ -1,163 +1,179 @@ -#ifndef __itkFslShCoefficientImageConverter_cpp -#define __itkFslShCoefficientImageConverter_cpp +#ifndef __itkShCoefficientImageImporter_cpp +#define __itkShCoefficientImageImporter_cpp #include #include #include -#include "itkFslShCoefficientImageConverter.h" +#include "itkShCoefficientImageImporter.h" #include #include using namespace boost::math; namespace itk { template< class PixelType, int ShOrder > -FslShCoefficientImageConverter< PixelType, ShOrder >::FslShCoefficientImageConverter() +ShCoefficientImageImporter< PixelType, ShOrder >::ShCoefficientImageImporter() + : m_Toolkit(FSL) { m_ShBasis.set_size(QBALL_ODFSIZE, (ShOrder+1)*(ShOrder+2)/2); - CalcShBasis(); } template< class PixelType, int ShOrder > -void FslShCoefficientImageConverter< PixelType, ShOrder > +void ShCoefficientImageImporter< PixelType, ShOrder > ::GenerateData() { + CalcShBasis(); if (m_InputImage.IsNull()) return; Vector spacing4 = m_InputImage->GetSpacing(); Point origin4 = m_InputImage->GetOrigin(); Matrix direction4 = m_InputImage->GetDirection(); ImageRegion<4> imageRegion4 = m_InputImage->GetLargestPossibleRegion(); Vector spacing3; Point origin3; Matrix direction3; ImageRegion<3> imageRegion3; spacing3[0] = spacing4[0]; spacing3[1] = spacing4[1]; spacing3[2] = spacing4[2]; origin3[0] = origin4[0]; origin3[1] = origin4[1]; origin3[2] = origin4[2]; for (int r=0; r<3; r++) for (int c=0; c<3; c++) direction3[r][c] = direction4[r][c]; imageRegion3.SetSize(0, imageRegion4.GetSize()[0]); imageRegion3.SetSize(1, imageRegion4.GetSize()[1]); imageRegion3.SetSize(2, imageRegion4.GetSize()[2]); m_QballImage = QballImageType::New(); m_QballImage->SetSpacing( spacing3 ); m_QballImage->SetOrigin( origin3 ); m_QballImage->SetDirection( direction3 ); m_QballImage->SetRegions( imageRegion3 ); m_QballImage->Allocate(); Vector< PixelType, QBALL_ODFSIZE > nullVec1; nullVec1.Fill(0.0); m_QballImage->FillBuffer(nullVec1); m_CoefficientImage = CoefficientImageType::New(); m_CoefficientImage->SetSpacing( spacing3 ); m_CoefficientImage->SetOrigin( origin3 ); m_CoefficientImage->SetDirection( direction3 ); m_CoefficientImage->SetRegions( imageRegion3 ); m_CoefficientImage->Allocate(); Vector< PixelType, (ShOrder*ShOrder + ShOrder + 2)/2 + ShOrder > nullVec2; nullVec2.Fill(0.0); m_CoefficientImage->FillBuffer(nullVec2); typedef ImageRegionConstIterator< InputImageType > InputIteratorType; int x = imageRegion4.GetSize(0); int y = imageRegion4.GetSize(1); int z = imageRegion4.GetSize(2); int numCoeffs = imageRegion4.GetSize(3); for (int a=0; a coeffs((ShOrder*ShOrder + ShOrder + 2)/2 + ShOrder,1); typename InputImageType::IndexType index; index.SetElement(0,a); index.SetElement(1,b); index.SetElement(2,c); typename CoefficientImageType::PixelType pix; for (int d=0; dGetPixel(index); coeffs[d][0] = pix[d]; } typename CoefficientImageType::IndexType index2; index2.SetElement(0,a); index2.SetElement(1,b); index2.SetElement(2,c); m_CoefficientImage->SetPixel(index2, pix); - vnl_matrix odf = m_ShBasis*coeffs; typename QballImageType::PixelType pix2; + vnl_matrix odf = m_ShBasis*coeffs; for (int d=0; dSetPixel(index2,pix2); } - } // generate spherical harmonic values of the desired order for each input direction template< class PixelType, int ShOrder > -void FslShCoefficientImageConverter< PixelType, ShOrder > +void ShCoefficientImageImporter< PixelType, ShOrder > ::CalcShBasis() { vnl_matrix_fixed sphCoords = GetSphericalOdfDirections(); int j, m; double mag, plm; for (int p=0; p(l,abs(m),cos(sphCoords(0,p))); - mag = sqrt((double)(2*l+1)/(4.0*M_PI)*factorial(l-abs(m))/factorial(l+abs(m)))*plm; - - if (m<0) - m_ShBasis(p,j) = sqrt(2.0)*mag*cos(-m*sphCoords(1,p)); - else if (m==0) - m_ShBasis(p,j) = mag; - else - m_ShBasis(p,j) = pow(-1.0, m)*sqrt(2.0)*mag*sin(m*sphCoords(1,p)); + switch (m_Toolkit) + { + case FSL: + plm = legendre_p(l,abs(m),cos(sphCoords(0,p))); + mag = sqrt((double)(2*l+1)/(4.0*M_PI)*factorial(l-abs(m))/factorial(l+abs(m)))*plm; + if (m<0) + m_ShBasis(p,j) = sqrt(2.0)*mag*cos(-m*sphCoords(1,p)); + else if (m==0) + m_ShBasis(p,j) = mag; + else + m_ShBasis(p,j) = pow(-1.0, m)*sqrt(2.0)*mag*sin(m*sphCoords(1,p)); + break; + case MRTRIX: + plm = legendre_p(l,abs(m),-cos(sphCoords(0,p))); + mag = sqrt((double)(2*l+1)/(4.0*M_PI)*factorial(l-abs(m))/factorial(l+abs(m)))*plm; + if (m>0) + m_ShBasis(p,j) = mag*cos(m*sphCoords(1,p)); + else if (m==0) + m_ShBasis(p,j) = mag; + else + m_ShBasis(p,j) = mag*sin(-m*sphCoords(1,p)); + break; + } + j++; } } } // convert cartesian to spherical coordinates template< class PixelType, int ShOrder > -vnl_matrix_fixed FslShCoefficientImageConverter< PixelType, ShOrder > +vnl_matrix_fixed ShCoefficientImageImporter< PixelType, ShOrder > ::GetSphericalOdfDirections() { itk::OrientationDistributionFunction< PixelType, QBALL_ODFSIZE > odf; vnl_matrix_fixed* dir = odf.GetDirections(); vnl_matrix_fixed sphCoords; for (int i=0; iget_column(i).magnitude(); if( magget(2,i)/mag); // theta sphCoords(1,i) = atan2(dir->get(1,i), dir->get(0,i)); // phi } } return sphCoords; } } -#endif // __itkFslShCoefficientImageConverter_cpp +#endif // __itkShCoefficientImageImporter_cpp diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkShCoefficientImageImporter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkShCoefficientImageImporter.h new file mode 100644 index 0000000000..74673b4e87 --- /dev/null +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkShCoefficientImageImporter.h @@ -0,0 +1,94 @@ +/*========================================================================= + + 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. + +=========================================================================*/ +#ifndef __itkShCoefficientImageImporter_h_ +#define __itkShCoefficientImageImporter_h_ + +#include + +namespace itk{ +/** \class ShCoefficientImageImporter + Converts FSL reconstructions of diffusionweighted images (4D images containing the sh coefficients) to qball images or 3D sh-coefficient images. +*/ + +template< class PixelType, int ShOrder > +class ShCoefficientImageImporter : public ProcessObject +{ + +public: + + enum NormalizationMethods { + NO_NORM, + SINGLE_VEC_NORM, + SPACING_COMPENSATION + }; + + enum Toolkit { ///< SH coefficient convention (depends on toolkit) + FSL, + MRTRIX + }; + + typedef ShCoefficientImageImporter Self; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + typedef ProcessObject Superclass; + typedef itk::Image< float, 4 > InputImageType; + typedef Image< Vector< PixelType, (ShOrder*ShOrder + ShOrder + 2)/2 + ShOrder >, 3 > CoefficientImageType; + typedef Image< Vector< PixelType, QBALL_ODFSIZE >, 3 > QballImageType; + + /** Method for creation through the object factory. */ + itkNewMacro(Self) + + /** Runtime information support. */ + itkTypeMacro(ShCoefficientImageImporter, ProcessObject) + + // input + itkSetMacro( InputImage, InputImageType::Pointer) ///< sh coefficient image in FSL file format + + // output + itkGetMacro( CoefficientImage, typename CoefficientImageType::Pointer) ///< mitk style image containing the SH coefficients + itkGetMacro( QballImage, typename QballImageType::Pointer) ///< mitk Q-Ball image generated from the coefficients + + itkSetMacro( Toolkit, Toolkit) ///< define SH coefficient convention (depends on toolkit) + itkGetMacro( Toolkit, Toolkit) ///< SH coefficient convention (depends on toolkit) + + void GenerateData(); + +protected: + ShCoefficientImageImporter(); + ~ShCoefficientImageImporter(){} + + void CalcShBasis(); + vnl_matrix_fixed GetSphericalOdfDirections(); + + InputImageType::Pointer m_InputImage; + typename CoefficientImageType::Pointer m_CoefficientImage; ///< mitk style image containing the SH coefficients + typename QballImageType::Pointer m_QballImage; ///< mitk Q-Ball image generated from the coefficients + vnl_matrix m_ShBasis; + Toolkit m_Toolkit; + +private: + +}; + +} + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkShCoefficientImageImporter.cpp" +#endif + +#endif //__itkShCoefficientImageImporter_h_ + diff --git a/Modules/DiffusionImaging/DiffusionCore/files.cmake b/Modules/DiffusionImaging/DiffusionCore/files.cmake index 75d68664b5..5c047608b9 100644 --- a/Modules/DiffusionImaging/DiffusionCore/files.cmake +++ b/Modules/DiffusionImaging/DiffusionCore/files.cmake @@ -1,116 +1,116 @@ set(CPP_FILES # DicomImport DicomImport/mitkDicomDiffusionImageReader.cpp DicomImport/mitkGroupDiffusionHeadersFilter.cpp DicomImport/mitkDicomDiffusionImageHeaderReader.cpp DicomImport/mitkGEDicomDiffusionImageHeaderReader.cpp DicomImport/mitkPhilipsDicomDiffusionImageHeaderReader.cpp DicomImport/mitkSiemensDicomDiffusionImageHeaderReader.cpp DicomImport/mitkSiemensMosaicDicomDiffusionImageHeaderReader.cpp # DataStructures IODataStructures/mitkDiffusionCoreObjectFactory.cpp # DataStructures -> DWI IODataStructures/DiffusionWeightedImages/mitkDiffusionImageHeaderInformation.cpp IODataStructures/DiffusionWeightedImages/mitkDiffusionImageSource.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageReader.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageWriter.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageIOFactory.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageWriterFactory.cpp IODataStructures/DiffusionWeightedImages/mitkDiffusionImageSerializer.cpp # DataStructures -> QBall IODataStructures/QBallImages/mitkQBallImageSource.cpp IODataStructures/QBallImages/mitkNrrdQBallImageReader.cpp IODataStructures/QBallImages/mitkNrrdQBallImageWriter.cpp IODataStructures/QBallImages/mitkNrrdQBallImageIOFactory.cpp IODataStructures/QBallImages/mitkNrrdQBallImageWriterFactory.cpp IODataStructures/QBallImages/mitkQBallImage.cpp IODataStructures/QBallImages/mitkQBallImageSerializer.cpp # DataStructures -> Tensor IODataStructures/TensorImages/mitkTensorImageSource.cpp IODataStructures/TensorImages/mitkNrrdTensorImageReader.cpp IODataStructures/TensorImages/mitkNrrdTensorImageWriter.cpp IODataStructures/TensorImages/mitkNrrdTensorImageIOFactory.cpp IODataStructures/TensorImages/mitkNrrdTensorImageWriterFactory.cpp IODataStructures/TensorImages/mitkTensorImage.cpp IODataStructures/TensorImages/mitkTensorImageSerializer.cpp # Rendering Rendering/vtkMaskedProgrammableGlyphFilter.cpp Rendering/mitkCompositeMapper.cpp Rendering/mitkVectorImageVtkGlyphMapper3D.cpp Rendering/vtkOdfSource.cxx Rendering/vtkThickPlane.cxx Rendering/mitkOdfNormalizationMethodProperty.cpp Rendering/mitkOdfScaleByProperty.cpp Rendering/mitkPlanarFigureMapper3D.cpp # Algorithms Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.cpp Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.cpp # Function Collection mitkDiffusionFunctionCollection.cpp ) set(H_FILES # function Collection mitkDiffusionFunctionCollection.h # Rendering Rendering/mitkDiffusionImageMapper.h Rendering/mitkOdfVtkMapper2D.h Rendering/mitkPlanarFigureMapper3D.h # Reconstruction Algorithms/Reconstruction/itkDiffusionQballReconstructionImageFilter.h Algorithms/Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.h Algorithms/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.h Algorithms/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h Algorithms/Reconstruction/itkPointShell.h Algorithms/Reconstruction/itkOrientationDistributionFunction.h Algorithms/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h # IO Datastructures IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h # Algorithms Algorithms/itkDiffusionQballGeneralizedFaImageFilter.h Algorithms/itkDiffusionQballPrepareVisualizationImageFilter.h Algorithms/itkTensorDerivedMeasurementsFilter.h Algorithms/itkBrainMaskExtractionImageFilter.h Algorithms/itkB0ImageExtractionImageFilter.h Algorithms/itkB0ImageExtractionToSeparateImageFilter.h Algorithms/itkTensorImageToDiffusionImageFilter.h Algorithms/itkTensorToL2NormImageFilter.h Algorithms/itkGaussianInterpolateImageFunction.h Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.h Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.h Algorithms/itkDiffusionTensorPrincipalDirectionImageFilter.h Algorithms/itkCartesianToPolarVectorImageFilter.h Algorithms/itkPolarToCartesianVectorImageFilter.h Algorithms/itkDistanceMapFilter.h Algorithms/itkProjectionFilter.h Algorithms/itkResidualImageFilter.h Algorithms/itkExtractChannelFromRgbaImageFilter.h Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.h Algorithms/itkMergeDiffusionImagesFilter.h Algorithms/itkDwiPhantomGenerationFilter.h Algorithms/itkFiniteDiffOdfMaximaExtractionFilter.h Algorithms/itkMrtrixPeakImageConverter.h Algorithms/itkFslPeakImageConverter.h - Algorithms/itkFslShCoefficientImageConverter.h + Algorithms/itkShCoefficientImageImporter.h Algorithms/itkOdfMaximaExtractionFilter.h Algorithms/itkResampleDwiImageFilter.h ) set( TOOL_FILES ) diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkOdfMaximaExtractionView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkOdfMaximaExtractionView.cpp index daaff92d8c..dbf50eec80 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkOdfMaximaExtractionView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkOdfMaximaExtractionView.cpp @@ -1,723 +1,751 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ //misc #define _USE_MATH_DEFINES #include #include // Blueberry #include #include // Qmitk #include "QmitkOdfMaximaExtractionView.h" // MITK #include #include #include #include #include #include // ITK #include #include #include #include #include -#include +#include #include const std::string QmitkOdfMaximaExtractionView::VIEW_ID = "org.mitk.views.odfmaximaextractionview"; using namespace mitk; QmitkOdfMaximaExtractionView::QmitkOdfMaximaExtractionView() : QmitkFunctionality() , m_Controls( 0 ) , m_MultiWidget( NULL ) { } // Destructor QmitkOdfMaximaExtractionView::~QmitkOdfMaximaExtractionView() { } void QmitkOdfMaximaExtractionView::CreateQtPartControl( QWidget *parent ) { // build up qt view, unless already done if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkOdfMaximaExtractionViewControls; m_Controls->setupUi( parent ); connect((QObject*) m_Controls->m_StartTensor, SIGNAL(clicked()), (QObject*) this, SLOT(StartTensor())); connect((QObject*) m_Controls->m_StartFiniteDiff, SIGNAL(clicked()), (QObject*) this, SLOT(StartFiniteDiff())); connect((QObject*) m_Controls->m_GenerateImageButton, SIGNAL(clicked()), (QObject*) this, SLOT(GenerateImage())); - connect((QObject*) m_Controls->m_ConvertFromFsl, SIGNAL(clicked()), (QObject*) this, SLOT(ConvertPeaksFromFsl())); - connect((QObject*) m_Controls->m_ConvertShFromFsl, SIGNAL(clicked()), (QObject*) this, SLOT(ConvertShCoeffsFromFsl())); - connect((QObject*) m_Controls->m_MrtrixPeaks, SIGNAL(clicked()), (QObject*) this, SLOT(ConvertPeaksFromMrtrix())); + connect((QObject*) m_Controls->m_ImportPeaks, SIGNAL(clicked()), (QObject*) this, SLOT(ConvertPeaks())); + connect((QObject*) m_Controls->m_ImportShCoeffs, SIGNAL(clicked()), (QObject*) this, SLOT(ConvertShCoeffs())); } } void QmitkOdfMaximaExtractionView::UpdateGui() { m_Controls->m_GenerateImageButton->setEnabled(false); m_Controls->m_StartFiniteDiff->setEnabled(false); m_Controls->m_StartTensor->setEnabled(false); m_Controls->m_CoeffImageFrame->setEnabled(false); if (!m_ImageNodes.empty() || !m_TensorImageNodes.empty()) { m_Controls->m_InputData->setTitle("Input Data"); if (!m_TensorImageNodes.empty()) { m_Controls->m_DwiFibLabel->setText(m_TensorImageNodes.front()->GetName().c_str()); m_Controls->m_StartTensor->setEnabled(true); } else { m_Controls->m_DwiFibLabel->setText(m_ImageNodes.front()->GetName().c_str()); m_Controls->m_StartFiniteDiff->setEnabled(true); m_Controls->m_GenerateImageButton->setEnabled(true); m_Controls->m_CoeffImageFrame->setEnabled(true); m_Controls->m_ShOrderBox->setEnabled(true); m_Controls->m_MaxNumPeaksBox->setEnabled(true); m_Controls->m_PeakThresholdBox->setEnabled(true); m_Controls->m_AbsoluteThresholdBox->setEnabled(true); } } else m_Controls->m_DwiFibLabel->setText("mandatory"); if (m_ImageNodes.empty()) { - m_Controls->m_ConvertFromFsl->setEnabled(false); - m_Controls->m_ConvertShFromFsl->setEnabled(false); + m_Controls->m_ImportPeaks->setEnabled(false); + m_Controls->m_ImportShCoeffs->setEnabled(false); } else { - m_Controls->m_ConvertFromFsl->setEnabled(true); - m_Controls->m_ConvertShFromFsl->setEnabled(true); + m_Controls->m_ImportPeaks->setEnabled(true); + m_Controls->m_ImportShCoeffs->setEnabled(true); } if (!m_BinaryImageNodes.empty()) { m_Controls->m_MaskLabel->setText(m_BinaryImageNodes.front()->GetName().c_str()); } else { m_Controls->m_MaskLabel->setText("optional"); } } template -void QmitkOdfMaximaExtractionView::TemplatedConvertShCoeffsFromFsl(mitk::Image* mitkImg) +void QmitkOdfMaximaExtractionView::TemplatedConvertShCoeffs(mitk::Image* mitkImg) { - typedef itk::FslShCoefficientImageConverter< float, shOrder > FilterType; + typedef itk::ShCoefficientImageImporter< float, shOrder > FilterType; typedef mitk::ImageToItk< itk::Image< float, 4 > > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(mitkImg); caster->Update(); typename FilterType::Pointer filter = FilterType::New(); + + switch (m_Controls->m_ToolkitBox->currentIndex()) + { + case 0: + filter->SetToolkit(FilterType::FSL); + break; + case 1: + filter->SetToolkit(FilterType::MRTRIX); + break; + default: + filter->SetToolkit(FilterType::FSL); + } + filter->SetInputImage(caster->GetOutput()); filter->GenerateData(); typename FilterType::QballImageType::Pointer itkQbi = filter->GetQballImage(); typename FilterType::CoefficientImageType::Pointer itkCi = filter->GetCoefficientImage(); { mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk( itkCi.GetPointer() ); img->SetVolume( itkCi->GetBufferPointer() ); DataNode::Pointer node = DataNode::New(); node->SetData(img); - node->SetName("FSL_ShCoefficientImage"); + node->SetName("_ShCoefficientImage"); GetDataStorage()->Add(node); } { mitk::QBallImage::Pointer img = mitk::QBallImage::New(); img->InitializeByItk( itkQbi.GetPointer() ); img->SetVolume( itkQbi->GetBufferPointer() ); DataNode::Pointer node = DataNode::New(); node->SetData(img); - node->SetName("FSL_QballImage"); + node->SetName("_QballImage"); GetDataStorage()->Add(node); } } -void QmitkOdfMaximaExtractionView::ConvertShCoeffsFromFsl() +void QmitkOdfMaximaExtractionView::ConvertShCoeffs() { if (m_ImageNodes.empty()) return; mitk::Image::Pointer mitkImg = dynamic_cast(m_ImageNodes.at(0)->GetData()); if (mitkImg->GetDimension()!=4) { MITK_INFO << "wrong image type (need 4 dimensions)"; return; } int nrCoeffs = mitkImg->GetLargestPossibleRegion().GetSize()[3]; // solve bx² + cx + d = 0 = shOrder² + 2*shOrder + 2-2*neededCoeffs; int c=3, d=2-2*nrCoeffs; double D = c*c-4*d; int shOrder; if (D>0) { shOrder = (-c+sqrt(D))/2.0; if (shOrder<0) shOrder = (-c-sqrt(D))/2.0; } else if (D==0) shOrder = -c/2.0; MITK_INFO << "using SH-order " << shOrder; switch (shOrder) { case 4: - TemplatedConvertShCoeffsFromFsl<4>(mitkImg); + TemplatedConvertShCoeffs<4>(mitkImg); break; case 6: - TemplatedConvertShCoeffsFromFsl<6>(mitkImg); + TemplatedConvertShCoeffs<6>(mitkImg); break; case 8: - TemplatedConvertShCoeffsFromFsl<8>(mitkImg); + TemplatedConvertShCoeffs<8>(mitkImg); break; case 10: - TemplatedConvertShCoeffsFromFsl<10>(mitkImg); + TemplatedConvertShCoeffs<10>(mitkImg); break; case 12: - TemplatedConvertShCoeffsFromFsl<12>(mitkImg); + TemplatedConvertShCoeffs<12>(mitkImg); break; default: MITK_INFO << "SH-order " << shOrder << " not supported"; } } -void QmitkOdfMaximaExtractionView::ConvertPeaksFromFsl() +void QmitkOdfMaximaExtractionView::ConvertPeaks() { if (m_ImageNodes.empty()) return; - typedef itk::Image< float, 4 > ItkImageType; - typedef itk::FslPeakImageConverter< float > FilterType; - FilterType::Pointer filter = FilterType::New(); - FilterType::InputType::Pointer inputVec = FilterType::InputType::New(); - mitk::Geometry3D::Pointer geom; - - for (int i=0; im_ToolkitBox->currentIndex()) { - mitk::Image::Pointer mitkImg = dynamic_cast(m_ImageNodes.at(i)->GetData()); - geom = mitkImg->GetGeometry(); - typedef mitk::ImageToItk< FilterType::InputImageType > CasterType; - CasterType::Pointer caster = CasterType::New(); - caster->SetInput(mitkImg); - caster->Update(); - FilterType::InputImageType::Pointer itkImg = caster->GetOutput(); - inputVec->InsertElement(inputVec->Size(), itkImg); - } + case 0: + { + typedef itk::Image< float, 4 > ItkImageType; + typedef itk::FslPeakImageConverter< float > FilterType; + FilterType::Pointer filter = FilterType::New(); + FilterType::InputType::Pointer inputVec = FilterType::InputType::New(); + mitk::Geometry3D::Pointer geom; - filter->SetInputImages(inputVec); - filter->GenerateData(); + for (int i=0; i(m_ImageNodes.at(i)->GetData()); + geom = mitkImg->GetGeometry(); + typedef mitk::ImageToItk< FilterType::InputImageType > CasterType; + CasterType::Pointer caster = CasterType::New(); + caster->SetInput(mitkImg); + caster->Update(); + FilterType::InputImageType::Pointer itkImg = caster->GetOutput(); + inputVec->InsertElement(inputVec->Size(), itkImg); + } - mitk::Vector3D outImageSpacing = geom->GetSpacing(); - float maxSpacing = 1; - if(outImageSpacing[0]>outImageSpacing[1] && outImageSpacing[0]>outImageSpacing[2]) - maxSpacing = outImageSpacing[0]; - else if (outImageSpacing[1] > outImageSpacing[2]) - maxSpacing = outImageSpacing[1]; - else - maxSpacing = outImageSpacing[2]; - - mitk::FiberBundleX::Pointer directions = filter->GetOutputFiberBundle(); - directions->SetGeometry(geom); - DataNode::Pointer node = DataNode::New(); - node->SetData(directions); - node->SetName("FSL_VectorField"); - node->SetProperty("Fiber2DSliceThickness", mitk::FloatProperty::New(maxSpacing)); - node->SetProperty("Fiber2DfadeEFX", mitk::BoolProperty::New(false)); - GetDataStorage()->Add(node); - - typedef FilterType::DirectionImageContainerType DirectionImageContainerType; - DirectionImageContainerType::Pointer container = filter->GetDirectionImageContainer(); - for (int i=0; iSize(); i++) - { - ItkDirectionImage3DType::Pointer itkImg = container->GetElement(i); - mitk::Image::Pointer img = mitk::Image::New(); - img->InitializeByItk( itkImg.GetPointer() ); - img->SetVolume( itkImg->GetBufferPointer() ); + filter->SetInputImages(inputVec); + filter->GenerateData(); + + mitk::Vector3D outImageSpacing = geom->GetSpacing(); + float maxSpacing = 1; + if(outImageSpacing[0]>outImageSpacing[1] && outImageSpacing[0]>outImageSpacing[2]) + maxSpacing = outImageSpacing[0]; + else if (outImageSpacing[1] > outImageSpacing[2]) + maxSpacing = outImageSpacing[1]; + else + maxSpacing = outImageSpacing[2]; + + mitk::FiberBundleX::Pointer directions = filter->GetOutputFiberBundle(); + directions->SetGeometry(geom); DataNode::Pointer node = DataNode::New(); - node->SetData(img); - QString name(m_ImageNodes.at(i)->GetName().c_str()); - name += "_Direction"; - name += QString::number(i+1); - node->SetName(name.toStdString().c_str()); + node->SetData(directions); + node->SetName("_VectorField"); + node->SetProperty("Fiber2DSliceThickness", mitk::FloatProperty::New(maxSpacing)); + node->SetProperty("Fiber2DfadeEFX", mitk::BoolProperty::New(false)); GetDataStorage()->Add(node); - } -} -void QmitkOdfMaximaExtractionView::ConvertPeaksFromMrtrix() -{ - if (m_ImageNodes.empty()) - return; + typedef FilterType::DirectionImageContainerType DirectionImageContainerType; + DirectionImageContainerType::Pointer container = filter->GetDirectionImageContainer(); + for (int i=0; iSize(); i++) + { + ItkDirectionImage3DType::Pointer itkImg = container->GetElement(i); + mitk::Image::Pointer img = mitk::Image::New(); + img->InitializeByItk( itkImg.GetPointer() ); + img->SetVolume( itkImg->GetBufferPointer() ); + DataNode::Pointer node = DataNode::New(); + node->SetData(img); + QString name(m_ImageNodes.at(i)->GetName().c_str()); + name += "_Direction"; + name += QString::number(i+1); + node->SetName(name.toStdString().c_str()); + GetDataStorage()->Add(node); + } + break; + } + case 1: + { + typedef itk::Image< float, 4 > ItkImageType; + typedef itk::MrtrixPeakImageConverter< float > FilterType; + FilterType::Pointer filter = FilterType::New(); - typedef itk::Image< float, 4 > ItkImageType; - typedef itk::MrtrixPeakImageConverter< float > FilterType; - FilterType::Pointer filter = FilterType::New(); + // cast to itk + mitk::Image::Pointer mitkImg = dynamic_cast(m_ImageNodes.at(0)->GetData()); + mitk::Geometry3D::Pointer geom = mitkImg->GetGeometry(); + typedef mitk::ImageToItk< FilterType::InputImageType > CasterType; + CasterType::Pointer caster = CasterType::New(); + caster->SetInput(mitkImg); + caster->Update(); + FilterType::InputImageType::Pointer itkImg = caster->GetOutput(); - // cast to itk - mitk::Image::Pointer mitkImg = dynamic_cast(m_ImageNodes.at(0)->GetData()); - mitk::Geometry3D::Pointer geom = mitkImg->GetGeometry(); - typedef mitk::ImageToItk< FilterType::InputImageType > CasterType; - CasterType::Pointer caster = CasterType::New(); - caster->SetInput(mitkImg); - caster->Update(); - FilterType::InputImageType::Pointer itkImg = caster->GetOutput(); + filter->SetInputImage(itkImg); + filter->GenerateData(); - filter->SetInputImage(itkImg); - filter->GenerateData(); + mitk::Vector3D outImageSpacing = geom->GetSpacing(); + float maxSpacing = 1; + if(outImageSpacing[0]>outImageSpacing[1] && outImageSpacing[0]>outImageSpacing[2]) + maxSpacing = outImageSpacing[0]; + else if (outImageSpacing[1] > outImageSpacing[2]) + maxSpacing = outImageSpacing[1]; + else + maxSpacing = outImageSpacing[2]; - mitk::Vector3D outImageSpacing = geom->GetSpacing(); - float maxSpacing = 1; - if(outImageSpacing[0]>outImageSpacing[1] && outImageSpacing[0]>outImageSpacing[2]) - maxSpacing = outImageSpacing[0]; - else if (outImageSpacing[1] > outImageSpacing[2]) - maxSpacing = outImageSpacing[1]; - else - maxSpacing = outImageSpacing[2]; - - mitk::FiberBundleX::Pointer directions = filter->GetOutputFiberBundle(); - directions->SetGeometry(geom); - DataNode::Pointer node = DataNode::New(); - node->SetData(directions); - QString name(m_ImageNodes.at(0)->GetName().c_str()); - name += "_VectorField"; - node->SetName(name.toStdString().c_str()); - node->SetProperty("Fiber2DSliceThickness", mitk::FloatProperty::New(maxSpacing)); - node->SetProperty("Fiber2DfadeEFX", mitk::BoolProperty::New(false)); - GetDataStorage()->Add(node); - - typedef FilterType::DirectionImageContainerType DirectionImageContainerType; - DirectionImageContainerType::Pointer container = filter->GetDirectionImageContainer(); - for (int i=0; iSize(); i++) - { - ItkDirectionImage3DType::Pointer itkImg = container->GetElement(i); - mitk::Image::Pointer img = mitk::Image::New(); - img->InitializeByItk( itkImg.GetPointer() ); - img->SetVolume( itkImg->GetBufferPointer() ); + mitk::FiberBundleX::Pointer directions = filter->GetOutputFiberBundle(); + directions->SetGeometry(geom); DataNode::Pointer node = DataNode::New(); - node->SetData(img); + node->SetData(directions); QString name(m_ImageNodes.at(0)->GetName().c_str()); - name += "_Direction"; - name += QString::number(i+1); + name += "_VectorField"; node->SetName(name.toStdString().c_str()); + node->SetProperty("Fiber2DSliceThickness", mitk::FloatProperty::New(maxSpacing)); + node->SetProperty("Fiber2DfadeEFX", mitk::BoolProperty::New(false)); GetDataStorage()->Add(node); + + typedef FilterType::DirectionImageContainerType DirectionImageContainerType; + DirectionImageContainerType::Pointer container = filter->GetDirectionImageContainer(); + for (int i=0; iSize(); i++) + { + ItkDirectionImage3DType::Pointer itkImg = container->GetElement(i); + mitk::Image::Pointer img = mitk::Image::New(); + img->InitializeByItk( itkImg.GetPointer() ); + img->SetVolume( itkImg->GetBufferPointer() ); + DataNode::Pointer node = DataNode::New(); + node->SetData(img); + QString name(m_ImageNodes.at(0)->GetName().c_str()); + name += "_Direction"; + name += QString::number(i+1); + node->SetName(name.toStdString().c_str()); + GetDataStorage()->Add(node); + } + break; + } } } void QmitkOdfMaximaExtractionView::GenerateImage() { if (!m_ImageNodes.empty()) GenerateDataFromDwi(); } void QmitkOdfMaximaExtractionView::StartTensor() { if (m_TensorImageNodes.empty()) return; typedef itk::DiffusionTensorPrincipalDirectionImageFilter< float, float > MaximaExtractionFilterType; MaximaExtractionFilterType::Pointer filter = MaximaExtractionFilterType::New(); mitk::Geometry3D::Pointer geometry; try{ TensorImage::Pointer img = dynamic_cast(m_TensorImageNodes.at(0)->GetData()); ItkTensorImage::Pointer itkImage = ItkTensorImage::New(); CastToItkImage(img, itkImage); filter->SetInput(itkImage); geometry = img->GetGeometry(); } catch(itk::ExceptionObject &e) { MITK_INFO << "wrong image type: " << e.what(); throw e; } if (!m_BinaryImageNodes.empty()) { ItkUcharImgType::Pointer itkMaskImage = ItkUcharImgType::New(); Image::Pointer mitkMaskImg = dynamic_cast(m_BinaryImageNodes.at(0)->GetData()); CastToItkImage(mitkMaskImg, itkMaskImage); filter->SetMaskImage(itkMaskImage); } if (m_Controls->m_NormalizationBox->currentIndex()==0) filter->SetNormalizeVectors(false); filter->Update(); if (m_Controls->m_OutputDirectionImagesBox->isChecked()) { MaximaExtractionFilterType::OutputImageType::Pointer itkImg = filter->GetOutput(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk( itkImg.GetPointer() ); img->SetVolume( itkImg->GetBufferPointer() ); DataNode::Pointer node = DataNode::New(); node->SetData(img); QString name(m_TensorImageNodes.at(0)->GetName().c_str()); name += "_PrincipalDirection"; node->SetName(name.toStdString().c_str()); GetDataStorage()->Add(node); } if (m_Controls->m_OutputNumDirectionsBox->isChecked()) { ItkUcharImgType::Pointer numDirImage = filter->GetNumDirectionsImage(); mitk::Image::Pointer image2 = mitk::Image::New(); image2->InitializeByItk( numDirImage.GetPointer() ); image2->SetVolume( numDirImage->GetBufferPointer() ); DataNode::Pointer node2 = DataNode::New(); node2->SetData(image2); QString name(m_TensorImageNodes.at(0)->GetName().c_str()); name += "_NumDirections"; node2->SetName(name.toStdString().c_str()); GetDataStorage()->Add(node2); } if (m_Controls->m_OutputVectorFieldBox->isChecked()) { mitk::Vector3D outImageSpacing = geometry->GetSpacing(); float minSpacing = 1; if(outImageSpacing[0]GetOutputFiberBundle(); directions->SetGeometry(geometry); DataNode::Pointer node = DataNode::New(); node->SetData(directions); QString name(m_TensorImageNodes.at(0)->GetName().c_str()); name += "_VectorField"; node->SetName(name.toStdString().c_str()); node->SetProperty("Fiber2DSliceThickness", mitk::FloatProperty::New(minSpacing)); node->SetProperty("Fiber2DfadeEFX", mitk::BoolProperty::New(false)); GetDataStorage()->Add(node); } } template void QmitkOdfMaximaExtractionView::StartMaximaExtraction() { typedef itk::FiniteDiffOdfMaximaExtractionFilter< float, shOrder, 20242 > MaximaExtractionFilterType; typename MaximaExtractionFilterType::Pointer filter = MaximaExtractionFilterType::New(); + switch (m_Controls->m_ToolkitBox->currentIndex()) + { + case 0: + filter->SetToolkit(MaximaExtractionFilterType::FSL); + break; + case 1: + filter->SetToolkit(MaximaExtractionFilterType::MRTRIX); + break; + default: + filter->SetToolkit(MaximaExtractionFilterType::FSL); + } + mitk::Geometry3D::Pointer geometry; try{ Image::Pointer img = dynamic_cast(m_ImageNodes.at(0)->GetData()); typedef ImageToItk< typename MaximaExtractionFilterType::CoefficientImageType > CasterType; typename CasterType::Pointer caster = CasterType::New(); caster->SetInput(img); caster->Update(); filter->SetInput(caster->GetOutput()); geometry = img->GetGeometry(); } catch(itk::ExceptionObject &e) { MITK_INFO << "wrong image type: " << e.what(); throw; } filter->SetAngularThreshold(cos((float)m_Controls->m_AngularThreshold->value()*M_PI/180)); filter->SetClusteringThreshold(cos((float)m_Controls->m_ClusteringAngleBox->value()*M_PI/180)); filter->SetMaxNumPeaks(m_Controls->m_MaxNumPeaksBox->value()); filter->SetPeakThreshold(m_Controls->m_PeakThresholdBox->value()); filter->SetAbsolutePeakThreshold(m_Controls->m_AbsoluteThresholdBox->value()); if (!m_BinaryImageNodes.empty()) { ItkUcharImgType::Pointer itkMaskImage = ItkUcharImgType::New(); Image::Pointer mitkMaskImg = dynamic_cast(m_BinaryImageNodes.at(0)->GetData()); CastToItkImage(mitkMaskImg, itkMaskImage); filter->SetMaskImage(itkMaskImage); } switch (m_Controls->m_NormalizationBox->currentIndex()) { case 0: filter->SetNormalizationMethod(MaximaExtractionFilterType::NO_NORM); break; case 1: filter->SetNormalizationMethod(MaximaExtractionFilterType::MAX_VEC_NORM); break; case 2: filter->SetNormalizationMethod(MaximaExtractionFilterType::SINGLE_VEC_NORM); break; } filter->Update(); if (m_Controls->m_OutputDirectionImagesBox->isChecked()) { typedef typename MaximaExtractionFilterType::ItkDirectionImageContainer ItkDirectionImageContainer; typename ItkDirectionImageContainer::Pointer container = filter->GetDirectionImageContainer(); for (int i=0; iSize(); i++) { typename MaximaExtractionFilterType::ItkDirectionImage::Pointer itkImg = container->GetElement(i); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk( itkImg.GetPointer() ); img->SetVolume( itkImg->GetBufferPointer() ); DataNode::Pointer node = DataNode::New(); node->SetData(img); QString name(m_ImageNodes.at(0)->GetName().c_str()); name += "_Direction"; name += QString::number(i+1); node->SetName(name.toStdString().c_str()); GetDataStorage()->Add(node); } } if (m_Controls->m_OutputNumDirectionsBox->isChecked()) { ItkUcharImgType::Pointer numDirImage = filter->GetNumDirectionsImage(); mitk::Image::Pointer image2 = mitk::Image::New(); image2->InitializeByItk( numDirImage.GetPointer() ); image2->SetVolume( numDirImage->GetBufferPointer() ); DataNode::Pointer node2 = DataNode::New(); node2->SetData(image2); QString name(m_ImageNodes.at(0)->GetName().c_str()); name += "_NumDirections"; node2->SetName(name.toStdString().c_str()); GetDataStorage()->Add(node2); } if (m_Controls->m_OutputVectorFieldBox->isChecked()) { mitk::Vector3D outImageSpacing = geometry->GetSpacing(); float minSpacing = 1; if(outImageSpacing[0]GetOutputFiberBundle(); directions->SetGeometry(geometry); DataNode::Pointer node = DataNode::New(); node->SetData(directions); QString name(m_ImageNodes.at(0)->GetName().c_str()); name += "_VectorField"; node->SetName(name.toStdString().c_str()); node->SetProperty("Fiber2DSliceThickness", mitk::FloatProperty::New(minSpacing)); node->SetProperty("Fiber2DfadeEFX", mitk::BoolProperty::New(false)); GetDataStorage()->Add(node); } } void QmitkOdfMaximaExtractionView::StartFiniteDiff() { if (m_ImageNodes.empty()) return; switch (m_Controls->m_ShOrderBox->currentIndex()) { case 0: StartMaximaExtraction<2>(); break; case 1: StartMaximaExtraction<4>(); break; case 2: StartMaximaExtraction<6>(); break; case 3: StartMaximaExtraction<8>(); break; case 4: StartMaximaExtraction<10>(); break; case 5: StartMaximaExtraction<12>(); break; } } void QmitkOdfMaximaExtractionView::GenerateDataFromDwi() { typedef itk::OdfMaximaExtractionFilter< float > MaximaExtractionFilterType; MaximaExtractionFilterType::Pointer filter = MaximaExtractionFilterType::New(); mitk::Geometry3D::Pointer geometry; if (!m_ImageNodes.empty()) { try{ Image::Pointer img = dynamic_cast(m_ImageNodes.at(0)->GetData()); typedef ImageToItk< MaximaExtractionFilterType::CoefficientImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(img); caster->Update(); filter->SetShCoeffImage(caster->GetOutput()); geometry = img->GetGeometry(); } catch(itk::ExceptionObject &e) { MITK_INFO << "wrong image type: " << e.what(); return; } } else return; filter->SetMaxNumPeaks(m_Controls->m_MaxNumPeaksBox->value()); filter->SetPeakThreshold(m_Controls->m_PeakThresholdBox->value()); if (!m_BinaryImageNodes.empty()) { ItkUcharImgType::Pointer itkMaskImage = ItkUcharImgType::New(); Image::Pointer mitkMaskImg = dynamic_cast(m_BinaryImageNodes.at(0)->GetData()); CastToItkImage(mitkMaskImg, itkMaskImage); filter->SetMaskImage(itkMaskImage); } switch (m_Controls->m_NormalizationBox->currentIndex()) { case 0: filter->SetNormalizationMethod(MaximaExtractionFilterType::NO_NORM); break; case 1: filter->SetNormalizationMethod(MaximaExtractionFilterType::MAX_VEC_NORM); break; case 2: filter->SetNormalizationMethod(MaximaExtractionFilterType::SINGLE_VEC_NORM); break; } filter->GenerateData(); ItkUcharImgType::Pointer numDirImage = filter->GetNumDirectionsImage(); if (m_Controls->m_OutputDirectionImagesBox->isChecked()) { typedef MaximaExtractionFilterType::ItkDirectionImageContainer ItkDirectionImageContainer; ItkDirectionImageContainer::Pointer container = filter->GetDirectionImageContainer(); for (int i=0; iSize(); i++) { MaximaExtractionFilterType::ItkDirectionImage::Pointer itkImg = container->GetElement(i); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk( itkImg.GetPointer() ); img->SetVolume( itkImg->GetBufferPointer() ); DataNode::Pointer node = DataNode::New(); node->SetData(img); QString name(m_ImageNodes.at(0)->GetName().c_str()); name += "_Direction"; name += QString::number(i+1); node->SetName(name.toStdString().c_str()); GetDataStorage()->Add(node); } } if (m_Controls->m_OutputNumDirectionsBox->isChecked()) { mitk::Image::Pointer image2 = mitk::Image::New(); image2->InitializeByItk( numDirImage.GetPointer() ); image2->SetVolume( numDirImage->GetBufferPointer() ); DataNode::Pointer node = DataNode::New(); node->SetData(image2); QString name(m_ImageNodes.at(0)->GetName().c_str()); name += "_NumDirections"; node->SetName(name.toStdString().c_str()); GetDataStorage()->Add(node); } if (m_Controls->m_OutputVectorFieldBox->isChecked()) { mitk::Vector3D outImageSpacing = geometry->GetSpacing(); float minSpacing = 1; if(outImageSpacing[0]GetOutputFiberBundle(); directions->SetGeometry(geometry); DataNode::Pointer node = DataNode::New(); node->SetData(directions); QString name(m_ImageNodes.at(0)->GetName().c_str()); name += "_VectorField"; node->SetName(name.toStdString().c_str()); node->SetProperty("Fiber2DSliceThickness", mitk::FloatProperty::New(minSpacing)); node->SetProperty("Fiber2DfadeEFX", mitk::BoolProperty::New(false)); GetDataStorage()->Add(node); } } void QmitkOdfMaximaExtractionView::StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget) { m_MultiWidget = &stdMultiWidget; } void QmitkOdfMaximaExtractionView::StdMultiWidgetNotAvailable() { m_MultiWidget = NULL; } void QmitkOdfMaximaExtractionView::OnSelectionChanged( std::vector nodes ) { m_Controls->m_InputData->setTitle("Please Select Input Data"); m_Controls->m_DwiFibLabel->setText("mandatory"); m_Controls->m_MaskLabel->setText("optional"); m_BinaryImageNodes.clear(); m_ImageNodes.clear(); m_TensorImageNodes.clear(); // iterate all selected objects, adjust warning visibility for( std::vector::iterator it = nodes.begin(); it != nodes.end(); ++it ) { mitk::DataNode::Pointer node = *it; if ( node.IsNotNull() && dynamic_cast(node->GetData()) ) { m_TensorImageNodes.push_back(node); } else if( node.IsNotNull() && dynamic_cast(node->GetData()) ) { bool isBinary = false; node->GetPropertyValue("binary", isBinary); if (isBinary) m_BinaryImageNodes.push_back(node); else m_ImageNodes.push_back(node); } } UpdateGui(); } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkOdfMaximaExtractionView.h b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkOdfMaximaExtractionView.h index eddc9844c9..0a0ad0a2c3 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkOdfMaximaExtractionView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkOdfMaximaExtractionView.h @@ -1,92 +1,91 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include "ui_QmitkOdfMaximaExtractionViewControls.h" #include #include #include /*! \brief View providing several methods to extract peaks from the spherical harmonic representation of ODFs or from tensors \sa QmitkFunctionality \ingroup Functionalities */ // Forward Qt class declarations class QmitkOdfMaximaExtractionView : public QmitkFunctionality { // this is needed for all Qt objects that should have a Qt meta-object // (everything that derives from QObject and wants to have signal/slots) Q_OBJECT public: static const std::string VIEW_ID; QmitkOdfMaximaExtractionView(); virtual ~QmitkOdfMaximaExtractionView(); virtual void CreateQtPartControl(QWidget *parent); virtual void StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget); virtual void StdMultiWidgetNotAvailable(); typedef itk::Image ItkUcharImgType; typedef itk::Image< itk::DiffusionTensor3D< float >, 3 > ItkTensorImage; typedef itk::Image< itk::Vector< float, 3>, 3 > ItkDirectionImage3DType; ///< contains a 3D vector in each voxel protected slots: - void ConvertShCoeffsFromFsl(); ///< convert fsl spherical harmonic coefficients to the according mitk datatype - void ConvertPeaksFromFsl(); ///< convert fsl peaks to the according mitk datatype - void ConvertPeaksFromMrtrix(); ///< convert mrtrix peaks to the according mitk datatype + void ConvertShCoeffs(); ///< convert spherical harmonic coefficients to the according mitk datatype + void ConvertPeaks(); ///< convert peak files from other toolkits to the according mitk datatype void GenerateImage(); ///< semicontinuous ODF maxima extraction void StartFiniteDiff(); ///< ODF maxima extraction using finite differences on the densely sampled sphere void StartTensor(); ///< extract principal eigenvectors from tensor image protected: /// \brief called by QmitkFunctionality when DataManager's selection has changed virtual void OnSelectionChanged( std::vector nodes ); Ui::QmitkOdfMaximaExtractionViewControls* m_Controls; QmitkStdMultiWidget* m_MultiWidget; std::vector< mitk::DataNode::Pointer > m_BinaryImageNodes; ///< mask images std::vector< mitk::DataNode::Pointer > m_ImageNodes; std::vector< mitk::DataNode::Pointer > m_TensorImageNodes; void UpdateGui(); ///< update button activity etc. dpending on current datamanager selection void GenerateDataFromDwi(); ///< semicontinuous ODF maxima extraction - template void TemplatedConvertShCoeffsFromFsl(mitk::Image* mitkImg); + template void TemplatedConvertShCoeffs(mitk::Image* mitkImg); template void StartMaximaExtraction(); ///< ODF maxima extraction using finite differences on the densely sampled sphere private: }; diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkOdfMaximaExtractionViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkOdfMaximaExtractionViewControls.ui index 30735a8fd0..d89adae693 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkOdfMaximaExtractionViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkOdfMaximaExtractionViewControls.ui @@ -1,469 +1,473 @@ QmitkOdfMaximaExtractionViewControls 0 0 392 - 769 + 761 Form false Extract ODF peaks using a semicontinuous method (Aganj et al. 2010). EXPERIMENTAL! Start Analytical Extraction (only SH order 4) Please Select Input Data Select a tensor image or a SH coefficient image (generate using Q-Ball reconstruction view). ShCoeff/DTI: Mask Image: <html><head/><body><p><span style=" color:#ff0000;">mandatory</span></p></body></html> <html><head/><body><p><span style=" color:#969696;">optional</span></p></body></html> Parameters QFrame::NoFrame QFrame::Raised 0 6 Vector normalization: <html><head/><body><p>The vector fields are always coorected for image spacing and using the lagest eigenvalue in case of the tensor peak extraction. This is done for visualizytion purposes. The output direction images are not affected.</p></body></html> 1 No Normalization MAX Normalize Single Vec Normalization false QFrame::NoFrame QFrame::Raised 0 6 SH Order: Absolute threshold: false Absolute peak threshold (only used for the finite differences method). The value is additionally scaled by 1/GFA. 3 0.000000000000000 1.000000000000000 0.001000000000000 0.010000000000000 Max. Peaks: false Peak threshold relative to the largest peak per voxel. 0.000000000000000 1.000000000000000 0.050000000000000 0.400000000000000 false Maximum number of peaks to extract. 1 1000 3 Relative threshold: false 1 2 4 6 8 10 12 Clustering angle: Cluster close directions. Define "close" here. 90 25 Angular threshold: Discard smaller peaks in the defined angle around the maximum peaks. 0 90 0 Qt::Vertical 20 259 Import From Other Tools - + false - Generate Q-Ball image and MITK compatible SH coefficient from the FSL qboot SH output. + Generate Q-Ball image and MITK compatible SH coefficient from other toolkits. - Convert FSL SH - Coefficients + Import SH - Coefficients - - - - false - + + - Generate vector field and direction images from the FSL qboot peak extraction output. - - - Convert FSL Peaks + Define SH coefficient convention (depends on toolkit) + + + FSL + + + + + MRtrix + + - - + + - true + false - + Generate vector field and direction images from the FSL qboot peak extraction output. - Convert MRtrix Peaks + Import Peak Image false Extract ODF peaks using finite differences on the densely sampled ODF surface. Start Finite Differences Extraction false Extract principal eigenvectors of input tensors. Start Tensor Principal Direction Extraction Output QFormLayout::AllNonFixedFieldsGrow Only for visualization purposes! The vectors are automatically corrected for image spacing and for the largest eigenvalue in case of the tensor peak extraction. Vector Field true Output unsigned char image containing the number of directions per voxel. #Directions per Voxel false Output one image per extracted direction containing the direction vecors as pixel values. Direction Images false