Medical Imaging Interaction Toolkit  2018.4.99-12ad79a3
Medical Imaging Interaction Toolkit
itkAdaptiveThresholdIterator.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 (DKFZ)
6 All rights reserved.
7 
8 Use of this source code is governed by a 3-clause BSD license that can be
9 found in the LICENSE file.
10 
11 ============================================================================*/
12 #ifndef _itkAdaptiveThresholdIterator_txx
13 #define _itkAdaptiveThresholdIterator_txx
14 
15 #include "itkAdaptiveThresholdIterator.h"
16 #include "mitkProgressBar.h"
17 #include <cmath>
18 
19 namespace itk
20 {
21  template <class TImage, class TFunction>
22  AdaptiveThresholdIterator<TImage, TFunction>::AdaptiveThresholdIterator(ImageType *imagePtr,
23  FunctionType *fnPtr,
24  IndexType startIndex)
25  : m_FineDetectionMode(false), m_DetectionStop(false)
26  {
27  this->m_OutputImage = imagePtr;
28  m_Function = fnPtr;
29  m_StartIndices.push_back(startIndex);
30 
31  // Set up the temporary image
32  this->InitializeIterator();
33  }
34 
35  template <class TImage, class TFunction>
36  AdaptiveThresholdIterator<TImage, TFunction>::AdaptiveThresholdIterator(ImageType *imagePtr,
37  FunctionType *fnPtr,
38  std::vector<IndexType> &startIndex)
39  : m_FineDetectionMode(false), m_DetectionStop(false)
40  {
41  this->m_OutputImage = imagePtr;
42  m_Function = fnPtr;
43  unsigned int i;
44  for (i = 0; i < startIndex.size(); i++)
45  {
46  m_StartIndices.push_back(startIndex[i]);
47  }
48 
49  // Set up the temporary image
50  this->InitializeIterator();
51  }
52 
53  template <class TImage, class TFunction>
54  AdaptiveThresholdIterator<TImage, TFunction>::AdaptiveThresholdIterator(ImageType *imagePtr, FunctionType *fnPtr)
55  : m_FineDetectionMode(false), m_DetectionStop(false)
56  {
57  this->m_OutputImage = imagePtr; // here we store the image, we have to wite the result to
58  m_Function = fnPtr;
59 
60  // Set up the temporary image
61  this->InitializeIterator();
62  }
63 
64  template <class TImage, class TFunction>
65  void AdaptiveThresholdIterator<TImage, TFunction>::InitializeIterator()
66  {
67  // Get the origin and spacing from the image in simple arrays
68  m_ImageOrigin = this->m_OutputImage->GetOrigin();
69  m_ImageSpacing = this->m_OutputImage->GetSpacing();
70  m_ImageRegion = this->m_OutputImage->GetBufferedRegion();
71 
72  this->InitRegionGrowingState();
73 
74  m_VoxelCounter = 0;
75  m_LastVoxelNumber = 0;
76  m_CurrentLeakageRatio = 0;
77  m_DetectedLeakagePoint = 0;
78  }
79 
80  template <class TImage, class TFunction>
81  void AdaptiveThresholdIterator<TImage, TFunction>::SetExpansionDirection(bool upwards)
82  {
83  m_UpwardsExpansion = upwards;
84  }
85 
86  template <class TImage, class TFunction>
87  void AdaptiveThresholdIterator<TImage, TFunction>::IncrementRegionGrowingState()
88  {
89  // make the progressbar go one step further
90  if (!m_FineDetectionMode)
91  mitk::ProgressBar::GetInstance()->Progress();
92 
93  // updating the thresholds
94  if (m_UpwardsExpansion)
95  {
96  this->ExpandThresholdUpwards();
97  }
98  else
99  {
100  this->ExpandThresholdDownwards();
101  }
102 
103  // leakage-detection
104  int criticalValue = 2000; // calculate a bit more "intelligent"
105 
106  if (!m_FineDetectionMode)
107  {
108  int diff = m_VoxelCounter - m_LastVoxelNumber;
109  if (diff > m_CurrentLeakageRatio)
110  {
111  m_CurrentLeakageRatio = diff;
112  m_DetectedLeakagePoint = m_RegionGrowingState;
113  }
114  m_LastVoxelNumber = m_VoxelCounter;
115  m_VoxelCounter = 0;
116  }
117  else // fine leakage detection
118  {
119  // counting voxels over interations; if above a critical value (to be extended) then set this to leakage
120  int diff = m_VoxelCounter - m_LastVoxelNumber;
121  // std::cout<<"diff is: "<<diff<<"\n";
122  if (diff <= criticalValue && (!m_DetectionStop))
123  {
124  // m_CurrentLeakageRatio = diff;
125  m_DetectedLeakagePoint = m_RegionGrowingState + 1; // TODO check why this is needed
126  // std::cout<<"setting CurrentLeakageRatio to: "<<diff<<" and leakagePoint to: "<<m_DetectedLeakagePoint<<"\n";
127  }
128  else
129  {
130  m_DetectionStop = true;
131  // std::cout<<"\n\n[ITERATOR] detection stop!!!\n";
132  }
133  m_LastVoxelNumber = m_VoxelCounter;
134  m_VoxelCounter = 0;
135  }
136  // increment the counter
137  m_RegionGrowingState++;
138  }
139 
140  template <class TImage, class TFunction>
141  void AdaptiveThresholdIterator<TImage, TFunction>::ExpandThresholdUpwards()
142  {
143  int upper = (int)m_Function->GetUpper();
144  upper++;
145  m_Function->ThresholdBetween(m_MinTH, upper);
146  }
147 
148  template <class TImage, class TFunction>
149  void AdaptiveThresholdIterator<TImage, TFunction>::ExpandThresholdDownwards()
150  {
151  int lower = (int)m_Function->GetLower();
152  lower--;
153  m_Function->ThresholdBetween(lower, m_MaxTH);
154  }
155 
156  template <class TImage, class TFunction>
157  void AdaptiveThresholdIterator<TImage, TFunction>::InitRegionGrowingState()
158  {
159  this->m_RegionGrowingState = 1;
160  }
161 
162  template <class TImage, class TFunction>
163  int AdaptiveThresholdIterator<TImage, TFunction>::EstimateDistance(IndexType tempIndex)
164  {
165  PixelType value = this->m_Function->GetInputImage()->GetPixel(tempIndex);
166  PixelType minPixel = PixelType(m_MinTH);
167  PixelType maxPixel = PixelType(m_MaxTH);
168 
169  if (value > maxPixel || value < minPixel)
170  {
171  return 0;
172  }
173  if (m_UpwardsExpansion)
174  {
175  return (int)(value - m_SeedPointValue);
176  }
177  else
178  {
179  return (int)(m_SeedPointValue - value);
180  }
181  }
182 
183  template <class TImage, class TFunction>
184  void AdaptiveThresholdIterator<TImage, TFunction>::SetMinTH(int min)
185  {
186  m_MinTH = min;
187  }
188 
189  template <class TImage, class TFunction>
190  void AdaptiveThresholdIterator<TImage, TFunction>::SetMaxTH(int max)
191  {
192  m_MaxTH = max;
193  }
194 
195  template <class TImage, class TFunction>
196  int AdaptiveThresholdIterator<TImage, TFunction>::GetSeedPointValue()
197  {
198  return this->m_SeedPointValue;
199  }
200 
201  template <class TImage, class TFunction>
202  void AdaptiveThresholdIterator<TImage, TFunction>::GoToBegin()
203  {
204  m_QueueMap.clear();
205 
206  m_SeedPointValue = 0;
207 
208  IndexType seedIndex = m_StartIndices[0];
209 
210  bool doAverage = false; // enable or disable manually!
211  if (doAverage)
212  {
213  // loops for creating the sum of the N27-neighborhood around the seedpoint
214  for (int i = -1; i <= 1; i++)
215  {
216  for (int j = -1; j <= 1; j++)
217  {
218  for (int k = -1; k <= 1; k++)
219  {
220  seedIndex[0] = seedIndex[0] + i;
221  seedIndex[1] = seedIndex[1] + j;
222  seedIndex[2] = seedIndex[2] + k;
223 
224  m_SeedPointValue += (int)m_Function->GetInputImage()->GetPixel(seedIndex);
225  }
226  }
227  }
228 
229  // value of seedpoint computed from mean of N27-neighborhood
230  m_SeedPointValue = m_SeedPointValue / 27;
231  }
232  else
233  {
234  m_SeedPointValue = (int)m_Function->GetInputImage()->GetPixel(seedIndex);
235  }
236 
237  this->CheckSeedPointValue();
238 
239  m_InitializeValue = (this->CalculateMaxRGS() + 1);
240  if (!m_FineDetectionMode)
241  mitk::ProgressBar::GetInstance()->AddStepsToDo(m_InitializeValue - 1);
242 
243  // only initialize with zeros for the first segmention (raw segmentation mode)
244  if (!m_FineDetectionMode)
245  {
246  this->m_OutputImage->FillBuffer((PixelType)0);
247  }
248  if (m_UpwardsExpansion)
249  {
250  m_Function->ThresholdBetween(m_MinTH, m_SeedPointValue);
251  }
252  else
253  {
254  m_Function->ThresholdBetween(m_SeedPointValue, m_MaxTH);
255  }
256 
257  this->m_IsAtEnd = true;
258 
259  seedIndex = m_StartIndices[0]; // warum noch mal? Steht doch schon in Zeile 224
260 
261  if (this->m_OutputImage->GetBufferedRegion().IsInside(seedIndex) && this->m_SeedPointValue > this->m_MinTH &&
262  this->m_SeedPointValue < this->m_MaxTH)
263  {
264  // Push the seed onto the queue
265  this->InsertIndexTypeIntoQueueMap(m_RegionGrowingState, seedIndex);
266 
267  // Obviously, we're at the beginning
268  this->m_IsAtEnd = false;
269 
270  this->m_OutputImage->SetPixel(seedIndex, (m_InitializeValue - m_RegionGrowingState));
271  }
272  }
273 
274  template <class TImage, class TFunction>
275  void AdaptiveThresholdIterator<TImage, TFunction>::CheckSeedPointValue()
276  {
277  // checks, if the value, that has been averaged over the N-27 neighborhood aorund the seedpoint is still inside
278  // the thresholds-ranges. if not, the actual value of the seedpoint (not averaged) is used
279  if (m_SeedPointValue < m_MinTH || m_SeedPointValue > m_MaxTH)
280  {
281  m_SeedPointValue = (int)m_Function->GetInputImage()->GetPixel(m_StartIndices[0]);
282  }
283  }
284 
285  template <class TImage, class TFunction>
286  unsigned int AdaptiveThresholdIterator<TImage, TFunction>::CalculateMaxRGS()
287  {
288  if (m_UpwardsExpansion)
289  {
290  return (m_MaxTH - m_SeedPointValue);
291  }
292  else
293  {
294  return (m_SeedPointValue - m_MinTH);
295  }
296  }
297 
298  template <class TImage, class TFunction>
299  bool AdaptiveThresholdIterator<TImage, TFunction>::IsPixelIncluded(const IndexType &index) const
300  {
301  // checks, if grayvalue of current voxel is inside the currently used thresholds
302  return this->m_Function->EvaluateAtIndex(index);
303  }
304 
305  template <class TImage, class TFunction>
306  void AdaptiveThresholdIterator<TImage, TFunction>::InsertIndexTypeIntoQueueMap(unsigned int key, IndexType index)
307  {
308  // first check if the key-specific queue already exists
309  if (m_QueueMap.count(key) == 0)
310  {
311  // if queue doesn´t exist, create it, push the IndexType onto it
312  // and insert it into the map
313 
314  IndexQueueType newQueue;
315  newQueue.push(index);
316 
317  typedef typename QueueMapType::value_type KeyIndexQueue;
318 
319  m_QueueMap.insert(KeyIndexQueue(key, newQueue));
320  }
321  else
322  {
323  // if queue already exists in the map, push IndexType onto its specific queue
324  (*(m_QueueMap.find(key))).second.push(index);
325  }
326  }
327 
328  template <class TImage, class TFunction>
329  void AdaptiveThresholdIterator<TImage, TFunction>::DoExtendedFloodStep()
330  {
331  // The index in the front of the queue should always be
332  // valid and be inside since this is what the iterator
333  // uses in the Set/Get methods. This is ensured by the
334  // GoToBegin() method.
335 
336  typename QueueMapType::iterator currentIt = m_QueueMap.find(m_RegionGrowingState);
337 
338  if (currentIt == m_QueueMap.end())
339  {
340  this->IncrementRegionGrowingState();
341  }
342  else
343  {
344  IndexQueueType *currentQueue = &(*currentIt).second;
345 
346  // Take the index in the front of the queue
347  const IndexType &topIndex = currentQueue->front();
348 
349  // Iterate through all possible dimensions
350  // NOTE: Replace this with a ShapeNeighborhoodIterator
351  for (unsigned int i = 0; i < NDimensions; i++)
352  {
353  // The j loop establishes either left or right neighbor (+-1)
354  for (int j = -1; j <= 1; j += 2)
355  {
356  IndexType tempIndex;
357 
358  // build the index of a neighbor
359  for (unsigned int k = 0; k < NDimensions; k++)
360  {
361  if (i != k)
362  {
363  tempIndex.m_Index[k] = topIndex[k];
364  }
365  else
366  {
367  tempIndex.m_Index[k] = topIndex[k] + j;
368  }
369  } // end build the index of a neighbor
370 
371  // If this is a valid index and have not been tested,
372  // then test it.
373  if (m_ImageRegion.IsInside(tempIndex))
374  {
375  // check if voxel hasn´t already been processed
376  if (this->m_OutputImage->GetPixel(tempIndex) == 0)
377  {
378  // if it is inside, push it into the queue
379  if (this->IsPixelIncluded(tempIndex))
380  {
381  // hier wird Voxel in momentan aktiven Stack und ins OutputImage geschrieben
382  this->InsertIndexTypeIntoQueueMap((m_RegionGrowingState), tempIndex);
383  this->m_OutputImage->SetPixel(tempIndex, (m_InitializeValue - m_RegionGrowingState));
384  }
385  else // If the pixel is not inside the current threshold
386  {
387  int distance = this->EstimateDistance(
388  tempIndex); // [!] sollte nicht estimateDistance sondern calculateDistance() heißen!
389  if (distance != 0)
390  {
391  // hier wird Voxel in entsprechenden Stack und ins OutputImage geschrieben
392  this->InsertIndexTypeIntoQueueMap((distance), tempIndex);
393  this->m_OutputImage->SetPixel(tempIndex, (m_InitializeValue - distance));
394  }
395  }
396  }
397  }
398  } // end left/right neighbor loop
399  } // end check all neighbors
400 
401  // Now that all the potential neighbors have been
402  // inserted we can get rid of the pixel in the front
403  currentQueue->pop();
404  m_VoxelCounter++;
405 
406  if (currentQueue->empty())
407  {
408  // if currently used queue is empty
409  this->IncrementRegionGrowingState();
410  }
411  }
412  if (m_RegionGrowingState >= (m_InitializeValue) || m_DetectionStop)
413  {
414  this->m_IsAtEnd = true;
415  // std::cout << "RegionGrowing finished !" << std::endl;
416  // std::cout << "Detected point of leakage: " << m_DetectedLeakagePoint << std::endl;
417  }
418  }
419 
420 } // end namespace itk
421 
422 #endif