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