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 ===================================================================*/
16 #ifndef _itkAdaptiveThresholdIterator_txx
17 #define _itkAdaptiveThresholdIterator_txx
19 #include "itkAdaptiveThresholdIterator.h"
20 #include "mitkProgressBar.h"
25 template <class TImage, class TFunction>
26 AdaptiveThresholdIterator<TImage, TFunction>::AdaptiveThresholdIterator(ImageType *imagePtr,
29 : m_FineDetectionMode(false), m_DetectionStop(false)
31 this->m_OutputImage = imagePtr;
33 m_StartIndices.push_back(startIndex);
35 // Set up the temporary image
36 this->InitializeIterator();
39 template <class TImage, class TFunction>
40 AdaptiveThresholdIterator<TImage, TFunction>::AdaptiveThresholdIterator(ImageType *imagePtr,
42 std::vector<IndexType> &startIndex)
43 : m_FineDetectionMode(false), m_DetectionStop(false)
45 this->m_OutputImage = imagePtr;
48 for (i = 0; i < startIndex.size(); i++)
50 m_StartIndices.push_back(startIndex[i]);
53 // Set up the temporary image
54 this->InitializeIterator();
57 template <class TImage, class TFunction>
58 AdaptiveThresholdIterator<TImage, TFunction>::AdaptiveThresholdIterator(ImageType *imagePtr, FunctionType *fnPtr)
59 : m_FineDetectionMode(false), m_DetectionStop(false)
61 this->m_OutputImage = imagePtr; // here we store the image, we have to wite the result to
64 // Set up the temporary image
65 this->InitializeIterator();
68 template <class TImage, class TFunction>
69 void AdaptiveThresholdIterator<TImage, TFunction>::InitializeIterator()
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();
76 this->InitRegionGrowingState();
79 m_LastVoxelNumber = 0;
80 m_CurrentLeakageRatio = 0;
81 m_DetectedLeakagePoint = 0;
84 template <class TImage, class TFunction>
85 void AdaptiveThresholdIterator<TImage, TFunction>::SetExpansionDirection(bool upwards)
87 m_UpwardsExpansion = upwards;
90 template <class TImage, class TFunction>
91 void AdaptiveThresholdIterator<TImage, TFunction>::IncrementRegionGrowingState()
93 // make the progressbar go one step further
94 if (!m_FineDetectionMode)
95 mitk::ProgressBar::GetInstance()->Progress();
97 // updating the thresholds
98 if (m_UpwardsExpansion)
100 this->ExpandThresholdUpwards();
104 this->ExpandThresholdDownwards();
108 int criticalValue = 2000; // calculate a bit more "intelligent"
110 if (!m_FineDetectionMode)
112 int diff = m_VoxelCounter - m_LastVoxelNumber;
113 if (diff > m_CurrentLeakageRatio)
115 m_CurrentLeakageRatio = diff;
116 m_DetectedLeakagePoint = m_RegionGrowingState;
118 m_LastVoxelNumber = m_VoxelCounter;
121 else // fine leakage detection
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))
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";
134 m_DetectionStop = true;
135 // std::cout<<"\n\n[ITERATOR] detection stop!!!\n";
137 m_LastVoxelNumber = m_VoxelCounter;
140 // increment the counter
141 m_RegionGrowingState++;
144 template <class TImage, class TFunction>
145 void AdaptiveThresholdIterator<TImage, TFunction>::ExpandThresholdUpwards()
147 int upper = (int)m_Function->GetUpper();
149 m_Function->ThresholdBetween(m_MinTH, upper);
152 template <class TImage, class TFunction>
153 void AdaptiveThresholdIterator<TImage, TFunction>::ExpandThresholdDownwards()
155 int lower = (int)m_Function->GetLower();
157 m_Function->ThresholdBetween(lower, m_MaxTH);
160 template <class TImage, class TFunction>
161 void AdaptiveThresholdIterator<TImage, TFunction>::InitRegionGrowingState()
163 this->m_RegionGrowingState = 1;
166 template <class TImage, class TFunction>
167 int AdaptiveThresholdIterator<TImage, TFunction>::EstimateDistance(IndexType tempIndex)
169 PixelType value = this->m_Function->GetInputImage()->GetPixel(tempIndex);
170 PixelType minPixel = PixelType(m_MinTH);
171 PixelType maxPixel = PixelType(m_MaxTH);
173 if (value > maxPixel || value < minPixel)
177 if (m_UpwardsExpansion)
179 return (int)(value - m_SeedPointValue);
183 return (int)(m_SeedPointValue - value);
187 template <class TImage, class TFunction>
188 void AdaptiveThresholdIterator<TImage, TFunction>::SetMinTH(int min)
193 template <class TImage, class TFunction>
194 void AdaptiveThresholdIterator<TImage, TFunction>::SetMaxTH(int max)
199 template <class TImage, class TFunction>
200 int AdaptiveThresholdIterator<TImage, TFunction>::GetSeedPointValue()
202 return this->m_SeedPointValue;
205 template <class TImage, class TFunction>
206 void AdaptiveThresholdIterator<TImage, TFunction>::GoToBegin()
210 m_SeedPointValue = 0;
212 IndexType seedIndex = m_StartIndices[0];
214 bool doAverage = false; // enable or disable manually!
217 // loops for creating the sum of the N27-neighborhood around the seedpoint
218 for (int i = -1; i <= 1; i++)
220 for (int j = -1; j <= 1; j++)
222 for (int k = -1; k <= 1; k++)
224 seedIndex[0] = seedIndex[0] + i;
225 seedIndex[1] = seedIndex[1] + j;
226 seedIndex[2] = seedIndex[2] + k;
228 m_SeedPointValue += (int)m_Function->GetInputImage()->GetPixel(seedIndex);
233 // value of seedpoint computed from mean of N27-neighborhood
234 m_SeedPointValue = m_SeedPointValue / 27;
238 m_SeedPointValue = (int)m_Function->GetInputImage()->GetPixel(seedIndex);
241 this->CheckSeedPointValue();
243 m_InitializeValue = (this->CalculateMaxRGS() + 1);
244 if (!m_FineDetectionMode)
245 mitk::ProgressBar::GetInstance()->AddStepsToDo(m_InitializeValue - 1);
247 // only initialize with zeros for the first segmention (raw segmentation mode)
248 if (!m_FineDetectionMode)
250 this->m_OutputImage->FillBuffer((PixelType)0);
252 if (m_UpwardsExpansion)
254 m_Function->ThresholdBetween(m_MinTH, m_SeedPointValue);
258 m_Function->ThresholdBetween(m_SeedPointValue, m_MaxTH);
261 this->m_IsAtEnd = true;
263 seedIndex = m_StartIndices[0]; // warum noch mal? Steht doch schon in Zeile 224
265 if (this->m_OutputImage->GetBufferedRegion().IsInside(seedIndex) && this->m_SeedPointValue > this->m_MinTH &&
266 this->m_SeedPointValue < this->m_MaxTH)
268 // Push the seed onto the queue
269 this->InsertIndexTypeIntoQueueMap(m_RegionGrowingState, seedIndex);
271 // Obviously, we're at the beginning
272 this->m_IsAtEnd = false;
274 this->m_OutputImage->SetPixel(seedIndex, (m_InitializeValue - m_RegionGrowingState));
278 template <class TImage, class TFunction>
279 void AdaptiveThresholdIterator<TImage, TFunction>::CheckSeedPointValue()
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)
285 m_SeedPointValue = (int)m_Function->GetInputImage()->GetPixel(m_StartIndices[0]);
289 template <class TImage, class TFunction>
290 unsigned int AdaptiveThresholdIterator<TImage, TFunction>::CalculateMaxRGS()
292 if (m_UpwardsExpansion)
294 return (m_MaxTH - m_SeedPointValue);
298 return (m_SeedPointValue - m_MinTH);
302 template <class TImage, class TFunction>
303 bool AdaptiveThresholdIterator<TImage, TFunction>::IsPixelIncluded(const IndexType &index) const
305 // checks, if grayvalue of current voxel is inside the currently used thresholds
306 return this->m_Function->EvaluateAtIndex(index);
309 template <class TImage, class TFunction>
310 void AdaptiveThresholdIterator<TImage, TFunction>::InsertIndexTypeIntoQueueMap(unsigned int key, IndexType index)
312 // first check if the key-specific queue already exists
313 if (m_QueueMap.count(key) == 0)
315 // if queue doesn´t exist, create it, push the IndexType onto it
316 // and insert it into the map
318 IndexQueueType newQueue;
319 newQueue.push(index);
321 typedef typename QueueMapType::value_type KeyIndexQueue;
323 m_QueueMap.insert(KeyIndexQueue(key, newQueue));
327 // if queue already exists in the map, push IndexType onto its specific queue
328 (*(m_QueueMap.find(key))).second.push(index);
332 template <class TImage, class TFunction>
333 void AdaptiveThresholdIterator<TImage, TFunction>::DoExtendedFloodStep()
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.
340 typename QueueMapType::iterator currentIt = m_QueueMap.find(m_RegionGrowingState);
342 if (currentIt == m_QueueMap.end())
344 this->IncrementRegionGrowingState();
348 IndexQueueType *currentQueue = &(*currentIt).second;
350 // Take the index in the front of the queue
351 const IndexType &topIndex = currentQueue->front();
353 // Iterate through all possible dimensions
354 // NOTE: Replace this with a ShapeNeighborhoodIterator
355 for (unsigned int i = 0; i < NDimensions; i++)
357 // The j loop establishes either left or right neighbor (+-1)
358 for (int j = -1; j <= 1; j += 2)
362 // build the index of a neighbor
363 for (unsigned int k = 0; k < NDimensions; k++)
367 tempIndex.m_Index[k] = topIndex[k];
371 tempIndex.m_Index[k] = topIndex[k] + j;
373 } // end build the index of a neighbor
375 // If this is a valid index and have not been tested,
377 if (m_ImageRegion.IsInside(tempIndex))
379 // check if voxel hasn´t already been processed
380 if (this->m_OutputImage->GetPixel(tempIndex) == 0)
382 // if it is inside, push it into the queue
383 if (this->IsPixelIncluded(tempIndex))
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));
389 else // If the pixel is not inside the current threshold
391 int distance = this->EstimateDistance(
392 tempIndex); // [!] sollte nicht estimateDistance sondern calculateDistance() heißen!
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));
402 } // end left/right neighbor loop
403 } // end check all neighbors
405 // Now that all the potential neighbors have been
406 // inserted we can get rid of the pixel in the front
410 if (currentQueue->empty())
412 // if currently used queue is empty
413 this->IncrementRegionGrowingState();
416 if (m_RegionGrowingState >= (m_InitializeValue) || m_DetectionStop)
418 this->m_IsAtEnd = true;
419 // std::cout << "RegionGrowing finished !" << std::endl;
420 // std::cout << "Detected point of leakage: " << m_DetectedLeakagePoint << std::endl;
424 } // end namespace itk