1 /*===================================================================
3 The Medical Imaging Interaction Toolkit (MITK)
5 Copyright (c) German Cancer Research Center,
6 Division of Medical and Biological Informatics.
9 This software is distributed WITHOUT ANY WARRANTY; without
10 even the implied warranty of MERCHANTABILITY or FITNESS FOR
13 See LICENSE.txt or http://www.mitk.org for details.
15 ===================================================================*/
17 #ifndef __itkShortestPathCostFunctionLiveWire_txx
18 #define __itkShortestPathCostFunctionLiveWire_txx
20 #include "itkShortestPathCostFunctionLiveWire.h"
24 #include <itkCannyEdgeDetectionImageFilter.h>
25 #include <itkCastImageFilter.h>
26 #include <itkGradientImageFilter.h>
27 #include <itkGradientMagnitudeImageFilter.h>
28 #include <itkLaplacianImageFilter.h>
29 #include <itkStatisticsImageFilter.h>
30 #include <itkZeroCrossingImageFilter.h>
35 template <class TInputImageType>
36 ShortestPathCostFunctionLiveWire<TInputImageType>::ShortestPathCostFunctionLiveWire()
38 m_UseRepulsivePoints = false;
40 m_Initialized = false;
45 template <class TInputImageType>
46 void ShortestPathCostFunctionLiveWire<TInputImageType>::AddRepulsivePoint(const IndexType &index)
48 this->m_MaskImage->SetPixel(index, 255);
49 m_UseRepulsivePoints = true;
52 template <class TInputImageType>
53 void ShortestPathCostFunctionLiveWire<TInputImageType>::RemoveRepulsivePoint(const IndexType &index)
55 this->m_MaskImage->SetPixel(index, 0);
58 template <class TInputImageType>
59 void ShortestPathCostFunctionLiveWire<TInputImageType>::SetImage(const TInputImageType *_arg)
61 if (this->m_Image != _arg)
65 // initialize mask image
66 this->m_MaskImage = UnsignedCharImageType::New();
67 this->m_MaskImage->SetRegions(this->m_Image->GetLargestPossibleRegion());
68 this->m_MaskImage->SetOrigin(this->m_Image->GetOrigin());
69 this->m_MaskImage->SetSpacing(this->m_Image->GetSpacing());
70 this->m_MaskImage->SetDirection(this->m_Image->GetDirection());
71 this->m_MaskImage->Allocate();
72 this->m_MaskImage->FillBuffer(0);
75 this->m_Initialized = false;
79 template <class TInputImageType>
80 void ShortestPathCostFunctionLiveWire<TInputImageType>::ClearRepulsivePoints()
82 m_UseRepulsivePoints = false;
83 this->m_MaskImage->FillBuffer(0);
86 template <class TInputImageType>
87 double ShortestPathCostFunctionLiveWire<TInputImageType>::GetCost(IndexType p1, IndexType p2)
89 // local component costs
96 // if we are on the mask, return asap
97 if (m_UseRepulsivePoints)
99 if ((this->m_MaskImage->GetPixel(p1) != 0) || (this->m_MaskImage->GetPixel(p2) != 0))
103 double gradientX, gradientY;
104 gradientX = gradientY = 0.0;
108 double gradientMagnitude;
110 // Gradient Magnitude costs
111 gradientMagnitude = this->m_GradientMagnitudeImage->GetPixel(p2);
112 gradientX = m_GradientImage->GetPixel(p2)[0];
113 gradientY = m_GradientImage->GetPixel(p2)[1];
115 if (m_UseCostMap && !m_CostMap.empty())
117 std::map<int, int>::iterator end = m_CostMap.end();
118 std::map<int, int>::iterator last = --(m_CostMap.end());
121 std::map<int, int>::iterator x;
122 // std::map< int, int >::key_type keyOfX = static_cast<std::map< int, int >::key_type>(gradientMagnitude * 1000);
123 int keyOfX = static_cast<int>(gradientMagnitude /* ShortestPathCostFunctionLiveWire::MAPSCALEFACTOR*/);
124 x = m_CostMap.find(keyOfX);
126 std::map<int, int>::iterator left2;
127 std::map<int, int>::iterator left1;
128 std::map<int, int>::iterator right1;
129 std::map<int, int>::iterator right2;
132 { // x can also be == end if the key is not in the map but between two other keys
133 // search next key within map from x upwards
134 right1 = m_CostMap.lower_bound(keyOfX);
141 if (right1 == end || right1 == last)
145 else //( right1 != (end-1) )
148 right2 = ++right1; // rght1 + 1
152 if (right1 == m_CostMap.begin())
157 else if (right1 == (++(m_CostMap.begin())))
160 left1 = --right1; // rght1 - 1
167 left1 = --right1; // rght1 - 1
168 left2 = --right1; // rght1 - 2
172 double partRight1, partRight2, partLeft1, partLeft2;
173 partRight1 = partRight2 = partLeft1 = partLeft2 = 0.0;
176 f(x) = v(bin) * e^ ( -1/2 * (|x-k(bin)| / sigma)^2 )
178 gaussian approximation
181 v(bin) is the value in the map
187 partLeft2 = ShortestPathCostFunctionLiveWire<TInputImageType>::Gaussian(keyOfX, left2->first, left2->second);
192 partLeft1 = ShortestPathCostFunctionLiveWire<TInputImageType>::Gaussian(keyOfX, left1->first, left1->second);
197 partRight1 = ShortestPathCostFunctionLiveWire<TInputImageType>::Gaussian(keyOfX, right1->first, right1->second);
202 partRight2 = ShortestPathCostFunctionLiveWire<TInputImageType>::Gaussian(keyOfX, right2->first, right2->second);
205 if (m_MaxMapCosts > 0.0)
207 gradientCost = 1.0 - ((partRight1 + partRight2 + partLeft1 + partLeft2) / m_MaxMapCosts);
210 { // use linear mapping
211 gradientCost = 1.0 - (gradientMagnitude / m_GradientMax);
215 { // use linear mapping
216 // value between 0 (good) and 1 (bad)
217 gradientCost = 1.0 - (gradientMagnitude / m_GradientMax);
220 // Laplacian zero crossing costs
221 // f(p) = 0; if I(p)=0
223 double laplacianCost;
224 typename Superclass::PixelType laplaceImageValue;
226 laplaceImageValue = m_EdgeImage->GetPixel(p2);
228 if (laplaceImageValue < 0 || laplaceImageValue > 0)
237 // gradient vector at p1
238 double nGradientAtP1[2];
239 nGradientAtP1[0] = gradientX; // previously computed for gradient magnitude
240 nGradientAtP1[1] = gradientY;
242 // gradient direction unit vector of p1
243 nGradientAtP1[0] /= gradientMagnitude;
244 nGradientAtP1[1] /= gradientMagnitude;
247 // gradient vector at p1
248 double nGradientAtP2[2];
250 nGradientAtP2[0] = m_GradientImage->GetPixel(p2)[0];
251 nGradientAtP2[1] = m_GradientImage->GetPixel(p2)[1];
253 nGradientAtP2[0] /= m_GradientMagnitudeImage->GetPixel(p2);
254 nGradientAtP2[1] /= m_GradientMagnitudeImage->GetPixel(p2);
256 double scalarProduct = (nGradientAtP1[0] * nGradientAtP2[0]) + (nGradientAtP1[1] * nGradientAtP2[1]);
257 if (abs(scalarProduct) >= 1.0)
259 // this should probably not happen; make sure the input for acos is valid
260 scalarProduct = 0.999999999;
263 double gradientDirectionCost = acos(scalarProduct) / 3.14159265;
265 if (this->m_UseCostMap)
277 costs = w1 * laplacianCost + w2 * gradientCost + w3 * gradientDirectionCost;
279 // scale by euclidian distance
281 if (p1[0] == p2[0] || p1[1] == p2[1])
283 // horizontal or vertical neighbor
289 costScale = sqrt(2.0);
297 template <class TInputImageType>
298 double ShortestPathCostFunctionLiveWire<TInputImageType>::GetMinCost()
303 template <class TInputImageType>
304 void ShortestPathCostFunctionLiveWire<TInputImageType>::Initialize()
308 typedef itk::CastImageFilter<TInputImageType, FloatImageType> CastFilterType;
309 typename CastFilterType::Pointer castFilter = CastFilterType::New();
310 castFilter->SetInput(this->m_Image);
312 // init gradient magnitude image
313 typedef itk::GradientMagnitudeImageFilter<FloatImageType, FloatImageType> GradientMagnitudeFilterType;
314 typename GradientMagnitudeFilterType::Pointer gradientFilter = GradientMagnitudeFilterType::New();
315 gradientFilter->SetInput(castFilter->GetOutput());
316 // gradientFilter->SetNumberOfThreads(4);
317 // gradientFilter->GetOutput()->SetRequestedRegion(m_RequestedRegion);
319 gradientFilter->Update();
320 this->m_GradientMagnitudeImage = gradientFilter->GetOutput();
322 typedef itk::StatisticsImageFilter<FloatImageType> StatisticsImageFilterType;
323 typename StatisticsImageFilterType::Pointer statisticsImageFilter = StatisticsImageFilterType::New();
324 statisticsImageFilter->SetInput(this->m_GradientMagnitudeImage);
325 statisticsImageFilter->Update();
327 m_GradientMax = statisticsImageFilter->GetMaximum();
329 typedef itk::GradientImageFilter<FloatImageType> GradientFilterType;
331 typename GradientFilterType::Pointer filter = GradientFilterType::New();
332 // sigma is specified in millimeters
333 // filter->SetSigma( 1.5 );
334 filter->SetInput(castFilter->GetOutput());
337 m_GradientImage = filter->GetOutput();
339 // init zero crossings
340 // typedef itk::ZeroCrossingImageFilter< TInputImageType, UnsignedCharImageType > ZeroCrossingImageFilterType;
341 // ZeroCrossingImageFilterType::Pointer zeroCrossingImageFilter = ZeroCrossingImageFilterType::New();
342 // zeroCrossingImageFilter->SetInput(this->m_Image);
343 // zeroCrossingImageFilter->SetBackgroundValue(1);
344 // zeroCrossingImageFilter->SetForegroundValue(0);
345 // zeroCrossingImageFilter->SetNumberOfThreads(4);
346 // zeroCrossingImageFilter->Update();
348 // m_EdgeImage = zeroCrossingImageFilter->GetOutput();
350 // cast image to float to apply canny edge dection filter
351 /*typedef itk::CastImageFilter< TInputImageType, FloatImageType > CastFilterType;
352 CastFilterType::Pointer castFilter = CastFilterType::New();
353 castFilter->SetInput(this->m_Image);*/
355 // typedef itk::LaplacianImageFilter<FloatImageType, FloatImageType > filterType;
356 // filterType::Pointer laplacianFilter = filterType::New();
357 // laplacianFilter->SetInput( castFilter->GetOutput() ); // NOTE: input image type must be double or float
358 // laplacianFilter->Update();
360 // m_EdgeImage = laplacianFilter->GetOutput();
362 // init canny edge detection
363 typedef itk::CannyEdgeDetectionImageFilter<FloatImageType, FloatImageType> CannyEdgeDetectionImageFilterType;
364 typename CannyEdgeDetectionImageFilterType::Pointer cannyEdgeDetectionfilter =
365 CannyEdgeDetectionImageFilterType::New();
366 cannyEdgeDetectionfilter->SetInput(castFilter->GetOutput());
367 cannyEdgeDetectionfilter->SetUpperThreshold(30);
368 cannyEdgeDetectionfilter->SetLowerThreshold(15);
369 cannyEdgeDetectionfilter->SetVariance(4);
370 cannyEdgeDetectionfilter->SetMaximumError(.01f);
372 cannyEdgeDetectionfilter->Update();
373 m_EdgeImage = cannyEdgeDetectionfilter->GetOutput();
376 minCosts = 0.0; // The lower, the more thouroughly! 0 = dijkstra. If estimate costs are lower than actual costs
377 // everything is fine. If estimation is higher than actual costs, you might not get the shortest
378 // but a different path.
380 m_Initialized = true;
383 // check start/end point value
384 startValue = this->m_Image->GetPixel(this->m_StartIndex);
385 endValue = this->m_Image->GetPixel(this->m_EndIndex);
388 template <class TInputImageType>
389 double ShortestPathCostFunctionLiveWire<TInputImageType>::SigmoidFunction(
390 double I, double max, double min, double alpha, double beta)
392 // Using the SIgmoid formula from ITK Software Guide 6.3.2 Non Linear Mappings
393 double Exponent = -1 * ((I - beta) / alpha);
394 double Factor = 1 / (1 + exp(Exponent));
395 double newI = (max - min) * Factor + min;
400 template <class TInputImageType>
401 double ShortestPathCostFunctionLiveWire<TInputImageType>::Gaussian(double x, double xOfGaussian, double yOfGaussian)
403 return yOfGaussian * exp(-0.5 * pow((x - xOfGaussian), 2));
406 } // end namespace itk
408 #endif // __itkShortestPathCostFunctionLiveWire_txx