Medical Imaging Interaction Toolkit  2018.4.99-07c45cb1
Medical Imaging Interaction Toolkit
mitkShowSegmentationAsSmoothedSurface.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 
15 #include "mitkImageCast.h"
16 #include "mitkImageToItk.h"
17 #include <itkAddImageFilter.h>
18 #include <itkBinaryMedianImageFilter.h>
19 #include <itkBinaryThresholdImageFilter.h>
20 #include <itkConnectedThresholdImageFilter.h>
21 #include <itkConstantPadImageFilter.h>
22 #include <itkDiscreteGaussianImageFilter.h>
23 #include <itkImageRegionIteratorWithIndex.h>
24 #include <itkMultiplyImageFilter.h>
25 #include <itkRegionOfInterestImageFilter.h>
26 #include <mitkGeometry3D.h>
27 #include <mitkImageTimeSelector.h>
29 #include <mitkProgressBar.h>
30 #include <mitkStatusBar.h>
31 #include <mitkUIDGenerator.h>
33 #include <vtkCleanPolyData.h>
34 #include <vtkPolyDataNormals.h>
35 #include <vtkQuadricDecimation.h>
36 
37 using namespace mitk;
38 using namespace std;
39 
40 ShowSegmentationAsSmoothedSurface::ShowSegmentationAsSmoothedSurface()
41 {
42 }
43 
44 ShowSegmentationAsSmoothedSurface::~ShowSegmentationAsSmoothedSurface()
45 {
46 }
47 
49 {
50  Superclass::Initialize(other);
51 
52  bool syncVisibility = false;
53 
54  if (other != nullptr)
55  other->GetParameter("Sync visibility", syncVisibility);
56 
57  SetParameter("Sync visibility", syncVisibility);
58  SetParameter("Wireframe", false);
59 
60  // The Smoothing value is used as variance for a Gauss filter.
61  // A reasonable default value equals the image spacing in mm.
62  SetParameter("Smoothing", 1.0);
63 
64  // Valid range for decimation value is [0, 1). High values
65  // increase decimation, especially when very close to 1.
66  // A value of 0 disables decimation.
67  SetParameter("Decimation", 0.5);
68 
69  // Valid range for closing value is [0, 1]. Higher values
70  // increase closing. A value of 0 disables closing.
71  SetParameter("Closing", 0.0);
72 }
73 
75 {
76  try
77  {
79  GetPointerParameter("Input", image);
80 
81  return image.IsNotNull() && GetGroupNode();
82  }
83  catch (const invalid_argument &)
84  {
85  return false;
86  }
87 }
88 
90 {
92  GetPointerParameter("Input", image);
93 
94  double smoothing;
95  GetParameter("Smoothing", smoothing);
96 
97  double decimation;
98  GetParameter("Decimation", decimation);
99 
100  double closing;
101  GetParameter("Closing", closing);
102 
103  int timeNr = 0;
104  GetParameter("TimeNr", timeNr);
105 
106  if (image->GetDimension() == 4)
107  MITK_INFO << "CREATING SMOOTHED POLYGON MODEL (t = " << timeNr << ')';
108  else
109  MITK_INFO << "CREATING SMOOTHED POLYGON MODEL";
110 
111  MITK_INFO << " Smoothing = " << smoothing;
112  MITK_INFO << " Decimation = " << decimation;
113  MITK_INFO << " Closing = " << closing;
114 
115  Geometry3D::Pointer geometry = dynamic_cast<Geometry3D *>(image->GetGeometry()->Clone().GetPointer());
116 
117  // Make ITK image out of MITK image
118 
119  typedef itk::Image<unsigned char, 3> CharImageType;
120  typedef itk::Image<unsigned short, 3> ShortImageType;
121  typedef itk::Image<float, 3> FloatImageType;
122 
123  if (image->GetDimension() == 4)
124  {
126  imageTimeSelector->SetInput(image);
127  imageTimeSelector->SetTimeNr(timeNr);
128  imageTimeSelector->UpdateLargestPossibleRegion();
129  image = imageTimeSelector->GetOutput(0);
130  }
131 
133 
134  try
135  {
136  imageToItkFilter->SetInput(image);
137  }
138  catch (const itk::ExceptionObject &e)
139  {
140  // Most probably the input image type is wrong. Binary images are expected to be
141  // >unsigned< char images.
142  MITK_ERROR << e.GetDescription() << endl;
143  return false;
144  }
145 
146  imageToItkFilter->Update();
147 
148  CharImageType::Pointer itkImage = imageToItkFilter->GetOutput();
149 
150  // Get bounding box and relabel
151 
152  MITK_INFO << "Extracting VOI...";
153 
154  int imageLabel = 1;
155  bool roiFound = false;
156 
157  CharImageType::IndexType minIndex;
159 
160  CharImageType::IndexType maxIndex;
162 
163  itk::ImageRegionIteratorWithIndex<CharImageType> iter(itkImage, itkImage->GetLargestPossibleRegion());
164 
165  for (iter.GoToBegin(); !iter.IsAtEnd(); ++iter)
166  {
167  if (iter.Get() == imageLabel)
168  {
169  roiFound = true;
170  iter.Set(1);
171 
172  CharImageType::IndexType currentIndex = iter.GetIndex();
173 
174  for (unsigned int dim = 0; dim < 3; ++dim)
175  {
176  minIndex[dim] = min(currentIndex[dim], minIndex[dim]);
177  maxIndex[dim] = max(currentIndex[dim], maxIndex[dim]);
178  }
179  }
180  else
181  {
182  iter.Set(0);
183  }
184  }
185 
186  if (!roiFound)
187  {
189  MITK_ERROR << "Didn't found segmentation labeled with " << imageLabel << "!" << endl;
190  return false;
191  }
192 
194 
195  // Extract and pad bounding box
196 
197  typedef itk::RegionOfInterestImageFilter<CharImageType, CharImageType> ROIFilterType;
198 
199  ROIFilterType::Pointer roiFilter = ROIFilterType::New();
200  CharImageType::RegionType region;
201  CharImageType::SizeType size;
202 
203  for (unsigned int dim = 0; dim < 3; ++dim)
204  {
205  size[dim] = maxIndex[dim] - minIndex[dim] + 1;
206  }
207 
208  region.SetIndex(minIndex);
209  region.SetSize(size);
210 
211  roiFilter->SetInput(itkImage);
212  roiFilter->SetRegionOfInterest(region);
213  roiFilter->ReleaseDataFlagOn();
214  roiFilter->ReleaseDataBeforeUpdateFlagOn();
215 
216  typedef itk::ConstantPadImageFilter<CharImageType, CharImageType> PadFilterType;
217 
218  PadFilterType::Pointer padFilter = PadFilterType::New();
219  const PadFilterType::SizeValueType pad[3] = {10, 10, 10};
220 
221  padFilter->SetInput(roiFilter->GetOutput());
222  padFilter->SetConstant(0);
223  padFilter->SetPadLowerBound(pad);
224  padFilter->SetPadUpperBound(pad);
225  padFilter->ReleaseDataFlagOn();
226  padFilter->ReleaseDataBeforeUpdateFlagOn();
227  padFilter->Update();
228 
229  CharImageType::Pointer roiImage = padFilter->GetOutput();
230 
231  roiImage->DisconnectPipeline();
232  roiFilter = nullptr;
233  padFilter = nullptr;
234 
235  // Correct origin of real geometry (changed by cropping and padding)
236 
237  typedef Geometry3D::TransformType TransformType;
238 
239  TransformType::Pointer transform = TransformType::New();
240  TransformType::OutputVectorType translation;
241 
242  for (unsigned int dim = 0; dim < 3; ++dim)
243  translation[dim] = (int)minIndex[dim] - (int)pad[dim];
244 
245  transform->SetIdentity();
246  transform->Translate(translation);
247  geometry->Compose(transform, true);
248 
250 
251  // Median
252 
253  MITK_INFO << "Median...";
254 
255  typedef itk::BinaryMedianImageFilter<CharImageType, CharImageType> MedianFilterType;
256 
257  MedianFilterType::Pointer medianFilter = MedianFilterType::New();
258  CharImageType::SizeType radius = {{0}};
259 
260  medianFilter->SetRadius(radius);
261  medianFilter->SetBackgroundValue(0);
262  medianFilter->SetForegroundValue(1);
263  medianFilter->SetInput(roiImage);
264  medianFilter->ReleaseDataFlagOn();
265  medianFilter->ReleaseDataBeforeUpdateFlagOn();
266  medianFilter->Update();
267 
269 
270  // Intelligent closing
271 
272  MITK_INFO << "Intelligent closing...";
273 
274  auto surfaceRatio = (unsigned int)((1.0f - closing) * 100.0f);
275 
277 
278  ClosingFilterType::Pointer closingFilter = ClosingFilterType::New();
279 
280  closingFilter->SetInput(medianFilter->GetOutput());
281  closingFilter->ReleaseDataFlagOn();
282  closingFilter->ReleaseDataBeforeUpdateFlagOn();
283  closingFilter->SetSurfaceRatio(surfaceRatio);
284  closingFilter->Update();
285 
286  ShortImageType::Pointer closedImage = closingFilter->GetOutput();
287 
288  closedImage->DisconnectPipeline();
289  roiImage = nullptr;
290  medianFilter = nullptr;
291  closingFilter = nullptr;
292 
294 
295  // Gaussian blur
296 
297  MITK_INFO << "Gauss...";
298 
299  typedef itk::BinaryThresholdImageFilter<ShortImageType, FloatImageType> BinaryThresholdToFloatFilterType;
300 
301  BinaryThresholdToFloatFilterType::Pointer binThresToFloatFilter = BinaryThresholdToFloatFilterType::New();
302 
303  binThresToFloatFilter->SetInput(closedImage);
304  binThresToFloatFilter->SetLowerThreshold(1);
305  binThresToFloatFilter->SetUpperThreshold(1);
306  binThresToFloatFilter->SetInsideValue(100);
307  binThresToFloatFilter->SetOutsideValue(0);
308  binThresToFloatFilter->ReleaseDataFlagOn();
309  binThresToFloatFilter->ReleaseDataBeforeUpdateFlagOn();
310 
311  typedef itk::DiscreteGaussianImageFilter<FloatImageType, FloatImageType> GaussianFilterType;
312 
313  // From the following line on, IntelliSense (VS 2008) is broken. Any idea how to fix it?
314  GaussianFilterType::Pointer gaussFilter = GaussianFilterType::New();
315 
316  gaussFilter->SetInput(binThresToFloatFilter->GetOutput());
317  gaussFilter->SetUseImageSpacing(true);
318  gaussFilter->SetVariance(smoothing);
319  gaussFilter->ReleaseDataFlagOn();
320  gaussFilter->ReleaseDataBeforeUpdateFlagOn();
321 
322  typedef itk::BinaryThresholdImageFilter<FloatImageType, CharImageType> BinaryThresholdFromFloatFilterType;
323 
324  BinaryThresholdFromFloatFilterType::Pointer binThresFromFloatFilter = BinaryThresholdFromFloatFilterType::New();
325 
326  binThresFromFloatFilter->SetInput(gaussFilter->GetOutput());
327  binThresFromFloatFilter->SetLowerThreshold(50);
328  binThresFromFloatFilter->SetUpperThreshold(255);
329  binThresFromFloatFilter->SetInsideValue(1);
330  binThresFromFloatFilter->SetOutsideValue(0);
331  binThresFromFloatFilter->ReleaseDataFlagOn();
332  binThresFromFloatFilter->ReleaseDataBeforeUpdateFlagOn();
333  binThresFromFloatFilter->Update();
334 
335  CharImageType::Pointer blurredImage = binThresFromFloatFilter->GetOutput();
336 
337  blurredImage->DisconnectPipeline();
338  closedImage = nullptr;
339  binThresToFloatFilter = nullptr;
340  gaussFilter = nullptr;
341 
343 
344  // Fill holes
345 
346  MITK_INFO << "Filling cavities...";
347 
348  typedef itk::ConnectedThresholdImageFilter<CharImageType, CharImageType> ConnectedThresholdFilterType;
349 
350  ConnectedThresholdFilterType::Pointer connectedThresFilter = ConnectedThresholdFilterType::New();
351 
352  CharImageType::IndexType corner;
353 
354  corner[0] = 0;
355  corner[1] = 0;
356  corner[2] = 0;
357 
358  connectedThresFilter->SetInput(blurredImage);
359  connectedThresFilter->SetSeed(corner);
360  connectedThresFilter->SetLower(0);
361  connectedThresFilter->SetUpper(0);
362  connectedThresFilter->SetReplaceValue(2);
363  connectedThresFilter->ReleaseDataFlagOn();
364  connectedThresFilter->ReleaseDataBeforeUpdateFlagOn();
365 
366  typedef itk::BinaryThresholdImageFilter<CharImageType, CharImageType> BinaryThresholdFilterType;
367 
368  BinaryThresholdFilterType::Pointer binThresFilter = BinaryThresholdFilterType::New();
369 
370  binThresFilter->SetInput(connectedThresFilter->GetOutput());
371  binThresFilter->SetLowerThreshold(0);
372  binThresFilter->SetUpperThreshold(0);
373  binThresFilter->SetInsideValue(50);
374  binThresFilter->SetOutsideValue(0);
375  binThresFilter->ReleaseDataFlagOn();
376  binThresFilter->ReleaseDataBeforeUpdateFlagOn();
377 
378  typedef itk::AddImageFilter<CharImageType, CharImageType, CharImageType> AddFilterType;
379 
380  AddFilterType::Pointer addFilter = AddFilterType::New();
381 
382  addFilter->SetInput1(blurredImage);
383  addFilter->SetInput2(binThresFilter->GetOutput());
384  addFilter->ReleaseDataFlagOn();
385  addFilter->ReleaseDataBeforeUpdateFlagOn();
386  addFilter->Update();
387 
389 
390  // Surface extraction
391 
392  MITK_INFO << "Surface extraction...";
393 
394  Image::Pointer filteredImage = Image::New();
395  CastToMitkImage(addFilter->GetOutput(), filteredImage);
396 
397  filteredImage->SetGeometry(geometry);
398 
400 
401  imageToSurfaceFilter->SetInput(filteredImage);
402  imageToSurfaceFilter->SetThreshold(50);
403  imageToSurfaceFilter->SmoothOn();
404  imageToSurfaceFilter->SetDecimate(ImageToSurfaceFilter::NoDecimation);
405 
406  m_Surface = imageToSurfaceFilter->GetOutput(0);
407 
409 
410  // Mesh decimation
411 
412  if (decimation > 0.0f && decimation < 1.0f)
413  {
414  MITK_INFO << "Quadric mesh decimation...";
415 
416  vtkQuadricDecimation *quadricDecimation = vtkQuadricDecimation::New();
417  quadricDecimation->SetInputData(m_Surface->GetVtkPolyData());
418  quadricDecimation->SetTargetReduction(decimation);
419  quadricDecimation->AttributeErrorMetricOn();
420  quadricDecimation->GlobalWarningDisplayOff();
421  quadricDecimation->Update();
422 
423  vtkCleanPolyData *cleaner = vtkCleanPolyData::New();
424  cleaner->SetInputConnection(quadricDecimation->GetOutputPort());
425  cleaner->PieceInvariantOn();
426  cleaner->ConvertLinesToPointsOn();
427  cleaner->ConvertStripsToPolysOn();
428  cleaner->PointMergingOn();
429  cleaner->Update();
430 
431  m_Surface->SetVtkPolyData(cleaner->GetOutput());
432  }
433 
435 
436  // Compute Normals
437 
438  vtkPolyDataNormals *computeNormals = vtkPolyDataNormals::New();
439  computeNormals->SetInputData(m_Surface->GetVtkPolyData());
440  computeNormals->SetFeatureAngle(360.0f);
441  computeNormals->AutoOrientNormalsOn();
442  computeNormals->FlipNormalsOff();
443  computeNormals->Update();
444 
445  m_Surface->SetVtkPolyData(computeNormals->GetOutput());
446 
447  return true;
448 }
449 
451 {
453 
454  bool wireframe = false;
455  GetParameter("Wireframe", wireframe);
456 
457  if (wireframe)
458  {
459  auto *representation =
460  dynamic_cast<VtkRepresentationProperty *>(node->GetProperty("material.representation"));
461 
462  if (representation != nullptr)
463  representation->SetRepresentationToWireframe();
464  }
465 
466  node->SetProperty("opacity", FloatProperty::New(1.0));
467  node->SetProperty("line width", IntProperty::New(1));
468  node->SetProperty("scalar visibility", BoolProperty::New(false));
469 
470  std::string groupNodeName = "surface";
471  DataNode *groupNode = GetGroupNode();
472 
473  if (groupNode != nullptr)
474  groupNode->GetName(groupNodeName);
475 
476  node->SetProperty("name", StringProperty::New(groupNodeName));
477  node->SetData(m_Surface);
478 
479  BaseProperty *colorProperty = groupNode->GetProperty("color");
480 
481  if (colorProperty != nullptr)
482  node->ReplaceProperty("color", colorProperty->Clone());
483  else
484  node->SetProperty("color", ColorProperty::New(1.0f, 0.0f, 0.0f));
485 
486  bool showResult = true;
487  GetParameter("Show result", showResult);
488 
489  bool syncVisibility = false;
490  GetParameter("Sync visibility", syncVisibility);
491 
493  GetPointerParameter("Input", image);
494 
495  BaseProperty *organTypeProperty = image->GetProperty("organ type");
496 
497  if (organTypeProperty != nullptr)
498  m_Surface->SetProperty("organ type", organTypeProperty);
499 
500  BaseProperty *visibleProperty = groupNode->GetProperty("visible");
501 
502  if (visibleProperty != nullptr && syncVisibility)
503  node->ReplaceProperty("visible", visibleProperty->Clone());
504  else
505  node->SetProperty("visible", BoolProperty::New(showResult));
506 
507  InsertBelowGroupNode(node);
508 
509  Superclass::ThreadedUpdateSuccessful();
510 }
mitk::BaseProperty * GetProperty(const char *propertyKey, const mitk::BaseRenderer *renderer=nullptr, bool fallBackOnDataProperties=true) const
Get the property (instance of BaseProperty) with key propertyKey from the PropertyList of the rendere...
void Progress(unsigned int steps=1)
Sets the current amount of progress to current progress + steps.
Pointer Clone() const
Standard implementation of BaseGeometry.
#define MITK_INFO
Definition: mitkLogMacros.h:18
#define MITK_ERROR
Definition: mitkLogMacros.h:20
static Pointer New()
GeometryTransformHolder::TransformType TransformType
STL namespace.
DataCollection - Class to facilitate loading/accessing structured data.
itk::Image< double, 3 > FloatImageType
Definition: CLBrainMask.cpp:31
static ProgressBar * GetInstance()
static method to get the GUI dependent ProgressBar-instance so the methods for steps to do and progre...
void GetParameter(const char *parameter, T &value) const
static Pointer New()
Abstract base class for properties.
static Pointer New()
static T max(T x, T y)
Definition: svm.cpp:56
mitk::Image::Pointer image
static Pointer New()
static Pointer New()
static T min(T x, T y)
Definition: svm.cpp:53
static Pointer New()
static Pointer New()
void CastToMitkImage(const itk::SmartPointer< ItkOutputImageType > &itkimage, itk::SmartPointer< mitk::Image > &mitkoutputimage)
Cast an itk::Image (with a specific type) to an mitk::Image.
Definition: mitkImageCast.h:74
void Initialize(const NonBlockingAlgorithm *other=nullptr) override
static Pointer New()
Class for nodes of the DataTree.
Definition: mitkDataNode.h:57
static Pointer New()
bool GetName(std::string &nodeName, const mitk::BaseRenderer *renderer=nullptr, const char *propertyKey="name") const
Convenience access method for accessing the name of an object (instance of StringProperty with proper...
Definition: mitkDataNode.h:369