Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
mitkExtractDirectedPlaneImageFilter.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 
19 //#include "mitkImageMapperGL2D.h"
20 
22 #include <mitkDataNode.h>
23 #include <mitkProperties.h>
25 
26 #include <vtkGeneralTransform.h>
27 #include <vtkImageChangeInformation.h>
28 #include <vtkImageData.h>
29 #include <vtkPoints.h>
30 #include <vtkSmartPointer.h>
31 #include <vtkTransform.h>
32 #include <vtkTransform.h>
33 
35 {
36  MITK_WARN << "Class ExtractDirectedPlaneImageFilter is deprecated! Use ExtractSliceFilter instead.";
37 
39 
40  m_TargetTimestep = 0;
42  m_ResliceInterpolationProperty = nullptr; // VtkResliceInterpolationProperty::New(); //TODO initial with value
44  m_ThickSlicesNum = 1;
45 }
46 
48 {
49  if (m_ResliceInterpolationProperty != nullptr)
50  m_ResliceInterpolationProperty->Delete();
51  m_Reslicer->Delete();
52 }
53 
55 {
56  // A world geometry must be set...
57  if (m_WorldGeometry == nullptr)
58  {
59  itkWarningMacro(<< "No world geometry has been set. Returning.");
60  return;
61  }
62 
63  Image *input = const_cast<ImageToImageFilter::InputImageType *>(this->GetInput());
64  input->Update();
65 
66  if (input == nullptr)
67  {
68  itkWarningMacro(<< "No input set.");
69  return;
70  }
71 
72  const TimeGeometry *inputTimeGeometry = input->GetTimeGeometry();
73  if ((inputTimeGeometry == nullptr) || (inputTimeGeometry->CountTimeSteps() == 0))
74  {
75  itkWarningMacro(<< "Error reading input image geometry.");
76  return;
77  }
78 
79  // Get the target timestep; if none is set, use the lowest given.
80  unsigned int timestep = m_TargetTimestep;
81 
82  if (inputTimeGeometry->IsValidTimeStep(timestep) == false)
83  {
84  itkWarningMacro(<< "This is not a valid timestep: " << timestep);
85  return;
86  }
87 
88  // check if there is something to display.
89  if (!input->IsVolumeSet(timestep))
90  {
91  itkWarningMacro(<< "No volume data existent at given timestep " << timestep);
92  return;
93  }
94 
95  Image::RegionType requestedRegion = input->GetLargestPossibleRegion();
96  requestedRegion.SetIndex(3, timestep);
97  requestedRegion.SetSize(3, 1);
98  requestedRegion.SetSize(4, 1);
99  input->SetRequestedRegion(&requestedRegion);
100  input->Update();
101 
102  vtkImageData *inputData = input->GetVtkImageData(timestep);
103 
104  if (inputData == nullptr)
105  {
106  itkWarningMacro(<< "Could not extract vtk image data for given timestep" << timestep);
107  return;
108  }
109 
110  double spacing[3];
111  inputData->GetSpacing(spacing);
112 
113  // how big the area is in physical coordinates: widthInMM x heightInMM pixels
114  mitk::ScalarType widthInMM, heightInMM;
115 
116  // where we want to sample
117  Point3D origin;
118  Vector3D right, bottom, normal;
119  Vector3D rightInIndex, bottomInIndex;
120 
121  assert(input->GetTimeGeometry() == inputTimeGeometry);
122 
123  // take transform of input image into account
124  BaseGeometry *inputGeometry = inputTimeGeometry->GetGeometryForTimeStep(timestep);
125  if (inputGeometry == nullptr)
126  {
127  itkWarningMacro(<< "There is no Geometry3D at given timestep " << timestep);
128  return;
129  }
130 
131  ScalarType mmPerPixel[2];
132 
133  // Bounds information for reslicing (only required if reference geometry
134  // is present)
135  double bounds[6];
136  bool boundsInitialized = false;
137 
138  for (auto &bound : bounds)
139  {
140  bound = 0.0;
141  }
142 
143  Vector2D extent;
144  extent.Fill(0.0);
145 
146  // Do we have a simple PlaneGeometry?
147  if (dynamic_cast<const PlaneGeometry *>(m_WorldGeometry) != nullptr &&
148  dynamic_cast<const AbstractTransformGeometry *>(m_WorldGeometry) == nullptr)
149  {
150  const PlaneGeometry *planeGeometry = static_cast<const PlaneGeometry *>(m_WorldGeometry);
151  origin = planeGeometry->GetOrigin();
152  right = planeGeometry->GetAxisVector(0);
153  bottom = planeGeometry->GetAxisVector(1);
154  normal = planeGeometry->GetNormal();
155 
156  if (m_InPlaneResampleExtentByGeometry)
157  {
158  // Resampling grid corresponds to the current world geometry. This
159  // means that the spacing of the output 2D image depends on the
160  // currently selected world geometry, and *not* on the image itself.
161 
162  extent[0] = m_WorldGeometry->GetExtent(0);
163  extent[1] = m_WorldGeometry->GetExtent(1);
164  }
165  else
166  {
167  // Resampling grid corresponds to the input geometry. This means that
168  // the spacing of the output 2D image is directly derived from the
169  // associated input image, regardless of the currently selected world
170  // geometry.
171  inputGeometry->WorldToIndex(right, rightInIndex);
172  inputGeometry->WorldToIndex(bottom, bottomInIndex);
173  extent[0] = rightInIndex.GetNorm();
174  extent[1] = bottomInIndex.GetNorm();
175  }
176 
177  // Get the extent of the current world geometry and calculate resampling
178  // spacing therefrom.
179  widthInMM = m_WorldGeometry->GetExtentInMM(0);
180  heightInMM = m_WorldGeometry->GetExtentInMM(1);
181 
182  mmPerPixel[0] = widthInMM / extent[0];
183  mmPerPixel[1] = heightInMM / extent[1];
184 
185  right.Normalize();
186  bottom.Normalize();
187  normal.Normalize();
188 
189  // origin += right * ( mmPerPixel[0] * 0.5 );
190  // origin += bottom * ( mmPerPixel[1] * 0.5 );
191 
192  // widthInMM -= mmPerPixel[0];
193  // heightInMM -= mmPerPixel[1];
194 
195  // Use inverse transform of the input geometry for reslicing the 3D image
196  m_Reslicer->SetResliceTransform(inputGeometry->GetVtkTransform()->GetLinearInverse());
197 
198  // Set background level to TRANSLUCENT (see PlaneGeometryDataVtkMapper3D)
199  m_Reslicer->SetBackgroundLevel(-32768);
200 
201  // Check if a reference geometry does exist (as would usually be the case for
202  // PlaneGeometry).
203  // Note: this is currently not strictly required, but could facilitate
204  // correct plane clipping.
205  if (m_WorldGeometry->GetReferenceGeometry())
206  {
207  // Calculate the actual bounds of the transformed plane clipped by the
208  // dataset bounding box; this is required for drawing the texture at the
209  // correct position during 3D mapping.
210  boundsInitialized =
211  this->CalculateClippedPlaneBounds(m_WorldGeometry->GetReferenceGeometry(), planeGeometry, bounds);
212  }
213  }
214 
215  // Do we have an AbstractTransformGeometry?
216  else if (dynamic_cast<const AbstractTransformGeometry *>(m_WorldGeometry))
217  {
218  const mitk::AbstractTransformGeometry *abstractGeometry =
219  dynamic_cast<const AbstractTransformGeometry *>(m_WorldGeometry);
220 
221  extent[0] = abstractGeometry->GetParametricExtent(0);
222  extent[1] = abstractGeometry->GetParametricExtent(1);
223 
224  widthInMM = abstractGeometry->GetParametricExtentInMM(0);
225  heightInMM = abstractGeometry->GetParametricExtentInMM(1);
226 
227  mmPerPixel[0] = widthInMM / extent[0];
228  mmPerPixel[1] = heightInMM / extent[1];
229 
230  origin = abstractGeometry->GetPlane()->GetOrigin();
231 
232  right = abstractGeometry->GetPlane()->GetAxisVector(0);
233  right.Normalize();
234 
235  bottom = abstractGeometry->GetPlane()->GetAxisVector(1);
236  bottom.Normalize();
237 
238  normal = abstractGeometry->GetPlane()->GetNormal();
239  normal.Normalize();
240 
241  // Use a combination of the InputGeometry *and* the possible non-rigid
242  // AbstractTransformGeometry for reslicing the 3D Image
243  vtkGeneralTransform *composedResliceTransform = vtkGeneralTransform::New();
244  composedResliceTransform->Identity();
245  composedResliceTransform->Concatenate(inputGeometry->GetVtkTransform()->GetLinearInverse());
246  composedResliceTransform->Concatenate(abstractGeometry->GetVtkAbstractTransform());
247 
248  m_Reslicer->SetResliceTransform(composedResliceTransform);
249 
250  // Set background level to BLACK instead of translucent, to avoid
251  // boundary artifacts (see PlaneGeometryDataVtkMapper3D)
252  m_Reslicer->SetBackgroundLevel(-1023);
253  composedResliceTransform->Delete();
254  }
255  else
256  {
257  itkWarningMacro(<< "World Geometry has to be a PlaneGeometry or an AbstractTransformGeometry.");
258  return;
259  }
260 
261  // Make sure that the image to be resliced has a certain minimum size.
262  if ((extent[0] <= 2) && (extent[1] <= 2))
263  {
264  itkWarningMacro(<< "Image is too small to be resliced...");
265  return;
266  }
267 
268  vtkSmartPointer<vtkImageChangeInformation> unitSpacingImageFilter = vtkImageChangeInformation::New();
269  unitSpacingImageFilter->SetOutputSpacing(1.0, 1.0, 1.0);
270  unitSpacingImageFilter->SetInputData(inputData);
271 
272  m_Reslicer->SetInputConnection(unitSpacingImageFilter->GetOutputPort());
273 
274  // m_Reslicer->SetInput( inputData );
275 
276  m_Reslicer->SetOutputDimensionality(2);
277  m_Reslicer->SetOutputOrigin(0.0, 0.0, 0.0);
278 
279  Vector2D pixelsPerMM;
280  pixelsPerMM[0] = 1.0 / mmPerPixel[0];
281  pixelsPerMM[1] = 1.0 / mmPerPixel[1];
282 
283  // calulate the originArray and the orientations for the reslice-filter
284  double originArray[3];
285  itk2vtk(origin, originArray);
286 
287  m_Reslicer->SetResliceAxesOrigin(originArray);
288 
289  double cosines[9];
290 
291  // direction of the X-axis of the sampled result
292  vnl2vtk(right.GetVnlVector(), cosines);
293 
294  // direction of the Y-axis of the sampled result
295  vnl2vtk(bottom.GetVnlVector(), cosines + 3);
296 
297  // normal of the plane
298  vnl2vtk(normal.GetVnlVector(), cosines + 6);
299 
300  m_Reslicer->SetResliceAxesDirectionCosines(cosines);
301 
302  int xMin, xMax, yMin, yMax;
303  if (boundsInitialized)
304  {
305  xMin = static_cast<int>(bounds[0] / mmPerPixel[0]); //+ 0.5 );
306  xMax = static_cast<int>(bounds[1] / mmPerPixel[0]); //+ 0.5 );
307  yMin = static_cast<int>(bounds[2] / mmPerPixel[1]); //+ 0.5);
308  yMax = static_cast<int>(bounds[3] / mmPerPixel[1]); //+ 0.5 );
309  }
310  else
311  {
312  // If no reference geometry is available, we also don't know about the
313  // maximum plane size; so the overlap is just ignored
314 
315  xMin = yMin = 0;
316  xMax = static_cast<int>(extent[0] - pixelsPerMM[0]); //+ 0.5 );
317  yMax = static_cast<int>(extent[1] - pixelsPerMM[1]); //+ 0.5 );
318  }
319 
320  m_Reslicer->SetOutputSpacing(mmPerPixel[0], mmPerPixel[1], 1.0);
321  // xMax and yMax are meant exclusive until now, whereas
322  // SetOutputExtent wants an inclusive bound. Thus, we need
323  // to subtract 1.
324  m_Reslicer->SetOutputExtent(xMin, xMax - 1, yMin, yMax - 1, 0, 1);
325 
326  // Do the reslicing. Modified() is called to make sure that the reslicer is
327  // executed even though the input geometry information did not change; this
328  // is necessary when the input /em data, but not the /em geometry changes.
329  m_Reslicer->Modified();
330  m_Reslicer->ReleaseDataFlagOn();
331 
332  m_Reslicer->Update();
333 
334  // 1. Check the result
335  vtkImageData *reslicedImage = m_Reslicer->GetOutput();
336 
337  if ((reslicedImage == nullptr) || (reslicedImage->GetDataDimension() < 1))
338  {
339  itkWarningMacro(<< "Reslicer returned empty image");
340  return;
341  }
342 
343  unsigned int dimensions[2];
344  dimensions[0] = (unsigned int)extent[0];
345  dimensions[1] = (unsigned int)extent[1];
346  Vector3D spacingVector;
347  FillVector3D(spacingVector, mmPerPixel[0], mmPerPixel[1], 1.0);
348 
349  mitk::Image::Pointer resultImage = this->GetOutput();
350  resultImage->Initialize(input->GetPixelType(), 2, dimensions);
351  resultImage->SetSpacing(spacingVector);
352 }
353 
355 {
356  Superclass::GenerateOutputInformation();
357 }
358 
360  const PlaneGeometry *planeGeometry,
361  double *bounds)
362 {
363  // Clip the plane with the bounding geometry. To do so, the corner points
364  // of the bounding box are transformed by the inverse transformation
365  // matrix, and the transformed bounding box edges derived therefrom are
366  // clipped with the plane z=0. The resulting min/max values are taken as
367  // bounds for the image reslicer.
368  const BoundingBox *boundingBox = boundingGeometry->GetBoundingBox();
369 
370  BoundingBox::PointType bbMin = boundingBox->GetMinimum();
371  BoundingBox::PointType bbMax = boundingBox->GetMaximum();
372 
373  vtkPoints *points = vtkPoints::New();
374  if (boundingGeometry->GetImageGeometry())
375  {
376  points->InsertPoint(0, bbMin[0] - 0.5, bbMin[1] - 0.5, bbMin[2] - 0.5);
377  points->InsertPoint(1, bbMin[0] - 0.5, bbMin[1] - 0.5, bbMax[2] - 0.5);
378  points->InsertPoint(2, bbMin[0] - 0.5, bbMax[1] - 0.5, bbMax[2] - 0.5);
379  points->InsertPoint(3, bbMin[0] - 0.5, bbMax[1] - 0.5, bbMin[2] - 0.5);
380  points->InsertPoint(4, bbMax[0] - 0.5, bbMin[1] - 0.5, bbMin[2] - 0.5);
381  points->InsertPoint(5, bbMax[0] - 0.5, bbMin[1] - 0.5, bbMax[2] - 0.5);
382  points->InsertPoint(6, bbMax[0] - 0.5, bbMax[1] - 0.5, bbMax[2] - 0.5);
383  points->InsertPoint(7, bbMax[0] - 0.5, bbMax[1] - 0.5, bbMin[2] - 0.5);
384  }
385  else
386  {
387  points->InsertPoint(0, bbMin[0], bbMin[1], bbMin[2]);
388  points->InsertPoint(1, bbMin[0], bbMin[1], bbMax[2]);
389  points->InsertPoint(2, bbMin[0], bbMax[1], bbMax[2]);
390  points->InsertPoint(3, bbMin[0], bbMax[1], bbMin[2]);
391  points->InsertPoint(4, bbMax[0], bbMin[1], bbMin[2]);
392  points->InsertPoint(5, bbMax[0], bbMin[1], bbMax[2]);
393  points->InsertPoint(6, bbMax[0], bbMax[1], bbMax[2]);
394  points->InsertPoint(7, bbMax[0], bbMax[1], bbMin[2]);
395  }
396 
397  vtkPoints *newPoints = vtkPoints::New();
398 
399  vtkTransform *transform = vtkTransform::New();
400  transform->Identity();
401  transform->Concatenate(planeGeometry->GetVtkTransform()->GetLinearInverse());
402 
403  transform->Concatenate(boundingGeometry->GetVtkTransform());
404 
405  transform->TransformPoints(points, newPoints);
406  transform->Delete();
407 
408  bounds[0] = bounds[2] = 10000000.0;
409  bounds[1] = bounds[3] = -10000000.0;
410  bounds[4] = bounds[5] = 0.0;
411 
412  this->LineIntersectZero(newPoints, 0, 1, bounds);
413  this->LineIntersectZero(newPoints, 1, 2, bounds);
414  this->LineIntersectZero(newPoints, 2, 3, bounds);
415  this->LineIntersectZero(newPoints, 3, 0, bounds);
416  this->LineIntersectZero(newPoints, 0, 4, bounds);
417  this->LineIntersectZero(newPoints, 1, 5, bounds);
418  this->LineIntersectZero(newPoints, 2, 6, bounds);
419  this->LineIntersectZero(newPoints, 3, 7, bounds);
420  this->LineIntersectZero(newPoints, 4, 5, bounds);
421  this->LineIntersectZero(newPoints, 5, 6, bounds);
422  this->LineIntersectZero(newPoints, 6, 7, bounds);
423  this->LineIntersectZero(newPoints, 7, 4, bounds);
424 
425  // clean up vtk data
426  points->Delete();
427  newPoints->Delete();
428 
429  if ((bounds[0] > 9999999.0) || (bounds[2] > 9999999.0) || (bounds[1] < -9999999.0) || (bounds[3] < -9999999.0))
430  {
431  return false;
432  }
433  else
434  {
435  // The resulting bounds must be adjusted by the plane spacing, since we
436  // we have so far dealt with index coordinates
437  const mitk::Vector3D planeSpacing = planeGeometry->GetSpacing();
438  bounds[0] *= planeSpacing[0];
439  bounds[1] *= planeSpacing[0];
440  bounds[2] *= planeSpacing[1];
441  bounds[3] *= planeSpacing[1];
442  bounds[4] *= planeSpacing[2];
443  bounds[5] *= planeSpacing[2];
444  return true;
445  }
446 }
447 
448 bool mitk::ExtractDirectedPlaneImageFilter::LineIntersectZero(vtkPoints *points, int p1, int p2, double *bounds)
449 {
450  double point1[3];
451  double point2[3];
452  points->GetPoint(p1, point1);
453  points->GetPoint(p2, point2);
454 
455  if ((point1[2] * point2[2] <= 0.0) && (point1[2] != point2[2]))
456  {
457  double x, y;
458  x = (point1[0] * point2[2] - point1[2] * point2[0]) / (point2[2] - point1[2]);
459  y = (point1[1] * point2[2] - point1[2] * point2[1]) / (point2[2] - point1[2]);
460 
461  if (x < bounds[0])
462  {
463  bounds[0] = x;
464  }
465  if (x > bounds[1])
466  {
467  bounds[1] = x;
468  }
469  if (y < bounds[2])
470  {
471  bounds[2] = y;
472  }
473  if (y > bounds[3])
474  {
475  bounds[3] = y;
476  }
477  bounds[4] = bounds[5] = 0.0;
478  return true;
479  }
480  return false;
481 }
const Point3D GetOrigin() const
Get the origin, e.g. the upper-left corner of the plane.
mitk::Point3D PointType
itk::BoundingBox< unsigned long, 3, ScalarType > BoundingBox
Standard 3D-BoundingBox typedef.
bool CalculateClippedPlaneBounds(const BaseGeometry *boundingGeometry, const PlaneGeometry *planeGeometry, double *bounds)
double ScalarType
static bool CalculateClippedPlaneBounds(const BaseGeometry *boundingGeometry, const PlaneGeometry *planeGeometry, double *bounds)
Calculate the bounding box of the resliced image. This is necessary for arbitrarily rotated planes in...
virtual bool IsVolumeSet(int t=0, int n=0) const override
Check whether volume at time t in channel n is set.
Definition: mitkImage.cpp:606
virtual void GenerateData() override
A version of GenerateData() specific for image processing filters.
virtual void SetRequestedRegion(const itk::DataObject *data) override
const mitk::Vector3D GetSpacing() const
Get the spacing (size of a pixel).
virtual vtkImageData * GetVtkImageData(int t=0, int n=0)
Get a volume at a specific time t of channel n as a vtkImageData.
Definition: mitkImage.cpp:221
const mitk::TimeGeometry * GetTimeGeometry() const
Return the TimeGeometry of the data as const pointer.
Definition: mitkBaseData.h:52
void FillVector3D(Tout &out, mitk::ScalarType x, mitk::ScalarType y, mitk::ScalarType z)
Definition: mitkArray.h:110
Vector3D GetNormal() const
Normal of the plane.
virtual const PlaneGeometry * GetPlane()
Get the rectangular area that is used for transformation by m_VtkAbstractTransform and therewith defi...
VtkResliceInterpolationProperty * m_ResliceInterpolationProperty
Vector3D GetAxisVector(unsigned int direction) const
Get vector along bounding-box in the specified direction in mm.
#define MITK_WARN
Definition: mitkLogMacros.h:23
void vnl2vtk(const vnl_vector< Tin > &in, Tout *out)
itk::ImageRegion< RegionDimension > RegionType
Image class for storing images.
Definition: mitkImage.h:76
const RegionType & GetLargestPossibleRegion() const
Describes a geometry defined by an vtkAbstractTransform and a plane.
const mitk::PixelType GetPixelType(int n=0) const
Returns the PixelType of channel n.
Definition: mitkImage.cpp:105
bool LineIntersectZero(vtkPoints *points, int p1, int p2, double *bounds)
void itk2vtk(const Tin &in, Tout &out)
mitk::ScalarType GetParametricExtent(int direction) const
Get the parametric extent.
static bool LineIntersectZero(vtkPoints *points, int p1, int p2, double *bounds)
Internal helper method for intersection testing used only in CalculateClippedPlaneBounds() ...
virtual BaseGeometry::Pointer GetGeometryForTimeStep(TimeStepType timeStep) const =0
Returns the geometry which corresponds to the given time step.
virtual mitk::ScalarType GetParametricExtentInMM(int direction) const
virtual vtkAbstractTransform * GetVtkAbstractTransform() const
Get the vtkAbstractTransform (stored in m_VtkAbstractTransform)
Describes a two-dimensional, rectangular plane.
vtkLinearTransform * GetVtkTransform() const
Get the m_IndexToWorldTransform as a vtkLinearTransform.
virtual bool GetImageGeometry() const
Is this an ImageGeometry?
virtual bool IsValidTimeStep(TimeStepType timeStep) const =0
Test for the given time step if a geometry is availible.
BaseGeometry Describes the geometry of a data object.
virtual const BoundingBoxType * GetBoundingBox()
static itkEventMacro(BoundingShapeInteractionEvent, itk::AnyEvent) class MITKBOUNDINGSHAPE_EXPORT BoundingShapeInteractor Pointer New()
Basic interaction methods for mitk::GeometryData.