Medical Imaging Interaction Toolkit  2018.4.99-12ad79a3
Medical Imaging Interaction Toolkit
mitkExtractSliceFilter.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 
13 #include "mitkExtractSliceFilter.h"
14 
16 #include <mitkPlaneClipping.h>
17 
18 #include <vtkGeneralTransform.h>
19 #include <vtkImageChangeInformation.h>
20 #include <vtkImageData.h>
21 #include <vtkImageExtractComponents.h>
22 #include <vtkLinearTransform.h>
23 
24 mitk::ExtractSliceFilter::ExtractSliceFilter(vtkImageReslice *reslicer): m_XMin(0), m_XMax(0), m_YMin(0), m_YMax(0)
25 {
26  if (reslicer == nullptr)
27  {
28  m_Reslicer = vtkSmartPointer<vtkImageReslice>::New();
29  }
30  else
31  {
32  m_Reslicer = reslicer;
33  }
34 
35  m_TimeStep = 0;
36  m_Reslicer->ReleaseDataFlagOn();
38  m_ResliceTransform = nullptr;
42  m_ZSpacing = 1.0;
43  m_ZMin = 0;
44  m_ZMax = 0;
45  m_VtkOutputRequested = false;
46  m_BackgroundLevel = -32768.0;
47  m_Component = 0;
48 }
49 
51 {
52  m_ResliceTransform = nullptr;
53  m_WorldGeometry = nullptr;
54  delete[] m_OutPutSpacing;
55 }
56 
58 {
59  Image::ConstPointer input = this->GetInput();
60 
61  if (input.IsNull())
62  {
63  return;
64  }
65 
66  if (nullptr == m_WorldGeometry)
67  {
68  return;
69  }
70 
71  Vector3D right, bottom;
72  double widthInMM, heightInMM;
73  Vector2D extent;
74 
75  // set the geometry from current worldgeometry for the reusultimage
76  // this is needed that the image has the correct mitk geometry
77  // the sliceGeometry is the Geometry of the result slice
78  PlaneGeometry::Pointer sliceGeometry = m_WorldGeometry->Clone();
79 
80  sliceGeometry->GetIndexToWorldTransform()->SetMatrix(m_WorldGeometry->GetIndexToWorldTransform()->GetMatrix());
81 
82  // the origin of the worldGeometry is transformed to center based coordinates to be an imageGeometry
83  Point3D sliceOrigin = sliceGeometry->GetOrigin();
84 
85  auto abstractGeometry =
86  dynamic_cast<const AbstractTransformGeometry *>(m_WorldGeometry);
87 
88  auto planeGeometry = dynamic_cast<const PlaneGeometry *>(m_WorldGeometry);
89 
90  if (abstractGeometry != nullptr)
91  {
92  extent[0] = abstractGeometry->GetParametricExtent(0);
93  extent[1] = abstractGeometry->GetParametricExtent(1);
94 
95  widthInMM = abstractGeometry->GetParametricExtentInMM(0);
96  heightInMM = abstractGeometry->GetParametricExtentInMM(1);
97 
98  right = abstractGeometry->GetPlane()->GetAxisVector(0);
99  bottom = abstractGeometry->GetPlane()->GetAxisVector(1);
100  }
101  else
102  {
103  if (planeGeometry != nullptr)
104  {
105  // if the worldGeomatry is a PlaneGeometry everything is straight forward
106  right = planeGeometry->GetAxisVector(0);
107  bottom = planeGeometry->GetAxisVector(1);
108 
110  {
111  // Resampling grid corresponds to the current world geometry. This
112  // means that the spacing of the output 2D image depends on the
113  // currently selected world geometry, and *not* on the image itself.
114  extent[0] = m_WorldGeometry->GetExtent(0);
115  extent[1] = m_WorldGeometry->GetExtent(1);
116  }
117  else
118  {
119  const TimeGeometry *inputTimeGeometry = input->GetTimeGeometry();
120  if ((inputTimeGeometry == nullptr) || (inputTimeGeometry->CountTimeSteps() <= 0))
121  {
122  itkWarningMacro(<< "Error reading input image TimeGeometry.");
123  return;
124  }
125 
126  // Resampling grid corresponds to the input geometry. This means that
127  // the spacing of the output 2D image is directly derived from the
128  // associated input image, regardless of the currently selected world
129  // geometry.
130  Vector3D rightInIndex, bottomInIndex;
131  inputTimeGeometry->GetGeometryForTimeStep(m_TimeStep)->WorldToIndex(right, rightInIndex);
132  inputTimeGeometry->GetGeometryForTimeStep(m_TimeStep)->WorldToIndex(bottom, bottomInIndex);
133  extent[0] = rightInIndex.GetNorm();
134  extent[1] = bottomInIndex.GetNorm();
135  }
136 
137  // Get the extent of the current world geometry and calculate resampling
138  // spacing therefrom.
139  widthInMM = m_WorldGeometry->GetExtentInMM(0);
140  heightInMM = m_WorldGeometry->GetExtentInMM(1);
141  }
142  else
143  {
144  mitkThrow()<<"mitk::ExtractSliceFilter: No fitting geometry for reslice axis!";
145  return;
146  }
147  }
148 
149  right.Normalize();
150  bottom.Normalize();
151 
152  m_OutPutSpacing[0] = widthInMM / extent[0];
153  m_OutPutSpacing[1] = heightInMM / extent[1];
154 
155  /*========== BEGIN setup extent of the slice ==========*/
156  int xMin, xMax, yMin, yMax;
157 
158  xMin = yMin = 0;
159  xMax = static_cast<int>(extent[0]);
160  yMax = static_cast<int>(extent[1]);
161 
163  {
164  double sliceBounds[6];
165  for (auto &sliceBound : sliceBounds)
166  {
167  sliceBound = 0.0;
168  }
169 
170  if (this->GetClippedPlaneBounds(m_WorldGeometry->GetReferenceGeometry(), planeGeometry, sliceBounds))
171  {
172  // Calculate output extent (integer values)
173  xMin = static_cast<int>(sliceBounds[0] / m_OutPutSpacing[0] + 0.5);
174  xMax = static_cast<int>(sliceBounds[1] / m_OutPutSpacing[0] + 0.5);
175  yMin = static_cast<int>(sliceBounds[2] / m_OutPutSpacing[1] + 0.5);
176  yMax = static_cast<int>(sliceBounds[3] / m_OutPutSpacing[1] + 0.5);
177  } // ELSE we use the default values
178  }
179 
180 
181  sliceOrigin += right * (m_OutPutSpacing[0] * 0.5);
182  sliceOrigin += bottom * (m_OutPutSpacing[1] * 0.5);
183 
184  // a worldGeometry is no imageGeometry, thus it is manually set to true
185  sliceGeometry->ImageGeometryOn();
186 
187  /*At this point we have to adjust the geometry because the origin isn't correct.
188  The wrong origin is related to the rotation of the current world geometry plane.
189  This causes errors on transferring world to index coordinates. We just shift the
190  origin in each direction about the amount of the expanding (needed while rotating
191  the plane).
192  */
193  Vector3D axis0 = sliceGeometry->GetAxisVector(0);
194  Vector3D axis1 = sliceGeometry->GetAxisVector(1);
195  axis0.Normalize();
196  axis1.Normalize();
197 
198  // adapt the origin. Note that for orthogonal planes the minima are '0' and thus the origin stays the same.
199  sliceOrigin += (axis0 * (xMin * m_OutPutSpacing[0])) + (axis1 * (yMin * m_OutPutSpacing[1]));
200 
201  sliceGeometry->SetOrigin(sliceOrigin);
202 
203  /*the bounds as well as the extent of the worldGeometry are not adapted correctly during crosshair rotation.
204  This is only a quick fix and has to be evaluated.
205  The new bounds are set via the max values of the calculated slice extent. It will look like [ 0, x, 0, y, 0, 1].
206  */
208  boundsCopy[0] = boundsCopy[2] = boundsCopy[4] = 0;
209  boundsCopy[5] = 1;
210  boundsCopy[1] = xMax - xMin;
211  boundsCopy[3] = yMax - yMin;
212  sliceGeometry->SetBounds(boundsCopy);
213 
214  sliceGeometry->Modified();
215 
216  Image::Pointer output = this->GetOutput();
217  output->Initialize(input->GetPixelType(), 2, *sliceGeometry);
218 
219  m_XMin = xMin;
220  m_XMax = xMax;
221  m_YMin = yMin;
222  m_YMax = yMax;
223  m_Right = right;
224  m_Bottom = bottom;
225 }
226 
228 {
229  // As we want all pixel information fo the image in our plane, the requested region
230  // is set to the largest possible region in the image.
231  // This is needed because an oblique plane has a larger extent then the image
232  // and the in pipeline it is checked via PropagateResquestedRegion(). But the
233  // extent of the slice is actually fitting because it is oblique within the image.
235  input->SetRequestedRegionToLargestPossibleRegion();
236 }
237 
239 {
240  return m_OutPutSpacing;
241 }
242 
244 {
245  mitk::Image *input = this->GetInput();
246 
247  if (!input)
248  {
249  MITK_ERROR << "mitk::ExtractSliceFilter: No input image available. Please set the input!" << std::endl;
250  itkExceptionMacro("mitk::ExtractSliceFilter: No input image available. Please set the input!");
251  return;
252  }
253 
254  if (!m_WorldGeometry)
255  {
256  MITK_ERROR << "mitk::ExtractSliceFilter: No Geometry for reslicing available." << std::endl;
257  itkExceptionMacro("mitk::ExtractSliceFilter: No Geometry for reslicing available.");
258  return;
259  }
260 
261  const TimeGeometry *inputTimeGeometry = this->GetInput()->GetTimeGeometry();
262  if ((inputTimeGeometry == nullptr) || (inputTimeGeometry->CountTimeSteps() <= 0))
263  {
264  itkWarningMacro(<< "Error reading input image TimeGeometry.");
265  return;
266  }
267 
268  // is it a valid timeStep?
269  if (inputTimeGeometry->IsValidTimeStep(m_TimeStep) == false)
270  {
271  itkWarningMacro(<< "This is not a valid timestep: " << m_TimeStep);
272  return;
273  }
274 
275  // check if there is something to display.
276  if (!input->IsVolumeSet(m_TimeStep))
277  {
278  itkWarningMacro(<< "No volume data existent at given timestep " << m_TimeStep);
279  return;
280  }
281 
282  /*================#BEGIN setup vtkImageReslice properties================*/
283  Point3D origin;
284  Vector3D normal;
285 
286  const auto *planeGeometry = dynamic_cast<const PlaneGeometry *>(m_WorldGeometry);
287  // Code for curved planes, mostly taken 1:1 from imageVtkMapper2D and not tested yet.
288  // Do we have an AbstractTransformGeometry?
289  // This is the case for AbstractTransformGeometry's (e.g. a ThinPlateSplineCurvedGeometry )
290  const auto *abstractGeometry =
291  dynamic_cast<const AbstractTransformGeometry *>(m_WorldGeometry);
292 
293  if (abstractGeometry != nullptr)
294  {
295  m_ResliceTransform = abstractGeometry;
296 
297  origin = abstractGeometry->GetPlane()->GetOrigin();
298 
299  normal = abstractGeometry->GetPlane()->GetNormal();
300  normal.Normalize();
301 
302  // Use a combination of the InputGeometry *and* the possible non-rigid
303  // AbstractTransformGeometry for reslicing the 3D Image
304  vtkSmartPointer<vtkGeneralTransform> composedResliceTransform = vtkSmartPointer<vtkGeneralTransform>::New();
305  composedResliceTransform->Identity();
306  composedResliceTransform->Concatenate(
307  inputTimeGeometry->GetGeometryForTimeStep(m_TimeStep)->GetVtkTransform()->GetLinearInverse());
308  composedResliceTransform->Concatenate(abstractGeometry->GetVtkAbstractTransform());
309 
310  m_Reslicer->SetResliceTransform(composedResliceTransform);
311 
312  // Set background level to BLACK instead of translucent, to avoid
313  // boundary artifacts (see PlaneGeometryDataVtkMapper3D)
314  // Note: Backgroundlevel was hardcoded before to -1023
315  m_Reslicer->SetBackgroundLevel(m_BackgroundLevel);
316  }
317  else
318  {
319  if (planeGeometry != nullptr)
320  {
321  // if the worldGeomatry is a PlaneGeometry everything is straight forward
322 
323  origin = planeGeometry->GetOrigin();
324  normal = planeGeometry->GetNormal();
325  normal.Normalize();
326 
327  /*
328  * Transform the origin to center based coordinates.
329  * Note:
330  * This is needed besause vtk's origin is center based too (!!!) ( see 'The VTK book' page 88 )
331  * and the worldGeometry surrouding the image is no imageGeometry. So the worldGeometry
332  * has its origin at the corner of the voxel and needs to be transformed.
333  */
334  origin += m_Right * (m_OutPutSpacing[0] * 0.5);
335  origin += m_Bottom * (m_OutPutSpacing[1] * 0.5);
336 
337  // set the tranform for reslicing.
338  // Use inverse transform of the input geometry for reslicing the 3D image.
339  // This is needed if the image volume already transformed
340  if (m_ResliceTransform.IsNotNull())
341  m_Reslicer->SetResliceTransform(m_ResliceTransform->GetVtkTransform()->GetLinearInverse());
342 
343  // Set background level to TRANSLUCENT (see PlaneGeometryDataVtkMapper3D),
344  // else the background of the image turns out gray
345  // Note: Backgroundlevel was hardcoded to -32768
346  m_Reslicer->SetBackgroundLevel(m_BackgroundLevel);
347  }
348  else
349  {
350  itkExceptionMacro("mitk::ExtractSliceFilter: No fitting geometry for reslice axis!");
351  return;
352  }
353  }
354 
355  if (m_ResliceTransform.IsNotNull())
356  {
357  // if the resliceTransform is set the reslice axis are recalculated.
358  // Thus the geometry information is not fitting. Therefor a unitSpacingFilter
359  // is used to set up a global spacing of 1 and compensate the transform.
360  vtkSmartPointer<vtkImageChangeInformation> unitSpacingImageFilter =
361  vtkSmartPointer<vtkImageChangeInformation>::New();
362  unitSpacingImageFilter->ReleaseDataFlagOn();
363 
364  unitSpacingImageFilter->SetOutputSpacing(1.0, 1.0, 1.0);
365  unitSpacingImageFilter->SetInputData(input->GetVtkImageData(m_TimeStep));
366 
367  m_Reslicer->SetInputConnection(unitSpacingImageFilter->GetOutputPort());
368  }
369  else
370  {
371  // if no transform is set the image can be used directly
372  m_Reslicer->SetInputData(input->GetVtkImageData(m_TimeStep));
373  }
374 
375  /*setup the plane where vktImageReslice extracts the slice*/
376 
377  // ResliceAxesOrigin is the anchor point of the plane
378  double originInVtk[3];
379  itk2vtk(origin, originInVtk);
380  m_Reslicer->SetResliceAxesOrigin(originInVtk);
381 
382  // the cosines define the plane: x and y are the direction vectors, n is the planes normal
383  // this specifies a matrix 3x3
384  // x1 y1 n1
385  // x2 y2 n2
386  // x3 y3 n3
387  double cosines[9];
388 
389  vnl2vtk(m_Right.GetVnlVector(), cosines); // x
390 
391  vnl2vtk(m_Bottom.GetVnlVector(), cosines + 3); // y
392 
393  vnl2vtk(normal.GetVnlVector(), cosines + 6); // n
394 
395  m_Reslicer->SetResliceAxesDirectionCosines(cosines);
396 
397  // we only have one slice, not a volume
398  m_Reslicer->SetOutputDimensionality(m_OutputDimension);
399 
400  // set the interpolation mode for slicing
401  switch (this->m_InterpolationMode)
402  {
403  case RESLICE_NEAREST:
404  m_Reslicer->SetInterpolationModeToNearestNeighbor();
405  break;
406  case RESLICE_LINEAR:
407  m_Reslicer->SetInterpolationModeToLinear();
408  break;
409  case RESLICE_CUBIC:
410  m_Reslicer->SetInterpolationModeToCubic();
411  break;
412  default:
413  // the default interpolation used by mitk
414  m_Reslicer->SetInterpolationModeToNearestNeighbor();
415  }
416 
417  /*========== BEGIN setup extent of the slice ==========*/
418  // Set the output extents! First included pixel index and last included pixel index
419  // xMax and yMax are one after the last pixel. so they have to be decremented by 1.
420  // In case we have a 2D image, xMax or yMax might be 0. in this case, do not decrement, but take 0.
421 
422  m_Reslicer->SetOutputExtent(m_XMin, std::max(0, m_XMax - 1), m_YMin, std::max(0, m_YMax - 1), m_ZMin, m_ZMax);
423  /*========== END setup extent of the slice ==========*/
424 
425  m_Reslicer->SetOutputOrigin(0.0, 0.0, 0.0);
426 
427  m_Reslicer->SetOutputSpacing(m_OutPutSpacing[0], m_OutPutSpacing[1], m_ZSpacing);
428 
429  // TODO check the following lines, they are responsible whether vtk error outputs appear or not
430  m_Reslicer->UpdateWholeExtent(); // this produces a bad allocation error for 2D images
431  // m_Reslicer->GetOutput()->UpdateInformation();
432  // m_Reslicer->GetOutput()->SetUpdateExtentToWholeExtent();
433 
434  // start the pipeline
435  m_Reslicer->Update();
436  /*================ #END setup vtkImageReslice properties================*/
437 
439  {
440  // no conversion to mitk
441  // no mitk geometry will be set, as the output is vtkImageData only!!!
442  // no image component will be extracted, as the caller might need the whole multi-component image as vtk output
443  return;
444  }
445  else
446  {
447  auto reslicedImage = vtkSmartPointer<vtkImageData>::New();
448  reslicedImage = m_Reslicer->GetOutput();
449 
450  if (nullptr == reslicedImage)
451  {
452  itkWarningMacro(<< "Reslicer returned empty image");
453  return;
454  }
455 
456  /*================ #BEGIN Extract component from image slice ================*/
457  int numberOfScalarComponent = reslicedImage->GetNumberOfScalarComponents();
458  if (numberOfScalarComponent > 1 && static_cast<unsigned int>(numberOfScalarComponent) >= m_Component)
459  {
460  // image has more than one component, extract the correct component information with the given 'component' parameter
461  auto vectorComponentExtractor = vtkSmartPointer<vtkImageExtractComponents>::New();
462  vectorComponentExtractor->SetInputData(reslicedImage);
463  vectorComponentExtractor->SetComponents(m_Component);
464  vectorComponentExtractor->Update();
465 
466  reslicedImage = vectorComponentExtractor->GetOutput();
467  }
468  /*================ #END Extract component from image slice ================*/
469 
470  /*================ #BEGIN Convert the slice to an mitk::Image ================*/
471  mitk::Image::Pointer resultImage = GetOutput();
472 
473  /*Temporary store the geometry that is already correct (set in GeneratOutputInformation())
474  but will be reset due to initialize.*/
475  mitk::BaseGeometry::Pointer resultGeometry = resultImage->GetGeometry();
476 
477  // initialize resultimage with the specs of the vtkImageData object returned from vtkImageReslice
478  if (reslicedImage->GetDataDimension() == 1)
479  {
480  // If original image was 2D, the slice might have an y extent of 0.
481  // Still i want to ensure here that Image is 2D
482  resultImage->Initialize(reslicedImage, 1, -1, -1, 1);
483  }
484  else
485  {
486  resultImage->Initialize(reslicedImage);
487  }
488 
489  // transfer the voxel data
490  resultImage->SetVolume(reslicedImage->GetScalarPointer());
491 
492  /*================ #END Convert the slice to an mitk::Image ================*/
493 
494  resultImage->SetGeometry(resultGeometry);
495  }
496 }
497 
499 {
500  if (!m_WorldGeometry || !this->GetInput())
501  return false;
502 
503  return this->GetClippedPlaneBounds(
504  m_WorldGeometry->GetReferenceGeometry(), dynamic_cast<const PlaneGeometry *>(m_WorldGeometry), bounds);
505 }
506 
508  const PlaneGeometry *planeGeometry,
509  double *bounds)
510 {
511  bool b = mitk::PlaneClipping::CalculateClippedPlaneBounds(boundingGeometry, planeGeometry, bounds);
512 
513  return b;
514 }
ResliceInterpolation m_InterpolationMode
#define MITK_ERROR
Definition: mitkLogMacros.h:20
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...
bool IsVolumeSet(int t=0, int n=0) const override
Check whether volume at time t in channel n is set.
Definition: mitkImage.cpp:602
ScalarType GetExtent(unsigned int direction) const
Set the time bounds (in ms)
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:217
Pointer Clone() const
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.
const mitk::TimeGeometry * GetTimeGeometry() const
Return the TimeGeometry of the data as const pointer.
Definition: mitkBaseData.h:61
mitk::ScalarType * GetOutputSpacing()
Get the spacing of the slice. returns mitk::ScalarType[2].
const PlaneGeometry * m_WorldGeometry
void vnl2vtk(const vnl_vector< Tin > &in, Tout *out)
ExtractSliceFilter(vtkImageReslice *reslicer=nullptr)
#define mitkThrow()
Image class for storing images.
Definition: mitkImage.h:72
Describes a geometry defined by an vtkAbstractTransform and a plane.
static T max(T x, T y)
Definition: svm.cpp:56
void itk2vtk(const Tin &in, Tout &out)
InputImageType * GetInput(void)
mitk::ScalarType * m_OutPutSpacing
virtual BaseGeometry::Pointer GetGeometryForTimeStep(TimeStepType timeStep) const =0
Returns the geometry which corresponds to the given time step.
bool GetClippedPlaneBounds(double bounds[6])
Get the bounding box of the slice [xMin, xMax, yMin, yMax, zMin, zMax] The method uses the input of t...
OutputType * GetOutput()
Get the output data of this image source object.
Describes a two-dimensional, rectangular plane.
void GenerateInputRequestedRegion() override
vtkSmartPointer< vtkImageReslice > m_Reslicer
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.
mitk::AffineTransform3D * GetIndexToWorldTransform()
Get the transformation used to convert from index to world coordinates.
BoundingBoxType::BoundsArrayType BoundsArrayType