Medical Imaging Interaction Toolkit  2018.4.99-12ad79a3
Medical Imaging Interaction Toolkit
mitkAutoCropImageFilter.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 
15 #include "mitkGeometry3D.h"
16 #include "mitkImageAccessByItk.h"
17 #include "mitkImageCast.h"
18 #include "mitkImageReadAccessor.h"
19 #include "mitkPlaneGeometry.h"
20 #include "mitkStatusBar.h"
21 
22 #include <itkImageRegionConstIterator.h>
23 #include <itkRegionOfInterestImageFilter.h>
25 
27  : m_BackgroundValue(0), m_MarginFactor(1.0), m_TimeSelector(nullptr), m_OverrideCroppingRegion(false)
28 {
29 }
30 
32 {
33 }
34 
35 template <typename TPixel, unsigned int VImageDimension>
36 void mitk::AutoCropImageFilter::ITKCrop3DImage(itk::Image<TPixel, VImageDimension> *inputItkImage,
37  unsigned int timestep)
38 {
39  if (inputItkImage == nullptr)
40  {
42  "An internal error occurred. Can't convert Image. Please report to bugs@mitk.org");
43  MITK_ERROR << "image is nullptr...returning" << std::endl;
44  return;
45  }
46 
47  typedef itk::Image<TPixel, VImageDimension> InternalImageType;
48  typedef typename InternalImageType::Pointer InternalImagePointer;
49 
50  typedef itk::RegionOfInterestImageFilter<InternalImageType, InternalImageType> ROIFilterType;
51  typedef typename itk::RegionOfInterestImageFilter<InternalImageType, InternalImageType>::Pointer ROIFilterPointer;
52 
53  InternalImagePointer outputItk = InternalImageType::New();
54 
55  ROIFilterPointer roiFilter = ROIFilterType::New();
56  roiFilter->SetInput(0, inputItkImage);
57  roiFilter->SetRegionOfInterest(this->GetCroppingRegion());
58  roiFilter->Update();
59  outputItk = roiFilter->GetOutput();
60  outputItk->DisconnectPipeline();
61 
62  mitk::Image::Pointer newMitkImage = mitk::Image::New();
63  mitk::CastToMitkImage(outputItk, newMitkImage);
64  MITK_INFO << "Crop-Output dimension: " << (newMitkImage->GetDimension() == 3)
65  << " Filter-Output dimension: " << this->GetOutput()->GetDimension() << " Timestep: " << timestep;
66 
67  mitk::ImageReadAccessor newMitkImgAcc(newMitkImage);
68  this->GetOutput()->SetVolume(newMitkImgAcc.GetData(), timestep);
69 }
70 
72 {
73  mitk::Image::Pointer input = const_cast<mitk::Image *>(this->GetInput());
74  mitk::Image::Pointer output = this->GetOutput();
75 
76  if (input->GetDimension() <= 2)
77  {
78  MITK_ERROR << "Only 3D any 4D images are supported." << std::endl;
79  return;
80  }
81 
83 
84  if ((output->IsInitialized()) && (output->GetPipelineMTime() <= m_TimeOfHeaderInitialization.GetMTime()))
85  return;
86 
87  itkDebugMacro(<< "GenerateOutputInformation()");
88 
89  // PART I: initialize input requested region. We do this already here (and not
90  // later when GenerateInputRequestedRegion() is called), because we
91  // also need the information to setup the output.
92 
93  // pre-initialize input-requested-region to largest-possible-region
94  // and correct time-region; spatial part will be cropped by
95  // bounding-box of bounding-object below
96  m_InputRequestedRegion = input->GetLargestPossibleRegion();
97 
98  // build region out of index and size calculated in ComputeNewImageBounds()
99 
101  index[0] = m_RegionIndex[0];
102  index[1] = m_RegionIndex[1];
103  index[2] = m_RegionIndex[2];
104  index[3] = m_InputRequestedRegion.GetIndex()[3];
105  index[4] = m_InputRequestedRegion.GetIndex()[4];
106 
108  size[0] = m_RegionSize[0];
109  size[1] = m_RegionSize[1];
110  size[2] = m_RegionSize[2];
111  size[3] = m_InputRequestedRegion.GetSize()[3];
112  size[4] = m_InputRequestedRegion.GetSize()[4];
113 
114  mitk::SlicedData::RegionType cropRegion(index, size);
115 
116  // crop input-requested-region with cropping region computed from the image data
117  if (m_InputRequestedRegion.Crop(cropRegion) == false)
118  {
119  // crop not possible => do nothing: set time size to 0.
120  size.Fill(0);
121  m_InputRequestedRegion.SetSize(size);
122  return;
123  }
124 
125  // set input-requested-region, because we access it later in
126  // GenerateInputRequestedRegion (there we just set the time)
127  input->SetRequestedRegion(&m_InputRequestedRegion);
128 
129  // PART II: initialize output image
130  unsigned int dimension = input->GetDimension();
131  auto dimensions = new unsigned int[dimension];
132  itk2vtk(m_InputRequestedRegion.GetSize(), dimensions);
133  if (dimension > 3)
134  memcpy(dimensions + 3, input->GetDimensions() + 3, (dimension - 3) * sizeof(unsigned int));
135 
136  // create basic slicedGeometry that will be initialized below
137  output->Initialize(mitk::PixelType(GetOutputPixelType()), dimension, dimensions);
138  delete[] dimensions;
139 
140  // clone the IndexToWorldTransform from the input, otherwise we will overwrite it, when adjusting the origin of the
141  // output image!!
142  itk::ScalableAffineTransform<mitk::ScalarType, 3>::Pointer cloneTransform =
143  itk::ScalableAffineTransform<mitk::ScalarType, 3>::New();
144  cloneTransform->Compose(input->GetGeometry()->GetIndexToWorldTransform());
145  output->GetGeometry()->SetIndexToWorldTransform(cloneTransform.GetPointer());
146 
147  // Position the output Image to match the corresponding region of the input image
148  mitk::SlicedGeometry3D *slicedGeometry = output->GetSlicedGeometry();
149  mitk::SlicedGeometry3D::Pointer inputGeometry = input->GetSlicedGeometry();
150  const mitk::SlicedData::IndexType &start = m_InputRequestedRegion.GetIndex();
151  mitk::Point3D origin;
152  vtk2itk(start, origin);
153  input->GetSlicedGeometry()->IndexToWorld(origin, origin);
154  slicedGeometry->SetOrigin(origin);
155 
156  // get the PlaneGeometry for the first slice of the original image
158  dynamic_cast<mitk::PlaneGeometry *>(inputGeometry->GetPlaneGeometry(0)->Clone().GetPointer());
159  assert(plane);
160 
161  // re-initialize the plane according to the new requirements:
162  // dimensions of the cropped image
163  // right- and down-vector as well as spacing do not change, so use the ones from
164  // input image
165  ScalarType dimX = output->GetDimensions()[0];
166  ScalarType dimY = output->GetDimensions()[1];
167  mitk::Vector3D right = plane->GetAxisVector(0);
168  mitk::Vector3D down = plane->GetAxisVector(1);
169  mitk::Vector3D spacing = plane->GetSpacing();
170  plane->InitializeStandardPlane(dimX, dimY, right, down, &spacing);
171  // set the new origin on the PlaneGeometry as well
172  plane->SetOrigin(origin);
173 
174  // re-initialize the slicedGeometry with the correct planeGeometry
175  // in order to get a fully initialized SlicedGeometry3D
176  slicedGeometry->InitializeEvenlySpaced(
177  plane, inputGeometry->GetSpacing()[2], output->GetSlicedGeometry()->GetSlices());
178 
179  mitk::TimeGeometry *timeSlicedGeometry = output->GetTimeGeometry();
180  auto *propTimeGeometry = dynamic_cast<ProportionalTimeGeometry *>(timeSlicedGeometry);
181  propTimeGeometry->Initialize(slicedGeometry, output->GetDimension(3));
182 
183  m_TimeOfHeaderInitialization.Modified();
184 
185  output->SetPropertyList(input->GetPropertyList()->Clone());
186 }
187 
189 {
190  mitk::Image::ConstPointer input = this->GetInput();
191  mitk::Image::Pointer output = this->GetOutput();
192 
193  if (input.IsNull())
194  return;
195 
196  if (input->GetDimension() <= 2)
197  {
198  MITK_ERROR << "Only 3D and 4D images supported";
199  return;
200  }
201 
202  if ((output->IsInitialized() == false))
203  return;
204 
205  if (m_TimeSelector.IsNull())
207 
208  m_TimeSelector->SetInput(input);
209 
210  mitk::SlicedData::RegionType outputRegion = input->GetRequestedRegion();
211 
212  int tstart = outputRegion.GetIndex(3);
213  int tmax = tstart + outputRegion.GetSize(3);
214 
215  for (int timestep = tstart; timestep < tmax; ++timestep)
216  {
217  m_TimeSelector->SetTimeNr(timestep);
218  m_TimeSelector->UpdateLargestPossibleRegion();
219 
220  AccessFixedDimensionByItk_1(m_TimeSelector->GetOutput(), ITKCrop3DImage, 3, timestep);
221  }
222 
223  // this->GetOutput()->Update(); // Not sure if this is necessary...
224 
225  m_TimeOfHeaderInitialization.Modified();
226 }
227 
229 {
230  mitk::Image::ConstPointer inputMitk = this->GetInput();
231 
233  {
234  for (unsigned int i = 0; i < 3; ++i)
235  {
236  m_RegionIndex[i] = m_CroppingRegion.GetIndex()[i];
237  m_RegionSize[i] = m_CroppingRegion.GetSize()[i];
238 
239  if (m_RegionIndex[i] >= static_cast<RegionType::IndexValueType>(inputMitk->GetDimension(i)))
240  {
241  itkExceptionMacro("Cropping index is not inside the image. " << std::endl
242  << "Index:"
243  << std::endl
244  << m_CroppingRegion.GetIndex()
245  << std::endl
246  << "Size:"
247  << std::endl
248  << m_CroppingRegion.GetSize());
249  }
250 
251  if (m_RegionIndex[i] + m_RegionSize[i] >= inputMitk->GetDimension(i))
252  {
253  m_RegionSize[i] = inputMitk->GetDimension(i) - m_RegionIndex[i];
254  }
255  }
256 
257  for (unsigned int i = 0; i < 3; ++i)
258  {
259  m_RegionIndex[i] = m_CroppingRegion.GetIndex()[i];
260  m_RegionSize[i] = m_CroppingRegion.GetSize()[i];
261  }
262  }
263  else
264  {
265  // Check if a 3D or 4D image is present
266  unsigned int timeSteps = 1;
267  if (inputMitk->GetDimension() == 4)
268  timeSteps = inputMitk->GetDimension(3);
269 
270  ImageType::IndexType minima, maxima;
271 
272  if (inputMitk->GetDimension() == 4)
273  {
274  // initialize with time step 0
276  m_TimeSelector->SetInput(inputMitk);
277  m_TimeSelector->SetTimeNr(0);
278  m_TimeSelector->UpdateLargestPossibleRegion();
279  inputMitk = m_TimeSelector->GetOutput();
280  }
281 
282  ImagePointer inputItk = ImageType::New();
283  mitk::CastToItkImage(inputMitk, inputItk);
284 
285  // it is assumed that all volumes in a time series have the same 3D dimensions
286  ImageType::RegionType origRegion = inputItk->GetLargestPossibleRegion();
287 
288  // Initialize min and max on the first (or only) time step
289  maxima = inputItk->GetLargestPossibleRegion().GetIndex();
290  minima[0] = inputItk->GetLargestPossibleRegion().GetSize()[0];
291  minima[1] = inputItk->GetLargestPossibleRegion().GetSize()[1];
292  minima[2] = inputItk->GetLargestPossibleRegion().GetSize()[2];
293 
294  typedef itk::ImageRegionConstIterator<ImageType> ConstIteratorType;
295 
296  for (unsigned int idx = 0; idx < timeSteps; ++idx)
297  {
298  // if 4D image, update time step and itk image
299  if (idx > 0)
300  {
301  m_TimeSelector->SetTimeNr(idx);
302  m_TimeSelector->UpdateLargestPossibleRegion();
303  inputMitk = m_TimeSelector->GetOutput();
304  mitk::CastToItkImage(inputMitk, inputItk);
305  }
306 
307  ConstIteratorType inIt(inputItk, origRegion);
308 
309  for (inIt.GoToBegin(); !inIt.IsAtEnd(); ++inIt)
310  {
311  float pix_val = inIt.Get();
312  if (fabs(pix_val - m_BackgroundValue) > mitk::eps)
313  {
314  for (int i = 0; i < 3; i++)
315  {
316  minima[i] = vnl_math_min((int)minima[i], (int)(inIt.GetIndex()[i]));
317  maxima[i] = vnl_math_max((int)maxima[i], (int)(inIt.GetIndex()[i]));
318  }
319  }
320  }
321  }
322 
323  typedef ImageType::RegionType::SizeType::SizeValueType SizeValueType;
324 
325  m_RegionSize[0] = (SizeValueType)(m_MarginFactor * (maxima[0] - minima[0] + 1));
326  m_RegionSize[1] = (SizeValueType)(m_MarginFactor * (maxima[1] - minima[1] + 1));
327  m_RegionSize[2] = (SizeValueType)(m_MarginFactor * (maxima[2] - minima[2] + 1));
328  m_RegionIndex = minima;
329 
330  m_RegionIndex[0] -= (m_RegionSize[0] - maxima[0] + minima[0] - 1) / 2;
331  m_RegionIndex[1] -= (m_RegionSize[1] - maxima[1] + minima[1] - 1) / 2;
332  m_RegionIndex[2] -= (m_RegionSize[2] - maxima[2] + minima[2] - 1) / 2;
333 
334  ImageType::RegionType cropRegion(m_RegionIndex, m_RegionSize);
335  origRegion.Crop(cropRegion);
336 
337  m_RegionSize[0] = origRegion.GetSize()[0];
338  m_RegionSize[1] = origRegion.GetSize()[1];
339  m_RegionSize[2] = origRegion.GetSize()[2];
340 
341  m_RegionIndex[0] = origRegion.GetIndex()[0];
342  m_RegionIndex[1] = origRegion.GetIndex()[1];
343  m_RegionIndex[2] = origRegion.GetIndex()[2];
344 
345  m_CroppingRegion = origRegion;
346  }
347 }
348 
350 {
351 }
352 
354 {
355  return this->GetInput()->GetPixelType();
356 }
357 
359 {
360  m_CroppingRegion = overrideRegion;
362 }
mitk::SlicedData::RegionType m_InputRequestedRegion
#define MITK_INFO
Definition: mitkLogMacros.h:18
virtual const PixelType GetOutputPixelType()
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
RegionType::IndexType m_RegionIndex
virtual void InitializeEvenlySpaced(mitk::PlaneGeometry *geometry2D, unsigned int slices)
Completely initialize this instance as evenly-spaced with slices parallel to the provided PlaneGeomet...
itk::Size< RegionDimension > SizeType
void DisplayErrorText(const char *t)
void SetCroppingRegion(RegionType overrideRegion)
itk::ImageRegion< RegionDimension > RegionType
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 vtk2itk(const Tin &in, Tout &out)
static Pointer New()
Describes the geometry of a data object consisting of slices.
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!
void itk2vtk(const Tin &in, Tout &out)
InputImageType * GetInput(void)
#define AccessFixedDimensionByItk_1(mitkImage, itkImageTypeFunction, dimension, arg1)
void CastToMitkImage(const itk::SmartPointer< ItkOutputImageType > &itkimage, itk::SmartPointer< mitk::Image > &mitkoutputimage)
Cast an itk::Image (with a specific type) to an mitk::Image.
Definition: mitkImageCast.h:74
void MITKCORE_EXPORT CastToItkImage(const mitk::Image *mitkImage, itk::SmartPointer< ItkOutputImageType > &itkOutputImage)
Cast an mitk::Image to an itk::Image with a specific type.
MITKCORE_EXPORT const ScalarType eps
void ITKCrop3DImage(itk::Image< TPixel, VImageDimension > *inputItkImage, unsigned int timestep)
OutputType * GetOutput()
Get the output data of this image source object.
Describes a two-dimensional, rectangular plane.
itk::Image< mitk::ScalarType, 3 > InternalImageType
virtual RegionType GetCroppingRegion()
ImageReadAccessor class to get locked read access for a particular image part.
void Initialize() override
Initilizes a new object with one time steps which contains an empty geometry.
const void * GetData() const
Gives const access to the data.
mitk::ImageTimeSelector::Pointer m_TimeSelector
static Pointer New()
Class for defining the data type of pixels.
Definition: mitkPixelType.h:51
itk::ImageRegion< 3 > RegionType