diff --git a/Modules/GraphAlgorithms/itkShortestPathCostFunction.h b/Modules/GraphAlgorithms/itkShortestPathCostFunction.h index 986d75d4a2..34595b3f81 100644 --- a/Modules/GraphAlgorithms/itkShortestPathCostFunction.h +++ b/Modules/GraphAlgorithms/itkShortestPathCostFunction.h @@ -1,97 +1,90 @@ /*=================================================================== 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. ===================================================================*/ #ifndef __itkShortestPathCostFunction_h #define __itkShortestPathCostFunction_h #include "itkObject.h" #include "itkObjectFactory.h" #include "itkShapedNeighborhoodIterator.h" -#include namespace itk { - // \brief this is a pure virtual superclass for all cost function for the itkShortestPathImageFilter + // \brief this is a pure virtual superclass for all cost functions used in itkShortestPathImageFilter template class ShortestPathCostFunction : public Object { public: /** Standard class typedefs. */ typedef ShortestPathCostFunction Self; typedef Object Superclass; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; - typedef ShapedNeighborhoodIterator< TInputImageType > itkShapedNeighborhoodIteratorType; - + typedef ShapedNeighborhoodIterator< TInputImageType > ShapedNeighborhoodIteratorType; /** Run-time type information (and related methods). */ itkTypeMacro(ShortestPathCostFunction, Object); /** Type definition for the input image. */ typedef TInputImageType ImageType; // More typdefs for convenience typedef typename TInputImageType::Pointer ImagePointer; typedef typename TInputImageType::ConstPointer ImageConstPointer; - - typedef typename TInputImageType::PixelType PixelType; - typedef typename TInputImageType::IndexType IndexType; - /** Set the input image. */ + // \brief Set the input image. itkSetConstObjectMacro(Image,TInputImageType); - // \brief calculates the costs for going from pixel1 to pixel2 + // \brief Calculate the cost for going from pixel p1 to pixel p2 virtual double GetCost( IndexType p1, IndexType p2) = 0; - // \brief returns the minimal costs possible (needed for A*) + // \brief Return the minimal possible cost (needed for A*) virtual double GetMinCost() = 0; // \brief Initialize the metric virtual void Initialize () = 0; - // \brief Set Starpoint for Path - void SetStartIndex (const IndexType & StartIndex); - - // \brief Set Endpoint for Path - void SetEndIndex(const IndexType & EndIndex); - + // \brief Set the starting index of a path + void SetStartIndex (const IndexType & index); - ShortestPathCostFunction(); + // \brief Set the ending index of a path + void SetEndIndex(const IndexType & index); protected: + ShortestPathCostFunction() {}; virtual ~ShortestPathCostFunction() {}; void PrintSelf(std::ostream& os, Indent indent) const; ImageConstPointer m_Image; IndexType m_StartIndex, m_EndIndex; private: ShortestPathCostFunction(const Self&); //purposely not implemented void operator=(const Self&); //purposely not implemented }; } // end namespace itk #include "itkShortestPathCostFunction.txx" #endif /* __itkShortestPathCostFunction_h */ diff --git a/Modules/GraphAlgorithms/itkShortestPathCostFunction.txx b/Modules/GraphAlgorithms/itkShortestPathCostFunction.txx index d620e89be6..e7b4141b06 100644 --- a/Modules/GraphAlgorithms/itkShortestPathCostFunction.txx +++ b/Modules/GraphAlgorithms/itkShortestPathCostFunction.txx @@ -1,64 +1,57 @@ /*=================================================================== 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. ===================================================================*/ -#ifndef __itkShortestPathCostFunction_cpp -#define __itkShortestPathCostFunction_cpp +#ifndef __itkShortestPathCostFunction_txx +#define __itkShortestPathCostFunction_txx +#include "itkShortestPathCostFunction.h" namespace itk { - // Constructor - template - ShortestPathCostFunction - ::ShortestPathCostFunction() - { - } - - template void ShortestPathCostFunction ::PrintSelf( std::ostream& os, Indent indent ) const { Superclass::PrintSelf(os,indent); } template void ShortestPathCostFunction:: - SetStartIndex (const typename TInputImageType::IndexType &StartIndex) + SetStartIndex (const typename TInputImageType::IndexType &index) { for (unsigned int i=0;i void ShortestPathCostFunction:: - SetEndIndex (const typename TInputImageType::IndexType &EndIndex) + SetEndIndex (const typename TInputImageType::IndexType &index) { for (unsigned int i=0;i - +#include "itkShortestPathCostFunction.h" +#include "itkImageRegionConstIterator.h" namespace itk { /** \brief Cost function for LiveWire purposes. Specific features are considered to calculate cummulative costs of a link between two pixels. These are: - Gradient Magnitude - Gradient Direction - Laplacian Zero Crossing - By default the Gradient Magnitude is mapped linear to costs + By default the Gradient Magnitude is mapped linearly to cost values between 0 (good) and 1 (bad). Via SetDynamicCostMap( std::map< int, int > &costMap) a cost map can be set to dynamically map Gradient Magnitude (non - linear). Thus lower values can be considered with lower costs + linear). Thus, lower values can be considered with lower costs than higher values of gradient magnitudes. To compute the costs of the gradient magnitude dynamically - a iverted map of the histogram of gradient magnitude image is used. + an iverted map of the histogram of gradient magnitude image is used. */ template class ITK_EXPORT ShortestPathCostFunctionLiveWire : public ShortestPathCostFunction { public: /** Standard class typedefs. */ typedef ShortestPathCostFunctionLiveWire Self; typedef ShortestPathCostFunction Superclass; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef itk::ImageRegionConstIterator ConstIteratorType; /** Method for creation through the object factory. */ itkNewMacro(Self); /** Run-time type information (and related methods). */ itkTypeMacro(ShortestPathCostFunctionLiveWire, ShortestPathCostFunction); - typedef itk::Image UnsignedCharImageType; - typedef itk::Image FloatImageType; + typedef itk::Image UnsignedCharImageType; + typedef itk::Image FloatImageType; typedef float ComponentType; - typedef itk::CovariantVector< ComponentType, 2 > OutputPixelType; - typedef itk::Image< OutputPixelType, 2 > VectorOutputImageType; - - typedef typename TInputImageType::IndexType IndexType; - typedef TInputImageType ImageType; - typedef itk::ImageRegion<2> RegionType; + typedef itk::CovariantVector< ComponentType, 2 > OutputPixelType; + typedef itk::Image< OutputPixelType, 2 > VectorOutputImageType; + typedef typename TInputImageType::IndexType IndexType; + typedef TInputImageType ImageType; + typedef itk::ImageRegion<2> RegionType; /** \brief calculates the costs for going from p1 to p2*/ virtual double GetCost(IndexType p1, IndexType p2); /** \brief returns the minimal costs possible (needed for A*)*/ virtual double GetMinCost(); /** \brief Initialize the metric*/ virtual void Initialize (); - - /** \brief Set repulsive path*/ - virtual void AddRepulsivePoint( itk::Index<3> ); + /** \brief Add repulsive point*/ + virtual void AddRepulsivePoint( IndexType index ); /** \brief Clear repulsive path*/ virtual void ClearRepulsivePoints( ); - - - ShortestPathCostFunctionLiveWire(); - - itkSetMacro (RequestedRegion, RegionType); itkGetMacro (RequestedRegion, RegionType); // Set/Get function for sigma parameter itkSetMacro (UseApproximateGradient, bool); itkGetMacro (UseApproximateGradient, bool); virtual void SetImage(const TInputImageType* _arg) { if (this->m_Image != _arg) { this->m_Image = _arg; this->Modified(); this->m_Initialized = false; } } void SetDynamicCostMap( std::map< int, int > &costMap) { this->m_CostMap = costMap; this->m_UseCostMap = true; this->m_MaxMapCosts = -1; this->Modified(); } void SetUseCostMap(bool useCostMap) { this->m_UseCostMap = useCostMap; } /** \brief Set the maximum of the dynamic cost map to save computation time. */ void SetCostMapMaximum(double max) { this->m_MaxMapCosts = max; } enum Constants{ MAPSCALEFACTOR = 10 }; /** \brief Returns the y value of gaussian with given offset and amplitude gaussian approximation f(x) = v(bin) * e^ ( -1/2 * (|x-k(bin)| / sigma)^2 ) \param x \param xOfGaussian - offset \param yOfGaussian - amplitude */ static double Gaussian(double x, double xOfGaussian, double yOfGaussian); protected: - virtual ~ShortestPathCostFunctionLiveWire() {}; + ShortestPathCostFunctionLiveWire(); + virtual ~ShortestPathCostFunctionLiveWire() {}; typename ImageType::Pointer m_GradientMagnImage; typename UnsignedCharImageType::Pointer m_ZeroCrossingsImage; typename FloatImageType::Pointer m_EdgeImage; typename VectorOutputImageType::Pointer m_GradientImage; double minCosts; - bool m_UseRepulsivePoint; + bool m_UseRepulsivePoints; - std::vector< itk::Index<3> > m_RepulsivePoints; + std::vector< IndexType > m_RepulsivePoints; typename Superclass::PixelType val; typename Superclass::PixelType startValue; typename Superclass::PixelType endValue; double m_GradientMax; RegionType m_RequestedRegion; bool m_UseApproximateGradient; bool m_Initialized; std::map< int, int > m_CostMap; bool m_UseCostMap; double m_MaxMapCosts; private: double SigmoidFunction(double I, double max, double min, double alpha, double beta); }; } // end namespace itk #ifndef ITK_MANUAL_INSTANTIATION #include "itkShortestPathCostFunctionLiveWire.txx" #endif #endif /* __itkShortestPathCostFunctionLiveWire_h */ diff --git a/Modules/GraphAlgorithms/itkShortestPathCostFunctionLiveWire.txx b/Modules/GraphAlgorithms/itkShortestPathCostFunctionLiveWire.txx index a9ad943a75..361facd306 100644 --- a/Modules/GraphAlgorithms/itkShortestPathCostFunctionLiveWire.txx +++ b/Modules/GraphAlgorithms/itkShortestPathCostFunctionLiveWire.txx @@ -1,441 +1,424 @@ /*=================================================================== 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. ===================================================================*/ #ifndef __itkShortestPathCostFunctionLiveWire_txx #define __itkShortestPathCostFunctionLiveWire_txx #include "itkShortestPathCostFunctionLiveWire.h" #include #include #include #include #include #include -#include #include -#include #include namespace itk { // Constructor template ShortestPathCostFunctionLiveWire ::ShortestPathCostFunctionLiveWire() { - m_UseRepulsivePoint = false; + m_UseRepulsivePoints = false; m_GradientMax = 0.0; m_Initialized = false; m_UseCostMap = false; m_MaxMapCosts = -1.0; } template void ShortestPathCostFunctionLiveWire - ::AddRepulsivePoint( itk::Index<3> c ) + ::AddRepulsivePoint( IndexType index ) { - m_RepulsivePoints.push_back(c); - m_UseRepulsivePoint = true; + m_RepulsivePoints.push_back(index); + m_UseRepulsivePoints = true; } - - template void ShortestPathCostFunctionLiveWire ::ClearRepulsivePoints() { m_RepulsivePoints.clear(); - m_UseRepulsivePoint = false; + m_UseRepulsivePoints = false; } template double ShortestPathCostFunctionLiveWire ::GetCost(IndexType p1 ,IndexType p2) { unsigned long xMAX = this->m_Image->GetLargestPossibleRegion().GetSize()[0]; unsigned long yMAX = this->m_Image->GetLargestPossibleRegion().GetSize()[1]; double gradientX, gradientY; gradientX = gradientY = 0.0; double gradientCost; double gradientMagnitude; - /* ++++++++++++++++++++ GradientMagnitude costs ++++++++++++++++++++++++++*/ - + // Gradient Magnitude costs gradientMagnitude = this->m_GradientMagnImage->GetPixel(p2); gradientX = m_GradientImage->GetPixel(p2)[0]; gradientY = m_GradientImage->GetPixel(p2)[1]; - - if(m_UseCostMap && !m_CostMap.empty()) { std::map< int, int >::iterator end = m_CostMap.end(); std::map< int, int >::iterator last = --(m_CostMap.end()); //current position std::map< int, int >::iterator x; //std::map< int, int >::key_type keyOfX = static_cast::key_type>(gradientMagnitude * 1000); int keyOfX = static_cast(gradientMagnitude /* ShortestPathCostFunctionLiveWire::MAPSCALEFACTOR*/); x = m_CostMap.find( keyOfX ); std::map< int, int >::iterator left2; std::map< int, int >::iterator left1; std::map< int, int >::iterator right1; std::map< int, int >::iterator right2; if( x == end ) {//x can also be == end if the key is not in the map but between two other keys //search next key within map from x upwards right1 = m_CostMap.lower_bound( keyOfX ); } else { right1 = x; } if(right1 == end || right1 == last ) { right2 = end; } else//( right1 != (end-1) ) { std::map< int, int >::iterator temp = right1; right2 = ++right1;//rght1 + 1 right1 = temp; } - if( right1 == m_CostMap.begin() ) { left1 = end; left2 = end; } else if( right1 == (++(m_CostMap.begin())) ) { std::map< int, int >::iterator temp = right1; left1 = --right1;//rght1 - 1 right1 = temp; left2 = end; } else { std::map< int, int >::iterator temp = right1; left1 = --right1;//rght1 - 1 left2 = --right1;//rght1 - 2 right1 = temp; } double partRight1, partRight2, partLeft1, partLeft2; partRight1 = partRight2 = partLeft1 = partLeft2 = 0.0; - - /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + /* f(x) = v(bin) * e^ ( -1/2 * (|x-k(bin)| / sigma)^2 ) gaussian approximation where v(bin) is the value in the map k(bin) is the key */ + if( left2 != end ) { partLeft2 = ShortestPathCostFunctionLiveWire::Gaussian(keyOfX, left2->first, left2->second); } if( left1 != end ) { partLeft1 = ShortestPathCostFunctionLiveWire::Gaussian(keyOfX, left1->first, left1->second); } if( right1 != end ) { partRight1 = ShortestPathCostFunctionLiveWire::Gaussian(keyOfX, right1->first, right1->second); } if( right2 != end ) { partRight2 = ShortestPathCostFunctionLiveWire::Gaussian(keyOfX, right2->first, right2->second); } - /*----------------------------------------------------------------------------*/ - if( m_MaxMapCosts > 0.0 ) { gradientCost = 1.0 - ( (partRight1 + partRight2 + partLeft1 + partLeft2) / m_MaxMapCosts ); } else {//use linear mapping gradientCost = 1.0 - (gradientMagnitude / m_GradientMax); } } else {//use linear mapping //value between 0 (good) and 1 (bad) gradientCost = 1.0 - (gradientMagnitude / m_GradientMax); - } - /* -----------------------------------------------------------------------------*/ - - - /* ++++++++++++++++++++ Laplacian zero crossing costs ++++++++++++++++++++++++++*/ + // Laplacian zero crossing costs // f(p) = 0; if I(p)=0 // or 1; if I(p)!=0 double laplacianCost; typename Superclass::PixelType laplaceImageValue; - laplaceImageValue = m_EdgeImage->GetPixel(p2); if(laplaceImageValue < 0 || laplaceImageValue > 0) { laplacianCost = 1.0; } else { laplacianCost = 0.0; } - /* -----------------------------------------------------------------------------*/ - - - - /* ++++++++++++++++++++ Gradient direction costs ++++++++++++++++++++++++++*/ + // Gradient direction costs //vector q-p i.e. p2-p1 double vQP[2]; vQP[0] = p2[0] - p1[0]; vQP[1] = p2[1] - p1[1]; //------- //vector p-q i.e. p1-p2 double vPQ[2]; vPQ[0] = p1[0] - p2[0]; vPQ[1] = p1[1] - p2[1]; //------- // gradient vector at p1 double nGradientAtP1[2]; nGradientAtP1[0] = gradientX;//previously computed for gradient magnitude nGradientAtP1[1] = gradientY; //gradient direction unit vector of p1 nGradientAtP1[0] /= gradientMagnitude; nGradientAtP1[1] /= gradientMagnitude; //------- // gradient vector at p1 double nGradientAtP2[2]; - nGradientAtP2[0] = m_GradientImage->GetPixel(p2)[0]; nGradientAtP2[1] = m_GradientImage->GetPixel(p2)[1]; nGradientAtP2[0] /= m_GradientMagnImage->GetPixel(p2); nGradientAtP2[1] /= m_GradientMagnImage->GetPixel(p2); double scalarProduct = (nGradientAtP1[0] * nGradientAtP2[0]) + (nGradientAtP1[1] * nGradientAtP2[1]); if( abs(scalarProduct) >= 1.0) { //this should probably not happen; make sure the input for acos is valid scalarProduct = 0.999999999; } double gradientDirectionCost = acos( scalarProduct ) / 3.14159265; - /*------------------------------------------------------------------------*/ - - - - /*+++++++++++++++++++++ local component costs +++++++++++++++++++++++++++*/ - /*weights*/ + // local component costs + // weights double w1; double w2; double w3; double costs = 0.0; if (this->m_UseCostMap){ w1 = 0.43; w2= 0.43; w3 = 0.14; }else{ w1 = 0.10; w2= 0.85; w3 = 0.05; } costs = w1 * laplacianCost + w2 * gradientCost + w3 * gradientDirectionCost; - - //scale by euclidian distance - double costScale; + //scale by euclidian distance + double costScale; if( p1[0] == p2[0] || p1[1] == p2[1]) { //horizontal or vertical neighbor costScale = 1.0; } else { //diagonal neighbor costScale = sqrt(2.0); } + // use repulsive points + + if (m_UseRepulsivePoints) + { + std::vector< IndexType >::const_iterator iter = m_RepulsivePoints.begin(); + + for (;iter != m_RepulsivePoints.end(); iter++) + { + if ( p1 == (*iter) || p2 == (*iter) ) + costs = 1.0; + } + } + costs *= costScale; return costs; } template double ShortestPathCostFunctionLiveWire ::GetMinCost() { return minCosts; } template void ShortestPathCostFunctionLiveWire ::Initialize() { if(!m_Initialized) { typedef itk::CastImageFilter< TInputImageType, FloatImageType > CastFilterType; typename CastFilterType::Pointer castFilter = CastFilterType::New(); castFilter->SetInput(this->m_Image); // init gradient magnitude image typedef itk::GradientMagnitudeImageFilter< FloatImageType, FloatImageType> GradientMagnitudeFilterType; typename GradientMagnitudeFilterType::Pointer gradientFilter = GradientMagnitudeFilterType::New(); gradientFilter->SetInput(castFilter->GetOutput()); //gradientFilter->SetNumberOfThreads(4); //gradientFilter->GetOutput()->SetRequestedRegion(m_RequestedRegion); gradientFilter->Update(); m_GradientMagnImage = gradientFilter->GetOutput(); typedef itk::StatisticsImageFilter StatisticsImageFilterType; typename StatisticsImageFilterType::Pointer statisticsImageFilter = StatisticsImageFilterType::New(); statisticsImageFilter->SetInput(this->m_GradientMagnImage); statisticsImageFilter->Update(); m_GradientMax = statisticsImageFilter->GetMaximum(); - - - //Filter class is instantiated - /*typedef itk::GradientRecursiveGaussianImageFilter GradientFilterType;*/ - typedef itk::GradientImageFilter< FloatImageType > GradientFilterType; typename GradientFilterType::Pointer filter = GradientFilterType::New(); //sigma is specified in millimeters //filter->SetSigma( 1.5 ); filter->SetInput(castFilter->GetOutput()); filter->Update(); m_GradientImage = filter->GetOutput(); // init zero crossings //typedef itk::ZeroCrossingImageFilter< TInputImageType, UnsignedCharImageType > ZeroCrossingImageFilterType; //ZeroCrossingImageFilterType::Pointer zeroCrossingImageFilter = ZeroCrossingImageFilterType::New(); //zeroCrossingImageFilter->SetInput(this->m_Image); //zeroCrossingImageFilter->SetBackgroundValue(1); //zeroCrossingImageFilter->SetForegroundValue(0); //zeroCrossingImageFilter->SetNumberOfThreads(4); //zeroCrossingImageFilter->Update(); //m_EdgeImage = zeroCrossingImageFilter->GetOutput(); //cast image to float to apply canny edge dection filter /*typedef itk::CastImageFilter< TInputImageType, FloatImageType > CastFilterType; CastFilterType::Pointer castFilter = CastFilterType::New(); castFilter->SetInput(this->m_Image);*/ //typedef itk::LaplacianImageFilter filterType; //filterType::Pointer laplacianFilter = filterType::New(); //laplacianFilter->SetInput( castFilter->GetOutput() ); // NOTE: input image type must be double or float //laplacianFilter->Update(); //m_EdgeImage = laplacianFilter->GetOutput(); //init canny edge detection typedef itk::CannyEdgeDetectionImageFilter CannyEdgeDetectionImageFilterType; typename CannyEdgeDetectionImageFilterType::Pointer cannyEdgeDetectionfilter = CannyEdgeDetectionImageFilterType::New(); cannyEdgeDetectionfilter->SetInput(castFilter->GetOutput()); cannyEdgeDetectionfilter->SetUpperThreshold(30); cannyEdgeDetectionfilter->SetLowerThreshold(15); cannyEdgeDetectionfilter->SetVariance(4); cannyEdgeDetectionfilter->SetMaximumError(.01f); cannyEdgeDetectionfilter->Update(); m_EdgeImage = cannyEdgeDetectionfilter->GetOutput(); // set minCosts minCosts = 0.0; // The lower, the more thouroughly! 0 = dijkstra. If estimate costs are lower than actual costs everything is fine. If estimation is higher than actual costs, you might not get the shortest but a different path. m_Initialized = true; } // check start/end point value startValue= this->m_Image->GetPixel(this->m_StartIndex); endValue= this->m_Image->GetPixel(this->m_EndIndex); } template double ShortestPathCostFunctionLiveWire::SigmoidFunction(double I, double max, double min, double alpha, double beta) { // Using the SIgmoid formula from ITK Software Guide 6.3.2 Non Linear Mappings double Exponent = -1 * ((I - beta) / alpha); double Factor = 1 / (1 + exp(Exponent)); double newI = (max - min) * Factor + min; return newI; } template double ShortestPathCostFunctionLiveWire::Gaussian(double x, double xOfGaussian, double yOfGaussian) { return yOfGaussian * exp( -0.5 * pow( (x - xOfGaussian), 2) ); } } // end namespace itk #endif // __itkShortestPathCostFunctionLiveWire_txx