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