diff --git a/Modules/GraphAlgorithms/itkShortestPathCostFunctionLiveWire.h b/Modules/GraphAlgorithms/itkShortestPathCostFunctionLiveWire.h index caee20bc4f..6cb5032ceb 100644 --- a/Modules/GraphAlgorithms/itkShortestPathCostFunctionLiveWire.h +++ b/Modules/GraphAlgorithms/itkShortestPathCostFunctionLiveWire.h @@ -1,194 +1,197 @@ /*=================================================================== 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_h #define __itkShortestPathCostFunctionLiveWire_h #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 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 than higher values of gradient magnitudes. To compute the costs of the gradient magnitude dynamically 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 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; /** \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 Add repulsive point*/ virtual void AddRepulsivePoint( IndexType index ); + virtual int GetNumberOfRepulsivePoints() + { return this->m_RepulsivePoints.size(); } + /** \brief Clear repulsive path*/ virtual void ClearRepulsivePoints( ); 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: 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_UseRepulsivePoints; 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 153aebb45a..b37d1f5e4a 100644 --- a/Modules/GraphAlgorithms/itkShortestPathCostFunctionLiveWire.txx +++ b/Modules/GraphAlgorithms/itkShortestPathCostFunctionLiveWire.txx @@ -1,425 +1,425 @@ /*=================================================================== 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 namespace itk { // Constructor template ShortestPathCostFunctionLiveWire ::ShortestPathCostFunctionLiveWire() { m_UseRepulsivePoints = false; m_GradientMax = 0.0; m_Initialized = false; m_UseCostMap = false; m_MaxMapCosts = -1.0; } template void ShortestPathCostFunctionLiveWire ::AddRepulsivePoint( IndexType index ) { m_RepulsivePoints.push_back(index); m_UseRepulsivePoints = true; } template void ShortestPathCostFunctionLiveWire ::ClearRepulsivePoints() { m_RepulsivePoints.clear(); m_UseRepulsivePoints = false; } template double ShortestPathCostFunctionLiveWire ::GetCost(IndexType p1 ,IndexType p2) { // local component costs // weights double w1; double w2; double w3; double costs = 0.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 = 1000; return costs; } } } 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; // 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 // 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 //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; 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; if( p1[0] == p2[0] || p1[1] == p2[1]) { //horizontal or vertical neighbor costScale = 1.0; } else { //diagonal neighbor costScale = sqrt(2.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(); 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