Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
mitkPlanarFigureMaskGenerator.cpp
Go to the documentation of this file.
2 #include <mitkBaseGeometry.h>
3 #include <mitkITKImageImport.h>
4 #include "mitkImageAccessByItk.h"
8 #include <mitkIOUtil.h>
9 
10 #include <itkCastImageFilter.h>
11 #include <itkVTKImageExport.h>
12 #include <itkVTKImageImport.h>
13 #include <itkImageDuplicator.h>
14 #include <itkExceptionObject.h>
15 
16 #include <vtkPoints.h>
17 #include <vtkImageStencil.h>
18 #include <vtkImageImport.h>
19 #include <vtkImageExport.h>
20 #include <vtkLassoStencilSource.h>
21 #include <vtkSmartPointer.h>
22 
23 
24 
25 
26 namespace mitk
27 {
28 
30 {
31  if ( planarFigure.IsNull() )
32  {
33  throw std::runtime_error( "Error: planar figure empty!" );
34  }
35  if ( !planarFigure->IsClosed() )
36  {
37  throw std::runtime_error( "Masking not possible for non-closed figures" );
38  }
39 
40  const PlaneGeometry *planarFigurePlaneGeometry = planarFigure->GetPlaneGeometry();
41  if ( planarFigurePlaneGeometry == nullptr )
42  {
43  throw std::runtime_error( "Planar-Figure not yet initialized!" );
44  }
45 
46  const PlaneGeometry *planarFigureGeometry =
47  dynamic_cast< const PlaneGeometry * >( planarFigurePlaneGeometry );
48  if ( planarFigureGeometry == nullptr )
49  {
50  throw std::runtime_error( "Non-planar planar figures not supported!" );
51  }
52 
53  if (planarFigure != m_PlanarFigure)
54  {
55  this->Modified();
56  m_PlanarFigure = planarFigure;
57  }
58 
59 }
60 
62 {
63  if (IsUpdateRequired())
64  {
65  this->CalculateMask();
66  }
67  return m_ReferenceImage;
68 }
69 
70 template < typename TPixel, unsigned int VImageDimension >
71 void PlanarFigureMaskGenerator::InternalCalculateMaskFromPlanarFigure(
72  const itk::Image< TPixel, VImageDimension > *image, unsigned int axis )
73 {
74  typedef itk::Image< TPixel, VImageDimension > ImageType;
75  typedef itk::Image< unsigned short, 2 > MaskImage2DType;
76 
77  typename MaskImage2DType::Pointer maskImage = MaskImage2DType::New();
78  maskImage->SetOrigin(image->GetOrigin());
79  maskImage->SetSpacing(image->GetSpacing());
80  maskImage->SetLargestPossibleRegion(image->GetLargestPossibleRegion());
81  maskImage->SetBufferedRegion(image->GetBufferedRegion());
82  maskImage->SetDirection(image->GetDirection());
83  maskImage->SetNumberOfComponentsPerPixel(image->GetNumberOfComponentsPerPixel());
84  maskImage->Allocate();
85  maskImage->FillBuffer(1);
86 
87  // all PolylinePoints of the PlanarFigure are stored in a vtkPoints object.
88  // These points are used by the vtkLassoStencilSource to create
89  // a vtkImageStencil.
90  const mitk::PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry();
91  const typename PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 );
92  const mitk::BaseGeometry *imageGeometry3D = m_inputImage->GetGeometry( 0 );
93  // If there is a second poly line in a closed planar figure, treat it as a hole.
94  PlanarFigure::PolyLineType planarFigureHolePolyline;
95 
96  if (m_PlanarFigure->GetPolyLinesSize() == 2)
97  planarFigureHolePolyline = m_PlanarFigure->GetPolyLine(1);
98 
99 
100  // Determine x- and y-dimensions depending on principal axis
101  // TODO use plane geometry normal to determine that automatically, then check whether the PF is aligned with one of the three principal axis
102  int i0, i1;
103  switch ( axis )
104  {
105  case 0:
106  i0 = 1;
107  i1 = 2;
108  break;
109 
110  case 1:
111  i0 = 0;
112  i1 = 2;
113  break;
114 
115  case 2:
116  default:
117  i0 = 0;
118  i1 = 1;
119  break;
120  }
121 
122  // store the polyline contour as vtkPoints object
123  bool outOfBounds = false;
124  vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
125  typename PlanarFigure::PolyLineType::const_iterator it;
126  for ( it = planarFigurePolyline.begin();
127  it != planarFigurePolyline.end();
128  ++it )
129  {
130  Point3D point3D;
131 
132  // Convert 2D point back to the local index coordinates of the selected
133  // image
134  // Fabian: From PlaneGeometry documentation:
135  // 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).
136  // 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.
137  planarFigurePlaneGeometry->Map( *it, point3D );
138 
139  // Polygons (partially) outside of the image bounds can not be processed
140  // further due to a bug in vtkPolyDataToImageStencil
141  if ( !imageGeometry3D->IsInside( point3D ) )
142  {
143  outOfBounds = true;
144  }
145 
146  imageGeometry3D->WorldToIndex( point3D, point3D );
147 
148  points->InsertNextPoint( point3D[i0], point3D[i1], 0 );
149  }
150 
151  vtkSmartPointer<vtkPoints> holePoints = nullptr;
152 
153  if (!planarFigureHolePolyline.empty())
154  {
155  holePoints = vtkSmartPointer<vtkPoints>::New();
156 
157  Point3D point3D;
158  PlanarFigure::PolyLineType::const_iterator end = planarFigureHolePolyline.end();
159 
160  for (it = planarFigureHolePolyline.begin(); it != end; ++it)
161  {
162  // Fabian: same as above
163  planarFigurePlaneGeometry->Map(*it, point3D);
164  imageGeometry3D->WorldToIndex(point3D, point3D);
165  holePoints->InsertNextPoint(point3D[i0], point3D[i1], 0);
166  }
167  }
168 
169  // mark a malformed 2D planar figure ( i.e. area = 0 ) as out of bounds
170  // this can happen when all control points of a rectangle lie on the same line = two of the three extents are zero
171  double bounds[6] = {0, 0, 0, 0, 0, 0};
172  points->GetBounds( bounds );
173  bool extent_x = (fabs(bounds[0] - bounds[1])) < mitk::eps;
174  bool extent_y = (fabs(bounds[2] - bounds[3])) < mitk::eps;
175  bool extent_z = (fabs(bounds[4] - bounds[5])) < mitk::eps;
176 
177  // throw an exception if a closed planar figure is deformed, i.e. has only one non-zero extent
178  if ( m_PlanarFigure->IsClosed() &&
179  ((extent_x && extent_y) || (extent_x && extent_z) || (extent_y && extent_z)))
180  {
181  mitkThrow() << "Figure has a zero area and cannot be used for masking.";
182  }
183 
184  if ( outOfBounds )
185  {
186  throw std::runtime_error( "Figure at least partially outside of image bounds!" );
187  }
188 
189  // create a vtkLassoStencilSource and set the points of the Polygon
190  vtkSmartPointer<vtkLassoStencilSource> lassoStencil = vtkSmartPointer<vtkLassoStencilSource>::New();
191  lassoStencil->SetShapeToPolygon();
192  lassoStencil->SetPoints( points );
193 
194  vtkSmartPointer<vtkLassoStencilSource> holeLassoStencil = nullptr;
195 
196  if (holePoints.GetPointer() != nullptr)
197  {
198  holeLassoStencil = vtkSmartPointer<vtkLassoStencilSource>::New();
199  holeLassoStencil->SetShapeToPolygon();
200  holeLassoStencil->SetPoints(holePoints);
201  }
202 
203  // Export from ITK to VTK (to use a VTK filter)
204  typedef itk::VTKImageImport< MaskImage2DType > ImageImportType;
205  typedef itk::VTKImageExport< MaskImage2DType > ImageExportType;
206 
207  typename ImageExportType::Pointer itkExporter = ImageExportType::New();
208  itkExporter->SetInput( maskImage );
209 // itkExporter->SetInput( castFilter->GetOutput() );
210 
211  vtkSmartPointer<vtkImageImport> vtkImporter = vtkSmartPointer<vtkImageImport>::New();
212  this->ConnectPipelines( itkExporter, vtkImporter );
213 
214  // Apply the generated image stencil to the input image
215  vtkSmartPointer<vtkImageStencil> imageStencilFilter = vtkSmartPointer<vtkImageStencil>::New();
216  imageStencilFilter->SetInputConnection( vtkImporter->GetOutputPort() );
217  imageStencilFilter->SetStencilConnection(lassoStencil->GetOutputPort());
218  imageStencilFilter->ReverseStencilOff();
219  imageStencilFilter->SetBackgroundValue( 0 );
220  imageStencilFilter->Update();
221 
222  vtkSmartPointer<vtkImageStencil> holeStencilFilter = nullptr;
223 
224  if (holeLassoStencil.GetPointer() != nullptr)
225  {
226  holeStencilFilter = vtkSmartPointer<vtkImageStencil>::New();
227  holeStencilFilter->SetInputConnection(imageStencilFilter->GetOutputPort());
228  holeStencilFilter->SetStencilConnection(holeLassoStencil->GetOutputPort());
229  holeStencilFilter->ReverseStencilOn();
230  holeStencilFilter->SetBackgroundValue(0);
231  holeStencilFilter->Update();
232  }
233 
234  // Export from VTK back to ITK
235  vtkSmartPointer<vtkImageExport> vtkExporter = vtkSmartPointer<vtkImageExport>::New();
236  vtkExporter->SetInputConnection( holeStencilFilter.GetPointer() == nullptr
237  ? imageStencilFilter->GetOutputPort()
238  : holeStencilFilter->GetOutputPort());
239  vtkExporter->Update();
240 
241  typename ImageImportType::Pointer itkImporter = ImageImportType::New();
242  this->ConnectPipelines( vtkExporter, itkImporter );
243  itkImporter->Update();
244 
245  typedef itk::ImageDuplicator< ImageImportType::OutputImageType > DuplicatorType;
247  duplicator->SetInputImage( itkImporter->GetOutput() );
248  duplicator->Update();
249 
250  // Store mask
251  m_InternalITKImageMask2D = duplicator->GetOutput();
252 }
253 
254 bool PlanarFigureMaskGenerator::GetPrincipalAxis(
255  const BaseGeometry *geometry, Vector3D vector,
256  unsigned int &axis )
257 {
258  vector.Normalize();
259  for ( unsigned int i = 0; i < 3; ++i )
260  {
261  Vector3D axisVector = geometry->GetAxisVector( i );
262  axisVector.Normalize();
263 
264  if ( fabs( fabs( axisVector * vector ) - 1.0) < mitk::eps )
265  {
266  axis = i;
267  return true;
268  }
269  }
270 
271  return false;
272 }
273 
274 void PlanarFigureMaskGenerator::CalculateMask()
275 {
276  if (m_inputImage.IsNull())
277  {
278  MITK_ERROR << "Image is not set.";
279  }
280 
281  if (m_PlanarFigure.IsNull())
282  {
283  MITK_ERROR << "PlanarFigure is not set.";
284  }
285 
286  if (m_TimeStep != 0)
287  {
288  MITK_WARN << "Multiple TimeSteps are not supported in PlanarFigureMaskGenerator (yet).";
289  }
290 
291  const BaseGeometry *imageGeometry = m_inputImage->GetGeometry();
292  if ( imageGeometry == nullptr )
293  {
294  throw std::runtime_error( "Image geometry invalid!" );
295  }
296 
297  if (m_inputImage->GetTimeSteps() > 0)
298  {
300  imgTimeSel->SetInput(m_inputImage);
301  imgTimeSel->SetTimeNr(m_TimeStep);
302  imgTimeSel->UpdateLargestPossibleRegion();
303  m_InternalTimeSliceImage = imgTimeSel->GetOutput();
304  }
305  else
306  {
307  m_InternalTimeSliceImage = m_inputImage;
308  }
309 
310  m_InternalITKImageMask2D = nullptr;
311  const PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry();
312  const PlaneGeometry *planarFigureGeometry = dynamic_cast< const PlaneGeometry * >( planarFigurePlaneGeometry );
313  //const BaseGeometry *imageGeometry = m_inputImage->GetGeometry();
314 
315  // Find principal direction of PlanarFigure in input image
316  unsigned int axis;
317  if ( !this->GetPrincipalAxis( imageGeometry,
318  planarFigureGeometry->GetNormal(), axis ) )
319  {
320  throw std::runtime_error( "Non-aligned planar figures not supported!" );
321  }
322  m_PlanarFigureAxis = axis;
323 
324  // Find slice number corresponding to PlanarFigure in input image
325  itk::Image< unsigned short, 3 >::IndexType index;
326  imageGeometry->WorldToIndex( planarFigureGeometry->GetOrigin(), index );
327 
328  unsigned int slice = index[axis];
329 
330  // extract image slice which corresponds to the planarFigure and store it in m_InternalImageSlice
331  mitk::Image::Pointer inputImageSlice = extract2DImageSlice(axis, slice);
332  //mitk::IOUtil::SaveImage(inputImageSlice, "/home/fabian/inputSliceImage.nrrd");
333  // Compute mask from PlanarFigure
334  AccessFixedDimensionByItk_1(inputImageSlice,
335  InternalCalculateMaskFromPlanarFigure,
336  2, axis)
337 
338  //convert itk mask to mitk::Image::Pointer and return it
339  mitk::Image::Pointer planarFigureMaskImage;
340  planarFigureMaskImage = mitk::GrabItkImageMemory(m_InternalITKImageMask2D);
341  //mitk::IOUtil::SaveImage(planarFigureMaskImage, "/home/fabian/planarFigureMaskImage.nrrd");
342 
343  //Convert2Dto3DImageFilter::Pointer sliceTo3DImageConverter = Convert2Dto3DImageFilter::New();
344  //sliceTo3DImageConverter->SetInput(planarFigureMaskImage);
345  //sliceTo3DImageConverter->Update();
346  //mitk::IOUtil::SaveImage(sliceTo3DImageConverter->GetOutput(), "/home/fabian/3DsliceImage.nrrd");
347 
348  m_ReferenceImage = inputImageSlice;
349  //mitk::IOUtil::SaveImage(m_ReferenceImage, "/home/fabian/referenceImage.nrrd");
350  m_InternalMask = planarFigureMaskImage;
351 }
352 
353 void PlanarFigureMaskGenerator::SetTimeStep(unsigned int timeStep)
354 {
355  if (timeStep != m_TimeStep)
356  {
357  m_TimeStep = timeStep;
358  }
359 }
360 
362 {
363  if (IsUpdateRequired())
364  {
365  this->CalculateMask();
366  this->Modified();
367  }
368 
369  m_InternalMaskUpdateTime = this->GetMTime();
370  return m_InternalMask;
371 }
372 
373 mitk::Image::Pointer PlanarFigureMaskGenerator::extract2DImageSlice(unsigned int axis, unsigned int slice)
374 {
375  // Extract slice with given position and direction from image
376  unsigned int dimension = m_InternalTimeSliceImage->GetDimension();
377  mitk::Image::Pointer imageSlice = mitk::Image::New();
378 
379  if (dimension == 3)
380  {
382  imageExtractor->SetInput( m_InternalTimeSliceImage );
383  imageExtractor->SetSliceDimension( axis );
384  imageExtractor->SetSliceIndex( slice );
385  imageExtractor->Update();
386  imageSlice = imageExtractor->GetOutput();
387  }
388  else if(dimension == 2)
389  {
390  imageSlice = m_InternalTimeSliceImage;
391  }
392  else
393  {
394  MITK_ERROR << "Unsupported image dimension. Dimension is: " << dimension << ". Only 2D and 3D images are supported.";
395  }
396 
397  return imageSlice;
398 }
399 
400 bool PlanarFigureMaskGenerator::IsUpdateRequired() const
401 {
402  unsigned long thisClassTimeStamp = this->GetMTime();
403  unsigned long internalMaskTimeStamp = m_InternalMask->GetMTime();
404  unsigned long planarFigureTimeStamp = m_PlanarFigure->GetMTime();
405  unsigned long inputImageTimeStamp = m_inputImage->GetMTime();
406 
407  if (thisClassTimeStamp > m_InternalMaskUpdateTime) // inputs have changed
408  {
409  return true;
410  }
411 
412  if (m_InternalMaskUpdateTime < planarFigureTimeStamp || m_InternalMaskUpdateTime < inputImageTimeStamp) // mask image has changed outside of this class
413  {
414  return true;
415  }
416 
417  if (internalMaskTimeStamp > m_InternalMaskUpdateTime) // internal mask has been changed outside of this class
418  {
419  return true;
420  }
421 
422  return false;
423 }
424 
425 }
426 
mitk::Image::Pointer m_inputImage
mitk::Image::Pointer GetReferenceImage()
GetReferenceImage per default returns the inputImage (as set by SetInputImage). If no input image is ...
itk::SmartPointer< Self > Pointer
void SetTimeStep(unsigned int timeStep)
SetTimeStep is used to set the time step for which the mask is to be generated.
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...
#define MITK_ERROR
Definition: mitkLogMacros.h:24
DataCollection - Class to facilitate loading/accessing structured data.
mitk::Image::Pointer GetMask()
GetMask Computes and returns the mask.
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)
map::core::discrete::Elements< 3 >::InternalImageType ImageType
class ITK_EXPORT Image
#define MITK_WARN
Definition: mitkLogMacros.h:23
Vector< ScalarType, 3 > Vector3D
Definition: mitkVector.h:134
#define mitkThrow()
mitk::Image::Pointer m_InternalMask
Point< ScalarType, 3 > Point3D
Definition: mitkPoint.h:99
static Pointer New()
#define AccessFixedDimensionByItk_1(mitkImage, itkImageTypeFunction, dimension, arg1)
Derived from MaskGenerator. This class is used to convert a mitk::PlanarFigure into a binary image ma...
MITKCORE_EXPORT const ScalarType eps
std::vector< PolyLineElement > PolyLineType
Describes a two-dimensional, rectangular plane.
bool IsInside(const mitk::Point3D &p) const
Test whether the point p (world coordinates in mm) is inside the bounding box.
BaseGeometry Describes the geometry of a data object.
static Pointer New()
void WorldToIndex(const mitk::Point3D &pt_mm, mitk::Point3D &pt_units) const
Convert world coordinates (in mm) of a point to (continuous!) index coordinates.
static itkEventMacro(BoundingShapeInteractionEvent, itk::AnyEvent) class MITKBOUNDINGSHAPE_EXPORT BoundingShapeInteractor Pointer New()
Basic interaction methods for mitk::GeometryData.