Medical Imaging Interaction Toolkit  2018.4.99-f51274ea
Medical Imaging Interaction Toolkit
mitkLabeledImageToSurfaceFilter.cpp
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 
14 
15 #include <vtkDecimatePro.h>
16 #include <vtkImageChangeInformation.h>
17 #include <vtkImageGaussianSmooth.h>
18 #include <vtkImageMarchingCubes.h>
19 #include <vtkImageThreshold.h>
20 #include <vtkLinearTransform.h>
21 #include <vtkMatrix4x4.h>
22 #include <vtkPolyData.h>
23 #include <vtkSmoothPolyDataFilter.h>
24 
25 #include <itkImageRegionIterator.h>
26 #include <itkNumericTraits.h>
27 #include <mitkImageAccessByItk.h>
29 
31  : m_GaussianStandardDeviation(1.5), m_GenerateAllLabels(true), m_Label(1), m_BackgroundLabel(0)
32 {
33 }
34 
36 {
37 }
38 
40 {
41  Superclass::GenerateOutputInformation();
42  //
43  // check which labels are available in the image
44  //
46  m_IdxToLabels.clear();
47 
48  //
49  // if we don't want to generate surfaces for all labels
50  // we have to remove all labels except m_Label and m_BackgroundLabel
51  // from the list of available labels
52  //
54  {
55  LabelMapType tmp;
56  LabelMapType::iterator it;
57  it = m_AvailableLabels.find(m_Label);
58  if (it != m_AvailableLabels.end())
59  tmp[m_Label] = it->second;
60  else
61  tmp[m_Label] = 0;
62 
64  if (it != m_AvailableLabels.end())
65  tmp[m_BackgroundLabel] = it->second;
66  else
67  tmp[m_BackgroundLabel] = 0;
68 
69  m_AvailableLabels = tmp;
70  }
71 
72  //
73  // check for the number of labels: if the whole image is filled, no
74  // background is available and thus the numberOfOutpus is equal to the
75  // number of available labels in the image (which is a special case).
76  // If we have background voxels, the number of outputs is one less than
77  // then number of available labels.
78  //
79  unsigned int numberOfOutputs = 0;
81  numberOfOutputs = m_AvailableLabels.size();
82  else
83  numberOfOutputs = m_AvailableLabels.size() - 1;
84  if (numberOfOutputs == 0)
85  {
86  itkWarningMacro("Number of outputs == 0");
87  }
88 
89  //
90  // determine the number of time steps of the input image
91  //
93 
94  unsigned int numberOfTimeSteps = image->GetTimeGeometry()->CountTimeSteps();
95 
96  //
97  // set the number of outputs to the number of labels used.
98  // initialize the output surfaces accordingly (incl. time steps)
99  //
100  this->SetNumberOfIndexedOutputs(numberOfOutputs);
101  this->SetNumberOfRequiredOutputs(numberOfOutputs);
102  for (unsigned int i = 0; i < numberOfOutputs; ++i)
103  {
104  if (!this->GetOutput(i))
105  {
106  mitk::Surface::Pointer output = static_cast<mitk::Surface *>(this->MakeOutput(0).GetPointer());
107  assert(output.IsNotNull());
108  output->Expand(numberOfTimeSteps);
109  this->SetNthOutput(i, output.GetPointer());
110  }
111  }
112 }
113 
115 {
117  if (image == nullptr)
118  {
119  itkWarningMacro("Image is nullptr");
120  return;
121  }
122 
123  mitk::Image::RegionType outputRegion = image->GetRequestedRegion();
124 
125  m_IdxToLabels.clear();
126 
127  if (this->GetNumberOfOutputs() == 0)
128  return;
129 
130  //
131  // traverse the known labels and create surfaces for them.
132  //
133  unsigned int currentOutputIndex = 0;
134  for (auto it = m_AvailableLabels.begin(); it != m_AvailableLabels.end(); ++it)
135  {
136  if (it->first == m_BackgroundLabel)
137  continue;
138  if ((it->second == 0) && m_GenerateAllLabels)
139  continue;
140 
141  assert(currentOutputIndex < this->GetNumberOfOutputs());
142  mitk::Surface::Pointer surface = this->GetOutput(currentOutputIndex);
143  assert(surface.IsNotNull());
144 
145  int tstart = outputRegion.GetIndex(3);
146  int tmax = tstart + outputRegion.GetSize(3); // GetSize()==1 - will aber 0 haben, wenn nicht zeitaufgeloet
147  int t;
148  for (t = tstart; t < tmax; ++t)
149  {
150  vtkImageData *vtkimagedata = image->GetVtkImageData(t);
151  CreateSurface(t, vtkimagedata, surface.GetPointer(), it->first);
152  }
153  m_IdxToLabels[currentOutputIndex] = it->first;
154  currentOutputIndex++;
155  }
156 }
157 
159  vtkImageData *vtkimage,
160  mitk::Surface *surface,
162 {
163  vtkImageChangeInformation *indexCoordinatesImageFilter = vtkImageChangeInformation::New();
164  indexCoordinatesImageFilter->SetInputData(vtkimage);
165  indexCoordinatesImageFilter->SetOutputOrigin(0.0, 0.0, 0.0);
166 
167  vtkImageThreshold *threshold = vtkImageThreshold::New();
168  threshold->SetInputConnection(indexCoordinatesImageFilter->GetOutputPort());
169  // indexCoordinatesImageFilter->Delete();
170  threshold->SetInValue(100);
171  threshold->SetOutValue(0);
172  threshold->ThresholdBetween(label, label);
173  threshold->SetOutputScalarTypeToUnsignedChar();
174  threshold->ReleaseDataFlagOn();
175 
176  vtkImageGaussianSmooth *gaussian = vtkImageGaussianSmooth::New();
177  gaussian->SetInputConnection(threshold->GetOutputPort());
178  // threshold->Delete();
179  gaussian->SetDimensionality(3);
180  gaussian->SetRadiusFactor(0.49);
181  gaussian->SetStandardDeviation(GetGaussianStandardDeviation());
182  gaussian->ReleaseDataFlagOn();
183  gaussian->UpdateInformation();
184  gaussian->Update();
185 
186  // MarchingCube -->create Surface
187  vtkMarchingCubes *skinExtractor = vtkMarchingCubes::New();
188  skinExtractor->ReleaseDataFlagOn();
189  skinExtractor->SetInputConnection(gaussian->GetOutputPort()); // RC++
190  indexCoordinatesImageFilter->Delete();
191  skinExtractor->SetValue(0, 50);
192 
193  vtkPolyData *polydata;
194  skinExtractor->Update();
195  polydata = skinExtractor->GetOutput();
196  polydata->Register(nullptr); // RC++
197  skinExtractor->Delete();
198 
199  if (m_Smooth)
200  {
201  vtkSmoothPolyDataFilter *smoother = vtkSmoothPolyDataFilter::New();
202  // read poly1 (poly1 can be the original polygon, or the decimated polygon)
203  smoother->SetInputData(polydata); // RC++
204  smoother->SetNumberOfIterations(m_SmoothIteration);
205  smoother->SetRelaxationFactor(m_SmoothRelaxation);
206  smoother->SetFeatureAngle(60);
207  smoother->FeatureEdgeSmoothingOff();
208  smoother->BoundarySmoothingOff();
209  smoother->SetConvergence(0);
210 
211  polydata->Delete(); // RC--
212  smoother->Update();
213  polydata = smoother->GetOutput();
214  polydata->Register(nullptr); // RC++
215  smoother->Delete();
216  }
217 
218  // decimate = to reduce number of polygons
219  if (m_Decimate == DecimatePro)
220  {
221  vtkDecimatePro *decimate = vtkDecimatePro::New();
222  decimate->SplittingOff();
223  decimate->SetErrorIsAbsolute(5);
224  decimate->SetFeatureAngle(30);
225  decimate->PreserveTopologyOn();
226  decimate->BoundaryVertexDeletionOff();
227  decimate->SetDegree(10); // std-value is 25!
228 
229  decimate->SetInputData(polydata); // RC++
230  decimate->SetTargetReduction(m_TargetReduction);
231  decimate->SetMaximumError(0.002);
232 
233  polydata->Delete(); // RC--
234  decimate->Update();
235  polydata = decimate->GetOutput();
236  polydata->Register(nullptr); // RC++
237  decimate->Delete();
238  }
239 
240  if (polydata->GetNumberOfPoints() > 0)
241  {
242  mitk::Vector3D spacing = GetInput()->GetGeometry(time)->GetSpacing();
243 
244  vtkPoints *points = polydata->GetPoints();
245  vtkMatrix4x4 *vtkmatrix = vtkMatrix4x4::New();
246  GetInput()->GetGeometry(time)->GetVtkTransform()->GetMatrix(vtkmatrix);
247  double(*matrix)[4] = vtkmatrix->Element;
248 
249  unsigned int i, j;
250  for (i = 0; i < 3; ++i)
251  for (j = 0; j < 3; ++j)
252  matrix[i][j] /= spacing[j];
253 
254  unsigned int n = points->GetNumberOfPoints();
255  double point[3];
256 
257  for (i = 0; i < n; i++)
258  {
259  points->GetPoint(i, point);
260  mitkVtkLinearTransformPoint(matrix, point, point);
261  points->SetPoint(i, point);
262  }
263  vtkmatrix->Delete();
264  }
265 
266  surface->SetVtkPolyData(polydata, time);
267  polydata->UnRegister(nullptr);
268 
269  gaussian->Delete();
270  threshold->Delete();
271 }
272 
273 template <typename TPixel, unsigned int VImageDimension>
274 void GetAvailableLabelsInternal(itk::Image<TPixel, VImageDimension> *image,
276 {
277  typedef itk::Image<TPixel, VImageDimension> ImageType;
278  typedef itk::ImageRegionIterator<ImageType> ImageRegionIteratorType;
279  availableLabels.clear();
280  ImageRegionIteratorType it(image, image->GetLargestPossibleRegion());
281  it.GoToBegin();
282  mitk::LabeledImageToSurfaceFilter::LabelMapType::iterator labelIt;
283  while (!it.IsAtEnd())
284  {
285  labelIt = availableLabels.find((mitk::LabeledImageToSurfaceFilter::LabelType)(it.Get()));
286  if (labelIt == availableLabels.end())
287  {
288  availableLabels[(mitk::LabeledImageToSurfaceFilter::LabelType)(it.Get())] = 1;
289  }
290  else
291  {
292  labelIt->second += 1;
293  }
294 
295  ++it;
296  }
297 }
298 
299 #define InstantiateAccessFunction_GetAvailableLabelsInternal(pixelType, dim) \
300  \
301 template void \
302  GetAvailableLabelsInternal(itk::Image<pixelType, dim> *, mitk::LabeledImageToSurfaceFilter::LabelMapType &);
303 
305 
307 {
309  LabelMapType availableLabels;
310  AccessFixedDimensionByItk_1(image, GetAvailableLabelsInternal, 3, availableLabels);
311  return availableLabels;
312 }
313 
315 {
316  itkWarningMacro("This function should never be called!");
317  assert(false);
318 }
319 
321  const unsigned int &idx)
322 {
323  auto it = m_IdxToLabels.find(idx);
324  if (it != m_IdxToLabels.end())
325  {
326  return it->second;
327  }
328  else
329  {
330  itkWarningMacro("Unknown index encountered: " << idx << ". There are " << this->GetNumberOfOutputs()
331  << " outputs available.");
333  }
334 }
335 
337 {
339 }
340 
343 {
344  // get the image spacing
346  const mitk::Vector3D spacing = image->GetSlicedGeometry()->GetSpacing();
347 
348  // get the number of voxels encountered for the given label,
349  // calculate the volume and return it.
350  auto it = m_AvailableLabels.find(label);
351  if (it != m_AvailableLabels.end())
352  {
353  return static_cast<float>(it->second) * (spacing[0] * spacing[1] * spacing[2] / 1000.0f);
354  }
355  else
356  {
357  itkWarningMacro("Unknown label encountered: " << label);
358  return 0.0;
359  }
360 }
virtual void CreateSurface(int time, vtkImageData *vtkimage, mitk::Surface *surface, LabelType label)
Class for storing surfaces (vtkPolyData).
Definition: mitkSurface.h:28
mitk::ScalarType GetVolumeForNthOutput(const unsigned int &i)
virtual TimeStepType CountTimeSteps() const =0
Returns the number of time steps.
itk::Image< unsigned char, 3 > ImageType
double ScalarType
InstantiateAccessFunctionForFixedDimension(GetAvailableLabelsInternal, 3)
vtkLinearTransform * GetVtkTransform() const
Get the m_IndexToWorldTransform as a vtkLinearTransform.
OutputType * GetOutput()
const mitk::Image * GetInput(void)
virtual vtkImageData * GetVtkImageData(int t=0, int n=0)
Get a volume at a specific time t of channel n as a vtkImageData.
Definition: mitkImage.cpp:217
LabelType GetLabelForNthOutput(const unsigned int &i)
void mitkVtkLinearTransformPoint(T1 matrix[4][4], T2 in[3], T3 out[3])
const mitk::TimeGeometry * GetTimeGeometry() const
Return the TimeGeometry of the data as const pointer.
Definition: mitkBaseData.h:66
itk::ImageRegion< RegionDimension > RegionType
Image class for storing images.
Definition: mitkImage.h:72
SlicedGeometry3D * GetSlicedGeometry(unsigned int t=0) const
Convenience access method for the geometry, which is of type SlicedGeometry3D (or a sub-class of it)...
virtual const RegionType & GetRequestedRegion() const
virtual void SetVtkPolyData(vtkPolyData *polydata, unsigned int t=0)
static T max(T x, T y)
Definition: svm.cpp:56
mitk::Image::Pointer image
void GetAvailableLabelsInternal(itk::Image< TPixel, VImageDimension > *image, mitk::LabeledImageToSurfaceFilter::LabelMapType &availableLabels)
mitk::ScalarType GetVolumeForLabel(const LabelType &label)
#define AccessFixedDimensionByItk_1(mitkImage, itkImageTypeFunction, dimension, arg1)
const mitk::Vector3D GetSpacing() const
Get the spacing (size of a pixel).
itk::DataObject::Pointer MakeOutput(DataObjectPointerArraySizeType idx) override
virtual double GetGaussianStandardDeviation()
std::map< LabelType, unsigned long > LabelMapType
mitk::BaseGeometry * GetGeometry(int t=0) const
Return the geometry, which is a TimeGeometry, of the data as non-const pointer.
Definition: mitkBaseData.h:143