Medical Imaging Interaction Toolkit  2018.4.99-08619e4f
Medical Imaging Interaction Toolkit
itkIntelligentBinaryClosingFilter.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 _itkIntelligentBinaryClosingFilter_txx
13 #define _itkIntelligentBinaryClosingFilter_txx
14 
15 #include "itkIntelligentBinaryClosingFilter.h"
16 #include <vector>
17 
18 namespace itk
19 {
20  template <class TInputImage, class TOutputImage>
21  IntelligentBinaryClosingFilter<TInputImage, TOutputImage>::IntelligentBinaryClosingFilter(): m_ClosingRadius(4.0), m_SurfaceRatio(70)
22  {
23  m_ErodeImageFilter = BinaryErodeImageFilterType::New();
24  m_DilateImageFilter = BinaryDilateImageFilterType::New();
25  m_SubtractImageFilter = SubtractImageFilterType::New();
26  m_ConnectedComponentImageFilter = ConnectedComponentImageFilterType::New();
27  m_RelabelComponentImageFilter = RelabelComponentImageFilterType::New();
28  m_RelabelComponentImageFilter->SetInPlace(true);
29  m_BorderDetectionDilateFilter = DilateComponentImageFilterType::New();
30  }
31 
32  template <class TInputImage, class TOutputImage>
33  void IntelligentBinaryClosingFilter<TInputImage, TOutputImage>::GenerateData()
34  {
35  const InputImageType *input = this->GetInput();
36  // Allocate the output image.
37  typename OutputImageType::Pointer output = this->GetOutput();
38  output->SetRequestedRegion(input->GetRequestedRegion());
39  output->SetBufferedRegion(input->GetBufferedRegion());
40  output->SetLargestPossibleRegion(input->GetLargestPossibleRegion());
41  output->Allocate();
42 
43  // set up structuring element for closing
44  StructuringElementType seClosing;
45  itk::SizeValueType radius[ImageDimension];
46  const typename InputImageType::SpacingType spacing = input->GetSpacing();
47  for (unsigned int d = 0; d < ImageDimension; d++)
48  { // closing works in voxel coordinates, so use spacing (and add 0.5 for correct rounding - cast just truncates)
49  radius[d] = (unsigned long)(m_ClosingRadius / spacing[d] + 0.5);
50  }
51  MITK_INFO << " Closing kernel size = [" << radius[0] << ", " << radius[1] << ", " << radius[2] << "]"
52  << std::endl;
53  seClosing.SetRadius(radius);
54  seClosing.CreateStructuringElement();
55 
56  // closing
57  m_DilateImageFilter->SetInput(this->GetInput());
58  m_DilateImageFilter->SetKernel(seClosing);
59  m_DilateImageFilter->SetDilateValue(1);
60  m_ErodeImageFilter->SetInput(m_DilateImageFilter->GetOutput());
61  m_ErodeImageFilter->SetKernel(seClosing);
62  m_ErodeImageFilter->SetErodeValue(1);
63 
64  // subtraction
65  m_SubtractImageFilter->SetInput1(m_ErodeImageFilter->GetOutput());
66  m_SubtractImageFilter->SetInput2(this->GetInput());
67 
68  // connected components
69  m_ConnectedComponentImageFilter->SetInput(m_SubtractImageFilter->GetOutput());
70  m_RelabelComponentImageFilter->SetInput(m_ConnectedComponentImageFilter->GetOutput());
71  m_RelabelComponentImageFilter->Update();
72 
73  // set up structuring element for border voxel detection
74  StructuringElementType seBorder;
75  for (auto &radiu : radius)
76  {
77  radiu = 1;
78  }
79  seBorder.SetRadius(radius);
80  seBorder.CreateStructuringElement();
81 
82  // dilate all components to detect border voxels
83  m_BorderDetectionDilateFilter->SetInput(m_RelabelComponentImageFilter->GetOutput());
84  m_BorderDetectionDilateFilter->SetKernel(seBorder);
85  m_BorderDetectionDilateFilter->Update();
86 
87  // count volumes and border voxels for all components
88  OutputIteratorType itComp(m_RelabelComponentImageFilter->GetOutput(),
89  m_RelabelComponentImageFilter->GetOutput()->GetLargestPossibleRegion());
90  OutputIteratorType itBorder(m_BorderDetectionDilateFilter->GetOutput(),
91  m_BorderDetectionDilateFilter->GetOutput()->GetLargestPossibleRegion());
92  ConstInputIteratorType itIn(input, input->GetLargestPossibleRegion());
93 
94  std::vector<unsigned int> volume(m_RelabelComponentImageFilter->GetNumberOfObjects() + 1, 0);
95  std::vector<unsigned int> border(m_RelabelComponentImageFilter->GetNumberOfObjects() + 1, 0);
96  std::vector<unsigned int> adjacent(m_RelabelComponentImageFilter->GetNumberOfObjects() + 1, 0);
97  typename OutputImageType::ValueType borderId, compId;
98  for (itComp.GoToBegin(), itBorder.GoToBegin(), itIn.GoToBegin(); !itComp.IsAtEnd(); ++itComp, ++itBorder, ++itIn)
99  {
100  borderId = itBorder.Get();
101  if (borderId != 0)
102  {
103  compId = itComp.Get();
104  if (compId != 0)
105  {
106  volume[compId]++;
107  }
108  else
109  {
110  // this is border country
111  border[borderId]++;
112  if (itIn.Get() != 0)
113  adjacent[borderId]++;
114  }
115  }
116  }
117 
118  // calculate ratios
119  std::vector<float> ratio(m_RelabelComponentImageFilter->GetNumberOfObjects() + 1, 0);
120  MITK_INFO << " " << m_RelabelComponentImageFilter->GetNumberOfObjects() << " components found" << std::endl;
121  for (unsigned int i = 0; i < ratio.size(); i++)
122  {
123  if (border[i] != 0)
124  ratio[i] = 100.0 * (float)(adjacent[i]) / (float)(border[i]);
125  }
126 
127  // fill output
128  OutputIteratorType itOut(output, output->GetLargestPossibleRegion());
129  for (itOut.GoToBegin(), itIn.GoToBegin(), itComp.GoToBegin(); !itOut.IsAtEnd(); ++itOut, ++itIn, ++itComp)
130  {
131  if (itIn.Get() != 0)
132  {
133  itOut.Set(1);
134  }
135  else
136  {
137  compId = itComp.Get();
138  if (ratio[compId] > m_SurfaceRatio)
139  itOut.Set(1);
140  else
141  itOut.Set(0);
142  }
143  }
144  }
145 
146  template <class TInputImage, class TOutputImage>
147  void IntelligentBinaryClosingFilter<TInputImage, TOutputImage>::PrintSelf(std::ostream &os, Indent indent) const
148  {
149  Superclass::PrintSelf(os, indent);
150  }
151 
152 } // end namespace itk
153 
154 #endif