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