Medical Imaging Interaction Toolkit  2018.4.99-12ad79a3
Medical Imaging Interaction Toolkit
mitkExtrudedContour.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 
13 #include "mitkExtrudedContour.h"
14 #include "mitkBaseProcess.h"
15 #include "mitkNumericTypes.h"
17 
18 #include <vtkCellArray.h>
19 #include <vtkClipPolyData.h>
20 #include <vtkLinearExtrusionFilter.h>
21 #include <vtkPlanes.h>
22 #include <vtkPoints.h>
23 #include <vtkPolyData.h>
24 
25 #include <vtkDoubleArray.h>
26 #include <vtkFloatArray.h>
27 #include <vtkPlane.h>
28 #include <vtkPolygon.h>
29 
30 // vtkButterflySubdivisionFilter * subdivs;
31 #include <vtkCubeSource.h>
32 #include <vtkImageData.h>
33 #include <vtkLinearSubdivisionFilter.h>
34 #include <vtkSampleFunction.h>
35 #include <vtkTriangleFilter.h>
36 
38  : m_Contour(nullptr), m_ClippingGeometry(nullptr), m_AutomaticVectorGeneration(false), m_Decimate(nullptr)
39 {
41  timeGeometry->Initialize(1);
42  SetTimeGeometry(timeGeometry);
43 
44  FillVector3D(m_Vector, 0.0, 0.0, 1.0);
45  m_RightVector.Fill(0.0);
46 
47  m_ExtrusionFilter = vtkLinearExtrusionFilter::New();
48 
49  m_ExtrusionFilter->CappingOff();
50  m_ExtrusionFilter->SetExtrusionTypeToVectorExtrusion();
51 
52  double vtkvector[3] = {0, 0, 1};
53 
54  // set extrusion vector
55  m_ExtrusionFilter->SetVector(vtkvector);
56 
57  m_TriangleFilter = vtkTriangleFilter::New();
58  m_TriangleFilter->SetInputConnection(m_ExtrusionFilter->GetOutputPort());
59 
60  m_SubdivisionFilter = vtkLinearSubdivisionFilter::New();
61  m_SubdivisionFilter->SetInputConnection(m_TriangleFilter->GetOutputPort());
62  m_SubdivisionFilter->SetNumberOfSubdivisions(4);
63 
64  m_ClippingBox = vtkPlanes::New();
65  m_ClipPolyDataFilter = vtkClipPolyData::New();
66  m_ClipPolyDataFilter->SetInputConnection(m_SubdivisionFilter->GetOutputPort());
67  m_ClipPolyDataFilter->SetClipFunction(m_ClippingBox);
68  m_ClipPolyDataFilter->InsideOutOn();
69 
70  m_Polygon = vtkPolygon::New();
71 
73 }
74 
76 {
77  m_ClipPolyDataFilter->Delete();
78  m_ClippingBox->Delete();
79  m_SubdivisionFilter->Delete();
80  m_TriangleFilter->Delete();
81  m_ExtrusionFilter->Delete();
82  m_Polygon->Delete();
83 }
84 
85 bool mitk::ExtrudedContour::IsInside(const Point3D &worldPoint) const
86 {
87  static double polygonNormal[3] = {0.0, 0.0, 1.0};
88 
89  // project point onto plane
90  float xt[3];
91  itk2vtk(worldPoint, xt);
92 
93  xt[0] = worldPoint[0] - m_Origin[0];
94  xt[1] = worldPoint[1] - m_Origin[1];
95  xt[2] = worldPoint[2] - m_Origin[2];
96 
97  float dist = xt[0] * m_Normal[0] + xt[1] * m_Normal[1] + xt[2] * m_Normal[2];
98  xt[0] -= dist * m_Normal[0];
99  xt[1] -= dist * m_Normal[1];
100  xt[2] -= dist * m_Normal[2];
101 
102  double x[3];
103 
104  x[0] = xt[0] * m_Right[0] + xt[1] * m_Right[1] + xt[2] * m_Right[2];
105  x[1] = xt[0] * m_Down[0] + xt[1] * m_Down[1] + xt[2] * m_Down[2];
106  x[2] = 0;
107 
108  // determine whether it's in the selection loop and then evaluate point
109  // in polygon only if absolutely necessary.
110  if (x[0] >= this->m_ProjectedContourBounds[0] && x[0] <= this->m_ProjectedContourBounds[1] &&
111  x[1] >= this->m_ProjectedContourBounds[2] && x[1] <= this->m_ProjectedContourBounds[3] &&
112  this->m_Polygon->PointInPolygon(x,
113  m_Polygon->Points->GetNumberOfPoints(),
114  ((vtkDoubleArray *)this->m_Polygon->Points->GetData())->GetPointer(0),
115  (double *)this->m_ProjectedContourBounds,
116  polygonNormal) == 1)
117  return true;
118  else
119  return false;
120 }
121 
123 {
124  return -1.0;
125 }
126 
128 {
129  if (this->GetSource())
130  {
131  this->GetSource()->UpdateOutputInformation();
132  }
134  {
135  BuildGeometry();
136  BuildSurface();
137  }
138  // if ( ( m_CalculateBoundingBox ) && ( m_PolyDataSeries.size() > 0 ) )
139  // CalculateBoundingBox();
140 }
141 
143 {
144  if (m_Contour.IsNull())
145  {
146  SetVtkPolyData(nullptr);
147  return;
148  }
149 
150  // set extrusion contour
151  vtkPolyData *polyData = vtkPolyData::New();
152  vtkCellArray *polys = vtkCellArray::New();
153 
154  polys->InsertNextCell(m_Polygon->GetPointIds());
155  polyData->SetPoints(m_Polygon->GetPoints());
156 
157  // float vtkpoint[3];
158  // unsigned int i, numPts = m_Polygon->GetNumberOfPoints();
159  // for(i=0; i<numPts; ++i)
160  //{
161  // float * vtkpoint = this->m_Polygon->Points->GetPoint(i);
162 
163  // pointids[i]=loopPoints->InsertNextPoint(vtkpoint);
164  //}
165  // polys->InsertNextCell( i, pointids );
166  // delete [] pointids;
167 
168  // polyData->SetPoints( loopPoints );
169  polyData->SetPolys(polys);
170  polys->Delete();
171 
172  m_ExtrusionFilter->SetInputData(polyData);
173 
174  polyData->Delete();
175 
176  // set extrusion scale factor
177  m_ExtrusionFilter->SetScaleFactor(GetGeometry()->GetExtentInMM(2));
178  SetVtkPolyData(m_SubdivisionFilter->GetOutput());
179  // if(m_ClippingGeometry.IsNull())
180  //{
181  // SetVtkPolyData(m_SubdivisionFilter->GetOutput());
182  //}
183  // else
184  //{
185  // m_ClipPolyDataFilter->SetInput(m_SubdivisionFilter->GetOutput());
186  // mitk::BoundingBox::BoundsArrayType bounds=m_ClippingGeometry->GetBounds();
187 
188  // m_ClippingBox->SetBounds(bounds[0], bounds[1], bounds[2], bounds[3], bounds[4], bounds[5]);
189  // m_ClippingBox->SetTransform(GetGeometry()->GetVtkTransform());
190  // m_ClipPolyDataFilter->SetClipFunction(m_ClippingBox);
191  // m_ClipPolyDataFilter->SetValue(0);
192  // SetVtkPolyData(m_ClipPolyDataFilter->GetOutput());
193  //}
194 
195  m_LastCalculateExtrusionTime.Modified();
196 }
197 
199 {
200  if (m_Contour.IsNull())
201  return;
202 
203  // Initialize(1);
204 
205  Vector3D nullvector;
206  nullvector.Fill(0.0);
207  float xProj[3];
208  unsigned int i;
209  unsigned int numPts = 20; // m_Contour->GetNumberOfPoints();
210  mitk::Contour::PathPointer path = m_Contour->GetContourPath();
211  mitk::Contour::PathType::InputType cstart = path->StartOfInput();
212  mitk::Contour::PathType::InputType cend = path->EndOfInput();
213  mitk::Contour::PathType::InputType cstep = (cend - cstart) / numPts;
214  mitk::Contour::PathType::InputType ccur;
215 
216  // Part I: guarantee/calculate legal vectors
217 
218  m_Vector.Normalize();
220  // check m_Vector
222  {
223  if (m_AutomaticVectorGeneration == false)
224  itkWarningMacro("Extrusion vector is 0 (" << m_Vector << "); trying to use normal of polygon");
225 
226  vtkPoints *loopPoints = vtkPoints::New();
227  // mitk::Contour::PointsContainerIterator pointsIt = m_Contour->GetPoints()->Begin();
228  double vtkpoint[3];
229 
230  unsigned int i = 0;
231  for (i = 0, ccur = cstart; i < numPts; ++i, ccur += cstep)
232  {
233  itk2vtk(path->Evaluate(ccur), vtkpoint);
234  loopPoints->InsertNextPoint(vtkpoint);
235  }
236 
237  // Make sure points define a loop with a m_Normal
238  vtkPolygon::ComputeNormal(loopPoints, m_Normal);
239  loopPoints->Delete();
240 
242  if (mitk::Equal(m_Vector, nullvector))
243  {
244  itkExceptionMacro("Cannot calculate normal of polygon");
245  }
246  }
247  // check m_RightVector
248  if ((mitk::Equal(m_RightVector, nullvector)) || (mitk::Equal(m_RightVector * m_Vector, 0.0) == false))
249  {
250  if (mitk::Equal(m_RightVector, nullvector))
251  {
252  itkDebugMacro("Right vector is 0. Calculating.");
253  }
254  else
255  {
256  itkWarningMacro("Right vector (" << m_RightVector << ") not perpendicular to extrusion vector " << m_Vector
257  << ": "
258  << m_RightVector * m_Vector);
259  }
260  // calculate a legal m_RightVector
261  if (mitk::Equal(m_Vector[1], 0.0f) == false)
262  {
263  FillVector3D(m_RightVector, 1.0f, -m_Vector[0] / m_Vector[1], 0.0f);
264  m_RightVector.Normalize();
265  }
266  else
267  {
268  FillVector3D(m_RightVector, 0.0f, 1.0f, 0.0f);
269  }
270  }
271 
272  // calculate down-vector
273  VnlVector rightDV = m_RightVector.GetVnlVector();
274  rightDV.normalize();
275  vnl2vtk(rightDV, m_Right);
276  VnlVector downDV = vnl_cross_3d(m_Vector.GetVnlVector(), rightDV);
277  downDV.normalize();
278  vnl2vtk(downDV, m_Down);
279 
280  // Part II: calculate plane as base for extrusion, project the contour
281  // on this plane and store as polygon for IsInside test and BoundingBox calculation
282 
283  // initialize m_ProjectionPlane, yet with origin at 0
284  m_ProjectionPlane->InitializeStandardPlane(rightDV, downDV);
285 
286  // create vtkPolygon from contour and simultaneously determine 2D bounds of
287  // contour projected on m_ProjectionPlane
288  // mitk::Contour::PointsContainerIterator pointsIt = m_Contour->GetPoints()->Begin();
289  m_Polygon->Points->Reset();
290  m_Polygon->Points->SetNumberOfPoints(numPts);
291  m_Polygon->PointIds->Reset();
292  m_Polygon->PointIds->SetNumberOfIds(numPts);
293  mitk::Point2D pt2d;
294  mitk::Point3D pt3d;
298  xProj[2] = 0.0;
299  for (i = 0, ccur = cstart; i < numPts; ++i, ccur += cstep)
300  {
301  pt3d.CastFrom(path->Evaluate(ccur));
302  m_ProjectionPlane->Map(pt3d, pt2d);
303  xProj[0] = pt2d[0];
304  if (pt2d[0] < min[0])
305  min[0] = pt2d[0];
306  if (pt2d[0] > max[0])
307  max[0] = pt2d[0];
308  xProj[1] = pt2d[1];
309  if (pt2d[1] < min[1])
310  min[1] = pt2d[1];
311  if (pt2d[1] > max[1])
312  max[1] = pt2d[1];
313  m_Polygon->Points->SetPoint(i, xProj);
314  m_Polygon->PointIds->SetId(i, i);
315  }
316  // shift parametric origin to (0,0)
317  for (i = 0; i < numPts; ++i)
318  {
319  double *pt = this->m_Polygon->Points->GetPoint(i);
320 
321  pt[0] -= min[0];
322  pt[1] -= min[1];
323  itkDebugMacro(<< i << ": (" << pt[0] << "," << pt[1] << "," << pt[2] << ")");
324  }
325  this->m_Polygon->GetBounds(m_ProjectedContourBounds);
326  // m_ProjectedContourBounds[4]=-1.0; m_ProjectedContourBounds[5]=1.0;
327 
328  // calculate origin (except translation along the normal) and bounds
329  // of m_ProjectionPlane:
330  // origin is composed of the minimum x-/y-coordinates of the polygon,
331  // bounds from the extent of the polygon, both after projecting on the plane
332  mitk::Point3D origin;
333  m_ProjectionPlane->Map(min, origin);
334  ScalarType bounds[6] = {0, max[0] - min[0], 0, max[1] - min[1], 0, 1};
335  m_ProjectionPlane->SetBounds(bounds);
336  m_ProjectionPlane->SetOrigin(origin);
337 
338  // Part III: initialize geometry
339  if (m_ClippingGeometry.IsNotNull())
340  {
343  unsigned char i;
344  for (i = 0; i < 8; ++i)
345  {
346  dist = m_ProjectionPlane->SignedDistance(m_ClippingGeometry->GetCornerPoint(i));
347  if (dist < min_dist)
348  min_dist = dist;
349  if (dist > max_dist)
350  max_dist = dist;
351  }
352  // incorporate translation along the normal into origin
353  origin = origin + m_Vector * min_dist;
354  m_ProjectionPlane->SetOrigin(origin);
355  bounds[5] = max_dist - min_dist;
356  }
357  else
358  bounds[5] = 20;
359 
360  itk2vtk(origin, m_Origin);
361 
363  assert(g3d.IsNotNull());
364  g3d->SetBounds(bounds);
365  g3d->SetIndexToWorldTransform(m_ProjectionPlane->GetIndexToWorldTransform());
366 
368  timeGeometry->Initialize(g3d, 1);
369  SetTimeGeometry(timeGeometry);
370 }
371 
372 unsigned long mitk::ExtrudedContour::GetMTime() const
373 {
374  unsigned long latestTime = Superclass::GetMTime();
375  if (m_Contour.IsNotNull())
376  {
377  unsigned long localTime;
378  localTime = m_Contour->GetMTime();
379  if (localTime > latestTime)
380  latestTime = localTime;
381  }
382  return latestTime;
383 }
itk::SmartPointer< mitk::BaseDataSource > GetSource() const
Get the process object that generated this data object.
vnl_vector< ScalarType > VnlVector
Definition: mitkVector.h:134
double ScalarType
mitk::ScalarType GetVolume() override
vtkLinearSubdivisionFilter * m_SubdivisionFilter
virtual void SetTimeGeometry(TimeGeometry *geometry)
Set the TimeGeometry of the data, which will be referenced (not copied!).
vtkLinearExtrusionFilter * m_ExtrusionFilter
mitk::PlaneGeometry::Pointer m_ProjectionPlane
void FillVector3D(Tout &out, mitk::ScalarType x, mitk::ScalarType y, mitk::ScalarType z)
Definition: mitkArray.h:106
mitk::Vector3D m_RightVector
void vnl2vtk(const vnl_vector< Tin > &in, Tout *out)
mitk::BaseGeometry::Pointer m_ClippingGeometry
unsigned long GetMTime() const override
void vtk2itk(const Tin &in, Tout &out)
virtual void SetVtkPolyData(vtkPolyData *polydata, unsigned int t=0)
static T max(T x, T y)
Definition: svm.cpp:56
void UpdateOutputInformation() override
PathType::Pointer PathPointer
Definition: mitkContour.h:39
static T min(T x, T y)
Definition: svm.cpp:53
void itk2vtk(const Tin &in, Tout &out)
vtkClipPolyData * m_ClipPolyDataFilter
MITKNEWMODULE_EXPORT bool Equal(mitk::ExampleDataStructure *leftHandSide, mitk::ExampleDataStructure *rightHandSide, mitk::ScalarType eps, bool verbose)
Returns true if the example data structures are considered equal.
mitk::Contour::Pointer m_Contour
static Pointer New()
itk::TimeStamp m_LastCalculateExtrusionTime
bool IsInside(const Point3D &p) const override
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
vtkTriangleFilter * m_TriangleFilter
float m_Right[3]
For fast projection on plane.