Medical Imaging Interaction Toolkit  2018.4.99-389bf124
Medical Imaging Interaction Toolkit
itkConnectedAdaptiveThresholdImageFilter.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 
13 #ifndef _itkConnectedAdaptiveThresholdImageFilter_txx
14 #define _itkConnectedAdaptiveThresholdImageFilter_txx
15 
16 #include "itkAdaptiveThresholdIterator.h"
17 #include "itkBinaryThresholdImageFunction.h"
18 #include "itkConnectedAdaptiveThresholdImageFilter.h"
19 #include "itkMinimumMaximumImageFilter.h"
20 #include "itkThresholdImageFilter.h"
21 
22 namespace itk
23 {
24  /**
25  * Constructor
26  */
27  template <class TInputImage, class TOutputImage>
28  ConnectedAdaptiveThresholdImageFilter<TInputImage, TOutputImage>::ConnectedAdaptiveThresholdImageFilter()
29  : m_OutoutImageMaskFineSegmentation(nullptr),
30  m_GrowingDirectionIsUpwards(true),
31  m_SeedpointValue(0),
32  m_DetectedLeakagePoint(0),
33  m_InitValue(0),
34  m_AdjLowerTh(0),
35  m_AdjUpperTh(0),
36  m_FineDetectionMode(false),
37  m_DiscardLastPreview(false),
38  m_SegmentationCancelled(false)
39  {
40  }
41 
42  template <class TInputImage, class TOutputImage>
43  void ConnectedAdaptiveThresholdImageFilter<TInputImage, TOutputImage>::GenerateData()
44  {
45  typename ConnectedAdaptiveThresholdImageFilter::InputImageConstPointer inputImage = this->GetInput();
46  typename ConnectedAdaptiveThresholdImageFilter::OutputImagePointer outputImage = this->GetOutput();
47 
48  typename Superclass::InputPixelObjectType::Pointer lowerThreshold = this->GetLowerInput();
49  typename Superclass::InputPixelObjectType::Pointer upperThreshold = this->GetUpperInput();
50 
51  // kommt drauf, wie wir hier die Pipeline aufbauen
52  this->SetLower(lowerThreshold->Get());
53  this->SetUpper(upperThreshold->Get());
54  typedef BinaryThresholdImageFunction<InputImageType> FunctionType;
55  typedef AdaptiveThresholdIterator<OutputImageType, FunctionType> IteratorType;
56 
57  int initValue = IteratorType::CalculateInitializeValue((int)(this->GetLower()), (int)(this->GetUpper()));
58 
59  // Initialize the output according to the segmentation (fine or raw)
60  if (m_FineDetectionMode)
61  {
62  outputImage = this->m_OutoutImageMaskFineSegmentation;
63  }
64 
65  typename ConnectedAdaptiveThresholdImageFilter::OutputImageRegionType region = outputImage->GetRequestedRegion();
66  outputImage->SetBufferedRegion(region);
67  outputImage->Allocate();
68  if (!m_FineDetectionMode)
69  { // only initalize the output image if we are using the raw segmentation mode
70  outputImage->FillBuffer((typename ConnectedAdaptiveThresholdImageFilter::OutputImagePixelType)initValue);
71  }
72 
73  typename FunctionType::Pointer function = FunctionType::New();
74  function->SetInputImage(inputImage);
75 
76  typename Superclass::SeedContainerType seeds;
77  seeds = this->GetSeeds();
78 
79  // pass parameters needed for region growing to iterator
80  IteratorType it(outputImage, function, seeds);
81  it.SetFineDetectionMode(m_FineDetectionMode);
82  it.SetExpansionDirection(m_GrowingDirectionIsUpwards);
83  it.SetMinTH((int)(this->GetLower()));
84  it.SetMaxTH((int)(this->GetUpper()));
85  it.GoToBegin();
86  this->m_SeedpointValue = it.GetSeedPointValue();
87 
88  if ((this->GetLower()) > this->m_SeedpointValue || this->m_SeedpointValue > (this->GetUpper()))
89  {
90  // set m_SegmentationCancelled to true, so if it doesn't reach the point where it is set back to false
91  // we can asssume that there was an error
92  this->m_SegmentationCancelled = true;
93  return;
94  }
95 
96  // iterate through image until
97  while (!it.IsAtEnd())
98  {
99  // make iterator go one step further (calls method DoFloodStep())
100  ++it;
101  }
102  this->m_DetectedLeakagePoint = it.GetLeakagePoint();
103  this->m_SegmentationCancelled = false;
104  }
105 
106  template <class TInputImage, class TOutputImage>
107  TOutputImage *itk::ConnectedAdaptiveThresholdImageFilter<TInputImage, TOutputImage>::GetResultImage()
108  {
109  return m_OutoutImageMaskFineSegmentation;
110  }
111 
112  template <class TInputImage, class TOutputImage>
113  typename ConnectedAdaptiveThresholdImageFilter<TInputImage, TOutputImage>::IndexType
114  itk::ConnectedAdaptiveThresholdImageFilter<TInputImage, TOutputImage>::CorrectSeedPointPosition(
115  unsigned int sizeOfVolume, int lowerTh, int upperTh)
116  {
117  typename ConnectedAdaptiveThresholdImageFilter::InputImageConstPointer inputImage = this->GetInput();
118  typedef typename TInputImage::IndexType IndexType;
119  IndexType itkIntelligentSeedIndex;
120 
121  int seedPixelValue = inputImage->GetPixel(m_SeedPointIndex);
122 
123  // set new seed index to the voxel with the darkest value and shortest distance to original seed
124  if (seedPixelValue > upperTh || seedPixelValue < lowerTh)
125  {
126  // MITK_INFO << "seed pixel value [BEFORE] = " << seedPixelValue;
127  // ToDo crop region
128  itk::Index<3> workindex;
129  for (int i = 0; i < 3; i++)
130  {
131  workindex[i] = m_SeedPointIndex[i] - sizeOfVolume / 2;
132  if (workindex[i] < 0)
133  workindex[i] = 0;
134  }
135  itk::Size<3> worksize;
136  for (int i = 0; i < 3; i++)
137  {
138  worksize[i] = sizeOfVolume;
139  }
140 
141  itk::ImageRegion<3> workregion(workindex, worksize);
142  itk::ImageRegionIterator<TInputImage> regionIt(const_cast<TInputImage *>(inputImage.GetPointer()), workregion);
143 
144  // int darkestGrayValue=seedPixelValue;
145  int currentGrayValue;
146 
147  float distance = (float)(sizeOfVolume / 2);
148  float relativeDistance = 1; // between 0 and 1
149 
150  mitk::Vector3D seedVector, currentVector;
151  mitk::FillVector3D(seedVector, m_SeedPointIndex[0], m_SeedPointIndex[1], m_SeedPointIndex[2]);
152  currentVector = seedVector;
153 
154  float costValue = 0; // beware, Depending on seeking upper or lower value...
155 
156  for (regionIt.GoToBegin(); !regionIt.IsAtEnd(); ++regionIt)
157  {
158  // get current gray value
159  currentGrayValue = regionIt.Value();
160  // get current seed index
161  m_SeedPointIndex = regionIt.GetIndex();
162 
163  // fill current vector
164  mitk::FillVector3D(currentVector, m_SeedPointIndex[0], m_SeedPointIndex[1], m_SeedPointIndex[2]);
165  // calculate distance from original seed to new seed
166  mitk::Vector3D distVector = currentVector - seedVector;
167  distance = fabs(distVector.GetSquaredNorm());
168  relativeDistance = distance / (sizeOfVolume / 2);
169 
170  // calculate "cost function"
171  float currentCostValue = (1 - relativeDistance) * currentGrayValue;
172 
173  if (currentCostValue < costValue && currentGrayValue < upperTh)
174  {
175  itkIntelligentSeedIndex = regionIt.GetIndex();
176  costValue = currentCostValue;
177 
178  // MITK_INFO <<"cost value="<< costValue;
179  // MITK_INFO <<"darkest and closest Voxel ="<< currentGrayValue;
180  // MITK_INFO <<"m_UPPER="<< upperTh;
181  }
182  }
183  // MITK_INFO<< "seed pixel value [AFTER] =" << inputImage->GetPixel(itkIntelligentSeedIndex) <<"\n";
184  }
185  else
186  { // no correction of the seed point is needed, just pass the original seed
187  itkIntelligentSeedIndex = m_SeedPointIndex;
188  }
189 
190  return itkIntelligentSeedIndex;
191  }
192 
193  template <class TInputImage, class TOutputImage>
194  void itk::ConnectedAdaptiveThresholdImageFilter<TInputImage, TOutputImage>::CropMask(unsigned int croppingSize)
195  {
196  // initialize center point of the working region
197  itk::Index<3> workindex;
198  for (int i = 0; i < 3; i++)
199  {
200  workindex[i] = m_SeedPointIndex[i] - croppingSize / 2;
201  if (workindex[i] < 0)
202  workindex[i] = 0;
203  }
204  // initialize working volume
205  itk::Size<3> worksize;
206  for (int i = 0; i < 3; i++)
207  {
208  worksize[i] = croppingSize;
209  }
210  // set working region
211  itk::ImageRegion<3> workregion(workindex, worksize);
212 
213  // check if the entire region is inside the image
214  if (!(m_OutoutImageMaskFineSegmentation->GetLargestPossibleRegion().IsInside(workregion)))
215  {
216  // if not then crop to the intersection of the image (gemeinsame Schnittmenge Bild und workingRegion)
217  if (!(workregion.Crop(m_OutoutImageMaskFineSegmentation->GetLargestPossibleRegion())))
218  {
219  MITK_ERROR << "Cropping working region failed!";
220  return;
221  }
222  }
223 
224  // initialize region iterator
225  itk::ImageRegionIterator<TOutputImage> regionIt(m_OutoutImageMaskFineSegmentation, workregion);
226  for (regionIt.GoToBegin(); !regionIt.IsAtEnd(); ++regionIt)
227  {
228  // and set all voxel inside the working region to zero
229  regionIt.Set(0);
230  }
231  }
232 
233  template <class TInputImage, class TOutputImage>
234  unsigned int itk::ConnectedAdaptiveThresholdImageFilter<TInputImage, TOutputImage>::AdjustIteratorMask()
235  {
236  typedef itk::ThresholdImageFilter<TOutputImage> ThresholdFilterType;
237  typedef itk::MinimumMaximumImageFilter<TOutputImage> MaxFilterType;
238 
239  typename ThresholdFilterType::Pointer thresholdFilter = ThresholdFilterType::New();
240  typename MaxFilterType::Pointer maxFilter = MaxFilterType::New();
241 
242  unsigned int maxValue;
243 
244  if (!m_DiscardLastPreview)
245  {
246  // get the biggest value of the image
247  maxFilter->SetInput(m_OutoutImageMaskFineSegmentation);
248  maxFilter->UpdateLargestPossibleRegion();
249  maxValue = maxFilter->GetMaximum();
250  }
251  else
252  { // use the last biggest value in the preview. This was set in SetParameterForFineSegmentation(...adjLowerTh...) []
253  maxValue = m_AdjLowerTh;
254  }
255 
256  // set all values <lower && >upper to zero (thresouldOutside uses < and > NOT <= and >=)
257  thresholdFilter->SetInput(m_OutoutImageMaskFineSegmentation);
258  thresholdFilter->SetOutsideValue(0);
259  thresholdFilter->ThresholdOutside(m_AdjLowerTh, maxValue);
260  thresholdFilter->UpdateLargestPossibleRegion();
261 
262  // set all values in between lower and upper (>=lower && <=upper) to the highest value in the image
263  thresholdFilter->SetInput(thresholdFilter->GetOutput());
264  thresholdFilter->SetOutsideValue(maxValue);
265  thresholdFilter->ThresholdOutside(0, m_AdjLowerTh - 1);
266  thresholdFilter->UpdateLargestPossibleRegion();
267 
268  m_OutoutImageMaskFineSegmentation = thresholdFilter->GetOutput();
269 
270  return maxValue;
271  }
272 
273  template <class TInputImage, class TOutputImage>
274  void itk::ConnectedAdaptiveThresholdImageFilter<TInputImage, TOutputImage>::SetParameterForFineSegmentation(
275  TOutputImage *iteratorMaskForFineSegmentation,
276  unsigned int adjLowerTh,
277  unsigned int adjUpperTh,
278  itk::Index<3> seedPoint,
279  bool discardLeafSegmentation)
280  {
281  // just to make sure we´re in the right mode and the mask exsits
282  if (m_FineDetectionMode && iteratorMaskForFineSegmentation)
283  {
284  m_OutoutImageMaskFineSegmentation = iteratorMaskForFineSegmentation;
285  m_AdjLowerTh = adjLowerTh;
286  m_AdjUpperTh = adjUpperTh; // still needed?
287  m_SeedPointIndex = seedPoint;
288  m_DiscardLastPreview = discardLeafSegmentation;
289  }
290  else
291  {
292  if (!m_FineDetectionMode)
293  {
294  MITK_ERROR << "Fine-detection-segmentation mode not set!";
295  }
296  else
297  {
298  MITK_ERROR << "Iterator-mask-image not set!";
299  }
300  }
301  }
302 
303 } // end namespace itk
304 
305 #endif