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