Medical Imaging Interaction Toolkit  2018.4.99-36d69b77
Medical Imaging Interaction Toolkit
mitkPlanarFigureMaskGenerator.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 #include <mitkBaseGeometry.h>
15 #include <mitkITKImageImport.h>
16 #include "mitkImageAccessByItk.h"
17 #include <mitkExtractImageFilter.h>
19 #include <mitkImageTimeSelector.h>
20 #include <mitkIOUtil.h>
21 
22 #include <itkCastImageFilter.h>
23 #include <itkVTKImageExport.h>
24 #include <itkVTKImageImport.h>
25 #include <itkImageDuplicator.h>
26 #include <itkExceptionObject.h>
27 #include <itkLineIterator.h>
28 
29 #include <vtkPoints.h>
30 #include <vtkImageStencil.h>
31 #include <vtkImageImport.h>
32 #include <vtkImageExport.h>
33 #include <vtkLassoStencilSource.h>
34 #include <vtkSmartPointer.h>
35 
36 
37 
38 namespace mitk
39 {
40 
42 {
43  if ( planarFigure.IsNull() )
44  {
45  throw std::runtime_error( "Error: planar figure empty!" );
46  }
47 
48  const PlaneGeometry *planarFigurePlaneGeometry = planarFigure->GetPlaneGeometry();
49  if ( planarFigurePlaneGeometry == nullptr )
50  {
51  throw std::runtime_error( "Planar-Figure not yet initialized!" );
52  }
53 
54  const auto *planarFigureGeometry =
55  dynamic_cast< const PlaneGeometry * >( planarFigurePlaneGeometry );
56  if ( planarFigureGeometry == nullptr )
57  {
58  throw std::runtime_error( "Non-planar planar figures not supported!" );
59  }
60 
61  if (planarFigure != m_PlanarFigure)
62  {
63  this->Modified();
64  m_PlanarFigure = planarFigure;
65  }
66 
67 }
68 
70 {
71  if (IsUpdateRequired())
72  {
73  this->CalculateMask();
74  }
75  return m_ReferenceImage;
76 }
77 
78 template < typename TPixel, unsigned int VImageDimension >
79 void PlanarFigureMaskGenerator::InternalCalculateMaskFromPlanarFigure(
80  const itk::Image< TPixel, VImageDimension > *image, unsigned int axis )
81 {
82  typedef itk::Image< unsigned short, 2 > MaskImage2DType;
83 
84  typename MaskImage2DType::Pointer maskImage = MaskImage2DType::New();
85  maskImage->SetOrigin(image->GetOrigin());
86  maskImage->SetSpacing(image->GetSpacing());
87  maskImage->SetLargestPossibleRegion(image->GetLargestPossibleRegion());
88  maskImage->SetBufferedRegion(image->GetBufferedRegion());
89  maskImage->SetDirection(image->GetDirection());
90  maskImage->SetNumberOfComponentsPerPixel(image->GetNumberOfComponentsPerPixel());
91  maskImage->Allocate();
92  maskImage->FillBuffer(1);
93 
94  // all PolylinePoints of the PlanarFigure are stored in a vtkPoints object.
95  // These points are used by the vtkLassoStencilSource to create
96  // a vtkImageStencil.
97  const mitk::PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry();
98  const typename PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 );
99  const mitk::BaseGeometry *imageGeometry3D = m_inputImage->GetGeometry( 0 );
100  // If there is a second poly line in a closed planar figure, treat it as a hole.
101  PlanarFigure::PolyLineType planarFigureHolePolyline;
102 
103  if (m_PlanarFigure->GetPolyLinesSize() == 2)
104  planarFigureHolePolyline = m_PlanarFigure->GetPolyLine(1);
105 
106 
107  // Determine x- and y-dimensions depending on principal axis
108  // TODO use plane geometry normal to determine that automatically, then check whether the PF is aligned with one of the three principal axis
109  int i0, i1;
110  switch ( axis )
111  {
112  case 0:
113  i0 = 1;
114  i1 = 2;
115  break;
116 
117  case 1:
118  i0 = 0;
119  i1 = 2;
120  break;
121 
122  case 2:
123  default:
124  i0 = 0;
125  i1 = 1;
126  break;
127  }
128 
129  // store the polyline contour as vtkPoints object
130  vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
131  for (const auto& point : planarFigurePolyline)
132  {
133  Point3D point3D;
134 
135  // Convert 2D point back to the local index coordinates of the selected image
136  planarFigurePlaneGeometry->Map(point, point3D);
137  imageGeometry3D->WorldToIndex(point3D, point3D);
138 
139  points->InsertNextPoint(point3D[i0], point3D[i1], 0);
140  }
141 
142  vtkSmartPointer<vtkPoints> holePoints;
143 
144  if (!planarFigureHolePolyline.empty())
145  {
146  holePoints = vtkSmartPointer<vtkPoints>::New();
147  Point3D point3D;
148 
149  for (const auto& point : planarFigureHolePolyline)
150  {
151  planarFigurePlaneGeometry->Map(point, point3D);
152  imageGeometry3D->WorldToIndex(point3D, point3D);
153  holePoints->InsertNextPoint(point3D[i0], point3D[i1], 0);
154  }
155  }
156 
157  // mark a malformed 2D planar figure ( i.e. area = 0 ) as out of bounds
158  // this can happen when all control points of a rectangle lie on the same line = two of the three extents are zero
159  double bounds[6] = {0};
160  points->GetBounds(bounds);
161  bool extent_x = (fabs(bounds[0] - bounds[1])) < mitk::eps;
162  bool extent_y = (fabs(bounds[2] - bounds[3])) < mitk::eps;
163  bool extent_z = (fabs(bounds[4] - bounds[5])) < mitk::eps;
164 
165  // throw an exception if a closed planar figure is deformed, i.e. has only one non-zero extent
166  if (m_PlanarFigure->IsClosed() && ((extent_x && extent_y) || (extent_x && extent_z) || (extent_y && extent_z)))
167  {
168  mitkThrow() << "Figure has a zero area and cannot be used for masking.";
169  }
170 
171  // create a vtkLassoStencilSource and set the points of the Polygon
172  vtkSmartPointer<vtkLassoStencilSource> lassoStencil = vtkSmartPointer<vtkLassoStencilSource>::New();
173  lassoStencil->SetShapeToPolygon();
174  lassoStencil->SetPoints(points);
175 
176  vtkSmartPointer<vtkLassoStencilSource> holeLassoStencil = nullptr;
177 
178  if (holePoints.GetPointer() != nullptr)
179  {
180  holeLassoStencil = vtkSmartPointer<vtkLassoStencilSource>::New();
181  holeLassoStencil->SetShapeToPolygon();
182  holeLassoStencil->SetPoints(holePoints);
183  }
184 
185  // Export from ITK to VTK (to use a VTK filter)
186  typedef itk::VTKImageImport< MaskImage2DType > ImageImportType;
187  typedef itk::VTKImageExport< MaskImage2DType > ImageExportType;
188 
189  typename ImageExportType::Pointer itkExporter = ImageExportType::New();
190  itkExporter->SetInput( maskImage );
191 // itkExporter->SetInput( castFilter->GetOutput() );
192 
193  vtkSmartPointer<vtkImageImport> vtkImporter = vtkSmartPointer<vtkImageImport>::New();
194  this->ConnectPipelines( itkExporter, vtkImporter );
195 
196  // Apply the generated image stencil to the input image
197  vtkSmartPointer<vtkImageStencil> imageStencilFilter = vtkSmartPointer<vtkImageStencil>::New();
198  imageStencilFilter->SetInputConnection( vtkImporter->GetOutputPort() );
199  imageStencilFilter->SetStencilConnection(lassoStencil->GetOutputPort());
200  imageStencilFilter->ReverseStencilOff();
201  imageStencilFilter->SetBackgroundValue( 0 );
202  imageStencilFilter->Update();
203 
204  vtkSmartPointer<vtkImageStencil> holeStencilFilter = nullptr;
205 
206  if (holeLassoStencil.GetPointer() != nullptr)
207  {
208  holeStencilFilter = vtkSmartPointer<vtkImageStencil>::New();
209  holeStencilFilter->SetInputConnection(imageStencilFilter->GetOutputPort());
210  holeStencilFilter->SetStencilConnection(holeLassoStencil->GetOutputPort());
211  holeStencilFilter->ReverseStencilOn();
212  holeStencilFilter->SetBackgroundValue(0);
213  holeStencilFilter->Update();
214  }
215 
216  // Export from VTK back to ITK
217  vtkSmartPointer<vtkImageExport> vtkExporter = vtkSmartPointer<vtkImageExport>::New();
218  vtkExporter->SetInputConnection( holeStencilFilter.GetPointer() == nullptr
219  ? imageStencilFilter->GetOutputPort()
220  : holeStencilFilter->GetOutputPort());
221  vtkExporter->Update();
222 
223  typename ImageImportType::Pointer itkImporter = ImageImportType::New();
224  this->ConnectPipelines( vtkExporter, itkImporter );
225  itkImporter->Update();
226 
227  typedef itk::ImageDuplicator< ImageImportType::OutputImageType > DuplicatorType;
228  DuplicatorType::Pointer duplicator = DuplicatorType::New();
229  duplicator->SetInputImage( itkImporter->GetOutput() );
230  duplicator->Update();
231 
232  // Store mask
233  m_InternalITKImageMask2D = duplicator->GetOutput();
234 }
235 
236 template < typename TPixel, unsigned int VImageDimension >
237 void PlanarFigureMaskGenerator::InternalCalculateMaskFromOpenPlanarFigure(
238  const itk::Image< TPixel, VImageDimension > *image, unsigned int axis )
239 {
240  typedef itk::Image< unsigned short, 2 > MaskImage2DType;
241  typedef itk::LineIterator< MaskImage2DType > LineIteratorType;
242  typedef MaskImage2DType::IndexType IndexType2D;
243  typedef std::vector< IndexType2D > IndexVecType;
244 
245  typename MaskImage2DType::Pointer maskImage = MaskImage2DType::New();
246  maskImage->SetOrigin(image->GetOrigin());
247  maskImage->SetSpacing(image->GetSpacing());
248  maskImage->SetLargestPossibleRegion(image->GetLargestPossibleRegion());
249  maskImage->SetBufferedRegion(image->GetBufferedRegion());
250  maskImage->SetDirection(image->GetDirection());
251  maskImage->SetNumberOfComponentsPerPixel(image->GetNumberOfComponentsPerPixel());
252  maskImage->Allocate();
253  maskImage->FillBuffer(0);
254 
255  // all PolylinePoints of the PlanarFigure are stored in a vtkPoints object.
256  const mitk::PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry();
257  const typename PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 );
258  const mitk::BaseGeometry *imageGeometry3D = m_inputImage->GetGeometry( 0 );
259 
260  // Determine x- and y-dimensions depending on principal axis
261  // TODO use plane geometry normal to determine that automatically, then check whether the PF is aligned with one of the three principal axis
262  int i0, i1;
263  switch ( axis )
264  {
265  case 0:
266  i0 = 1;
267  i1 = 2;
268  break;
269 
270  case 1:
271  i0 = 0;
272  i1 = 2;
273  break;
274 
275  case 2:
276  default:
277  i0 = 0;
278  i1 = 1;
279  break;
280  }
281 
282  int numPolyLines = m_PlanarFigure->GetPolyLinesSize();
283  for ( int lineId = 0; lineId < numPolyLines; ++lineId )
284  {
285  // store the polyline contour as vtkPoints object
286  IndexVecType pointIndices;
287  for(const auto& point : planarFigurePolyline)
288  {
289  Point3D point3D;
290 
291  planarFigurePlaneGeometry->Map(point, point3D);
292  imageGeometry3D->WorldToIndex(point3D, point3D);
293 
294  IndexType2D index2D;
295  index2D[0] = point3D[i0];
296  index2D[1] = point3D[i1];
297 
298  pointIndices.push_back( index2D );
299  }
300 
301  size_t numLineSegments = pointIndices.size() - 1;
302  for (size_t i = 0; i < numLineSegments; ++i)
303  {
304  LineIteratorType lineIt(maskImage, pointIndices[i], pointIndices[i+1]);
305  while (!lineIt.IsAtEnd())
306  {
307  lineIt.Set(1);
308  ++lineIt;
309  }
310  }
311  }
312 
313  // Store mask
314  m_InternalITKImageMask2D = maskImage;
315 }
316 
318 {
319  if (!planarGeometry) return false;
320  if (!geometry) return false;
321 
322  unsigned int axis;
323  return GetPrincipalAxis(geometry,planarGeometry->GetNormal(), axis);
324 }
325 
326 bool PlanarFigureMaskGenerator::GetPrincipalAxis(
327  const BaseGeometry *geometry, Vector3D vector,
328  unsigned int &axis )
329 {
330  vector.Normalize();
331  for ( unsigned int i = 0; i < 3; ++i )
332  {
333  Vector3D axisVector = geometry->GetAxisVector( i );
334  axisVector.Normalize();
335 
336  //normal mitk::eps is to pedantic for this check. See e.g. T27122
337  //therefore choose a larger epsilon. The value was set a) as small as
338  //possible but b) still allowing to datasets like in (T27122) to pass
339  //when floating rounding errors sum up.
340  const double epsilon = 5e-5;
341  if ( fabs( fabs( axisVector * vector ) - 1.0) < epsilon)
342  {
343  axis = i;
344  return true;
345  }
346  }
347 
348  return false;
349 }
350 
351 void PlanarFigureMaskGenerator::CalculateMask()
352 {
353  if (m_inputImage.IsNull())
354  {
355  MITK_ERROR << "Image is not set.";
356  }
357 
358  if (m_PlanarFigure.IsNull())
359  {
360  MITK_ERROR << "PlanarFigure is not set.";
361  }
362 
363  if (m_TimeStep != 0)
364  {
365  MITK_WARN << "Multiple TimeSteps are not supported in PlanarFigureMaskGenerator (yet).";
366  }
367 
368  const BaseGeometry *imageGeometry = m_inputImage->GetGeometry();
369  if ( imageGeometry == nullptr )
370  {
371  throw std::runtime_error( "Image geometry invalid!" );
372  }
373 
374  if (m_inputImage->GetTimeSteps() > 0)
375  {
377  imgTimeSel->SetInput(m_inputImage);
378  imgTimeSel->SetTimeNr(m_TimeStep);
379  imgTimeSel->UpdateLargestPossibleRegion();
380  m_InternalTimeSliceImage = imgTimeSel->GetOutput();
381  }
382  else
383  {
384  m_InternalTimeSliceImage = m_inputImage;
385  }
386 
387  m_InternalITKImageMask2D = nullptr;
388  const PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry();
389  const auto *planarFigureGeometry = dynamic_cast< const PlaneGeometry * >( planarFigurePlaneGeometry );
390  //const BaseGeometry *imageGeometry = m_inputImage->GetGeometry();
391 
392  // Find principal direction of PlanarFigure in input image
393  unsigned int axis;
394  if ( !this->GetPrincipalAxis( imageGeometry,
395  planarFigureGeometry->GetNormal(), axis ) )
396  {
397  throw std::runtime_error( "Non-aligned planar figures not supported!" );
398  }
399  m_PlanarFigureAxis = axis;
400 
401  // Find slice number corresponding to PlanarFigure in input image
402  itk::Image< unsigned short, 3 >::IndexType index;
403  imageGeometry->WorldToIndex( planarFigureGeometry->GetOrigin(), index );
404 
405  unsigned int slice = index[axis];
406  m_PlanarFigureSlice = slice;
407 
408  // extract image slice which corresponds to the planarFigure and store it in m_InternalImageSlice
409  mitk::Image::ConstPointer inputImageSlice = extract2DImageSlice(axis, slice);
410  //mitk::IOUtil::Save(inputImageSlice, "/home/fabian/inputSliceImage.nrrd");
411  // Compute mask from PlanarFigure
412  // rastering for open planar figure:
413  if ( !m_PlanarFigure->IsClosed() )
414  {
415  AccessFixedDimensionByItk_1(inputImageSlice,
416  InternalCalculateMaskFromOpenPlanarFigure,
417  2, axis)
418  }
419  else//for closed planar figure
420  {
421  AccessFixedDimensionByItk_1(inputImageSlice,
422  InternalCalculateMaskFromPlanarFigure,
423  2, axis)
424  }
425 
426  //convert itk mask to mitk::Image::Pointer and return it
427  mitk::Image::Pointer planarFigureMaskImage;
428  planarFigureMaskImage = mitk::GrabItkImageMemory(m_InternalITKImageMask2D);
429  //mitk::IOUtil::Save(planarFigureMaskImage, "/home/fabian/planarFigureMaskImage.nrrd");
430 
431  //Convert2Dto3DImageFilter::Pointer sliceTo3DImageConverter = Convert2Dto3DImageFilter::New();
432  //sliceTo3DImageConverter->SetInput(planarFigureMaskImage);
433  //sliceTo3DImageConverter->Update();
434  //mitk::IOUtil::Save(sliceTo3DImageConverter->GetOutput(), "/home/fabian/3DsliceImage.nrrd");
435 
436  m_ReferenceImage = inputImageSlice;
437  //mitk::IOUtil::Save(m_ReferenceImage, "/home/fabian/referenceImage.nrrd");
438  m_InternalMask = planarFigureMaskImage;
439 }
440 
441 void PlanarFigureMaskGenerator::SetTimeStep(unsigned int timeStep)
442 {
443  if (timeStep != m_TimeStep)
444  {
445  m_TimeStep = timeStep;
446  }
447 }
448 
450 {
451  if (IsUpdateRequired())
452  {
453  this->CalculateMask();
454  this->Modified();
455  }
456 
457  m_InternalMaskUpdateTime = this->GetMTime();
458  return m_InternalMask;
459 }
460 
461 mitk::Image::ConstPointer PlanarFigureMaskGenerator::extract2DImageSlice(unsigned int axis, unsigned int slice)
462 {
463  // Extract slice with given position and direction from image
464  unsigned int dimension = m_InternalTimeSliceImage->GetDimension();
465 
466  if (dimension == 3)
467  {
469  imageExtractor->SetInput( m_InternalTimeSliceImage );
470  imageExtractor->SetSliceDimension( axis );
471  imageExtractor->SetSliceIndex( slice );
472  imageExtractor->Update();
473  return imageExtractor->GetOutput();
474  }
475  else if(dimension == 2)
476  {
477  return m_InternalTimeSliceImage;
478  }
479  else
480  {
481  MITK_ERROR << "Unsupported image dimension. Dimension is: " << dimension << ". Only 2D and 3D images are supported.";
482  return nullptr;
483  }
484 }
485 
486 bool PlanarFigureMaskGenerator::IsUpdateRequired() const
487 {
488  unsigned long thisClassTimeStamp = this->GetMTime();
489  unsigned long internalMaskTimeStamp = m_InternalMask->GetMTime();
490  unsigned long planarFigureTimeStamp = m_PlanarFigure->GetMTime();
491  unsigned long inputImageTimeStamp = m_inputImage->GetMTime();
492 
493  if (thisClassTimeStamp > m_InternalMaskUpdateTime) // inputs have changed
494  {
495  return true;
496  }
497 
498  if (m_InternalMaskUpdateTime < planarFigureTimeStamp || m_InternalMaskUpdateTime < inputImageTimeStamp) // mask image has changed outside of this class
499  {
500  return true;
501  }
502 
503  if (internalMaskTimeStamp > m_InternalMaskUpdateTime) // internal mask has been changed outside of this class
504  {
505  return true;
506  }
507 
508  return false;
509 }
510 
511 }
512 
mitk::Image::ConstPointer m_inputImage
Vector3D GetAxisVector(unsigned int direction) const
Get vector along bounding-box in the specified direction in mm.
#define MITK_ERROR
Definition: mitkLogMacros.h:20
DataCollection - Class to facilitate loading/accessing structured data.
virtual bool Map(const mitk::Point3D &pt3d_mm, mitk::Point2D &pt2d_mm) const
Project a 3D point given in mm (pt3d_mm) onto the 2D geometry. The result is a 2D point in mm (pt2d_m...
static Pointer New()
Image::Pointer GrabItkImageMemory(itk::SmartPointer< ItkOutputImageType > &itkimage, mitk::Image *mitkImage=nullptr, const BaseGeometry *geometry=nullptr, bool update=true)
Grabs the memory of an itk::Image (with a specific type) and puts it into an mitk::Image.The memory is managed by the mitk::Image after calling this function. The itk::Image remains valid until the mitk::Image decides to free the memory.
void SetPlanarFigure(mitk::PlanarFigure::Pointer planarFigure)
#define MITK_WARN
Definition: mitkLogMacros.h:19
void SetTimeStep(unsigned int timeStep) override
SetTimeStep is used to set the time step for which the mask is to be generated.
#define mitkThrow()
mitk::Image::Pointer m_InternalMask
mitk::Image::Pointer image
void WorldToIndex(const mitk::Point3D &pt_mm, mitk::Point3D &pt_units) const
Convert world coordinates (in mm) of a point to (continuous!) index coordinates.
mitk::Image::Pointer GetMask() override
GetMask Computes and returns the mask.
Vector3D GetNormal() const
Normal of the plane.
#define AccessFixedDimensionByItk_1(mitkImage, itkImageTypeFunction, dimension, arg1)
MITKCORE_EXPORT const ScalarType eps
std::vector< PolyLineElement > PolyLineType
Describes a two-dimensional, rectangular plane.
static bool CheckPlanarFigureIsNotTilted(const PlaneGeometry *planarGeometry, const BaseGeometry *geometry)
mitk::Image::ConstPointer GetReferenceImage() override
GetReferenceImage per default returns the inputImage (as set by SetInputImage). If no input image is ...
BaseGeometry Describes the geometry of a data object.
static Pointer New()