diff --git a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkFiniteDiffOdfMaximaExtractionFilter.h b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkFiniteDiffOdfMaximaExtractionFilter.h index 63ce806a2d..a11ed68a6b 100644 --- a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkFiniteDiffOdfMaximaExtractionFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkFiniteDiffOdfMaximaExtractionFilter.h @@ -1,149 +1,149 @@ /*========================================================================= 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< unsigned char, 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) + MAX_VEC_NORM ///< normalize all peaks according to their length in comparison to the largest peak in the respective voxel (0-1) }; typedef FiniteDiffOdfMaximaExtractionFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ImageToImageFilter< Image< Vector< PixelType, (ShOrder*ShOrder + ShOrder + 2)/2 + ShOrder >, 3 >, Image< unsigned char, 3 > > Superclass; /** Method for creation through the object factory. */ itkFactorylessNewMacro(Self) itkCloneMacro(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 typename Superclass::InputImageRegionType InputImageRegionType; typedef Image< float, 4 > PeakImageType; typedef OrientationDistributionFunction OdfType; typedef itk::Image ItkUcharImgType; typedef vnl_vector_fixed< double, 3 > DirectionType; // 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 itkSetMacro( FlipX, bool) ///< flip peaks in x direction itkSetMacro( FlipY, bool) ///< flip peaks in y direction itkSetMacro( FlipZ, bool) ///< flip peaks in z direction itkSetMacro( ApplyDirectionMatrix, bool) // output itkGetMacro( NumDirectionsImage, ItkUcharImgType::Pointer ) itkGetMacro( PeakImage, PeakImageType::Pointer ) 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, ThreadIdType 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 PeakImageType::Pointer m_PeakImage; ItkUcharImgType::Pointer m_NumDirectionsImage; ///< number of peaks per voxel ItkUcharImgType::Pointer m_MaskImage; ///< only voxels inside the binary mask are processed Toolkit m_Toolkit; bool m_FlipX; bool m_FlipY; bool m_FlipZ; bool m_ApplyDirectionMatrix; }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkFiniteDiffOdfMaximaExtractionFilter.cpp" #endif #endif //__itkFiniteDiffOdfMaximaExtractionFilter_h_ diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.cpp index a5213178b1..962cb70c17 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.cpp @@ -1,386 +1,412 @@ #include "itkTractsToVectorImageFilter.h" // VTK #include #include #include // ITK #include #include // misc #define _USE_MATH_DEFINES #include #include namespace itk{ static inline bool CompareVectorLengths(const vnl_vector_fixed< double, 3 >& v1, const vnl_vector_fixed< double, 3 >& v2) { - return (v1.magnitude()>v2.magnitude()); + return (v1.magnitude()>v2.magnitude()); } template< class PixelType > TractsToVectorImageFilter< PixelType >::TractsToVectorImageFilter(): - m_AngularThreshold(0.7), - m_Epsilon(0.999), - m_NormalizeVectors(false), - m_UseWorkingCopy(true), - m_MaxNumDirections(3), - m_SizeThreshold(0.3) + m_NormalizationMethod(GLOBAL_MAX), + m_AngularThreshold(0.7), + m_Epsilon(0.999), + m_UseWorkingCopy(true), + m_MaxNumDirections(3), + m_SizeThreshold(0.3) { - this->SetNumberOfRequiredOutputs(1); + this->SetNumberOfRequiredOutputs(1); } template< class PixelType > TractsToVectorImageFilter< PixelType >::~TractsToVectorImageFilter() { } template< class PixelType > vnl_vector_fixed TractsToVectorImageFilter< PixelType >::GetVnlVector(double point[]) { - vnl_vector_fixed vnlVector; - vnlVector[0] = point[0]; - vnlVector[1] = point[1]; - vnlVector[2] = point[2]; - return vnlVector; + vnl_vector_fixed vnlVector; + vnlVector[0] = point[0]; + vnlVector[1] = point[1]; + vnlVector[2] = point[2]; + return vnlVector; } template< class PixelType > itk::Point TractsToVectorImageFilter< PixelType >::GetItkPoint(double point[]) { - itk::Point itkPoint; - itkPoint[0] = point[0]; - itkPoint[1] = point[1]; - itkPoint[2] = point[2]; - return itkPoint; + itk::Point itkPoint; + itkPoint[0] = point[0]; + itkPoint[1] = point[1]; + itkPoint[2] = point[2]; + return itkPoint; } template< class PixelType > void TractsToVectorImageFilter< PixelType >::GenerateData() { - mitk::BaseGeometry::Pointer geometry = m_FiberBundle->GetGeometry(); - - // calculate new image parameters - itk::Vector spacing3; - itk::Point origin3; - itk::Matrix direction3; - ImageRegion<3> imageRegion3; - if (!m_MaskImage.IsNull()) + mitk::BaseGeometry::Pointer geometry = m_FiberBundle->GetGeometry(); + + // calculate new image parameters + itk::Vector spacing3; + itk::Point origin3; + itk::Matrix direction3; + ImageRegion<3> imageRegion3; + if (!m_MaskImage.IsNull()) + { + spacing3 = m_MaskImage->GetSpacing(); + imageRegion3 = m_MaskImage->GetLargestPossibleRegion(); + origin3 = m_MaskImage->GetOrigin(); + direction3 = m_MaskImage->GetDirection(); + } + else + { + spacing3 = geometry->GetSpacing(); + origin3 = geometry->GetOrigin(); + mitk::BaseGeometry::BoundsArrayType bounds = geometry->GetBounds(); + origin3[0] += bounds.GetElement(0); + origin3[1] += bounds.GetElement(2); + origin3[2] += bounds.GetElement(4); + + for (int i=0; i<3; i++) + for (int j=0; j<3; j++) + direction3[j][i] = geometry->GetMatrixColumn(i)[j]; + imageRegion3.SetSize(0, geometry->GetExtent(0)+1); + imageRegion3.SetSize(1, geometry->GetExtent(1)+1); + imageRegion3.SetSize(2, geometry->GetExtent(2)+1); + + + m_MaskImage = ItkUcharImgType::New(); + m_MaskImage->SetSpacing( spacing3 ); + m_MaskImage->SetOrigin( origin3 ); + m_MaskImage->SetDirection( direction3 ); + m_MaskImage->SetRegions( imageRegion3 ); + m_MaskImage->Allocate(); + m_MaskImage->FillBuffer(1); + } + OutputImageType::RegionType::SizeType outImageSize = imageRegion3.GetSize(); + m_OutImageSpacing = m_MaskImage->GetSpacing(); + m_ClusteredDirectionsContainer = ContainerType::New(); + + // initialize num directions image + 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); + + itk::Vector spacing4; + itk::Point origin4; + itk::Matrix direction4; + itk::ImageRegion<4> imageRegion4; + + spacing4[0] = spacing3[0]; spacing4[1] = spacing3[1]; spacing4[2] = spacing3[2]; spacing4[3] = 1; + origin4[0] = origin3[0]; origin4[1] = origin3[1]; origin4[2] = origin3[2]; origin3[3] = 0; + for (int r=0; r<3; r++) + for (int c=0; c<3; c++) + direction4[r][c] = direction3[r][c]; + direction4[3][3] = 1; + imageRegion4.SetSize(0, imageRegion3.GetSize()[0]); + imageRegion4.SetSize(1, imageRegion3.GetSize()[1]); + imageRegion4.SetSize(2, imageRegion3.GetSize()[2]); + imageRegion4.SetSize(3, m_MaxNumDirections*3); + + m_DirectionImage = ItkDirectionImageType::New(); + m_DirectionImage->SetSpacing( spacing4 ); + m_DirectionImage->SetOrigin( origin4 ); + m_DirectionImage->SetDirection( direction4 ); + m_DirectionImage->SetRegions( imageRegion4 ); + m_DirectionImage->Allocate(); + m_DirectionImage->FillBuffer(0.0); + + // resample fiber bundle + double minSpacing = 1; + if(m_OutImageSpacing[0]GetDeepCopy(); + + // resample fiber bundle for sufficient voxel coverage + m_FiberBundle->ResampleLinear(minSpacing/10); + + // iterate over all fibers + vtkSmartPointer fiberPolyData = m_FiberBundle->GetFiberPolyData(); + int numFibers = m_FiberBundle->GetNumFibers(); + m_DirectionsContainer = ContainerType::New(); + + VectorContainer< unsigned int, std::vector< double > >::Pointer peakLengths = VectorContainer< unsigned int, std::vector< double > >::New(); + + MITK_INFO << "Generating directions from tractogram"; + boost::progress_display disp(numFibers); + for( int i=0; iGetCell(i); + int numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); + if (numPoints<2) + continue; + + vnl_vector_fixed dir; + itk::Point worldPos; + vnl_vector v; + + + float fiberWeight = m_FiberBundle->GetFiberWeight(i); + + for( int j=0; jGetSpacing(); - imageRegion3 = m_MaskImage->GetLargestPossibleRegion(); - origin3 = m_MaskImage->GetOrigin(); - direction3 = m_MaskImage->GetDirection(); - } - else - { - spacing3 = geometry->GetSpacing(); - origin3 = geometry->GetOrigin(); - mitk::BaseGeometry::BoundsArrayType bounds = geometry->GetBounds(); - origin3[0] += bounds.GetElement(0); - origin3[1] += bounds.GetElement(2); - origin3[2] += bounds.GetElement(4); - - for (int i=0; i<3; i++) - for (int j=0; j<3; j++) - direction3[j][i] = geometry->GetMatrixColumn(i)[j]; - imageRegion3.SetSize(0, geometry->GetExtent(0)+1); - imageRegion3.SetSize(1, geometry->GetExtent(1)+1); - imageRegion3.SetSize(2, geometry->GetExtent(2)+1); - - - m_MaskImage = ItkUcharImgType::New(); - m_MaskImage->SetSpacing( spacing3 ); - m_MaskImage->SetOrigin( origin3 ); - m_MaskImage->SetDirection( direction3 ); - m_MaskImage->SetRegions( imageRegion3 ); - m_MaskImage->Allocate(); - m_MaskImage->FillBuffer(1); + // get current position along fiber in world coordinates + double* temp = points->GetPoint(j); + worldPos = GetItkPoint(temp); + itk::Index<3> index; + m_MaskImage->TransformPhysicalPointToIndex(worldPos, index); + if (!m_MaskImage->GetLargestPossibleRegion().IsInside(index) || m_MaskImage->GetPixel(index)==0) + continue; + + // get fiber tangent direction at this position + v = GetVnlVector(temp); + dir = GetVnlVector(points->GetPoint(j+1))-v; + if (dir.is_zero()) + continue; + dir.normalize(); + + // add direction to container + unsigned int idx = index[0] + outImageSize[0]*(index[1] + outImageSize[1]*index[2]); + DirectionContainerType::Pointer dirCont; + if (m_DirectionsContainer->IndexExists(idx)) + { + peakLengths->ElementAt(idx).push_back(fiberWeight); + + dirCont = m_DirectionsContainer->GetElement(idx); + if (dirCont.IsNull()) + { + dirCont = DirectionContainerType::New(); + dirCont->push_back(dir); + m_DirectionsContainer->InsertElement(idx, dirCont); + } + else + dirCont->push_back(dir); + } + else + { + dirCont = DirectionContainerType::New(); + dirCont->push_back(dir); + m_DirectionsContainer->InsertElement(idx, dirCont); + + std::vector< double > lengths; lengths.push_back(fiberWeight); + peakLengths->InsertElement(idx, lengths); + } } - OutputImageType::RegionType::SizeType outImageSize = imageRegion3.GetSize(); - m_OutImageSpacing = m_MaskImage->GetSpacing(); - m_ClusteredDirectionsContainer = ContainerType::New(); - - // initialize num directions image - 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); - - itk::Vector spacing4; - itk::Point origin4; - itk::Matrix direction4; - itk::ImageRegion<4> imageRegion4; - - spacing4[0] = spacing3[0]; spacing4[1] = spacing3[1]; spacing4[2] = spacing3[2]; spacing4[3] = 1; - origin4[0] = origin3[0]; origin4[1] = origin3[1]; origin4[2] = origin3[2]; origin3[3] = 0; - for (int r=0; r<3; r++) - for (int c=0; c<3; c++) - direction4[r][c] = direction3[r][c]; - direction4[3][3] = 1; - imageRegion4.SetSize(0, imageRegion3.GetSize()[0]); - imageRegion4.SetSize(1, imageRegion3.GetSize()[1]); - imageRegion4.SetSize(2, imageRegion3.GetSize()[2]); - imageRegion4.SetSize(3, m_MaxNumDirections*3); - - m_DirectionImage = ItkDirectionImageType::New(); - m_DirectionImage->SetSpacing( spacing4 ); - m_DirectionImage->SetOrigin( origin4 ); - m_DirectionImage->SetDirection( direction4 ); - m_DirectionImage->SetRegions( imageRegion4 ); - m_DirectionImage->Allocate(); - m_DirectionImage->FillBuffer(0.0); - - // resample fiber bundle - double minSpacing = 1; - if(m_OutImageSpacing[0]GetDeepCopy(); + } - // resample fiber bundle for sufficient voxel coverage - m_FiberBundle->ResampleLinear(minSpacing/10); + itk::ImageRegionIterator dirIt(m_NumDirectionsImage, m_NumDirectionsImage->GetLargestPossibleRegion()); - // iterate over all fibers - vtkSmartPointer fiberPolyData = m_FiberBundle->GetFiberPolyData(); - int numFibers = m_FiberBundle->GetNumFibers(); - m_DirectionsContainer = ContainerType::New(); + MITK_INFO << "Clustering directions"; + float max_dir_mag = 0; + boost::progress_display disp2(outImageSize[0]*outImageSize[1]*outImageSize[2]); + while(!dirIt.IsAtEnd()) + { + ++disp2; + OutputImageType::IndexType idx3 = dirIt.GetIndex(); + int idx_lin = idx3[0]+(idx3[1]+idx3[2]*outImageSize[1])*outImageSize[0]; - VectorContainer< unsigned int, std::vector< double > >::Pointer peakLengths = VectorContainer< unsigned int, std::vector< double > >::New(); + itk::Index<4> idx4; idx4[0] = idx3[0]; idx4[1] = idx3[1]; idx4[2] = idx3[2]; - MITK_INFO << "Generating directions from tractogram"; - boost::progress_display disp(numFibers); - for( int i=0; iIndexExists(idx_lin)) { - ++disp; - vtkCell* cell = fiberPolyData->GetCell(i); - int numPoints = cell->GetNumberOfPoints(); - vtkPoints* points = cell->GetPoints(); - if (numPoints<2) - continue; - - vnl_vector_fixed dir; - itk::Point worldPos; - vnl_vector v; - - - float fiberWeight = m_FiberBundle->GetFiberWeight(i); - - for( int j=0; jGetPoint(j); - worldPos = GetItkPoint(temp); - itk::Index<3> index; - m_MaskImage->TransformPhysicalPointToIndex(worldPos, index); - if (!m_MaskImage->GetLargestPossibleRegion().IsInside(index) || m_MaskImage->GetPixel(index)==0) - continue; - - // get fiber tangent direction at this position - v = GetVnlVector(temp); - dir = GetVnlVector(points->GetPoint(j+1))-v; - if (dir.is_zero()) - continue; - dir.normalize(); - - // add direction to container - unsigned int idx = index[0] + outImageSize[0]*(index[1] + outImageSize[1]*index[2]); - DirectionContainerType::Pointer dirCont; - if (m_DirectionsContainer->IndexExists(idx)) - { - peakLengths->ElementAt(idx).push_back(fiberWeight); - - dirCont = m_DirectionsContainer->GetElement(idx); - if (dirCont.IsNull()) - { - dirCont = DirectionContainerType::New(); - dirCont->push_back(dir); - m_DirectionsContainer->InsertElement(idx, dirCont); - } - else - dirCont->push_back(dir); - } - else - { - dirCont = DirectionContainerType::New(); - dirCont->push_back(dir); - m_DirectionsContainer->InsertElement(idx, dirCont); - - std::vector< double > lengths; lengths.push_back(fiberWeight); - peakLengths->InsertElement(idx, lengths); - } - } + ++dirIt; + continue; + } + DirectionContainerType::Pointer dirCont = m_DirectionsContainer->GetElement(idx_lin); + if (dirCont.IsNull() || dirCont->empty()) + { + ++dirIt; + continue; } - vtkSmartPointer m_VtkCellArray = vtkSmartPointer::New(); - vtkSmartPointer m_VtkPoints = vtkSmartPointer::New(); - - itk::ImageRegionIterator dirIt(m_NumDirectionsImage, m_NumDirectionsImage->GetLargestPossibleRegion()); - - MITK_INFO << "Clustering directions"; - boost::progress_display disp2(outImageSize[0]*outImageSize[1]*outImageSize[2]); - while(!dirIt.IsAtEnd()) + DirectionContainerType::Pointer directions; + if (m_MaxNumDirections>0) { - ++disp2; - OutputImageType::IndexType idx3 = dirIt.GetIndex(); - int idx_lin = idx3[0]+(idx3[1]+idx3[2]*outImageSize[1])*outImageSize[0]; + directions = FastClustering(dirCont, peakLengths->GetElement(idx_lin)); + std::sort( directions->begin(), directions->end(), CompareVectorLengths ); + } + else + directions = dirCont; - itk::Index<4> idx4; idx4[0] = idx3[0]; idx4[1] = idx3[1]; idx4[2] = idx3[2]; + unsigned int numDir = directions->size(); + if (m_MaxNumDirections>0 && numDir>m_MaxNumDirections) + numDir = m_MaxNumDirections; - if (!m_DirectionsContainer->IndexExists(idx_lin)) - { - ++dirIt; - continue; - } - DirectionContainerType::Pointer dirCont = m_DirectionsContainer->GetElement(idx_lin); - if (dirCont.IsNull() || dirCont->empty()) - { - ++dirIt; - continue; - } - -// std::vector< double > lengths; lengths.resize(dirCont->size(), 1); // all peaks have size 1 - DirectionContainerType::Pointer directions; - if (m_MaxNumDirections>0) - { - directions = FastClustering(dirCont, peakLengths->GetElement(idx_lin)); - std::sort( directions->begin(), directions->end(), CompareVectorLengths ); - } - else - directions = dirCont; + float voxel_max_mag = 0; + for (unsigned int i=0; iat(i); + float mag = dir.magnitude(); - unsigned int numDir = directions->size(); - if (m_MaxNumDirections>0 && numDir>m_MaxNumDirections) - numDir = m_MaxNumDirections; + if (mag>voxel_max_mag) + voxel_max_mag = mag; + if (mag>max_dir_mag) + max_dir_mag = mag; + } - int count = 0; - for (unsigned int i=0; iat(i); - - if (dir.magnitude()SetPixel(idx4, dir[j]); - } - } - dirIt.Set(count); - ++dirIt; + int count = 0; + for (unsigned int i=0; iat(i); + count++; + + float mag = dir.magnitude(); + if (m_NormalizationMethod==MAX_VEC_NORM && voxel_max_mag>mitk::eps) + dir /= voxel_max_mag; + else if (m_NormalizationMethod==SINGLE_VEC_NORM && mag>mitk::eps) + dir.normalize(); + + for (unsigned int j = 0; j<3; j++) + { + idx4[3] = i*3 + j; + m_DirectionImage->SetPixel(idx4, dir[j]); + } } + dirIt.Set(count); + ++dirIt; + } + + if (m_NormalizationMethod==GLOBAL_MAX && max_dir_mag>0) + { + itk::ImageRegionIterator dirImgIt(m_DirectionImage, m_DirectionImage->GetLargestPossibleRegion()); + while(!dirImgIt.IsAtEnd()) + { + dirImgIt.Set(dirImgIt.Get()/max_dir_mag); + ++dirImgIt; + } + } } template< class PixelType > TractsToVectorImageFilter< PixelType >::DirectionContainerType::Pointer TractsToVectorImageFilter< PixelType >::FastClustering(DirectionContainerType::Pointer inDirs, std::vector< double > lengths) { - DirectionContainerType::Pointer outDirs = DirectionContainerType::New(); - if (inDirs->size()<2) - return inDirs; - - DirectionType oldMean, currentMean; - std::vector< int > touched; - - // initialize - touched.resize(inDirs->size(), 0); - bool free = true; - currentMean = inDirs->at(0); // initialize first seed - currentMean.normalize(); - double length = lengths.at(0); - touched[0] = 1; - std::vector< double > newLengths; - bool meanChanged = false; - - double max = 0; - while (free) + DirectionContainerType::Pointer outDirs = DirectionContainerType::New(); + if (inDirs->size()<2) + { + if (inDirs->size()==1) + inDirs->SetElement(0, inDirs->at(0)*lengths.at(0)); + return inDirs; + } + + DirectionType oldMean, currentMean; + std::vector< int > touched; + + // initialize + touched.resize(inDirs->size(), 0); + bool free = true; + currentMean = inDirs->at(0); // initialize first seed + currentMean.normalize(); + double length = lengths.at(0); + touched[0] = 1; + std::vector< double > newLengths; + bool meanChanged = false; + + double max = 0; + while (free) + { + oldMean.fill(0.0); + + // start mean-shift clustering + double angle = 0; + + while (fabs(dot_product(currentMean, oldMean))<0.99) { - oldMean.fill(0.0); - - // start mean-shift clustering - double angle = 0; - - while (fabs(dot_product(currentMean, oldMean))<0.99) + oldMean = currentMean; + currentMean.fill(0.0); + for (unsigned int i=0; isize(); i++) + { + angle = dot_product(oldMean, inDirs->at(i)); + if (angle>=m_AngularThreshold) { - oldMean = currentMean; - currentMean.fill(0.0); - for (unsigned int i=0; isize(); i++) - { - angle = dot_product(oldMean, inDirs->at(i)); - if (angle>=m_AngularThreshold) - { - currentMean += inDirs->at(i); - if (meanChanged) - length += lengths.at(i); - touched[i] = 1; - meanChanged = true; - } - else if (-angle>=m_AngularThreshold) - { - currentMean -= inDirs->at(i); - if (meanChanged) - length += lengths.at(i); - touched[i] = 1; - meanChanged = true; - } - } - if(!meanChanged) - currentMean = oldMean; - else - currentMean.normalize(); + currentMean += inDirs->at(i); + if (meanChanged) + length += lengths.at(i); + touched[i] = 1; + meanChanged = true; } - - // found stable mean - outDirs->push_back(currentMean); - newLengths.push_back(length); - if (length>max) - max = length; - - // find next unused seed - free = false; - for (unsigned int i=0; iat(i); - free = true; - meanChanged = false; - length = lengths.at(i); - touched[i] = 1; - break; - } + else if (-angle>=m_AngularThreshold) + { + currentMean -= inDirs->at(i); + if (meanChanged) + length += lengths.at(i); + touched[i] = 1; + meanChanged = true; + } + } + if(!meanChanged) + currentMean = oldMean; + else + currentMean.normalize(); } - if (inDirs->size()==outDirs->size()) + // found stable mean + outDirs->push_back(currentMean); + newLengths.push_back(length); + if (length>max) + max = length; + + // find next unused seed + free = false; + for (unsigned int i=0; iat(i); + free = true; + meanChanged = false; + length = lengths.at(i); + touched[i] = 1; + break; + } + } + + if (inDirs->size()==outDirs->size()) + { + if (max>0) { - if (max>0) - for (unsigned int i=0; isize(); i++) - outDirs->SetElement(i, outDirs->at(i)*newLengths.at(i)/max); - return outDirs; + for (unsigned int i=0; isize(); i++) + outDirs->SetElement(i, outDirs->at(i)*newLengths.at(i)); } - else - return FastClustering(outDirs, newLengths); + return outDirs; + } + else + return FastClustering(outDirs, newLengths); } } diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.h index 60a27eb591..e4ed2d9fab 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.h @@ -1,100 +1,106 @@ #ifndef __itkTractsToVectorImageFilter_h__ #define __itkTractsToVectorImageFilter_h__ // MITK #include // ITK #include #include // VTK #include #include #include #include #include namespace itk{ /** * \brief Extracts the voxel-wise main directions of the input fiber bundle. */ template< class PixelType > class TractsToVectorImageFilter : public ImageSource< Image< PixelType, 4 > > { public: - typedef TractsToVectorImageFilter Self; - typedef ProcessObject Superclass; - typedef SmartPointer< Self > Pointer; - typedef SmartPointer< const Self > ConstPointer; - - typedef itk::Vector OutputVectorType; - typedef itk::Image OutputImageType; - typedef std::vector< OutputImageType::Pointer > OutputImageContainerType; - - typedef vnl_vector_fixed< double, 3 > DirectionType; - typedef VectorContainer< unsigned int, DirectionType > DirectionContainerType; - typedef VectorContainer< unsigned int, DirectionContainerType::Pointer > ContainerType; - typedef Image< PixelType, 4 > ItkDirectionImageType; - typedef itk::Image ItkUcharImgType; - typedef itk::Image ItkDoubleImgType; - - itkFactorylessNewMacro(Self) - itkCloneMacro(Self) - itkTypeMacro( TractsToVectorImageFilter, ImageSource ) - - itkSetMacro( SizeThreshold, float) - itkGetMacro( SizeThreshold, float) - itkSetMacro( AngularThreshold, float) ///< cluster directions that are closer together than the specified threshold - itkGetMacro( AngularThreshold, float) ///< cluster directions that are closer together than the specified threshold - itkSetMacro( NormalizeVectors, bool) ///< Normalize vectors to length 1 - itkGetMacro( NormalizeVectors, bool) ///< Normalize vectors to length 1 - itkSetMacro( UseWorkingCopy, bool) ///< Do not modify input fiber bundle. Use a copy. - itkGetMacro( UseWorkingCopy, bool) ///< Do not modify input fiber bundle. Use a copy. - itkSetMacro( MaxNumDirections, unsigned long) ///< If more directions are extracted, only the largest are kept. - itkGetMacro( MaxNumDirections, unsigned long) ///< If more directions are extracted, only the largest are kept. - itkSetMacro( MaskImage, ItkUcharImgType::Pointer) ///< only process voxels inside mask - itkSetMacro( FiberBundle, mitk::FiberBundle::Pointer) ///< input fiber bundle - itkGetMacro( ClusteredDirectionsContainer, ContainerType::Pointer) ///< output directions - itkGetMacro( NumDirectionsImage, ItkUcharImgType::Pointer) ///< number of directions per voxel - itkGetMacro( DirectionImage, typename ItkDirectionImageType::Pointer) ///< output directions - - void GenerateData() override; + + enum NormalizationMethods { + GLOBAL_MAX, ///< global maximum normalization + 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 in the voxel (0-1) + }; + + typedef TractsToVectorImageFilter Self; + typedef ProcessObject Superclass; + typedef SmartPointer< Self > Pointer; + typedef SmartPointer< const Self > ConstPointer; + + typedef itk::Vector OutputVectorType; + typedef itk::Image OutputImageType; + typedef std::vector< OutputImageType::Pointer > OutputImageContainerType; + + typedef vnl_vector_fixed< double, 3 > DirectionType; + typedef VectorContainer< unsigned int, DirectionType > DirectionContainerType; + typedef VectorContainer< unsigned int, DirectionContainerType::Pointer > ContainerType; + typedef Image< PixelType, 4 > ItkDirectionImageType; + typedef itk::Image ItkUcharImgType; + typedef itk::Image ItkDoubleImgType; + + itkFactorylessNewMacro(Self) + itkCloneMacro(Self) + itkTypeMacro( TractsToVectorImageFilter, ImageSource ) + + itkSetMacro( SizeThreshold, float) + itkGetMacro( SizeThreshold, float) + itkSetMacro( AngularThreshold, float) ///< cluster directions that are closer together than the specified threshold + itkGetMacro( AngularThreshold, float) ///< cluster directions that are closer together than the specified threshold + itkSetMacro( NormalizationMethod, NormalizationMethods) ///< normalization method of peaks + itkSetMacro( UseWorkingCopy, bool) ///< Do not modify input fiber bundle. Use a copy. + itkGetMacro( UseWorkingCopy, bool) ///< Do not modify input fiber bundle. Use a copy. + itkSetMacro( MaxNumDirections, unsigned long) ///< If more directions are extracted, only the largest are kept. + itkGetMacro( MaxNumDirections, unsigned long) ///< If more directions are extracted, only the largest are kept. + itkSetMacro( MaskImage, ItkUcharImgType::Pointer) ///< only process voxels inside mask + itkSetMacro( FiberBundle, mitk::FiberBundle::Pointer) ///< input fiber bundle + itkGetMacro( ClusteredDirectionsContainer, ContainerType::Pointer) ///< output directions + itkGetMacro( NumDirectionsImage, ItkUcharImgType::Pointer) ///< number of directions per voxel + itkGetMacro( DirectionImage, typename ItkDirectionImageType::Pointer) ///< output directions + + void GenerateData() override; protected: - DirectionContainerType::Pointer FastClustering(DirectionContainerType::Pointer inDirs, std::vector< double > lengths); ///< cluster fiber directions + DirectionContainerType::Pointer FastClustering(DirectionContainerType::Pointer inDirs, std::vector< double > lengths); ///< cluster fiber directions - vnl_vector_fixed GetVnlVector(double point[3]); - itk::Point GetItkPoint(double point[3]); + vnl_vector_fixed GetVnlVector(double point[3]); + itk::Point GetItkPoint(double point[3]); - TractsToVectorImageFilter(); - virtual ~TractsToVectorImageFilter(); + TractsToVectorImageFilter(); + virtual ~TractsToVectorImageFilter(); - mitk::FiberBundle::Pointer m_FiberBundle; ///< input fiber bundle - float m_AngularThreshold; ///< cluster directions that are closer together than the specified threshold - float m_Epsilon; ///< epsilon for vector equality check - ItkUcharImgType::Pointer m_MaskImage; ///< only voxels inside the binary mask are processed - bool m_NormalizeVectors; ///< normalize vectors to length 1 - itk::Vector m_OutImageSpacing; ///< spacing of output image - ContainerType::Pointer m_DirectionsContainer; ///< container for fiber directions - bool m_UseWorkingCopy; ///< do not modify input fiber bundle but work on copy - unsigned long m_MaxNumDirections; ///< if more directions per voxel are extracted, only the largest are kept - float m_SizeThreshold; + NormalizationMethods m_NormalizationMethod; ///< normalization method of peaks + mitk::FiberBundle::Pointer m_FiberBundle; ///< input fiber bundle + float m_AngularThreshold; ///< cluster directions that are closer together than the specified threshold + float m_Epsilon; ///< epsilon for vector equality check + ItkUcharImgType::Pointer m_MaskImage; ///< only voxels inside the binary mask are processed + itk::Vector m_OutImageSpacing; ///< spacing of output image + ContainerType::Pointer m_DirectionsContainer; ///< container for fiber directions + bool m_UseWorkingCopy; ///< do not modify input fiber bundle but work on copy + unsigned long m_MaxNumDirections; ///< if more directions per voxel are extracted, only the largest are kept + float m_SizeThreshold; - // output datastructures - typename ItkDirectionImageType::Pointer m_DirectionImage; - ContainerType::Pointer m_ClusteredDirectionsContainer; ///< contains direction vectors for each voxel - ItkUcharImgType::Pointer m_NumDirectionsImage; ///< shows number of fibers per voxel + // output datastructures + typename ItkDirectionImageType::Pointer m_DirectionImage; + ContainerType::Pointer m_ClusteredDirectionsContainer; ///< contains direction vectors for each voxel + ItkUcharImgType::Pointer m_NumDirectionsImage; ///< shows number of fibers per voxel }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkTractsToVectorImageFilter.cpp" #endif #endif // __itkTractsToVectorImageFilter_h__ diff --git a/Modules/DiffusionImaging/FiberTracking/Testing/mitkLocalFiberPlausibilityTest.cpp b/Modules/DiffusionImaging/FiberTracking/Testing/mitkLocalFiberPlausibilityTest.cpp index da55d9400b..0dd0aa5f74 100755 --- a/Modules/DiffusionImaging/FiberTracking/Testing/mitkLocalFiberPlausibilityTest.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Testing/mitkLocalFiberPlausibilityTest.cpp @@ -1,161 +1,161 @@ /*=================================================================== 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 #include #include #include #include #include #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include using namespace std; int mitkLocalFiberPlausibilityTest(int argc, char* argv[]) { omp_set_num_threads(1); MITK_TEST_BEGIN("mitkLocalFiberPlausibilityTest"); MITK_TEST_CONDITION_REQUIRED(argc==8,"check for input data") string fibFile = argv[1]; vector< string > referenceImages; referenceImages.push_back(argv[2]); referenceImages.push_back(argv[3]); string LDFP_ERROR_IMAGE = argv[4]; string LDFP_NUM_DIRECTIONS = argv[5]; string LDFP_VECTOR_FIELD = argv[6]; string LDFP_ERROR_IMAGE_IGNORE = argv[7]; float angularThreshold = 30; try { typedef itk::Image ItkUcharImgType; typedef itk::Image< itk::Vector< float, 3>, 3 > ItkDirectionImage3DType; typedef itk::VectorContainer< unsigned int, ItkDirectionImage3DType::Pointer > ItkDirectionImageContainerType; typedef itk::EvaluateDirectionImagesFilter< float > EvaluationFilterType; // load fiber bundle mitk::FiberBundle::Pointer inputTractogram = dynamic_cast(mitk::IOUtil::Load(fibFile)[0].GetPointer()); // load reference directions ItkDirectionImageContainerType::Pointer referenceImageContainer = ItkDirectionImageContainerType::New(); for (unsigned int i=0; i(mitk::IOUtil::Load(referenceImages.at(i))[0].GetPointer()); typedef mitk::ImageToItk< ItkDirectionImage3DType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(img); caster->Update(); ItkDirectionImage3DType::Pointer itkImg = caster->GetOutput(); referenceImageContainer->InsertElement(referenceImageContainer->Size(),itkImg); } catch(...){ MITK_INFO << "could not load: " << referenceImages.at(i); } } ItkUcharImgType::Pointer itkMaskImage = ItkUcharImgType::New(); ItkDirectionImage3DType::Pointer dirImg = referenceImageContainer->GetElement(0); itkMaskImage->SetSpacing( dirImg->GetSpacing() ); itkMaskImage->SetOrigin( dirImg->GetOrigin() ); itkMaskImage->SetDirection( dirImg->GetDirection() ); itkMaskImage->SetLargestPossibleRegion( dirImg->GetLargestPossibleRegion() ); itkMaskImage->SetBufferedRegion( dirImg->GetLargestPossibleRegion() ); itkMaskImage->SetRequestedRegion( dirImg->GetLargestPossibleRegion() ); itkMaskImage->Allocate(); itkMaskImage->FillBuffer(1); // extract directions from fiber bundle itk::TractsToVectorImageFilter::Pointer fOdfFilter = itk::TractsToVectorImageFilter::New(); fOdfFilter->SetFiberBundle(inputTractogram); fOdfFilter->SetMaskImage(itkMaskImage); fOdfFilter->SetAngularThreshold(cos(angularThreshold*M_PI/180)); - fOdfFilter->SetNormalizeVectors(true); + fOdfFilter->SetNormalizationMethod(itk::TractsToVectorImageFilter::NormalizationMethods::SINGLE_VEC_NORM); fOdfFilter->SetMaxNumDirections(3); fOdfFilter->SetSizeThreshold(0.3); fOdfFilter->SetUseWorkingCopy(false); fOdfFilter->SetNumberOfThreads(1); fOdfFilter->Update(); itk::TractsToVectorImageFilter::ItkDirectionImageType::Pointer direction_image = fOdfFilter->GetDirectionImage(); // Get directions and num directions image ItkUcharImgType::Pointer numDirImage = fOdfFilter->GetNumDirectionsImage(); mitk::Image::Pointer mitkNumDirImage = mitk::Image::New(); mitkNumDirImage->InitializeByItk( numDirImage.GetPointer() ); mitkNumDirImage->SetVolume( numDirImage->GetBufferPointer() ); // mitk::FiberBundle::Pointer testDirections = fOdfFilter->GetOutputFiberBundle(); // evaluate directions with missing directions EvaluationFilterType::Pointer evaluationFilter = EvaluationFilterType::New(); // evaluationFilter->SetImageSet(directionImageContainer); evaluationFilter->SetReferenceImageSet(referenceImageContainer); evaluationFilter->SetMaskImage(itkMaskImage); evaluationFilter->SetIgnoreMissingDirections(false); evaluationFilter->Update(); EvaluationFilterType::OutputImageType::Pointer angularErrorImage = evaluationFilter->GetOutput(0); mitk::Image::Pointer mitkAngularErrorImage = mitk::Image::New(); mitkAngularErrorImage->InitializeByItk( angularErrorImage.GetPointer() ); mitkAngularErrorImage->SetVolume( angularErrorImage->GetBufferPointer() ); // evaluate directions without missing directions evaluationFilter->SetIgnoreMissingDirections(true); evaluationFilter->Update(); EvaluationFilterType::OutputImageType::Pointer angularErrorImageIgnore = evaluationFilter->GetOutput(0); mitk::Image::Pointer mitkAngularErrorImageIgnore = mitk::Image::New(); mitkAngularErrorImageIgnore->InitializeByItk( angularErrorImageIgnore.GetPointer() ); mitkAngularErrorImageIgnore->SetVolume( angularErrorImageIgnore->GetBufferPointer() ); mitk::Image::Pointer gtAngularErrorImageIgnore = dynamic_cast(mitk::IOUtil::Load(LDFP_ERROR_IMAGE_IGNORE)[0].GetPointer()); mitk::Image::Pointer gtAngularErrorImage = dynamic_cast(mitk::IOUtil::Load(LDFP_ERROR_IMAGE)[0].GetPointer()); mitk::Image::Pointer gtNumTestDirImage = dynamic_cast(mitk::IOUtil::Load(LDFP_NUM_DIRECTIONS)[0].GetPointer()); MITK_ASSERT_EQUAL(gtAngularErrorImageIgnore, mitkAngularErrorImageIgnore, "Check if error images are equal (ignored missing directions)."); MITK_ASSERT_EQUAL(gtAngularErrorImage, mitkAngularErrorImage, "Check if error images are equal."); MITK_ASSERT_EQUAL(gtNumTestDirImage, mitkNumDirImage, "Check if num direction images are equal."); } catch (itk::ExceptionObject e) { MITK_INFO << e; return EXIT_FAILURE; } catch (std::exception e) { MITK_INFO << e.what(); return EXIT_FAILURE; } catch (...) { MITK_INFO << "ERROR!?!"; return EXIT_FAILURE; } MITK_TEST_END(); } diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberDirectionExtraction.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberDirectionExtraction.cpp index f2357dac75..9867591755 100755 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberDirectionExtraction.cpp +++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberDirectionExtraction.cpp @@ -1,169 +1,180 @@ /*=================================================================== 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 #include "mitkCommandLineParser.h" #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include using namespace std; /*! \brief Extract principal fiber directions from a tractogram */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Fiber Direction Extraction"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setDescription("Extract principal fiber directions from a tractogram"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("input", "i", mitkCommandLineParser::InputFile, "Input:", "input tractogram (.fib/.trk)", us::Any(), false); parser.addArgument("out", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output root", us::Any(), false); parser.addArgument("mask", "m", mitkCommandLineParser::InputFile, "Mask:", "mask image"); parser.addArgument("athresh", "a", mitkCommandLineParser::Float, "Angular threshold:", "angular threshold in degrees. closer fiber directions are regarded as one direction and clustered together.", 25, true); parser.addArgument("peakthresh", "t", mitkCommandLineParser::Float, "Peak size threshold:", "peak size threshold relative to largest peak in voxel", 0.2, true); parser.addArgument("verbose", "v", mitkCommandLineParser::Bool, "Verbose:", "output optional and intermediate calculation results"); parser.addArgument("numdirs", "d", mitkCommandLineParser::Int, "Max. num. directions:", "maximum number of fibers per voxel", 3, true); - parser.addArgument("normalize", "n", mitkCommandLineParser::Bool, "Normalize:", "normalize vectors"); + parser.addArgument("normalization", "n", mitkCommandLineParser::Int, "Normalization method:", "1=global maximum, 2=single vector, 3=voxel-wise maximum", 1); parser.addArgument("file_ending", "f", mitkCommandLineParser::String, "Image type:", ".nrrd, .nii, .nii.gz"); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; string fibFile = us::any_cast(parsedArgs["input"]); string maskImage(""); if (parsedArgs.count("mask")) maskImage = us::any_cast(parsedArgs["mask"]); float peakThreshold = 0.2; if (parsedArgs.count("peakthresh")) peakThreshold = us::any_cast(parsedArgs["peakthresh"]); float angularThreshold = 25; if (parsedArgs.count("athresh")) angularThreshold = us::any_cast(parsedArgs["athresh"]); string outRoot = us::any_cast(parsedArgs["out"]); bool verbose = false; if (parsedArgs.count("verbose")) verbose = us::any_cast(parsedArgs["verbose"]); int maxNumDirs = 3; if (parsedArgs.count("numdirs")) maxNumDirs = us::any_cast(parsedArgs["numdirs"]); - bool normalize = false; - if (parsedArgs.count("normalize")) - normalize = us::any_cast(parsedArgs["normalize"]); + int normalization = 1; + if (parsedArgs.count("normalization")) + normalization = us::any_cast(parsedArgs["normalization"]); std::string file_ending = ".nrrd"; if (parsedArgs.count("file_ending")) file_ending = us::any_cast(parsedArgs["file_ending"]); try { typedef itk::Image ItkUcharImgType; // load fiber bundle mitk::FiberBundle::Pointer inputTractogram = dynamic_cast(mitk::IOUtil::Load(fibFile)[0].GetPointer()); // load/create mask image ItkUcharImgType::Pointer itkMaskImage = nullptr; if (maskImage.compare("")!=0) { std::cout << "Using mask image"; itkMaskImage = ItkUcharImgType::New(); mitk::Image::Pointer mitkMaskImage = dynamic_cast(mitk::IOUtil::Load(maskImage)[0].GetPointer()); mitk::CastToItkImage(mitkMaskImage, itkMaskImage); } // extract directions from fiber bundle itk::TractsToVectorImageFilter::Pointer fOdfFilter = itk::TractsToVectorImageFilter::New(); fOdfFilter->SetFiberBundle(inputTractogram); fOdfFilter->SetMaskImage(itkMaskImage); fOdfFilter->SetAngularThreshold(cos(angularThreshold*M_PI/180)); - fOdfFilter->SetNormalizeVectors(normalize); + switch (normalization) + { + case 1: + fOdfFilter->SetNormalizationMethod(itk::TractsToVectorImageFilter::NormalizationMethods::GLOBAL_MAX); + break; + case 2: + fOdfFilter->SetNormalizationMethod(itk::TractsToVectorImageFilter::NormalizationMethods::SINGLE_VEC_NORM); + break; + case 3: + fOdfFilter->SetNormalizationMethod(itk::TractsToVectorImageFilter::NormalizationMethods::MAX_VEC_NORM); + break; + } fOdfFilter->SetUseWorkingCopy(false); fOdfFilter->SetSizeThreshold(peakThreshold); fOdfFilter->SetMaxNumDirections(maxNumDirs); fOdfFilter->Update(); { itk::TractsToVectorImageFilter::ItkDirectionImageType::Pointer itkImg = fOdfFilter->GetDirectionImage(); typedef itk::ImageFileWriter< itk::TractsToVectorImageFilter::ItkDirectionImageType > WriterType; WriterType::Pointer writer = WriterType::New(); string outfilename = outRoot; outfilename.append("_DIRECTIONS"); outfilename.append(file_ending); writer->SetFileName(outfilename.c_str()); writer->SetInput(itkImg); writer->Update(); } if (verbose) { // write num direction image ItkUcharImgType::Pointer numDirImage = fOdfFilter->GetNumDirectionsImage(); typedef itk::ImageFileWriter< ItkUcharImgType > WriterType; WriterType::Pointer writer = WriterType::New(); string outfilename = outRoot; outfilename.append("_NUM_DIRECTIONS"); outfilename.append(file_ending); writer->SetFileName(outfilename.c_str()); writer->SetInput(numDirImage); writer->Update(); } } catch (itk::ExceptionObject e) { std::cout << e; return EXIT_FAILURE; } catch (std::exception e) { std::cout << e.what(); return EXIT_FAILURE; } catch (...) { std::cout << "ERROR!?!"; return EXIT_FAILURE; } return EXIT_SUCCESS; } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationView.cpp index d786fd9f77..4ff98294ba 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationView.cpp @@ -1,444 +1,450 @@ /*=================================================================== 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. ===================================================================*/ // Blueberry #include #include // Qmitk #include "QmitkFiberQuantificationView.h" // Qt #include // MITK #include #include #include #include #include #include #include #include #include #include #include // ITK #include #include #include #include #include #include #include #include #include #include const std::string QmitkFiberQuantificationView::VIEW_ID = "org.mitk.views.fiberquantification"; const std::string id_DataManager = "org.mitk.views.datamanager"; using namespace mitk; QmitkFiberQuantificationView::QmitkFiberQuantificationView() : QmitkAbstractView() , m_Controls( 0 ) , m_UpsamplingFactor(5) { } // Destructor QmitkFiberQuantificationView::~QmitkFiberQuantificationView() { } void QmitkFiberQuantificationView::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::QmitkFiberQuantificationViewControls; m_Controls->setupUi( parent ); connect( m_Controls->m_ProcessFiberBundleButton, SIGNAL(clicked()), this, SLOT(ProcessSelectedBundles()) ); connect( m_Controls->m_ExtractFiberPeaks, SIGNAL(clicked()), this, SLOT(CalculateFiberDirections()) ); } } void QmitkFiberQuantificationView::SetFocus() { m_Controls->m_ProcessFiberBundleButton->setFocus(); } void QmitkFiberQuantificationView::CalculateFiberDirections() { typedef itk::Image ItkUcharImgType; // load fiber bundle mitk::FiberBundle::Pointer inputTractogram = dynamic_cast(m_SelectedFB.back()->GetData()); itk::TractsToVectorImageFilter::Pointer fOdfFilter = itk::TractsToVectorImageFilter::New(); if (m_SelectedImage.IsNotNull()) { ItkUcharImgType::Pointer itkMaskImage = ItkUcharImgType::New(); mitk::CastToItkImage(m_SelectedImage, itkMaskImage); fOdfFilter->SetMaskImage(itkMaskImage); } // extract directions from fiber bundle fOdfFilter->SetFiberBundle(inputTractogram); fOdfFilter->SetAngularThreshold(cos(m_Controls->m_AngularThreshold->value()*M_PI/180)); - fOdfFilter->SetNormalizeVectors(m_Controls->m_NormalizeDirectionsBox->isChecked()); + switch (m_Controls->m_FiberDirNormBox->currentIndex()) + { + case 0: + fOdfFilter->SetNormalizationMethod(itk::TractsToVectorImageFilter::NormalizationMethods::GLOBAL_MAX); + break; + case 1: + fOdfFilter->SetNormalizationMethod(itk::TractsToVectorImageFilter::NormalizationMethods::SINGLE_VEC_NORM); + break; + case 2: + fOdfFilter->SetNormalizationMethod(itk::TractsToVectorImageFilter::NormalizationMethods::MAX_VEC_NORM); + break; + } fOdfFilter->SetUseWorkingCopy(true); fOdfFilter->SetSizeThreshold(m_Controls->m_PeakThreshold->value()); fOdfFilter->SetMaxNumDirections(m_Controls->m_MaxNumDirections->value()); fOdfFilter->Update(); QString name = m_SelectedFB.back()->GetName().c_str(); if (m_Controls->m_NumDirectionsBox->isChecked()) { mitk::Image::Pointer mitkImage = mitk::Image::New(); mitkImage->InitializeByItk( fOdfFilter->GetNumDirectionsImage().GetPointer() ); mitkImage->SetVolume( fOdfFilter->GetNumDirectionsImage()->GetBufferPointer() ); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(mitkImage); node->SetName((name+"_NUM_DIRECTIONS").toStdString().c_str()); GetDataStorage()->Add(node, m_SelectedFB.back()); } - itk::TractsToVectorImageFilter::ItkDirectionImageType::Pointer itkImg = fOdfFilter->GetDirectionImage(); - - if (itkImg.IsNull()) - return; - - mitk::Image::Pointer mitkImage = dynamic_cast(PeakImage::New().GetPointer()); - mitkImage->InitializeByItk( itkImg.GetPointer() ); - mitkImage->SetVolume( itkImg->GetBufferPointer() ); + Image::Pointer mitkImage = dynamic_cast(PeakImage::New().GetPointer()); + mitk::CastToMitkImage(fOdfFilter->GetDirectionImage(), mitkImage); + mitkImage->SetVolume(fOdfFilter->GetDirectionImage()->GetBufferPointer()); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(mitkImage); node->SetName( (name+"_DIRECTIONS").toStdString().c_str()); GetDataStorage()->Add(node, m_SelectedFB.back()); } void QmitkFiberQuantificationView::UpdateGui() { m_Controls->m_ProcessFiberBundleButton->setEnabled(!m_SelectedFB.empty()); m_Controls->m_ExtractFiberPeaks->setEnabled(!m_SelectedFB.empty()); } void QmitkFiberQuantificationView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList& nodes) { //reset existing Vectors containing FiberBundles and PlanarFigures from a previous selection m_SelectedFB.clear(); m_SelectedSurfaces.clear(); m_SelectedImage = nullptr; for (mitk::DataNode::Pointer node: nodes) { if ( dynamic_cast(node->GetData()) ) { m_SelectedFB.push_back(node); } else if (dynamic_cast(node->GetData())) m_SelectedImage = dynamic_cast(node->GetData()); else if (dynamic_cast(node->GetData())) { m_SelectedSurfaces.push_back(dynamic_cast(node->GetData())); } } UpdateGui(); GenerateStats(); } void QmitkFiberQuantificationView::GenerateStats() { if ( m_SelectedFB.empty() ) return; QString stats(""); for( unsigned int i=0; i(node->GetData())) { if (i>0) stats += "\n-----------------------------\n"; stats += QString(node->GetName().c_str()) + "\n"; mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); stats += "Number of fibers: "+ QString::number(fib->GetNumFibers()) + "\n"; stats += "Number of points: "+ QString::number(fib->GetNumberOfPoints()) + "\n"; stats += "Min. length: "+ QString::number(fib->GetMinFiberLength(),'f',1) + " mm\n"; stats += "Max. length: "+ QString::number(fib->GetMaxFiberLength(),'f',1) + " mm\n"; stats += "Mean length: "+ QString::number(fib->GetMeanFiberLength(),'f',1) + " mm\n"; stats += "Median length: "+ QString::number(fib->GetMedianFiberLength(),'f',1) + " mm\n"; stats += "Standard deviation: "+ QString::number(fib->GetLengthStDev(),'f',1) + " mm\n"; vtkSmartPointer weights = fib->GetFiberWeights(); if (weights!=nullptr) { float weight = -1; int c = 0; for (int i=0; iGetSize(); i++) if (!mitk::Equal(weights->GetValue(i),weight,0.0001)) { weight = weights->GetValue(i); c++; if (c>1) break; } if (c>1) stats += "Detected fiber weights. Fibers are not weighted uniformly.\n"; else stats += "Fibers are weighted equally.\n"; } else stats += "No fiber weight array found.\n"; } } this->m_Controls->m_StatsTextEdit->setText(stats); } void QmitkFiberQuantificationView::ProcessSelectedBundles() { if ( m_SelectedFB.empty() ){ QMessageBox::information( nullptr, "Warning", "No fibe bundle selected!"); MITK_WARN("QmitkFiberQuantificationView") << "no fibe bundle selected"; return; } int generationMethod = m_Controls->m_GenerationBox->currentIndex(); for( unsigned int i=0; i(node->GetData())) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); QString name(node->GetName().c_str()); DataNode::Pointer newNode = nullptr; switch(generationMethod){ case 0: newNode = GenerateTractDensityImage(fib, false, true); name += "_TDI"; break; case 1: newNode = GenerateTractDensityImage(fib, false, false); name += "_TDI"; break; case 2: newNode = GenerateTractDensityImage(fib, true, false); name += "_envelope"; break; case 3: newNode = GenerateColorHeatmap(fib); break; case 4: newNode = GenerateFiberEndingsImage(fib); name += "_fiber_endings"; break; case 5: newNode = GenerateFiberEndingsPointSet(fib); name += "_fiber_endings"; break; } if (newNode.IsNotNull()) { newNode->SetName(name.toStdString()); GetDataStorage()->Add(newNode); } } } } // generate pointset displaying the fiber endings mitk::DataNode::Pointer QmitkFiberQuantificationView::GenerateFiberEndingsPointSet(mitk::FiberBundle::Pointer fib) { mitk::PointSet::Pointer pointSet = mitk::PointSet::New(); vtkSmartPointer fiberPolyData = fib->GetFiberPolyData(); int count = 0; int numFibers = fib->GetNumFibers(); for( int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (numPoints>0) { double* point = points->GetPoint(0); itk::Point itkPoint; itkPoint[0] = point[0]; itkPoint[1] = point[1]; itkPoint[2] = point[2]; pointSet->InsertPoint(count, itkPoint); count++; } if (numPoints>2) { double* point = points->GetPoint(numPoints-1); itk::Point itkPoint; itkPoint[0] = point[0]; itkPoint[1] = point[1]; itkPoint[2] = point[2]; pointSet->InsertPoint(count, itkPoint); count++; } } mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( pointSet ); return node; } // generate image displaying the fiber endings mitk::DataNode::Pointer QmitkFiberQuantificationView::GenerateFiberEndingsImage(mitk::FiberBundle::Pointer fib) { typedef unsigned int OutPixType; typedef itk::Image OutImageType; typedef itk::TractsToFiberEndingsImageFilter< OutImageType > ImageGeneratorType; ImageGeneratorType::Pointer generator = ImageGeneratorType::New(); generator->SetFiberBundle(fib); generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value()); if (m_SelectedImage.IsNotNull()) { OutImageType::Pointer itkImage = OutImageType::New(); CastToItkImage(m_SelectedImage, itkImage); generator->SetInputImage(itkImage); generator->SetUseImageGeometry(true); } generator->Update(); // get output image OutImageType::Pointer outImg = generator->GetOutput(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); // init data node mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(img); return node; } // generate rgba heatmap from fiber bundle mitk::DataNode::Pointer QmitkFiberQuantificationView::GenerateColorHeatmap(mitk::FiberBundle::Pointer fib) { typedef itk::RGBAPixel OutPixType; typedef itk::Image OutImageType; typedef itk::TractsToRgbaImageFilter< OutImageType > ImageGeneratorType; ImageGeneratorType::Pointer generator = ImageGeneratorType::New(); generator->SetFiberBundle(fib); generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value()); if (m_SelectedImage.IsNotNull()) { itk::Image::Pointer itkImage = itk::Image::New(); CastToItkImage(m_SelectedImage, itkImage); generator->SetInputImage(itkImage); generator->SetUseImageGeometry(true); } generator->Update(); // get output image typedef itk::Image OutType; OutType::Pointer outImg = generator->GetOutput(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); // init data node mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(img); return node; } // generate tract density image from fiber bundle mitk::DataNode::Pointer QmitkFiberQuantificationView::GenerateTractDensityImage(mitk::FiberBundle::Pointer fib, bool binary, bool absolute) { mitk::DataNode::Pointer node = mitk::DataNode::New(); if (binary) { typedef unsigned char OutPixType; typedef itk::Image OutImageType; itk::TractDensityImageFilter< OutImageType >::Pointer generator = itk::TractDensityImageFilter< OutImageType >::New(); generator->SetFiberBundle(fib); generator->SetBinaryOutput(binary); generator->SetOutputAbsoluteValues(absolute); generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value()); if (m_SelectedImage.IsNotNull()) { OutImageType::Pointer itkImage = OutImageType::New(); CastToItkImage(m_SelectedImage, itkImage); generator->SetInputImage(itkImage); generator->SetUseImageGeometry(true); } generator->Update(); // get output image typedef itk::Image OutType; OutType::Pointer outImg = generator->GetOutput(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); // init data node node->SetData(img); } else { typedef float OutPixType; typedef itk::Image OutImageType; itk::TractDensityImageFilter< OutImageType >::Pointer generator = itk::TractDensityImageFilter< OutImageType >::New(); generator->SetFiberBundle(fib); generator->SetBinaryOutput(binary); generator->SetOutputAbsoluteValues(absolute); generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value()); if (m_SelectedImage.IsNotNull()) { OutImageType::Pointer itkImage = OutImageType::New(); CastToItkImage(m_SelectedImage, itkImage); generator->SetInputImage(itkImage); generator->SetUseImageGeometry(true); } //generator->SetDoFiberResampling(false); generator->Update(); // get output image typedef itk::Image OutType; OutType::Pointer outImg = generator->GetOutput(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); // init data node node->SetData(img); } return node; } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationViewControls.ui index 300b989135..399ba8db30 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationViewControls.ui @@ -1,344 +1,357 @@ QmitkFiberQuantificationViewControls 0 0 484 574 Form Qt::Vertical 20 40 Fiber-derived images QFormLayout::AllNonFixedFieldsGrow false 0 0 200 16777215 11 Perform selected operation on all selected fiber bundles. Generate Image 0 0 QFrame::NoFrame QFrame::Raised 0 0 0 0 0 0 0 0 Tract Density Image (TDI) Normalized TDI Binary Envelope Fiber Bundle Image Fiber Endings Image Fiber Endings Pointset 0 0 Upsampling factor 1 0.100000000000000 10.000000000000000 0.100000000000000 1.000000000000000 Fiber Statistics Courier 10 Pitch false true Principal Fiber Directions Maximum number of fiber directions per voxel. 100 3 Image containing the number of distinct fiber clusters per voxel. Output #Directions per Voxel true false 0 0 200 16777215 11 Start Max. clusters: Size Threshold: - - - - <html><head/><body><p>Normalize output directions to length 1 otherwise directions are max normalized (voxel-wise).</p></body></html> - - - Normalize directions - - - true - - - Angular Threshold: + + + + <html><head/><body><p>Directions shorter than the defined threshold are discarded.</p></body></html> + + + 3 + + + 1.000000000000000 + + + 0.100000000000000 + + + 0.300000000000000 + + + Fiber directions with an angle smaller than the defined threshold are clustered. 2 0.000000000000000 90.000000000000000 1.000000000000000 30.000000000000000 - - - - <html><head/><body><p>Directions shorter than the defined threshold are discarded.</p></body></html> - - - 3 - - - 1.000000000000000 - - - 0.100000000000000 - - - 0.300000000000000 + + + + Normalization: + + + + + Global maximum + + + + + Single vector + + + + + Voxel-wise maximum + + + +