Medical Imaging Interaction Toolkit  2018.4.99-a3d2e8fb
Medical Imaging Interaction Toolkit
mitkBoundingShapeCropper.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 "mitkGeometry3D.h"
15 #include "mitkImageAccessByItk.h"
16 #include "mitkImageCast.h"
17 #include "mitkImageToItk.h"
18 #include "mitkStatusBar.h"
19 #include "mitkTimeHelper.h"
20 
21 #include <cmath>
22 
23 #include "vtkMatrix4x4.h"
24 #include "vtkSmartPointer.h"
25 #include "vtkTransform.h"
26 
27 #include "itkImageRegionIteratorWithIndex.h"
28 #include <itkImageIOBase.h>
29 #include <itkImageRegionConstIterator.h>
30 #include <itkRGBAPixel.h>
31 #include <itkRGBPixel.h>
32 
33 namespace mitk
34 {
36  : m_Geometry(nullptr),
37  m_OutsideValue(0),
38  m_UseCropTimeStepOnly(false),
39  m_CurrentTimeStep(0),
40  m_UseWholeInputRegion(false),
41  m_InputTimeSelector(mitk::ImageTimeSelector::New()),
42  m_OutputTimeSelector(mitk::ImageTimeSelector::New())
43  {
44  this->SetNumberOfIndexedInputs(2);
45  this->SetNumberOfRequiredInputs(2);
46  }
47 
49 
50  template <typename TPixel, unsigned int VImageDimension>
51  void BoundingShapeCropper::CutImage(itk::Image<TPixel, VImageDimension> *inputItkImage, int timeStep)
52  {
53  MITK_INFO << "Scalar Pixeltype" << std::endl;
54 
55  typedef TPixel TOutputPixel;
56  typedef itk::Image<TPixel, VImageDimension> ItkInputImageType;
57  typedef itk::Image<TOutputPixel, VImageDimension> ItkOutputImageType;
58  typedef typename itk::ImageBase<VImageDimension>::RegionType ItkRegionType;
59  typedef itk::ImageRegionIteratorWithIndex<ItkInputImageType> ItkInputImageIteratorType;
60  typedef itk::ImageRegionIteratorWithIndex<ItkOutputImageType> ItkOutputImageIteratorType;
61 
62  TOutputPixel outsideValue = this->GetOutsideValue();
63  // currently 0 if not set in advance
64  // TODO: change default value to itk::NumericTraits<TOutputPixel>::min();
65 
66  if (this->m_Geometry.IsNull())
67  return;
68 
69  if (inputItkImage == nullptr)
70  {
72  "An internal error occurred. Can't convert Image. Please report to bugs@mitk.org");
73  std::cout << " image is nullptr...returning" << std::endl;
74  return;
75  }
76 
77  // first convert the index
78  typename ItkRegionType::IndexType::IndexValueType tmpIndex[3];
79  itk2vtk(this->m_InputRequestedRegion.GetIndex(), tmpIndex);
80  typename ItkRegionType::IndexType index;
81  index.SetIndex(tmpIndex);
82 
83  // then convert the size
84  typename ItkRegionType::SizeType::SizeValueType tmpSize[3];
85  itk2vtk(this->m_InputRequestedRegion.GetSize(), tmpSize);
86  typename ItkRegionType::SizeType size;
87  size.SetSize(tmpSize);
88 
89  // create the ITK-image-region out of index and size
90  ItkRegionType inputRegionOfInterest(index, size);
91 
92  // Get access to the MITK output image via an ITK image
93  typename mitk::ImageToItk<ItkOutputImageType>::Pointer outputimagetoitk =
95  outputimagetoitk->SetInput(this->m_OutputTimeSelector->GetOutput());
96  outputimagetoitk->Update();
97  typename ItkOutputImageType::Pointer outputItkImage = outputimagetoitk->GetOutput();
98 
99  // create the iterators
100  ItkInputImageIteratorType inputIt(inputItkImage, inputRegionOfInterest);
101  ItkOutputImageIteratorType outputIt(outputItkImage, outputItkImage->GetLargestPossibleRegion());
102 
103  // Cut the boundingbox out of the image by iterating through all images
104  // TODO: use more efficient method by using the contour instead off all single pixels
105  mitk::Point3D p;
106  mitk::BaseGeometry *inputGeometry = this->GetInput()->GetGeometry(timeStep);
107 
108  // calculates translation based on offset+extent not on the transformation matrix
109  // NOTE: center of the box is
110  vtkSmartPointer<vtkMatrix4x4> imageTransform = this->m_Geometry->GetGeometry()->GetVtkTransform()->GetMatrix();
111  Point3D center = this->m_Geometry->GetGeometry()->GetCenter();
112  auto translation = vtkSmartPointer<vtkTransform>::New();
113  translation->Translate(center[0] - imageTransform->GetElement(0, 3),
114  center[1] - imageTransform->GetElement(1, 3),
115  center[2] - imageTransform->GetElement(2, 3));
116  auto transform = vtkSmartPointer<vtkTransform>::New();
117  transform->SetMatrix(imageTransform);
118  transform->PostMultiply();
119  transform->Concatenate(translation);
120  transform->Update();
121 
122  mitk::Vector3D extent;
123  for (unsigned int i = 0; i < 3; ++i)
124  extent[i] = (this->m_Geometry->GetGeometry()->GetExtent(i));
125 
126  for (inputIt.GoToBegin(), outputIt.GoToBegin(); !inputIt.IsAtEnd(); ++inputIt, ++outputIt)
127  {
128  vtk2itk(inputIt.GetIndex(), p);
129  inputGeometry->IndexToWorld(p, p);
130  ScalarType p2[4];
131  p2[0] = p[0];
132  p2[1] = p[1];
133  p2[2] = p[2];
134  p2[3] = 1;
135  // transform point from world to object coordinates
136  transform->GetInverse()->TransformPoint(p2, p2);
137  // check if the world point is within bounds
138  bool isInside = (p2[0] >= (-extent[0] / 2.0)) && (p2[0] <= (extent[0] / 2.0)) && (p2[1] >= (-extent[1] / 2.0)) &&
139  (p2[1] <= (extent[1] / 2.0)) && (p2[2] >= (-extent[2] / 2.0)) && (p2[2] <= (extent[2] / 2.0));
140 
141  if ((!this->m_UseCropTimeStepOnly && isInside) ||
142  (this->m_UseCropTimeStepOnly && timeStep == this->m_CurrentTimeStep && isInside))
143  {
144  outputIt.Set((TOutputPixel)inputIt.Value());
145  }
146  else
147  {
148  outputIt.Set(outsideValue);
149  }
150  }
151  }
152 
154  {
155  m_Geometry = const_cast<mitk::GeometryData *>(geometry);
156  // Process object is not const-correct so the const_cast is required here
157  this->ProcessObject::SetNthInput(1, const_cast<mitk::GeometryData *>(geometry));
158  }
159 
160  // const mitk::GeometryData* BoundingShapeCropper::GetGeometryData() const
161  //{
162  // return m_Geometry.GetPointer();
163  //}
164 
167  {
168  mitk::Image *output = this->GetOutput();
169  if ((output->IsInitialized() == false) || (m_Geometry.IsNull()) ||
170  (m_Geometry->GetTimeGeometry()->CountTimeSteps() == 0))
171  return;
172 
173  GenerateTimeInInputRegion(output, this->GetInput());
174  }
175 
177  {
178  // Set Cropping region
179  mitk::Image::Pointer output = this->GetOutput();
180  if ((output->IsInitialized()) && (output->GetPipelineMTime() <= m_TimeOfHeaderInitialization.GetMTime()))
181  return;
182 
183  mitk::Image::Pointer input = this->GetInput();
184 
185  if (input.IsNull())
186  {
187  mitkThrow() << "Input is not a mitk::Image";
188  }
189  itkDebugMacro(<< "GenerateOutputInformation()");
190 
191  unsigned int dimension = input->GetDimension();
192  if (dimension < 3)
193  {
194  mitkThrow() << "ImageCropper cannot handle 1D or 2D Objects.";
195  }
196 
197  if ((m_Geometry.IsNull()) || (m_Geometry->GetTimeGeometry()->CountTimeSteps() == 0))
198  return;
199 
200  mitk::BaseGeometry *bsGeometry = m_Geometry->GetGeometry();
201  mitk::BaseGeometry *inputImageGeometry = input->GetSlicedGeometry();
202 
203  // calculate bounding box
204  mitk::BoundingBox::Pointer bsBoxRelativeToImage =
205  bsGeometry->CalculateBoundingBoxRelativeToTransform(inputImageGeometry->GetIndexToWorldTransform());
206 
207  // pre-initialize input-requested-region to largest-possible-region
208  m_InputRequestedRegion = input->GetLargestPossibleRegion();
209  // build region out of bounding-box of index and size of the bounding box
210  mitk::SlicedData::IndexType index = m_InputRequestedRegion.GetIndex(); // init times and channels
211  mitk::BoundingBox::PointType min = bsBoxRelativeToImage->GetMinimum();
212  mitk::SlicedData::SizeType size = m_InputRequestedRegion.GetSize(); // init times and channels
213  mitk::BoundingBox::PointType max = bsBoxRelativeToImage->GetMaximum();
215  maxCorrected[0] = max[0];
216  maxCorrected[1] = max[1];
217  maxCorrected[2] = max[2];
218  maxCorrected[3] = input->GetDimensions()[3];
219  maxCorrected[4] = 0;
220 
221  for (unsigned int i = 0; i < dimension; i++)
222  {
223  index[i] = (mitk::SlicedData::IndexType::IndexValueType)(std::ceil(min[i]));
224  size[i] = (mitk::SlicedData::SizeType::SizeValueType)(std::ceil(maxCorrected[i]) - index[i]);
225  }
226  mitk::SlicedData::RegionType bsRegion(index, size);
227 
228  if (m_UseWholeInputRegion == false)
229  {
230  // crop input-requested-region with region of bounding-object
231  if (m_InputRequestedRegion.Crop(bsRegion) == false)
232  {
233  // crop not possible => do nothing: set time size to 0.
234  size.Fill(0);
235  m_InputRequestedRegion.SetSize(size);
236  bsRegion.SetSize(size);
237  mitkThrow() << "No overlap of the image and the cropping object.";
238  }
239  }
240 
241  // initialize output image
242  auto dimensions = new unsigned int[dimension];
243  if (dimension > 3 && !this->GetUseCropTimeStepOnly())
244  memcpy(dimensions + 3, input->GetDimensions() + 3, (dimension - 3) * sizeof(unsigned int));
245  else
246  dimension = 3; // set timeStep to zero if GetUseCropTimeStepOnly is true
247 
248  itk2vtk(m_InputRequestedRegion.GetSize(), dimensions);
249 
250  output->Initialize(mitk::PixelType(GetOutputPixelType()), dimension, dimensions);
251  delete[] dimensions;
252 
253  // Apply transform of the input image to the new generated output image
254  mitk::BoundingShapeCropper::RegionType outputRegion = output->GetRequestedRegion();
255 
256  m_TimeOfHeaderInitialization.Modified();
257  }
258 
260  {
261  // examine dimension and pixeltype
262  if ((image == nullptr) || (image->GetDimension() > 4) || (image->GetDimension() <= 2))
263  {
264  MITK_ERROR << "Filter cannot handle dimensions less than 2 and greater than 4" << std::endl;
265  itkExceptionMacro("Filter cannot handle dimensions less than 2 and greater than 4");
266  return;
267  }
268 
269  AccessByItk_1(image, CutImage, boTimeStep);
270  }
271 
273  {
274  MITK_INFO << "Generate Data" << std::endl;
275  mitk::Image::ConstPointer input = this->GetInput();
276  mitk::Image::Pointer output = this->GetOutput();
277 
278  if (input.IsNull())
279  return;
280 
281  if ((output->IsInitialized() == false) || (m_Geometry.IsNull()) ||
282  (m_Geometry->GetTimeGeometry()->CountTimeSteps() == 0))
283  return;
284 
285  m_InputTimeSelector->SetInput(input);
286  m_OutputTimeSelector->SetInput(this->GetOutput());
287 
288  mitk::BoundingShapeCropper::RegionType outputRegion = output->GetRequestedRegion();
289  mitk::BaseGeometry *inputImageGeometry = input->GetSlicedGeometry();
290 
291  // iterate over all time steps and perform cropping or masking on all or a specific timestep (perviously specified
292  // by UseCurrentTimeStepOnly)
293  int tstart = outputRegion.GetIndex(3);
294  int tmax = tstart + outputRegion.GetSize(3);
295 
296  if (this->m_UseCropTimeStepOnly)
297  {
298  mitk::SlicedGeometry3D *slicedGeometry = output->GetSlicedGeometry(tstart);
299  auto indexToWorldTransform = AffineTransform3D::New();
300  indexToWorldTransform->SetParameters(input->GetSlicedGeometry(tstart)->GetIndexToWorldTransform()->GetParameters());
301  slicedGeometry->SetIndexToWorldTransform(indexToWorldTransform);
302  const mitk::SlicedData::IndexType &start = m_InputRequestedRegion.GetIndex();
303  mitk::Point3D origin;
304  vtk2itk(start, origin);
305  inputImageGeometry->IndexToWorld(origin, origin);
306  slicedGeometry->SetOrigin(origin);
307  m_InputTimeSelector->SetTimeNr(m_CurrentTimeStep);
308  m_InputTimeSelector->UpdateLargestPossibleRegion();
309  m_OutputTimeSelector->SetTimeNr(tstart);
310  m_OutputTimeSelector->UpdateLargestPossibleRegion();
311  ComputeData(m_InputTimeSelector->GetOutput(), m_CurrentTimeStep);
312  }
313  else
314  {
315  int t;
316  for (t = tstart; t < tmax; ++t)
317  {
318  mitk::SlicedGeometry3D *slicedGeometry = output->GetSlicedGeometry(t);
319  auto indexToWorldTransform = AffineTransform3D::New();
320  indexToWorldTransform->SetParameters(input->GetSlicedGeometry(t)->GetIndexToWorldTransform()->GetParameters());
321  slicedGeometry->SetIndexToWorldTransform(indexToWorldTransform);
322  const mitk::SlicedData::IndexType &start = m_InputRequestedRegion.GetIndex();
323  mitk::Point3D origin;
324  vtk2itk(start, origin);
325  inputImageGeometry->IndexToWorld(origin, origin);
326  slicedGeometry->SetOrigin(origin);
327  m_InputTimeSelector->SetTimeNr(t);
328  m_InputTimeSelector->UpdateLargestPossibleRegion();
329  m_OutputTimeSelector->SetTimeNr(t);
330  m_OutputTimeSelector->UpdateLargestPossibleRegion();
331  ComputeData(m_InputTimeSelector->GetOutput(), t);
332  }
333  }
334  m_InputTimeSelector->SetInput(nullptr);
335  m_OutputTimeSelector->SetInput(nullptr);
336  m_TimeOfHeaderInitialization.Modified();
337  }
338 } // of namespace mitk
void IndexToWorld(const mitk::Vector3D &vec_units, mitk::Vector3D &vec_mm) const
Convert (continuous or discrete) index coordinates of a vector vec_units to world coordinates (in mm)...
virtual ScalarType GetOutsideValue()
void SetIndexToWorldTransform(mitk::AffineTransform3D *transform)
#define MITK_INFO
Definition: mitkLogMacros.h:18
const mitk::PixelType GetPixelType(int n=0) const
Returns the PixelType of channel n.
Definition: mitkImage.cpp:101
#define MITK_ERROR
Definition: mitkLogMacros.h:20
double ScalarType
DataCollection - Class to facilitate loading/accessing structured data.
#define AccessByItk_1(mitkImage, itkImageTypeFunction, arg1)
itk::Size< RegionDimension > SizeType
void DisplayErrorText(const char *t)
virtual bool GetUseCropTimeStepOnly()
unsigned int GetDimension() const
Get dimension of the image.
Definition: mitkImage.cpp:106
itk::ImageRegion< RegionDimension > RegionType
#define mitkThrow()
void SetOrigin(const Point3D &origin)
Set the origin, i.e. the upper-left corner of the plane.
Image class for storing images.
Definition: mitkImage.h:72
void CutImage(itk::Image< TPixel, VImageDimension > *inputItkImage, int timeStep)
Template Function for cropping and masking images with scalar pixel type.
Provides access to a volume at a specific time of the input image.
void vtk2itk(const Tin &in, Tout &out)
virtual void ComputeData(mitk::Image *input3D, int boTimeStep)
Process the image and create the output.
static T max(T x, T y)
Definition: svm.cpp:56
Data class only having a BaseGeometry but not containing any specific data.
mitk::Image::Pointer image
void GenerateInputRequestedRegion() override
Reimplemented from ImageToImageFilter.
Describes the geometry of a data object consisting of slices.
void GenerateTimeInInputRegion(const mitk::TimeGeometry *outputTimeGeometry, const TOutputRegion &outputRegion, const mitk::TimeGeometry *inputTimeGeometry, TInputRegion &inputRegion)
static StatusBar * GetInstance()
static method to get the GUI dependent StatusBar-instance so the methods DisplayText, etc. can be called No reference counting, cause of decentral static use!
static T min(T x, T y)
Definition: svm.cpp:53
static Pointer New()
void itk2vtk(const Tin &in, Tout &out)
InputImageType * GetInput(void)
void SetGeometry(const mitk::GeometryData *geometry)
Set geometry of the bounding object.
virtual bool IsInitialized() const
Check whether the data has been initialized, i.e., at least the Geometry and other header data has be...
virtual const PixelType GetOutputPixelType()
mitk::BoundingBox::Pointer CalculateBoundingBoxRelativeToTransform(const mitk::AffineTransform3D *transform) const
Calculates a bounding-box around the geometry relative to a coordinate system defined by a transform...
OutputType * GetOutput()
Get the output data of this image source object.
mitk::BaseGeometry * GetGeometry(int t=0) const
Return the geometry, which is a TimeGeometry, of the data as non-const pointer.
Definition: mitkBaseData.h:138
BaseGeometry Describes the geometry of a data object.
Class for defining the data type of pixels.
Definition: mitkPixelType.h:51
mitk::AffineTransform3D * GetIndexToWorldTransform()
Get the transformation used to convert from index to world coordinates.