Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
mitkLabelSetImageToSurfaceFilter.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 <mitkImageAccessByItk.h>
20 #include <mitkImageCast.h>
21 
22 // itk
23 #include <itkAntiAliasBinaryImageFilter.h>
24 #include <itkAutoCropLabelMapFilter.h>
25 #include <itkBinaryThresholdImageFilter.h>
26 #include <itkLabelImageToLabelMapFilter.h>
27 #include <itkLabelMap.h>
28 #include <itkLabelMapToLabelImageFilter.h>
29 #include <itkLabelObject.h>
30 #include <itkNumericTraits.h>
31 #include <itkSmoothingRecursiveGaussianImageFilter.h>
32 
33 // vtk
34 #include <vtkCleanPolyData.h>
35 #include <vtkImageChangeInformation.h>
36 #include <vtkImageData.h>
37 #include <vtkLinearTransform.h>
38 #include <vtkMarchingCubes.h>
39 #include <vtkSmartPointer.h>
40 
42  : m_GenerateAllLabels(false), m_RequestedLabel(1), m_BackgroundLabel(0), m_UseSmoothing(0), m_Sigma(0.1)
43 {
44 }
45 
47 {
48 }
49 
51 {
52  // Process object is not const-correct so the const_cast is required here
53  this->ProcessObject::SetNthInput(0, const_cast<mitk::Image *>(image));
54 }
55 /*
56 void mitk::LabelSetImageToSurfaceFilter::SetObserver(mitk::ProcessObserver::Pointer observer)
57 {
58  m_Observer = observer;
59 }
60 */
62 {
63  if (this->GetNumberOfInputs() < 1)
64  {
65  return nullptr;
66  }
67 
68  return static_cast<const mitk::Image *>(this->ProcessObject::GetInput(0));
69 }
70 
72 {
73  itkDebugMacro(<< "GenerateOutputInformation()");
74 }
75 
77 {
78  Image::ConstPointer inputImage = this->GetInput();
79  if (inputImage.IsNull())
80  return;
81 
82  mitk::Surface *outputSurface = this->GetOutput();
83  if (!outputSurface)
84  return;
85 
86  AccessFixedDimensionByItk_1(inputImage, InternalProcessing, 3, outputSurface);
87 }
88 
89 template <typename TPixel, unsigned int VDimension>
90 void mitk::LabelSetImageToSurfaceFilter::InternalProcessing(const itk::Image<TPixel, VDimension> *input,
91  mitk::Surface * /*surface*/)
92 {
93  typedef itk::Image<TPixel, VDimension> ImageType;
94 
95  typedef itk::BinaryThresholdImageFilter<ImageType, ImageType> BinaryThresholdFilterType;
96  typedef itk::LabelObject<TPixel, VDimension> LabelObjectType;
97  typedef itk::LabelMap<LabelObjectType> LabelMapType;
98  typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> Image2LabelMapType;
99  typedef itk::AutoCropLabelMapFilter<LabelMapType> AutoCropType;
100  typedef itk::LabelMapToLabelImageFilter<LabelMapType, ImageType> LabelMap2ImageType;
101 
102  typedef itk::Image<float, VDimension> RealImageType;
103 
104  typedef itk::AntiAliasBinaryImageFilter<ImageType, RealImageType> AntiAliasFilterType;
105  typedef itk::SmoothingRecursiveGaussianImageFilter<RealImageType, RealImageType> GaussianFilterType;
106 
108  thresholdFilter->SetInput(input);
109  thresholdFilter->SetLowerThreshold(m_RequestedLabel);
110  thresholdFilter->SetUpperThreshold(m_RequestedLabel);
111  thresholdFilter->SetOutsideValue(0);
112  thresholdFilter->SetInsideValue(1);
113  // thresholdFilter->ReleaseDataFlagOn();
114  thresholdFilter->Update();
115 
116  typename Image2LabelMapType::Pointer image2label = Image2LabelMapType::New();
117  image2label->SetInput(thresholdFilter->GetOutput());
118 
119  typename AutoCropType::SizeType border;
120  border[0] = 3;
121  border[1] = 3;
122  border[2] = 3;
123 
124  typename AutoCropType::Pointer autoCropFilter = AutoCropType::New();
125  autoCropFilter->SetInput(image2label->GetOutput());
126  autoCropFilter->SetCropBorder(border);
127  autoCropFilter->InPlaceOn();
128 
129  typename LabelMap2ImageType::Pointer label2image = LabelMap2ImageType::New();
130  label2image->SetInput(autoCropFilter->GetOutput());
131 
132  label2image->Update();
133 
134  typename AntiAliasFilterType::Pointer antiAliasFilter = AntiAliasFilterType::New();
135  antiAliasFilter->SetInput(label2image->GetOutput());
136  antiAliasFilter->SetMaximumRMSError(0.001);
137  antiAliasFilter->SetNumberOfLayers(3);
138  antiAliasFilter->SetUseImageSpacing(false);
139  antiAliasFilter->SetNumberOfIterations(40);
140 
141  antiAliasFilter->Update();
142 
143  typename RealImageType::Pointer result;
144 
145  if (m_UseSmoothing)
146  {
147  typename GaussianFilterType::Pointer gaussianFilter = GaussianFilterType::New();
148  gaussianFilter->SetSigma(m_Sigma);
149  gaussianFilter->SetInput(antiAliasFilter->GetOutput());
150  gaussianFilter->Update();
151  result = gaussianFilter->GetOutput();
152  }
153  else
154  {
155  result = antiAliasFilter->GetOutput();
156  }
157 
158  result->DisconnectPipeline();
159 
160  typename ImageType::RegionType cropRegion;
161  cropRegion = autoCropFilter->GetOutput()->GetLargestPossibleRegion();
162 
163  const typename ImageType::IndexType &cropIndex = cropRegion.GetIndex();
164 
165  m_ResultImage = mitk::Image::New();
166  mitk::CastToMitkImage(result, m_ResultImage);
167 
168  mitk::BaseGeometry *newGeometry = m_ResultImage->GetSlicedGeometry();
169  mitk::Point3D origin;
170  vtk2itk(cropIndex, origin);
171  this->GetInput()->GetGeometry()->IndexToWorld(origin, origin);
172  newGeometry->SetOrigin(origin);
173 
174  vtkImageData *vtkimage = const_cast<vtkImageData *>(m_ResultImage->GetVtkImageData(0));
175 
176  vtkSmartPointer<vtkImageChangeInformation> indexCoordinatesImageFilter =
178  indexCoordinatesImageFilter->SetInputData(vtkimage);
179  indexCoordinatesImageFilter->SetOutputOrigin(0.0, 0.0, 0.0);
180 
181  vtkSmartPointer<vtkMarchingCubes> marching = vtkSmartPointer<vtkMarchingCubes>::New();
182  marching->ComputeScalarsOff();
183  marching->ComputeNormalsOn();
184  marching->ComputeGradientsOn();
185  marching->SetInputConnection(indexCoordinatesImageFilter->GetOutputPort());
186  marching->SetValue(0, 0.0);
187 
188  marching->Update();
189 
190  vtkPolyData *polydata = marching->GetOutput();
191 
192  if ((!polydata) || (!polydata->GetNumberOfPoints()))
193  throw itk::ExceptionObject(__FILE__, __LINE__, "marching cubes has failed.");
194 
195  mitk::Vector3D spacing = newGeometry->GetSpacing();
196 
197  vtkPoints *points = polydata->GetPoints();
198  vtkMatrix4x4 *vtkmatrix = vtkMatrix4x4::New();
199  newGeometry->GetVtkTransform()->GetMatrix(vtkmatrix);
200  double(*matrix)[4] = vtkmatrix->Element;
201 
202  for (int i = 0; i < 3; ++i)
203  for (int j = 0; j < 3; ++j)
204  matrix[i][j] /= spacing[j];
205 
206  unsigned int n = points->GetNumberOfPoints();
207  double point[3];
208 
209  for (unsigned int i = 0; i < n; i++)
210  {
211  points->GetPoint(i, point);
212  mitkVtkLinearTransformPoint(matrix, point, point);
213  points->SetPoint(i, point);
214  }
215  vtkmatrix->Delete();
216 
217  vtkSmartPointer<vtkCleanPolyData> cleanPolyDataFilter = vtkSmartPointer<vtkCleanPolyData>::New();
218  cleanPolyDataFilter->SetInputData(polydata);
219  cleanPolyDataFilter->PieceInvariantOff();
220  cleanPolyDataFilter->ConvertLinesToPointsOff();
221  cleanPolyDataFilter->ConvertPolysToLinesOff();
222  cleanPolyDataFilter->ConvertStripsToPolysOff();
223  cleanPolyDataFilter->PointMergingOn();
224  cleanPolyDataFilter->Update();
225 
226  mitk::Surface::Pointer output = this->GetOutput(0);
227  output->SetVtkPolyData(cleanPolyDataFilter->GetOutput(), 0);
228 }
Class for storing surfaces (vtkPolyData).
Definition: mitkSurface.h:32
itk::SmartPointer< Self > Pointer
void InternalProcessing(const itk::Image< TPixel, VImageDimension > *input, mitk::Surface *surface)
const mitk::Vector3D GetSpacing() const
Get the spacing (size of a pixel).
std::map< LabelType, unsigned long > LabelMapType
map::core::discrete::Elements< 3 >::InternalImageType ImageType
void SetOrigin(const Point3D &origin)
Set the origin, i.e. the upper-left corner of the plane.
Image class for storing images.
Definition: mitkImage.h:76
void vtk2itk(const Tin &in, Tout &out)
static Pointer New()
#define AccessFixedDimensionByItk_1(mitkImage, itkImageTypeFunction, dimension, arg1)
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:78
vtkLinearTransform * GetVtkTransform() const
Get the m_IndexToWorldTransform as a vtkLinearTransform.
virtual void SetInput(const mitk::Image *image)
BaseGeometry Describes the geometry of a data object.
static itkEventMacro(BoundingShapeInteractionEvent, itk::AnyEvent) class MITKBOUNDINGSHAPE_EXPORT BoundingShapeInteractor Pointer New()
Basic interaction methods for mitk::GeometryData.