Medical Imaging Interaction Toolkit  2018.4.99-12ad79a3
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  bool outOfBounds = false;
131  vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
132  typename PlanarFigure::PolyLineType::const_iterator it;
133  for ( it = planarFigurePolyline.begin();
134  it != planarFigurePolyline.end();
135  ++it )
136  {
137  Point3D point3D;
138 
139  // Convert 2D point back to the local index coordinates of the selected
140  // image
141  // Fabian: From PlaneGeometry documentation:
142  // Converts a 2D point given in mm (pt2d_mm) relative to the upper-left corner of the geometry into the corresponding world-coordinate (a 3D point in mm, pt3d_mm).
143  // To convert a 2D point given in units (e.g., pixels in case of an image) into a 2D point given in mm (as required by this method), use IndexToWorld.
144  planarFigurePlaneGeometry->Map( *it, point3D );
145 
146  // Polygons (partially) outside of the image bounds can not be processed
147  // further due to a bug in vtkPolyDataToImageStencil
148  if ( !imageGeometry3D->IsInside( point3D ) )
149  {
150  outOfBounds = true;
151  }
152 
153  imageGeometry3D->WorldToIndex( point3D, point3D );
154 
155  points->InsertNextPoint( point3D[i0], point3D[i1], 0 );
156  }
157 
158  vtkSmartPointer<vtkPoints> holePoints = nullptr;
159 
160  if (!planarFigureHolePolyline.empty())
161  {
162  holePoints = vtkSmartPointer<vtkPoints>::New();
163 
164  Point3D point3D;
165  PlanarFigure::PolyLineType::const_iterator end = planarFigureHolePolyline.end();
166 
167  for (it = planarFigureHolePolyline.begin(); it != end; ++it)
168  {
169  // Fabian: same as above
170  planarFigurePlaneGeometry->Map(*it, point3D);
171  imageGeometry3D->WorldToIndex(point3D, point3D);
172  holePoints->InsertNextPoint(point3D[i0], point3D[i1], 0);
173  }
174  }
175 
176  // mark a malformed 2D planar figure ( i.e. area = 0 ) as out of bounds
177  // this can happen when all control points of a rectangle lie on the same line = two of the three extents are zero
178  double bounds[6] = {0, 0, 0, 0, 0, 0};
179  points->GetBounds( bounds );
180  bool extent_x = (fabs(bounds[0] - bounds[1])) < mitk::eps;
181  bool extent_y = (fabs(bounds[2] - bounds[3])) < mitk::eps;
182  bool extent_z = (fabs(bounds[4] - bounds[5])) < mitk::eps;
183 
184  // throw an exception if a closed planar figure is deformed, i.e. has only one non-zero extent
185  if ( m_PlanarFigure->IsClosed() &&
186  ((extent_x && extent_y) || (extent_x && extent_z) || (extent_y && extent_z)))
187  {
188  mitkThrow() << "Figure has a zero area and cannot be used for masking.";
189  }
190 
191  if ( outOfBounds )
192  {
193  throw std::runtime_error( "Figure at least partially outside of image bounds!" );
194  }
195 
196  // create a vtkLassoStencilSource and set the points of the Polygon
197  vtkSmartPointer<vtkLassoStencilSource> lassoStencil = vtkSmartPointer<vtkLassoStencilSource>::New();
198  lassoStencil->SetShapeToPolygon();
199  lassoStencil->SetPoints( points );
200 
201  vtkSmartPointer<vtkLassoStencilSource> holeLassoStencil = nullptr;
202 
203  if (holePoints.GetPointer() != nullptr)
204  {
205  holeLassoStencil = vtkSmartPointer<vtkLassoStencilSource>::New();
206  holeLassoStencil->SetShapeToPolygon();
207  holeLassoStencil->SetPoints(holePoints);
208  }
209 
210  // Export from ITK to VTK (to use a VTK filter)
211  typedef itk::VTKImageImport< MaskImage2DType > ImageImportType;
212  typedef itk::VTKImageExport< MaskImage2DType > ImageExportType;
213 
214  typename ImageExportType::Pointer itkExporter = ImageExportType::New();
215  itkExporter->SetInput( maskImage );
216 // itkExporter->SetInput( castFilter->GetOutput() );
217 
218  vtkSmartPointer<vtkImageImport> vtkImporter = vtkSmartPointer<vtkImageImport>::New();
219  this->ConnectPipelines( itkExporter, vtkImporter );
220 
221  // Apply the generated image stencil to the input image
222  vtkSmartPointer<vtkImageStencil> imageStencilFilter = vtkSmartPointer<vtkImageStencil>::New();
223  imageStencilFilter->SetInputConnection( vtkImporter->GetOutputPort() );
224  imageStencilFilter->SetStencilConnection(lassoStencil->GetOutputPort());
225  imageStencilFilter->ReverseStencilOff();
226  imageStencilFilter->SetBackgroundValue( 0 );
227  imageStencilFilter->Update();
228 
229  vtkSmartPointer<vtkImageStencil> holeStencilFilter = nullptr;
230 
231  if (holeLassoStencil.GetPointer() != nullptr)
232  {
233  holeStencilFilter = vtkSmartPointer<vtkImageStencil>::New();
234  holeStencilFilter->SetInputConnection(imageStencilFilter->GetOutputPort());
235  holeStencilFilter->SetStencilConnection(holeLassoStencil->GetOutputPort());
236  holeStencilFilter->ReverseStencilOn();
237  holeStencilFilter->SetBackgroundValue(0);
238  holeStencilFilter->Update();
239  }
240 
241  // Export from VTK back to ITK
242  vtkSmartPointer<vtkImageExport> vtkExporter = vtkSmartPointer<vtkImageExport>::New();
243  vtkExporter->SetInputConnection( holeStencilFilter.GetPointer() == nullptr
244  ? imageStencilFilter->GetOutputPort()
245  : holeStencilFilter->GetOutputPort());
246  vtkExporter->Update();
247 
248  typename ImageImportType::Pointer itkImporter = ImageImportType::New();
249  this->ConnectPipelines( vtkExporter, itkImporter );
250  itkImporter->Update();
251 
252  typedef itk::ImageDuplicator< ImageImportType::OutputImageType > DuplicatorType;
253  DuplicatorType::Pointer duplicator = DuplicatorType::New();
254  duplicator->SetInputImage( itkImporter->GetOutput() );
255  duplicator->Update();
256 
257  // Store mask
258  m_InternalITKImageMask2D = duplicator->GetOutput();
259 }
260 
261 template < typename TPixel, unsigned int VImageDimension >
262 void PlanarFigureMaskGenerator::InternalCalculateMaskFromOpenPlanarFigure(
263  const itk::Image< TPixel, VImageDimension > *image, unsigned int axis )
264 {
265  typedef itk::Image< unsigned short, 2 > MaskImage2DType;
266  typedef itk::LineIterator< MaskImage2DType > LineIteratorType;
267  typedef MaskImage2DType::IndexType IndexType2D;
268  typedef std::vector< IndexType2D > IndexVecType;
269 
270  typename MaskImage2DType::Pointer maskImage = MaskImage2DType::New();
271  maskImage->SetOrigin(image->GetOrigin());
272  maskImage->SetSpacing(image->GetSpacing());
273  maskImage->SetLargestPossibleRegion(image->GetLargestPossibleRegion());
274  maskImage->SetBufferedRegion(image->GetBufferedRegion());
275  maskImage->SetDirection(image->GetDirection());
276  maskImage->SetNumberOfComponentsPerPixel(image->GetNumberOfComponentsPerPixel());
277  maskImage->Allocate();
278  maskImage->FillBuffer(0);
279 
280  // all PolylinePoints of the PlanarFigure are stored in a vtkPoints object.
281  const mitk::PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry();
282  const typename PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 );
283  const mitk::BaseGeometry *imageGeometry3D = m_inputImage->GetGeometry( 0 );
284 
285  // Determine x- and y-dimensions depending on principal axis
286  // TODO use plane geometry normal to determine that automatically, then check whether the PF is aligned with one of the three principal axis
287  int i0, i1;
288  switch ( axis )
289  {
290  case 0:
291  i0 = 1;
292  i1 = 2;
293  break;
294 
295  case 1:
296  i0 = 0;
297  i1 = 2;
298  break;
299 
300  case 2:
301  default:
302  i0 = 0;
303  i1 = 1;
304  break;
305  }
306 
307  int numPolyLines = m_PlanarFigure->GetPolyLinesSize();
308  for ( int lineId = 0; lineId < numPolyLines; ++lineId )
309  {
310  // store the polyline contour as vtkPoints object
311  bool outOfBounds = false;
312  IndexVecType pointIndices;
313  typename PlanarFigure::PolyLineType::const_iterator it;
314  for ( it = planarFigurePolyline.begin();
315  it != planarFigurePolyline.end();
316  ++it )
317  {
318  Point3D point3D;
319 
320  planarFigurePlaneGeometry->Map( *it, point3D );
321 
322  if ( !imageGeometry3D->IsInside( point3D ) )
323  {
324  outOfBounds = true;
325  }
326 
327  imageGeometry3D->WorldToIndex( point3D, point3D );
328 
329  IndexType2D index2D;
330  index2D[0] = point3D[i0];
331  index2D[1] = point3D[i1];
332 
333  pointIndices.push_back( index2D );
334  }
335 
336  if ( outOfBounds )
337  {
338  throw std::runtime_error( "Figure at least partially outside of image bounds!" );
339  }
340 
341  for ( IndexVecType::const_iterator it = pointIndices.begin(); it != pointIndices.end()-1; ++it )
342  {
343  IndexType2D ind1 = *it;
344  IndexType2D ind2 = *(it+1);
345 
346  LineIteratorType lineIt( maskImage, ind1, ind2 );
347  while ( !lineIt.IsAtEnd() )
348  {
349  lineIt.Set( 1 );
350  ++lineIt;
351  }
352  }
353  }
354 
355  // Store mask
356  m_InternalITKImageMask2D = maskImage;
357 }
358 
359 bool PlanarFigureMaskGenerator::GetPrincipalAxis(
360  const BaseGeometry *geometry, Vector3D vector,
361  unsigned int &axis )
362 {
363  vector.Normalize();
364  for ( unsigned int i = 0; i < 3; ++i )
365  {
366  Vector3D axisVector = geometry->GetAxisVector( i );
367  axisVector.Normalize();
368 
369  if ( fabs( fabs( axisVector * vector ) - 1.0) < mitk::eps )
370  {
371  axis = i;
372  return true;
373  }
374  }
375 
376  return false;
377 }
378 
379 void PlanarFigureMaskGenerator::CalculateMask()
380 {
381  if (m_inputImage.IsNull())
382  {
383  MITK_ERROR << "Image is not set.";
384  }
385 
386  if (m_PlanarFigure.IsNull())
387  {
388  MITK_ERROR << "PlanarFigure is not set.";
389  }
390 
391  if (m_TimeStep != 0)
392  {
393  MITK_WARN << "Multiple TimeSteps are not supported in PlanarFigureMaskGenerator (yet).";
394  }
395 
396  const BaseGeometry *imageGeometry = m_inputImage->GetGeometry();
397  if ( imageGeometry == nullptr )
398  {
399  throw std::runtime_error( "Image geometry invalid!" );
400  }
401 
402  if (m_inputImage->GetTimeSteps() > 0)
403  {
405  imgTimeSel->SetInput(m_inputImage);
406  imgTimeSel->SetTimeNr(m_TimeStep);
407  imgTimeSel->UpdateLargestPossibleRegion();
408  m_InternalTimeSliceImage = imgTimeSel->GetOutput();
409  }
410  else
411  {
412  m_InternalTimeSliceImage = m_inputImage;
413  }
414 
415  m_InternalITKImageMask2D = nullptr;
416  const PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry();
417  const auto *planarFigureGeometry = dynamic_cast< const PlaneGeometry * >( planarFigurePlaneGeometry );
418  //const BaseGeometry *imageGeometry = m_inputImage->GetGeometry();
419 
420  // Find principal direction of PlanarFigure in input image
421  unsigned int axis;
422  if ( !this->GetPrincipalAxis( imageGeometry,
423  planarFigureGeometry->GetNormal(), axis ) )
424  {
425  throw std::runtime_error( "Non-aligned planar figures not supported!" );
426  }
427  m_PlanarFigureAxis = axis;
428 
429  // Find slice number corresponding to PlanarFigure in input image
430  itk::Image< unsigned short, 3 >::IndexType index;
431  imageGeometry->WorldToIndex( planarFigureGeometry->GetOrigin(), index );
432 
433  unsigned int slice = index[axis];
434  m_PlanarFigureSlice = slice;
435 
436  // extract image slice which corresponds to the planarFigure and store it in m_InternalImageSlice
437  mitk::Image::ConstPointer inputImageSlice = extract2DImageSlice(axis, slice);
438  //mitk::IOUtil::Save(inputImageSlice, "/home/fabian/inputSliceImage.nrrd");
439  // Compute mask from PlanarFigure
440  // rastering for open planar figure:
441  if ( !m_PlanarFigure->IsClosed() )
442  {
443  AccessFixedDimensionByItk_1(inputImageSlice,
444  InternalCalculateMaskFromOpenPlanarFigure,
445  2, axis)
446  }
447  else//for closed planar figure
448  {
449  AccessFixedDimensionByItk_1(inputImageSlice,
450  InternalCalculateMaskFromPlanarFigure,
451  2, axis)
452  }
453 
454  //convert itk mask to mitk::Image::Pointer and return it
455  mitk::Image::Pointer planarFigureMaskImage;
456  planarFigureMaskImage = mitk::GrabItkImageMemory(m_InternalITKImageMask2D);
457  //mitk::IOUtil::Save(planarFigureMaskImage, "/home/fabian/planarFigureMaskImage.nrrd");
458 
459  //Convert2Dto3DImageFilter::Pointer sliceTo3DImageConverter = Convert2Dto3DImageFilter::New();
460  //sliceTo3DImageConverter->SetInput(planarFigureMaskImage);
461  //sliceTo3DImageConverter->Update();
462  //mitk::IOUtil::Save(sliceTo3DImageConverter->GetOutput(), "/home/fabian/3DsliceImage.nrrd");
463 
464  m_ReferenceImage = inputImageSlice;
465  //mitk::IOUtil::Save(m_ReferenceImage, "/home/fabian/referenceImage.nrrd");
466  m_InternalMask = planarFigureMaskImage;
467 }
468 
469 void PlanarFigureMaskGenerator::SetTimeStep(unsigned int timeStep)
470 {
471  if (timeStep != m_TimeStep)
472  {
473  m_TimeStep = timeStep;
474  }
475 }
476 
478 {
479  if (IsUpdateRequired())
480  {
481  this->CalculateMask();
482  this->Modified();
483  }
484 
485  m_InternalMaskUpdateTime = this->GetMTime();
486  return m_InternalMask;
487 }
488 
489 mitk::Image::ConstPointer PlanarFigureMaskGenerator::extract2DImageSlice(unsigned int axis, unsigned int slice)
490 {
491  // Extract slice with given position and direction from image
492  unsigned int dimension = m_InternalTimeSliceImage->GetDimension();
493 
494  if (dimension == 3)
495  {
497  imageExtractor->SetInput( m_InternalTimeSliceImage );
498  imageExtractor->SetSliceDimension( axis );
499  imageExtractor->SetSliceIndex( slice );
500  imageExtractor->Update();
501  return imageExtractor->GetOutput();
502  }
503  else if(dimension == 2)
504  {
505  return m_InternalTimeSliceImage;
506  }
507  else
508  {
509  MITK_ERROR << "Unsupported image dimension. Dimension is: " << dimension << ". Only 2D and 3D images are supported.";
510  return nullptr;
511  }
512 }
513 
514 bool PlanarFigureMaskGenerator::IsUpdateRequired() const
515 {
516  unsigned long thisClassTimeStamp = this->GetMTime();
517  unsigned long internalMaskTimeStamp = m_InternalMask->GetMTime();
518  unsigned long planarFigureTimeStamp = m_PlanarFigure->GetMTime();
519  unsigned long inputImageTimeStamp = m_inputImage->GetMTime();
520 
521  if (thisClassTimeStamp > m_InternalMaskUpdateTime) // inputs have changed
522  {
523  return true;
524  }
525 
526  if (m_InternalMaskUpdateTime < planarFigureTimeStamp || m_InternalMaskUpdateTime < inputImageTimeStamp) // mask image has changed outside of this class
527  {
528  return true;
529  }
530 
531  if (internalMaskTimeStamp > m_InternalMaskUpdateTime) // internal mask has been changed outside of this class
532  {
533  return true;
534  }
535 
536  return false;
537 }
538 
539 }
540 
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
bool IsInside(const mitk::Point3D &p) const
Test whether the point p (world coordinates in mm) is inside the bounding box.
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.
#define AccessFixedDimensionByItk_1(mitkImage, itkImageTypeFunction, dimension, arg1)
MITKCORE_EXPORT const ScalarType eps
std::vector< PolyLineElement > PolyLineType
Describes a two-dimensional, rectangular plane.
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()