Medical Imaging Interaction Toolkit  2018.4.99-3e3f1a6e
Medical Imaging Interaction Toolkit
mitkSurfaceStampImageFilter.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 "mitkImageAccessByItk.h"
15 #include "mitkImageWriteAccessor.h"
16 #include "mitkTimeHelper.h"
17 #include <mitkImageCast.h>
18 
19 #include <vtkCell.h>
20 #include <vtkLinearTransform.h>
21 #include <vtkPolyData.h>
22 #include <vtkSmartPointer.h>
23 #include <vtkTransform.h>
24 #include <vtkTransformPolyDataFilter.h>
25 
26 #include <itkTriangleMeshToBinaryImageFilter.h>
27 
29  : m_MakeOutputBinary(false), m_OverwriteBackground(false), m_BackgroundValue(0.0), m_ForegroundValue(1.0)
30 {
31 }
32 
34 {
35 }
36 
38 {
39  mitk::Image *outputImage = this->GetOutput();
40  if ((outputImage->IsInitialized() == false))
41  return;
42 
43  GenerateTimeInInputRegion(outputImage, this->GetInput());
44 }
45 
47 {
48  mitk::Image::ConstPointer inputImage = this->GetInput();
49  mitk::Image::Pointer outputImage = this->GetOutput();
50 
51  itkDebugMacro(<< "GenerateOutputInformation()");
52 
53  if (inputImage.IsNull() || (inputImage->IsInitialized() == false) || (inputImage->GetTimeGeometry() == nullptr))
54  return;
55 
57  outputImage->Initialize(mitk::MakeScalarPixelType<unsigned char>(), *inputImage->GetTimeGeometry());
58  else
59  outputImage->Initialize(inputImage->GetPixelType(), *inputImage->GetTimeGeometry());
60 
61  outputImage->SetPropertyList(inputImage->GetPropertyList()->Clone());
62 }
63 
65 {
66  m_Surface = surface;
67 }
68 
70 {
71  mitk::Image::ConstPointer inputImage = this->GetInput();
72 
73  if (inputImage.IsNull())
74  return;
75 
76  mitk::Image::Pointer outputImage = this->GetOutput();
77  if (outputImage->IsInitialized() == false)
78  return;
79 
80  if (m_Surface.IsNull())
81  return;
82 
83  mitk::Image::RegionType outputRegion = outputImage->GetRequestedRegion();
84 
85  int tstart = outputRegion.GetIndex(3);
86  int tmax = tstart + outputRegion.GetSize(3);
87 
88  if (tmax > 0)
89  {
90  int t;
91  for (t = tstart; t < tmax; ++t)
92  {
93  this->SurfaceStamp(t);
94  }
95  }
96  else
97  {
98  this->SurfaceStamp(0);
99  }
100 }
101 
103 {
104  mitk::Image::Pointer inputImage = this->GetInput();
105 
106  const mitk::TimeGeometry *surfaceTimeGeometry = GetInput()->GetTimeGeometry();
107  const mitk::TimeGeometry *imageTimeGeometry = inputImage->GetTimeGeometry();
108 
109  // Convert time step from image time-frame to surface time-frame
110  mitk::TimePointType matchingTimePoint = imageTimeGeometry->TimeStepToTimePoint(time);
111  mitk::TimeStepType surfaceTimeStep = surfaceTimeGeometry->TimePointToTimeStep(matchingTimePoint);
112 
113  vtkPolyData *polydata = m_Surface->GetVtkPolyData(surfaceTimeStep);
114  if (!polydata)
115  mitkThrow() << "Polydata is null.";
116 
117  vtkSmartPointer<vtkTransformPolyDataFilter> transformFilter = vtkSmartPointer<vtkTransformPolyDataFilter>::New();
118  transformFilter->SetInputData(polydata);
119  // transformFilter->ReleaseDataFlagOn();
120 
121  vtkSmartPointer<vtkTransform> transform = vtkSmartPointer<vtkTransform>::New();
122  BaseGeometry::Pointer geometry = surfaceTimeGeometry->GetGeometryForTimeStep(surfaceTimeStep);
123 
124  transform->PostMultiply();
125  transform->Concatenate(geometry->GetVtkTransform()->GetMatrix());
126  // take image geometry into account. vtk-Image information will be changed to unit spacing and zero origin below.
127  BaseGeometry::Pointer imageGeometry = imageTimeGeometry->GetGeometryForTimeStep(time);
128 
129  transform->Concatenate(imageGeometry->GetVtkTransform()->GetLinearInverse());
130  transformFilter->SetTransform(transform);
131  transformFilter->Update();
132 
133  polydata = transformFilter->GetOutput();
134 
135  if (!polydata || !polydata->GetNumberOfPoints())
136  mitkThrow() << "Polydata retrieved from transformation is null or has no points.";
137 
138  MeshType::Pointer mesh = MeshType::New();
139  mesh->SetCellsAllocationMethod(MeshType::CellsAllocatedDynamicallyCellByCell);
140  unsigned int numberOfPoints = polydata->GetNumberOfPoints();
141  mesh->GetPoints()->Reserve(numberOfPoints);
142 
143  vtkPoints *points = polydata->GetPoints();
144 
145  MeshType::PointType point;
146  for (unsigned int i = 0; i < numberOfPoints; i++)
147  {
148  double *aux = points->GetPoint(i);
149  point[0] = aux[0];
150  point[1] = aux[1];
151  point[2] = aux[2];
152  mesh->SetPoint(i, point);
153  }
154 
155  // Load the polygons into the itk::Mesh
156  typedef MeshType::CellAutoPointer CellAutoPointerType;
157  typedef MeshType::CellType CellType;
158  typedef itk::TriangleCell<CellType> TriangleCellType;
159  typedef MeshType::PointIdentifier PointIdentifierType;
160  typedef MeshType::CellIdentifier CellIdentifierType;
161 
162  // Read the number of polygons
163  CellIdentifierType numberOfPolygons = 0;
164  numberOfPolygons = polydata->GetNumberOfPolys();
165 
166  PointIdentifierType numberOfCellPoints = 3;
167  CellIdentifierType i = 0;
168 
169  for (i = 0; i < numberOfPolygons; i++)
170  {
171  vtkIdList *cellIds;
172  vtkCell *vcell = polydata->GetCell(i);
173  cellIds = vcell->GetPointIds();
174 
175  CellAutoPointerType cell;
176  auto *triangleCell = new TriangleCellType;
177  PointIdentifierType k;
178  for (k = 0; k < numberOfCellPoints; k++)
179  {
180  triangleCell->SetPointId(k, cellIds->GetId(k));
181  }
182 
183  cell.TakeOwnership(triangleCell);
184  mesh->SetCell(i, cell);
185  }
186 
187  if (!mesh->GetNumberOfPoints())
188  mitkThrow() << "Generated itk mesh is empty.";
189 
190  if (m_MakeOutputBinary)
191  {
193  }
194  else
195  {
197  }
198 }
199 
201 {
202  auto *inputImage = this->GetInput();
203 
204  mitk::Image::Pointer outputImage = this->GetOutput();
205 
206  typedef itk::Image<unsigned char, 3> BinaryImageType;
207  BinaryImageType::Pointer itkInputImage;
208  mitk::CastToItkImage(inputImage, itkInputImage);
209 
210  typedef itk::TriangleMeshToBinaryImageFilter<MeshType, BinaryImageType> FilterType;
211 
212  FilterType::Pointer filter = FilterType::New();
213  filter->SetInput(mesh);
214  filter->SetInfoImage(itkInputImage);
215  filter->SetInsideValue(1);
216  filter->SetOutsideValue(0);
217  filter->Update();
218 
219  BinaryImageType::Pointer resultImage = filter->GetOutput();
220  resultImage->DisconnectPipeline();
221 
222  mitk::CastToMitkImage(resultImage, outputImage);
223 }
224 
225 template <typename TPixel>
226 void mitk::SurfaceStampImageFilter::SurfaceStampProcessing(itk::Image<TPixel, 3> *input, MeshType *mesh)
227 {
228  typedef itk::Image<TPixel, 3> ImageType;
229  typedef itk::Image<unsigned char, 3> BinaryImageType;
230 
231  typedef itk::TriangleMeshToBinaryImageFilter<mitk::SurfaceStampImageFilter::MeshType, BinaryImageType> FilterType;
232 
233  BinaryImageType::Pointer binaryInput = BinaryImageType::New();
234  binaryInput->SetSpacing(input->GetSpacing());
235  binaryInput->SetOrigin(input->GetOrigin());
236  binaryInput->SetDirection(input->GetDirection());
237  binaryInput->SetRegions(input->GetLargestPossibleRegion());
238  binaryInput->Allocate();
239  binaryInput->FillBuffer(0);
240 
241  FilterType::Pointer filter = FilterType::New();
242  filter->SetInput(mesh);
243  filter->SetInfoImage(binaryInput);
244  filter->SetInsideValue(1);
245  filter->SetOutsideValue(0);
246  filter->Update();
247 
248  BinaryImageType::Pointer resultImage = filter->GetOutput();
249  resultImage->DisconnectPipeline();
250 
251  mitk::Image::Pointer outputImage = this->GetOutput();
252  typename ImageType::Pointer itkOutputImage;
253  mitk::CastToItkImage(outputImage, itkOutputImage);
254 
255  typedef itk::ImageRegionConstIterator<BinaryImageType> BinaryIteratorType;
256  typedef itk::ImageRegionConstIterator<ImageType> InputIteratorType;
257  typedef itk::ImageRegionIterator<ImageType> OutputIteratorType;
258 
259  BinaryIteratorType sourceIter(resultImage, resultImage->GetLargestPossibleRegion());
260  sourceIter.GoToBegin();
261 
262  InputIteratorType inputIter(input, input->GetLargestPossibleRegion());
263  inputIter.GoToBegin();
264 
265  OutputIteratorType outputIter(itkOutputImage, itkOutputImage->GetLargestPossibleRegion());
266  outputIter.GoToBegin();
267 
268  typename ImageType::PixelType inputValue;
269  unsigned char sourceValue;
270 
271  auto fgValue = static_cast<typename ImageType::PixelType>(m_ForegroundValue);
272  auto bgValue = static_cast<typename ImageType::PixelType>(m_BackgroundValue);
273 
274  while (!sourceIter.IsAtEnd())
275  {
276  sourceValue = static_cast<unsigned char>(sourceIter.Get());
277  inputValue = static_cast<typename ImageType::PixelType>(inputIter.Get());
278 
279  if (sourceValue != 0)
280  outputIter.Set(fgValue);
281  else if (m_OverwriteBackground)
282  outputIter.Set(bgValue);
283  else
284  outputIter.Set(inputValue);
285 
286  ++sourceIter;
287  ++inputIter;
288  ++outputIter;
289  }
290 }
Class for storing surfaces (vtkPolyData).
Definition: mitkSurface.h:28
float k(1.0)
itk::Image< unsigned char, 3 > ImageType
const mitk::TimeGeometry * GetTimeGeometry() const
Return the TimeGeometry of the data as const pointer.
Definition: mitkBaseData.h:61
void SetSurface(mitk::Surface *surface)
itk::ImageRegion< RegionDimension > RegionType
virtual TimeStepType TimePointToTimeStep(TimePointType timePoint) const =0
Converts a time point to the corresponding time step.
#define mitkThrow()
void SurfaceStampProcessing(itk::Image< TPixel, 3 > *input, MeshType *mesh)
Image class for storing images.
Definition: mitkImage.h:72
virtual TimePointType TimeStepToTimePoint(TimeStepType timeStep) const =0
Converts a time step to a time point.
void SurfaceStampBinaryOutputProcessing(MeshType *mesh)
mitk::ScalarType TimePointType
std::vcl_size_t TimeStepType
void GenerateTimeInInputRegion(const mitk::TimeGeometry *outputTimeGeometry, const TOutputRegion &outputRegion, const mitk::TimeGeometry *inputTimeGeometry, TInputRegion &inputRegion)
itk::QuadEdgeMesh< double, 3 > MeshType
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.
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.