Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
itkShortestPathCostFunctionLiveWire.txx
Go to the documentation of this file.
1 /*===================================================================
2 
3 The Medical Imaging Interaction Toolkit (MITK)
4 
5 Copyright (c) German Cancer Research Center,
6 Division of Medical and Biological Informatics.
7 All rights reserved.
8 
9 This software is distributed WITHOUT ANY WARRANTY; without
10 even the implied warranty of MERCHANTABILITY or FITNESS FOR
11 A PARTICULAR PURPOSE.
12 
13 See LICENSE.txt or http://www.mitk.org for details.
14 
15 ===================================================================*/
16 
17 #ifndef __itkShortestPathCostFunctionLiveWire_txx
18 #define __itkShortestPathCostFunctionLiveWire_txx
19 
20 #include "itkShortestPathCostFunctionLiveWire.h"
21 
22 #include <math.h>
23 
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>
31 
32 namespace itk
33 {
34  // Constructor
35  template <class TInputImageType>
36  ShortestPathCostFunctionLiveWire<TInputImageType>::ShortestPathCostFunctionLiveWire()
37  {
38  m_UseRepulsivePoints = false;
39  m_GradientMax = 0.0;
40  m_Initialized = false;
41  m_UseCostMap = false;
42  m_MaxMapCosts = -1.0;
43  }
44 
45  template <class TInputImageType>
46  void ShortestPathCostFunctionLiveWire<TInputImageType>::AddRepulsivePoint(const IndexType &index)
47  {
48  this->m_MaskImage->SetPixel(index, 255);
49  m_UseRepulsivePoints = true;
50  }
51 
52  template <class TInputImageType>
53  void ShortestPathCostFunctionLiveWire<TInputImageType>::RemoveRepulsivePoint(const IndexType &index)
54  {
55  this->m_MaskImage->SetPixel(index, 0);
56  }
57 
58  template <class TInputImageType>
59  void ShortestPathCostFunctionLiveWire<TInputImageType>::SetImage(const TInputImageType *_arg)
60  {
61  if (this->m_Image != _arg)
62  {
63  this->m_Image = _arg;
64 
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);
73 
74  this->Modified();
75  this->m_Initialized = false;
76  }
77  }
78 
79  template <class TInputImageType>
80  void ShortestPathCostFunctionLiveWire<TInputImageType>::ClearRepulsivePoints()
81  {
82  m_UseRepulsivePoints = false;
83  this->m_MaskImage->FillBuffer(0);
84  }
85 
86  template <class TInputImageType>
87  double ShortestPathCostFunctionLiveWire<TInputImageType>::GetCost(IndexType p1, IndexType p2)
88  {
89  // local component costs
90  // weights
91  double w1;
92  double w2;
93  double w3;
94  double costs = 0.0;
95 
96  // if we are on the mask, return asap
97  if (m_UseRepulsivePoints)
98  {
99  if ((this->m_MaskImage->GetPixel(p1) != 0) || (this->m_MaskImage->GetPixel(p2) != 0))
100  return 1000;
101  }
102 
103  double gradientX, gradientY;
104  gradientX = gradientY = 0.0;
105 
106  double gradientCost;
107 
108  double gradientMagnitude;
109 
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];
114 
115  if (m_UseCostMap && !m_CostMap.empty())
116  {
117  std::map<int, int>::iterator end = m_CostMap.end();
118  std::map<int, int>::iterator last = --(m_CostMap.end());
119 
120  // current position
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);
125 
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;
130 
131  if (x == end)
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);
135  }
136  else
137  {
138  right1 = x;
139  }
140 
141  if (right1 == end || right1 == last)
142  {
143  right2 = end;
144  }
145  else //( right1 != (end-1) )
146  {
147  auto temp = right1;
148  right2 = ++right1; // rght1 + 1
149  right1 = temp;
150  }
151 
152  if (right1 == m_CostMap.begin())
153  {
154  left1 = end;
155  left2 = end;
156  }
157  else if (right1 == (++(m_CostMap.begin())))
158  {
159  auto temp = right1;
160  left1 = --right1; // rght1 - 1
161  right1 = temp;
162  left2 = end;
163  }
164  else
165  {
166  auto temp = right1;
167  left1 = --right1; // rght1 - 1
168  left2 = --right1; // rght1 - 2
169  right1 = temp;
170  }
171 
172  double partRight1, partRight2, partLeft1, partLeft2;
173  partRight1 = partRight2 = partLeft1 = partLeft2 = 0.0;
174 
175  /*
176  f(x) = v(bin) * e^ ( -1/2 * (|x-k(bin)| / sigma)^2 )
177 
178  gaussian approximation
179 
180  where
181  v(bin) is the value in the map
182  k(bin) is the key
183  */
184 
185  if (left2 != end)
186  {
187  partLeft2 = ShortestPathCostFunctionLiveWire<TInputImageType>::Gaussian(keyOfX, left2->first, left2->second);
188  }
189 
190  if (left1 != end)
191  {
192  partLeft1 = ShortestPathCostFunctionLiveWire<TInputImageType>::Gaussian(keyOfX, left1->first, left1->second);
193  }
194 
195  if (right1 != end)
196  {
197  partRight1 = ShortestPathCostFunctionLiveWire<TInputImageType>::Gaussian(keyOfX, right1->first, right1->second);
198  }
199 
200  if (right2 != end)
201  {
202  partRight2 = ShortestPathCostFunctionLiveWire<TInputImageType>::Gaussian(keyOfX, right2->first, right2->second);
203  }
204 
205  if (m_MaxMapCosts > 0.0)
206  {
207  gradientCost = 1.0 - ((partRight1 + partRight2 + partLeft1 + partLeft2) / m_MaxMapCosts);
208  }
209  else
210  { // use linear mapping
211  gradientCost = 1.0 - (gradientMagnitude / m_GradientMax);
212  }
213  }
214  else
215  { // use linear mapping
216  // value between 0 (good) and 1 (bad)
217  gradientCost = 1.0 - (gradientMagnitude / m_GradientMax);
218  }
219 
220  // Laplacian zero crossing costs
221  // f(p) = 0; if I(p)=0
222  // or 1; if I(p)!=0
223  double laplacianCost;
224  typename Superclass::PixelType laplaceImageValue;
225 
226  laplaceImageValue = m_EdgeImage->GetPixel(p2);
227 
228  if (laplaceImageValue < 0 || laplaceImageValue > 0)
229  {
230  laplacianCost = 1.0;
231  }
232  else
233  {
234  laplacianCost = 0.0;
235  }
236 
237  // gradient vector at p1
238  double nGradientAtP1[2];
239  nGradientAtP1[0] = gradientX; // previously computed for gradient magnitude
240  nGradientAtP1[1] = gradientY;
241 
242  // gradient direction unit vector of p1
243  nGradientAtP1[0] /= gradientMagnitude;
244  nGradientAtP1[1] /= gradientMagnitude;
245  //-------
246 
247  // gradient vector at p1
248  double nGradientAtP2[2];
249 
250  nGradientAtP2[0] = m_GradientImage->GetPixel(p2)[0];
251  nGradientAtP2[1] = m_GradientImage->GetPixel(p2)[1];
252 
253  nGradientAtP2[0] /= m_GradientMagnitudeImage->GetPixel(p2);
254  nGradientAtP2[1] /= m_GradientMagnitudeImage->GetPixel(p2);
255 
256  double scalarProduct = (nGradientAtP1[0] * nGradientAtP2[0]) + (nGradientAtP1[1] * nGradientAtP2[1]);
257  if (abs(scalarProduct) >= 1.0)
258  {
259  // this should probably not happen; make sure the input for acos is valid
260  scalarProduct = 0.999999999;
261  }
262 
263  double gradientDirectionCost = acos(scalarProduct) / 3.14159265;
264 
265  if (this->m_UseCostMap)
266  {
267  w1 = 0.43;
268  w2 = 0.43;
269  w3 = 0.14;
270  }
271  else
272  {
273  w1 = 0.10;
274  w2 = 0.85;
275  w3 = 0.05;
276  }
277  costs = w1 * laplacianCost + w2 * gradientCost + w3 * gradientDirectionCost;
278 
279  // scale by euclidian distance
280  double costScale;
281  if (p1[0] == p2[0] || p1[1] == p2[1])
282  {
283  // horizontal or vertical neighbor
284  costScale = 1.0;
285  }
286  else
287  {
288  // diagonal neighbor
289  costScale = sqrt(2.0);
290  }
291 
292  costs *= costScale;
293 
294  return costs;
295  }
296 
297  template <class TInputImageType>
298  double ShortestPathCostFunctionLiveWire<TInputImageType>::GetMinCost()
299  {
300  return minCosts;
301  }
302 
303  template <class TInputImageType>
304  void ShortestPathCostFunctionLiveWire<TInputImageType>::Initialize()
305  {
306  if (!m_Initialized)
307  {
308  typedef itk::CastImageFilter<TInputImageType, FloatImageType> CastFilterType;
309  typename CastFilterType::Pointer castFilter = CastFilterType::New();
310  castFilter->SetInput(this->m_Image);
311 
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);
318 
319  gradientFilter->Update();
320  this->m_GradientMagnitudeImage = gradientFilter->GetOutput();
321 
322  typedef itk::StatisticsImageFilter<FloatImageType> StatisticsImageFilterType;
323  typename StatisticsImageFilterType::Pointer statisticsImageFilter = StatisticsImageFilterType::New();
324  statisticsImageFilter->SetInput(this->m_GradientMagnitudeImage);
325  statisticsImageFilter->Update();
326 
327  m_GradientMax = statisticsImageFilter->GetMaximum();
328 
329  typedef itk::GradientImageFilter<FloatImageType> GradientFilterType;
330 
331  typename GradientFilterType::Pointer filter = GradientFilterType::New();
332  // sigma is specified in millimeters
333  // filter->SetSigma( 1.5 );
334  filter->SetInput(castFilter->GetOutput());
335  filter->Update();
336 
337  m_GradientImage = filter->GetOutput();
338 
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();
347 
348  // m_EdgeImage = zeroCrossingImageFilter->GetOutput();
349 
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);*/
354 
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();
359 
360  // m_EdgeImage = laplacianFilter->GetOutput();
361 
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);
371 
372  cannyEdgeDetectionfilter->Update();
373  m_EdgeImage = cannyEdgeDetectionfilter->GetOutput();
374 
375  // set minCosts
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.
379 
380  m_Initialized = true;
381  }
382 
383  // check start/end point value
384  startValue = this->m_Image->GetPixel(this->m_StartIndex);
385  endValue = this->m_Image->GetPixel(this->m_EndIndex);
386  }
387 
388  template <class TInputImageType>
389  double ShortestPathCostFunctionLiveWire<TInputImageType>::SigmoidFunction(
390  double I, double max, double min, double alpha, double beta)
391  {
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;
396 
397  return newI;
398  }
399 
400  template <class TInputImageType>
401  double ShortestPathCostFunctionLiveWire<TInputImageType>::Gaussian(double x, double xOfGaussian, double yOfGaussian)
402  {
403  return yOfGaussian * exp(-0.5 * pow((x - xOfGaussian), 2));
404  }
405 
406 } // end namespace itk
407 
408 #endif // __itkShortestPathCostFunctionLiveWire_txx