Medical Imaging Interaction Toolkit  2018.4.99-f51274ea
Medical Imaging Interaction Toolkit
mitkHeightFieldSurfaceClipImageFilter.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 "mitkImageTimeSelector.h"
15 #include "mitkProperties.h"
16 #include "mitkTimeHelper.h"
17 
18 #include "mitkImageAccessByItk.h"
19 #include "mitkImageToItk.h"
20 
21 #include <itkImageRegionConstIterator.h>
22 #include <itkImageRegionIteratorWithIndex.h>
23 #include <itkImageSliceConstIteratorWithIndex.h>
24 
25 #include <vtkCellLocator.h>
26 #include <vtkPolyData.h>
27 
28 #include <limits>
29 
30 namespace mitk
31 {
33  : m_ClippingMode(CLIPPING_MODE_CONSTANT),
34  m_ClippingConstant(0.0),
35  m_MultiplicationFactor(2.0),
36  m_MultiPlaneValue(2),
37  m_HeightFieldResolutionX(256),
38  m_HeightFieldResolutionY(256),
39  m_MaxHeight(1024.0)
40  {
41  this->SetNumberOfIndexedInputs(8);
42  this->SetNumberOfRequiredInputs(2);
43 
46  }
47 
50  {
51  this->SetNthInput(1, clippingSurface);
52  }
53 
55  {
56  if (planeList.size() > 7)
57  {
58  MITK_WARN << "Only 7 clipping planes are allowed!";
59  }
60 
61  for (unsigned int i = 0; i < planeList.size(); ++i)
62  {
63  this->SetNthInput(i + 1, planeList.at(i));
64  }
65  }
66 
68  {
69  return dynamic_cast<const Surface *>(itk::ProcessObject::GetInput(1));
70  }
71 
76  {
78  }
79 
81  {
83  }
84 
86  {
87  Image *outputImage = this->GetOutput();
88  Image *inputImage = this->GetInput(0);
89  const Surface *inputSurface = dynamic_cast<const Surface *>(this->GetInput(1));
90 
91  if (!outputImage->IsInitialized() || inputSurface == nullptr)
92  {
93  return;
94  }
95 
97 
98  GenerateTimeInInputRegion(outputImage, inputImage);
99  }
100 
102  {
103  const Image *inputImage = this->GetInput(0);
104  Image *outputImage = this->GetOutput();
105 
106  if (outputImage->IsInitialized() && (this->GetMTime() <= m_TimeOfHeaderInitialization.GetMTime()))
107  {
108  return;
109  }
110 
111  itkDebugMacro(<< "GenerateOutputInformation()");
112 
113  unsigned int i;
114  auto tmpDimensions = new unsigned int[inputImage->GetDimension()];
115 
116  for (i = 0; i < inputImage->GetDimension(); ++i)
117  {
118  tmpDimensions[i] = inputImage->GetDimension(i);
119  }
120 
121  outputImage->Initialize(
122  inputImage->GetPixelType(), inputImage->GetDimension(), tmpDimensions, inputImage->GetNumberOfChannels());
123 
124  delete[] tmpDimensions;
125 
126  outputImage->SetGeometry(static_cast<Geometry3D *>(inputImage->GetGeometry()->Clone().GetPointer()));
127 
128  outputImage->SetPropertyList(inputImage->GetPropertyList()->Clone());
129 
130  m_TimeOfHeaderInitialization.Modified();
131  }
132 
133  template <typename TPixel, unsigned int VImageDimension>
135  itk::Image<TPixel, VImageDimension> *inputItkImage,
136  HeightFieldSurfaceClipImageFilter *clipImageFilter,
137  vtkPolyData *clippingPolyData,
138  AffineTransform3D *imageToPlaneTransform)
139  {
140  typedef itk::Image<TPixel, VImageDimension> ItkInputImageType;
141  typedef itk::Image<TPixel, VImageDimension> ItkOutputImageType;
142 
143  typedef itk::ImageSliceConstIteratorWithIndex<ItkInputImageType> ItkInputImageIteratorType;
144  typedef itk::ImageRegionIteratorWithIndex<ItkOutputImageType> ItkOutputImageIteratorType;
145 
147  outputimagetoitk->SetInput(clipImageFilter->m_OutputTimeSelector->GetOutput());
148  outputimagetoitk->Update();
149 
150  typename ItkOutputImageType::Pointer outputItkImage = outputimagetoitk->GetOutput();
151  std::vector<double> test;
152 
153  // create the iterators
154  typename ItkInputImageType::RegionType inputRegionOfInterest = inputItkImage->GetLargestPossibleRegion();
155  ItkInputImageIteratorType inputIt(inputItkImage, inputRegionOfInterest);
156  ItkOutputImageIteratorType outputIt(outputItkImage, inputRegionOfInterest);
157 
158  // Get bounds of clipping data
159  clippingPolyData->ComputeBounds();
160  double *bounds = clippingPolyData->GetBounds();
161 
162  double xWidth = bounds[1] - bounds[0];
163  double yWidth = bounds[3] - bounds[2];
164 
165  // Create vtkCellLocator for clipping poly data
166  vtkCellLocator *cellLocator = vtkCellLocator::New();
167  cellLocator->SetDataSet(clippingPolyData);
168  cellLocator->CacheCellBoundsOn();
169  cellLocator->AutomaticOn();
170  cellLocator->BuildLocator();
171 
172  // Allocate memory for 2D image to hold the height field generated by
173  // projecting the clipping data onto the plane
174  auto heightField = new double[m_HeightFieldResolutionX * m_HeightFieldResolutionY];
175 
176  // Walk through height field and for each entry calculate height of the
177  // clipping poly data at this point by means of vtkCellLocator. The
178  // clipping data x/y bounds are used for converting from poly data space to
179  // image (height-field) space.
180  MITK_INFO << "Calculating Height Field..." << std::endl;
181  for (unsigned int y = 0; y < m_HeightFieldResolutionY; ++y)
182  {
183  for (unsigned int x = 0; x < m_HeightFieldResolutionX; ++x)
184  {
185  double p0[3], p1[3], surfacePoint[3], pcoords[3];
186  p0[0] = bounds[0] + xWidth * x / (double)m_HeightFieldResolutionX;
187  p0[1] = bounds[2] + yWidth * y / (double)m_HeightFieldResolutionY;
188  p0[2] = -m_MaxHeight;
189 
190  p1[0] = p0[0];
191  p1[1] = p0[1];
192  p1[2] = m_MaxHeight;
193 
194  double t, distance;
195  int subId;
196  if (cellLocator->IntersectWithLine(p0, p1, 0.1, t, surfacePoint, pcoords, subId))
197  {
198  distance = (2.0 * t - 1.0) * m_MaxHeight;
199  }
200  else
201  {
202  distance = -65536.0;
203  }
204  heightField[y * m_HeightFieldResolutionX + x] = distance;
205  itk::Image<double, 2>::IndexType index;
206  index[0] = x;
207  index[1] = y;
208  }
209  }
210 
211  // Walk through entire input image and for each point determine its distance
212  // from the x/y plane.
213  MITK_INFO << "Performing clipping..." << std::endl;
214 
215  auto factor = static_cast<TPixel>(clipImageFilter->m_MultiplicationFactor);
216  TPixel clippingConstant = clipImageFilter->m_ClippingConstant;
217 
218  inputIt.SetFirstDirection(0);
219  inputIt.SetSecondDirection(1);
220  // through all slices
221  for (inputIt.GoToBegin(), outputIt.GoToBegin(); !inputIt.IsAtEnd(); inputIt.NextSlice())
222  {
223  // through all lines of a slice
224  for (; !inputIt.IsAtEndOfSlice(); inputIt.NextLine())
225  {
226  // Transform the start(line) point from the image to the plane
227  Point3D imageP0, planeP0;
228  imageP0[0] = inputIt.GetIndex()[0];
229  imageP0[1] = inputIt.GetIndex()[1];
230  imageP0[2] = inputIt.GetIndex()[2];
231  planeP0 = imageToPlaneTransform->TransformPoint(imageP0);
232 
233  // Transform the end point (line) from the image to the plane
234  Point3D imageP1, planeP1;
235  imageP1[0] = imageP0[0] + inputRegionOfInterest.GetSize(0);
236  imageP1[1] = imageP0[1];
237  imageP1[2] = imageP0[2];
238  planeP1 = imageToPlaneTransform->TransformPoint(imageP1);
239 
240  // calculate the step size (if the plane is rotate, you go "crossway" through the image)
241  Vector3D step = (planeP1 - planeP0) / (double)inputRegionOfInterest.GetSize(0);
242 
243  // over all pixel
244  for (; !inputIt.IsAtEndOfLine(); ++inputIt, ++outputIt, planeP0 += step)
245  {
246  // Only ConstantMode: if image pixel value == constant mode value-->set output pixel value directly
247  if ((clipImageFilter->m_ClippingMode == CLIPPING_MODE_CONSTANT) &&
248  ((TPixel)inputIt.Get() == clippingConstant))
249  {
250  outputIt.Set(clippingConstant);
251  }
252 
253  else
254  {
255  auto x0 = (int)((double)(m_HeightFieldResolutionX) * (planeP0[0] - bounds[0]) / xWidth);
256  auto y0 = (int)((double)(m_HeightFieldResolutionY) * (planeP0[1] - bounds[2]) / yWidth);
257 
258  bool clip;
259 
260  // if the current point is outside of the plane region (RegionOfInterest)-->clip the pixel allways
261  if ((x0 < 0) || (x0 >= (int)m_HeightFieldResolutionX) || (y0 < 0) || (y0 >= (int)m_HeightFieldResolutionY))
262  {
263  clip = true;
264  }
265 
266  else
267  {
268  // Calculate bilinearly interpolated height field value at plane point
269  int x1 = x0 + 1;
270  int y1 = y0 + 1;
271  if (x1 >= (int)m_HeightFieldResolutionX)
272  {
273  x1 = x0;
274  }
275  if (y1 >= (int)m_HeightFieldResolutionY)
276  {
277  y1 = y0;
278  }
279 
280  // Get the neighbour points for the interpolation
281  ScalarType q00, q01, q10, q11;
282  q00 = heightField[y0 * m_HeightFieldResolutionX + x0];
283  q01 = heightField[y0 * m_HeightFieldResolutionX + x1];
284  q10 = heightField[y1 * m_HeightFieldResolutionX + x0];
285  q11 = heightField[y1 * m_HeightFieldResolutionX + x1];
286 
287  double p00 = ((double)(m_HeightFieldResolutionX) * (planeP0[0] - bounds[0]) / xWidth);
288  double p01 = ((double)(m_HeightFieldResolutionY) * (planeP0[1] - bounds[2]) / yWidth);
289 
290  ScalarType q =
291  q00 * ((double)x1 - p00) * ((double)y1 - p01) + q01 * (p00 - (double)x0) * ((double)y1 - p01) +
292  q10 * ((double)x1 - p00) * (p01 - (double)y0) + q11 * (p00 - (double)x0) * (p01 - (double)y0);
293 
294  if (q - planeP0[2] < 0)
295  {
296  clip = true;
297  }
298  else
299  {
300  clip = false;
301  }
302  }
303 
304  // different modes: differnt values for the clipped pixel
305  if (clip)
306  {
307  if (clipImageFilter->m_ClippingMode == CLIPPING_MODE_CONSTANT)
308  {
309  outputIt.Set(clipImageFilter->m_ClippingConstant);
310  }
311  else if (clipImageFilter->m_ClippingMode == CLIPPING_MODE_MULTIPLYBYFACTOR)
312  {
313  outputIt.Set(inputIt.Get() * factor);
314  }
315  else if (clipImageFilter->m_ClippingMode == CLIPPING_MODE_MULTIPLANE)
316  {
317  if (inputIt.Get() != 0)
318  outputIt.Set(inputIt.Get() + m_MultiPlaneValue);
319  else
320  outputIt.Set(inputIt.Get());
321  }
322  }
323  // the non-clipped pixel keeps his value
324  else
325  {
326  outputIt.Set(inputIt.Get());
327  }
328  }
329  }
330  }
331  }
332 
333  MITK_INFO << "DONE!" << std::endl;
334 
335  // Clean-up
336  cellLocator->Delete();
337  }
338 
340  {
341  const Image *inputImage = this->GetInput(0);
342 
343  const Image *outputImage = this->GetOutput();
344 
345  m_InputTimeSelector->SetInput(inputImage);
346  m_OutputTimeSelector->SetInput(outputImage);
347 
348  Image::RegionType outputRegion = outputImage->GetRequestedRegion();
349  const TimeGeometry *outputTimeGeometry = outputImage->GetTimeGeometry();
350  const TimeGeometry *inputTimeGeometry = inputImage->GetTimeGeometry();
351  ScalarType timeInMS;
352 
353  int timestep = 0;
354  int tstart = outputRegion.GetIndex(3);
355  int tmax = tstart + outputRegion.GetSize(3);
356 
357  for (unsigned int i = 1; i < this->GetNumberOfInputs(); ++i)
358  {
359  Surface *inputSurface = dynamic_cast<Surface *>(itk::ProcessObject::GetInput(i));
360 
361  if (!outputImage->IsInitialized() || inputSurface == nullptr)
362  return;
363 
364  MITK_INFO << "Plane: " << i;
365  MITK_INFO << "Clipping: Start\n";
366 
367  // const PlaneGeometry *clippingGeometryOfCurrentTimeStep = nullptr;
368 
369  int t;
370  for (t = tstart; t < tmax; ++t)
371  {
372  timeInMS = outputTimeGeometry->TimeStepToTimePoint(t);
373  timestep = inputTimeGeometry->TimePointToTimeStep(timeInMS);
374 
375  m_InputTimeSelector->SetTimeNr(timestep);
376  m_InputTimeSelector->UpdateLargestPossibleRegion();
377  m_OutputTimeSelector->SetTimeNr(t);
378  m_OutputTimeSelector->UpdateLargestPossibleRegion();
379 
380  // Compose IndexToWorld transform of image with WorldToIndexTransform of
381  // clipping data for conversion from image index space to plane index space
382  AffineTransform3D::Pointer planeWorldToIndexTransform = AffineTransform3D::New();
383  inputSurface->GetGeometry(t)->GetIndexToWorldTransform()->GetInverse(planeWorldToIndexTransform);
384 
385  AffineTransform3D::Pointer imageToPlaneTransform = AffineTransform3D::New();
386  imageToPlaneTransform->SetIdentity();
387 
388  imageToPlaneTransform->Compose(inputTimeGeometry->GetGeometryForTimeStep(t)->GetIndexToWorldTransform());
389  imageToPlaneTransform->Compose(planeWorldToIndexTransform);
390 
391  MITK_INFO << "Accessing ITK function...\n";
392  if (i == 1)
393  {
394  AccessByItk_3(m_InputTimeSelector->GetOutput(),
396  this,
397  inputSurface->GetVtkPolyData(t),
398  imageToPlaneTransform);
399  }
400  else
401  {
402  mitk::Image::Pointer extensionImage = m_OutputTimeSelector->GetOutput()->Clone();
404  extensionImage, _InternalComputeClippedImage, this, inputSurface->GetVtkPolyData(t), imageToPlaneTransform);
405  }
408  }
409  }
410 
411  m_TimeOfHeaderInitialization.Modified();
412  }
413 
414 } // namespace
void SetRequestedRegionToLargestPossibleRegion() override
Class for storing surfaces (vtkPolyData).
Definition: mitkSurface.h:28
#define AccessByItk_3(mitkImage, itkImageTypeFunction, arg1, arg2, arg3)
virtual vtkPolyData * GetVtkPolyData(unsigned int t=0) const
#define MITK_INFO
Definition: mitkLogMacros.h:18
virtual void Initialize(const mitk::PixelType &type, unsigned int dimension, const unsigned int *dimensions, unsigned int channels=1)
Definition: mitkImage.cpp:838
const mitk::PixelType GetPixelType(int n=0) const
Returns the PixelType of channel n.
Definition: mitkImage.cpp:101
double ScalarType
Follow Up Storage - Class to facilitate loading/accessing structured follow-up data.
Definition: testcase.h:28
Pointer Clone() const
DataCollection - Class to facilitate loading/accessing structured data.
void SetClippingMode(int mode)
Specifies whether clipped part of the image shall be replaced by a constant or multiplied by a user-s...
void _InternalComputeClippedImage(itk::Image< TPixel, VImageDimension > *itkImage, HeightFieldSurfaceClipImageFilter *clipImageFilter, vtkPolyData *clippingPolyData, AffineTransform3D *imageToPlaneTransform)
void SetGeometry(BaseGeometry *aGeometry3D) override
Sets a geometry to an image.
Definition: mitkImage.cpp:1321
void SetPropertyList(PropertyList *propertyList)
Set the data&#39;s property list.
const mitk::TimeGeometry * GetTimeGeometry() const
Return the TimeGeometry of the data as const pointer.
Definition: mitkBaseData.h:66
#define MITK_WARN
Definition: mitkLogMacros.h:19
unsigned int GetNumberOfChannels() const
Get the number of channels.
unsigned int GetDimension() const
Get dimension of the image.
Definition: mitkImage.cpp:106
void SetClippingSurface(Surface *clippingSurface)
Set/Get the surface defining a height field as a triangle mesh.
void SetClippingModeToMultiplyByFactor()
Specifies whether clipped part of the image shall be replaced by a constant or multiplied by a user-s...
void SetClippingModeToConstant()
Specifies whether clipped part of the image shall be replaced by a constant or multiplied by a user-s...
itk::ImageRegion< RegionDimension > RegionType
virtual TimeStepType TimePointToTimeStep(TimePointType timePoint) const =0
Converts a time point to the corresponding time step.
Image class for storing images.
Definition: mitkImage.h:72
virtual TimePointType TimeStepToTimePoint(TimeStepType timeStep) const =0
Converts a time step to a time point.
const Surface * GetClippingSurface() const
Set/Get the surface defining a height field as a triangle mesh.
virtual const RegionType & GetRequestedRegion() const
int GetClippingMode()
Specifies whether clipped part of the image shall be replaced by a constant or multiplied by a user-s...
itk::AffineGeometryFrame< ScalarType, 3 >::TransformType AffineTransform3D
void SetClippingSurfaces(ClippingPlaneList planeList)
Set/Get the surfaces defining a height field as a triangle mesh.
itk::TimeStamp m_TimeOfHeaderInitialization
Time when Header was last initialized.
void GenerateTimeInInputRegion(const mitk::TimeGeometry *outputTimeGeometry, const TOutputRegion &outputRegion, const mitk::TimeGeometry *inputTimeGeometry, TInputRegion &inputRegion)
mitk::PropertyList::Pointer GetPropertyList() const
Get the data&#39;s property list.
static Pointer New()
InputImageType * GetInput(void)
virtual bool IsInitialized() const
Check whether the data has been initialized, i.e., at least the Geometry and other header data has be...
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.
Filter for clipping an image with an height-field represented by an mitk::Surface.
mitk::BaseGeometry * GetGeometry(int t=0) const
Return the geometry, which is a TimeGeometry, of the data as non-const pointer.
Definition: mitkBaseData.h:143
static Pointer New()
mitk::AffineTransform3D * GetIndexToWorldTransform()
Get the transformation used to convert from index to world coordinates.