1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkPolyhedron.cxx
5
6 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7 All rights reserved.
8 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9
10 This software is distributed WITHOUT ANY WARRANTY; without even
11 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12 PURPOSE. See the above copyright notice for more information.
13
14 =========================================================================*/
15
16 // Hide VTK_DEPRECATED_IN_9_0_0() warnings for this class.
17 #define VTK_DEPRECATION_LEVEL 0
18
19 #include "vtkPolyhedron.h"
20 #include "vtkCellArray.h"
21 #include "vtkCellData.h"
22 #include "vtkCellLocator.h"
23 #include "vtkEdgeTable.h"
24 #include "vtkGenericCell.h"
25 #include "vtkIdTypeArray.h"
26 #include "vtkLine.h"
27 #include "vtkMath.h"
28 #include "vtkMeanValueCoordinatesInterpolator.h"
29 #include "vtkOrderedTriangulator.h"
30 #include "vtkPointData.h"
31 #include "vtkPointLocator.h"
32 #include "vtkPolyData.h"
33 #include "vtkPolygon.h"
34 #include "vtkQuad.h"
35 #include "vtkTetra.h"
36 #include "vtkTriangle.h"
37 #include "vtkVector.h"
38
39 #include <cmath>
40 #include <functional>
41 #include <map>
42 #include <set>
43 #include <unordered_map>
44 #include <unordered_set>
45 #include <vector>
46
47 vtkStandardNewMacro(vtkPolyhedron);
48
49 // Special typedef
50 typedef std::vector<vtkIdType> vtkIdVectorType;
51 class vtkPointIdMap : public std::map<vtkIdType, vtkIdType>
52 {
53 };
54
55 // an edge consists of two id's and their order
56 // is *not* important. To that end special hash and
57 // equals functions have been made
58 // typedef std::pair<vtkIdType, vtkIdType> Edge;
59
60 struct Edge : public std::pair<vtkIdType, vtkIdType>
61 {
62 public:
63 Edge() = default;
EdgeEdge64 Edge(vtkIdType a, vtkIdType b)
65 : std::pair<vtkIdType, vtkIdType>(a, b)
66 {
67 }
EdgeEdge68 Edge(vtkCell* edge)
69 : std::pair<vtkIdType, vtkIdType>(edge->GetPointId(0), edge->GetPointId(1))
70 {
71 }
operator <<(ostream & stream,const Edge & e)72 friend ostream& operator<<(ostream& stream, const Edge& e)
73 {
74 stream << e.first << " - " << e.second;
75 return stream;
76 }
77 };
78
79 struct hash_fn
80 {
operator ()hash_fn81 size_t operator()(Edge const& p) const
82 {
83 size_t i = (size_t)p.first;
84 size_t j = (size_t)p.second;
85
86 // first make order-independent, i.e. hash(i,j) == hash(j,i)
87 if (i < j)
88 {
89 size_t tmp = i;
90 i = j;
91 j = tmp;
92 }
93
94 // then XOR both together , multiplied by two primes to try to prevent collisions
95 return (17 * i) ^ (31 * j);
96 }
97 };
98
99 struct equal_fn
100 {
operator ()equal_fn101 bool operator()(Edge const& e1, Edge const& e2) const
102 {
103 return (e1.first == e2.first && e1.second == e2.second) ||
104 (e1.second == e2.first && e1.first == e2.second);
105 }
106 };
107
108 // these typedefs are for the contouoring code. There the order of two edges does not matter
109 // so we use the specially crafted equals and hash functions defined above.
110
111 typedef std::vector<Edge> EdgeVector;
112
113 typedef std::vector<EdgeVector> FaceEdgesVector;
114 typedef std::unordered_map<Edge, std::set<vtkIdType>, hash_fn, equal_fn> EdgeFaceSetMap;
115
116 typedef std::unordered_multimap<vtkIdType, Edge> PointIndexEdgeMultiMap;
117 typedef std::unordered_map<Edge, vtkIdType, hash_fn, equal_fn> EdgePointIndexMap;
118
119 typedef std::unordered_set<Edge, hash_fn, equal_fn> EdgeSet;
120
121 typedef vtkIdVectorType Face;
122 typedef std::vector<Face> FaceVector;
123
124 //// Special class for iterating through polyhedron faces
125 ////----------------------------------------------------------------------------
126 class vtkPolyhedronFaceIterator
127 {
128 public:
129 vtkIdType CurrentPolygonSize;
130 vtkIdType* Polygon;
131 vtkIdType* Current;
132 vtkIdType NumberOfPolygons;
133 vtkIdType Id;
134
vtkPolyhedronFaceIterator(vtkIdType numFaces,vtkIdType * t)135 vtkPolyhedronFaceIterator(vtkIdType numFaces, vtkIdType* t)
136 {
137 this->CurrentPolygonSize = t[0];
138 this->Polygon = t;
139 this->Current = t + 1;
140 this->NumberOfPolygons = numFaces;
141 this->Id = 0;
142 }
operator ++()143 vtkIdType* operator++()
144 {
145 this->Current += this->CurrentPolygonSize + 1;
146 this->Polygon = this->Current - 1;
147 this->Id++;
148 if (this->Id < this->NumberOfPolygons)
149 {
150 this->CurrentPolygonSize = this->Polygon[0];
151 }
152 else
153 {
154 this->CurrentPolygonSize = VTK_ID_MAX;
155 }
156 return this->Current;
157 }
158 };
159
160 //------------------------------------------------------------------------------
161 // Construct the hexahedron with eight points.
vtkPolyhedron()162 vtkPolyhedron::vtkPolyhedron()
163 {
164 this->Line = vtkLine::New();
165 this->Triangle = vtkTriangle::New();
166 this->Quad = vtkQuad::New();
167 this->Polygon = vtkPolygon::New();
168 this->Tetra = vtkTetra::New();
169 this->GlobalFaces = vtkIdTypeArray::New();
170 this->FaceLocations = vtkIdTypeArray::New();
171 this->PointIdMap = new vtkPointIdMap;
172
173 this->EdgesGenerated = 0;
174 this->EdgeTable = vtkEdgeTable::New();
175 this->Edges = vtkIdTypeArray::New();
176 this->Edges->SetNumberOfComponents(2);
177 this->EdgeFaces = vtkIdTypeArray::New();
178 this->EdgeFaces->SetNumberOfComponents(2);
179
180 this->FacesGenerated = 0;
181 this->Faces = vtkIdTypeArray::New();
182
183 this->BoundsComputed = 0;
184
185 this->PolyDataConstructed = 0;
186 this->PolyData = vtkPolyData::New();
187 this->Polys = vtkCellArray::New();
188 this->LocatorConstructed = 0;
189 this->CellLocator = vtkCellLocator::New();
190 this->CellIds = vtkIdList::New();
191 this->Cell = vtkGenericCell::New();
192
193 this->ValenceAtPoint = nullptr;
194 }
195
196 //------------------------------------------------------------------------------
~vtkPolyhedron()197 vtkPolyhedron::~vtkPolyhedron()
198 {
199 if (this->ValenceAtPoint != nullptr)
200 {
201 delete this->ValenceAtPoint;
202 for (int i = 0; i < this->GetNumberOfPoints(); i++)
203 {
204 delete this->PointToIncidentFaces[i];
205 }
206 delete this->PointToIncidentFaces;
207 }
208 this->Line->Delete();
209 this->Triangle->Delete();
210 this->Quad->Delete();
211 this->Polygon->Delete();
212 this->Tetra->Delete();
213 this->GlobalFaces->Delete();
214 this->FaceLocations->Delete();
215 delete this->PointIdMap;
216 this->EdgeTable->Delete();
217 this->Edges->Delete();
218 this->EdgeFaces->Delete();
219 this->Faces->Delete();
220 this->PolyData->Delete();
221 this->Polys->Delete();
222 this->CellLocator->Delete();
223 this->CellIds->Delete();
224 this->Cell->Delete();
225 }
226
227 //------------------------------------------------------------------------------
ComputeBounds()228 void vtkPolyhedron::ComputeBounds()
229 {
230 if (this->BoundsComputed)
231 {
232 return;
233 }
234
235 this->Superclass::GetBounds(); // stored in this->Bounds
236 this->BoundsComputed = 1;
237 }
238
239 //------------------------------------------------------------------------------
ConstructPolyData()240 void vtkPolyhedron::ConstructPolyData()
241 {
242 if (this->PolyDataConstructed)
243 {
244 return;
245 }
246
247 // Here's a trick, we're going to use the Faces array as the connectivity
248 // array. Note that the Faces have an added nfaces value at the beginning
249 // of the array. Other than that,it's a vtkCellArray. So we play games
250 // with the pointers.
251 this->GenerateFaces();
252
253 if (this->Faces->GetNumberOfTuples() == 0)
254 {
255 return;
256 }
257
258 const vtkIdType numCells = *this->Faces->GetPointer(0);
259 const vtkIdType connSize = this->Faces->GetNumberOfValues() - numCells - 1;
260 this->Polys->AllocateExact(numCells, connSize);
261 this->Polys->ImportLegacyFormat(this->Faces->GetPointer(1), this->Faces->GetNumberOfValues() - 1);
262
263 // Standard setup
264 this->PolyData->Initialize();
265 this->PolyData->SetPoints(this->Points);
266 this->PolyData->SetPolys(this->Polys);
267
268 this->PolyDataConstructed = 1;
269 }
270
GetPolyData()271 vtkPolyData* vtkPolyhedron::GetPolyData()
272 {
273 if (!this->PolyDataConstructed)
274 {
275 this->ConstructPolyData();
276 }
277
278 return this->PolyData;
279 }
280 //------------------------------------------------------------------------------
ConstructLocator()281 void vtkPolyhedron::ConstructLocator()
282 {
283 if (this->LocatorConstructed)
284 {
285 return;
286 }
287
288 this->ConstructPolyData();
289
290 // With the polydata set up, we can assign it to the locator
291 this->CellLocator->Initialize();
292 this->CellLocator->SetDataSet(this->PolyData);
293 this->CellLocator->BuildLocator();
294
295 this->LocatorConstructed = 1;
296 }
297
298 //------------------------------------------------------------------------------
ComputeParametricCoordinate(const double x[3],double pc[3])299 void vtkPolyhedron::ComputeParametricCoordinate(const double x[3], double pc[3])
300 {
301 this->ComputeBounds();
302 double* bounds = this->Bounds;
303
304 pc[0] = (x[0] - bounds[0]) / (bounds[1] - bounds[0]);
305 pc[1] = (x[1] - bounds[2]) / (bounds[3] - bounds[2]);
306 pc[2] = (x[2] - bounds[4]) / (bounds[5] - bounds[4]);
307 }
308
309 //------------------------------------------------------------------------------
ComputePositionFromParametricCoordinate(const double pc[3],double x[3])310 void vtkPolyhedron::ComputePositionFromParametricCoordinate(const double pc[3], double x[3])
311 {
312 this->ComputeBounds();
313 double* bounds = this->Bounds;
314 x[0] = (1 - pc[0]) * bounds[0] + pc[0] * bounds[1];
315 x[1] = (1 - pc[1]) * bounds[2] + pc[1] * bounds[3];
316 x[2] = (1 - pc[2]) * bounds[4] + pc[2] * bounds[5];
317 }
318
319 //------------------------------------------------------------------------------
320 // Should be called by GetCell() prior to any other method invocation and after the
321 // points, point ids, and faces have been loaded.
Initialize()322 void vtkPolyhedron::Initialize()
323 {
324 // Clear out any remaining memory.
325 this->PointIdMap->clear();
326
327 // We need to create a reverse map from the point ids to their canonical cell
328 // ids. This is a fancy way of saying that we have to be able to rapidly go
329 // from a PointId[i] to the location i in the cell.
330 vtkIdType i, id, numPointIds = this->PointIds->GetNumberOfIds();
331 for (i = 0; i < numPointIds; ++i)
332 {
333 id = this->PointIds->GetId(i);
334 (*this->PointIdMap)[id] = i;
335 }
336
337 // Edges have to be reset
338 this->EdgesGenerated = 0;
339 this->EdgeTable->Reset();
340 this->Edges->Reset();
341 this->EdgeFaces->Reset();
342 this->Faces->Reset();
343
344 // Polys have to be reset
345 this->Polys->Reset();
346
347 // Faces may need renumbering later. This means converting the face ids from
348 // global ids to local, canonical ids.
349 this->FacesGenerated = 0;
350
351 // No bounds have been computed as of yet.
352 this->BoundsComputed = 0;
353
354 // No supplemental geometric stuff created
355 this->PolyDataConstructed = 0;
356 this->LocatorConstructed = 0;
357 }
358
359 //------------------------------------------------------------------------------
GetNumberOfEdges()360 int vtkPolyhedron::GetNumberOfEdges()
361 {
362 // Make sure edges have been generated.
363 if (!this->EdgesGenerated)
364 {
365 this->GenerateEdges();
366 }
367
368 return static_cast<int>(this->Edges->GetNumberOfTuples());
369 }
370
371 //------------------------------------------------------------------------------
372 // This method requires that GenerateEdges() is invoked beforehand.
GetEdge(int edgeId)373 vtkCell* vtkPolyhedron::GetEdge(int edgeId)
374 {
375 // Make sure edges have been generated.
376 if (!this->EdgesGenerated)
377 {
378 this->GenerateEdges();
379 }
380
381 // Make sure requested edge is within range
382 vtkIdType numEdges = this->Edges->GetNumberOfTuples();
383
384 if (edgeId < 0 || edgeId >= numEdges)
385 {
386 return nullptr;
387 }
388
389 // Return the requested edge
390 vtkIdType edge[2];
391 this->Edges->GetTypedTuple(edgeId, edge);
392
393 // Recall that edge tuples are stored in canonical numbering
394 for (int i = 0; i < 2; i++)
395 {
396 this->Line->PointIds->SetId(i, this->PointIds->GetId(edge[i]));
397 this->Line->Points->SetPoint(i, this->Points->GetPoint(edge[i]));
398 }
399
400 return this->Line;
401 }
402
403 //------------------------------------------------------------------------------
GenerateEdges()404 int vtkPolyhedron::GenerateEdges()
405 {
406 if (this->EdgesGenerated)
407 {
408 return this->Edges->GetNumberOfTuples();
409 }
410
411 // check the number of faces and return if there aren't any
412 if (this->GlobalFaces->GetNumberOfTuples() == 0 || this->GlobalFaces->GetValue(0) <= 0)
413 {
414 return 0;
415 }
416
417 // Loop over all faces, inserting edges into the table
418 vtkIdType* faces = this->GlobalFaces->GetPointer(0);
419 vtkIdType nfaces = faces[0];
420 vtkIdType* face = faces + 1;
421 vtkIdType fid, i, edge[2], npts, edgeFaces[2], edgeId;
422 edgeFaces[1] = -1;
423
424 this->EdgeTable->InitEdgeInsertion(this->Points->GetNumberOfPoints(), 1);
425 for (fid = 0; fid < nfaces; ++fid)
426 {
427 npts = face[0];
428 for (i = 1; i <= npts; ++i)
429 {
430 edge[0] = (*this->PointIdMap)[face[i]];
431 edge[1] = (*this->PointIdMap)[(i != npts ? face[i + 1] : face[1])];
432 edgeFaces[0] = fid;
433 if ((edgeId = this->EdgeTable->IsEdge(edge[0], edge[1])) == (-1))
434 {
435 edgeId = this->EdgeTable->InsertEdge(edge[0], edge[1]);
436 this->Edges->InsertNextTypedTuple(edge);
437 this->EdgeFaces->InsertTypedTuple(edgeId, edgeFaces);
438 }
439 else
440 {
441 this->EdgeFaces->SetComponent(edgeId, 1, fid);
442 }
443 }
444 face += face[0] + 1;
445 } // for all faces
446
447 // Okay all done
448 this->EdgesGenerated = 1;
449 return this->Edges->GetNumberOfTuples();
450 }
451
452 //------------------------------------------------------------------------------
GetNumberOfFaces()453 int vtkPolyhedron::GetNumberOfFaces()
454 {
455 // Make sure faces have been generated.
456 if (!this->FacesGenerated)
457 {
458 this->GenerateFaces();
459 }
460
461 if (this->GlobalFaces->GetNumberOfTuples() == 0)
462 {
463 return 0;
464 }
465
466 return static_cast<int>(this->GlobalFaces->GetValue(0));
467 }
468
469 //------------------------------------------------------------------------------
GenerateFaces()470 void vtkPolyhedron::GenerateFaces()
471 {
472 if (this->FacesGenerated)
473 {
474 return;
475 }
476
477 if (this->GlobalFaces->GetNumberOfTuples() == 0)
478 {
479 return;
480 }
481
482 // Basically we just run through the faces and change the global ids to the
483 // canonical ids using the PointIdMap.
484 this->Faces->SetNumberOfTuples(this->GlobalFaces->GetNumberOfTuples());
485 vtkIdType* gFaces = this->GlobalFaces->GetPointer(0);
486 vtkIdType* faces = this->Faces->GetPointer(0);
487 vtkIdType nfaces = gFaces[0];
488 faces[0] = nfaces;
489 vtkIdType* gFace = gFaces + 1;
490 vtkIdType* face = faces + 1;
491 vtkIdType fid, i, id, npts;
492
493 for (fid = 0; fid < nfaces; ++fid)
494 {
495 npts = gFace[0];
496 face[0] = npts;
497 for (i = 1; i <= npts; ++i)
498 {
499 id = (*this->PointIdMap)[gFace[i]];
500 face[i] = id;
501 }
502 gFace += gFace[0] + 1;
503 face += face[0] + 1;
504 } // for all faces
505
506 // Okay we've done the deed
507 this->FacesGenerated = 1;
508 }
509
510 //------------------------------------------------------------------------------
GetFace(int faceId)511 vtkCell* vtkPolyhedron::GetFace(int faceId)
512 {
513 if (faceId < 0 || faceId >= this->GlobalFaces->GetValue(0))
514 {
515 return nullptr;
516 }
517
518 this->GenerateFaces();
519
520 // Okay load up the polygon
521 vtkIdType i, p, loc = this->FaceLocations->GetValue(faceId);
522 vtkIdType* face = this->GlobalFaces->GetPointer(loc);
523
524 this->Polygon->PointIds->SetNumberOfIds(face[0]);
525 this->Polygon->Points->SetNumberOfPoints(face[0]);
526
527 // grab faces in global id space
528 for (i = 0; i < face[0]; ++i)
529 {
530 this->Polygon->PointIds->SetId(i, face[i + 1]);
531 p = (*this->PointIdMap)[face[i + 1]];
532 this->Polygon->Points->SetPoint(i, this->Points->GetPoint(p));
533 }
534
535 return this->Polygon;
536 }
537
538 //------------------------------------------------------------------------------
539 // Specify the faces for this cell.
SetFaces(vtkIdType * faces)540 void vtkPolyhedron::SetFaces(vtkIdType* faces)
541 {
542 // Set up face structure
543 this->GlobalFaces->Reset();
544 this->FaceLocations->Reset();
545
546 if (!faces)
547 {
548 return;
549 }
550
551 vtkIdType nfaces = faces[0];
552 this->FaceLocations->SetNumberOfValues(nfaces);
553
554 this->GlobalFaces->InsertNextValue(nfaces);
555 vtkIdType* face = faces + 1;
556 vtkIdType faceLoc = 1;
557 vtkIdType i, fid, npts;
558
559 for (fid = 0; fid < nfaces; ++fid)
560 {
561 npts = face[0];
562 this->GlobalFaces->InsertNextValue(npts);
563 for (i = 1; i <= npts; ++i)
564 {
565 this->GlobalFaces->InsertNextValue(face[i]);
566 }
567 this->FaceLocations->SetValue(fid, faceLoc);
568
569 faceLoc += face[0] + 1;
570 face = faces + faceLoc;
571 } // for all faces
572 }
573
574 //------------------------------------------------------------------------------
575 // Return the list of faces for this cell.
GetFaces()576 vtkIdType* vtkPolyhedron::GetFaces()
577 {
578 if (!this->GlobalFaces->GetNumberOfTuples())
579 {
580 return nullptr;
581 }
582
583 return this->GlobalFaces->GetPointer(0);
584 }
585
586 //------------------------------------------------------------------------------
IntersectWithLine(const double p1[3],const double p2[3],double tol,double & tMin,double xMin[3],double pc[3],int & subId)587 int vtkPolyhedron::IntersectWithLine(const double p1[3], const double p2[3], double tol,
588 double& tMin, double xMin[3], double pc[3], int& subId)
589 {
590 // It's easiest if this is done in canonical space
591 this->GenerateFaces();
592
593 // Loop over all the faces, intersecting them in turn.
594 vtkIdType* face = this->Faces->GetPointer(0);
595 vtkIdType nfaces = *face++;
596 vtkIdType npts, i, fid, numHits = 0;
597 double t = VTK_FLOAT_MAX;
598 double x[3];
599
600 tMin = VTK_FLOAT_MAX;
601 for (fid = 0; fid < nfaces; ++fid)
602 {
603 npts = face[0];
604 vtkIdType hit = 0;
605 switch (npts)
606 {
607 case 3: // triangle
608 for (i = 0; i < 3; i++)
609 {
610 this->Triangle->Points->SetPoint(i, this->Points->GetPoint(face[i + 1]));
611 this->Triangle->PointIds->SetId(i, face[i + 1]);
612 }
613 hit = this->Triangle->IntersectWithLine(p1, p2, tol, t, x, pc, subId);
614 break;
615 case 4: // quad
616 for (i = 0; i < 4; i++)
617 {
618 this->Quad->Points->SetPoint(i, this->Points->GetPoint(face[i + 1]));
619 this->Quad->PointIds->SetId(i, face[i + 1]);
620 }
621 hit = this->Quad->IntersectWithLine(p1, p2, tol, t, x, pc, subId);
622 break;
623 default: // general polygon
624 this->Polygon->GetPoints()->SetNumberOfPoints(npts);
625 this->Polygon->GetPointIds()->SetNumberOfIds(npts);
626 for (i = 0; i < npts; i++)
627 {
628 this->Polygon->Points->SetPoint(i, this->Points->GetPoint(face[i + 1]));
629 this->Polygon->PointIds->SetId(i, face[i + 1]);
630 }
631 hit = this->Polygon->IntersectWithLine(p1, p2, tol, t, x, pc, subId);
632 break;
633 }
634
635 // Update minimum hit
636 if (hit)
637 {
638 numHits++;
639 if (t < tMin)
640 {
641 tMin = t;
642 xMin[0] = x[0];
643 xMin[1] = x[1];
644 xMin[2] = x[2];
645 }
646 }
647
648 face += face[0] + 1;
649 } // for all faces
650
651 // Compute parametric coordinates
652 this->ComputeParametricCoordinate(xMin, pc);
653
654 return (numHits > 0);
655 }
656
657 static const int VTK_MAX_ITER = 10; // Maximum iterations for ray-firing
658 static const int VTK_VOTE_THRESHOLD = 3;
659
660 //------------------------------------------------------------------------------
661 // Shoot random rays and count the number of intersections
IsInside(const double x[3],double tolerance)662 int vtkPolyhedron::IsInside(const double x[3], double tolerance)
663 {
664 // do a quick bounds check
665 this->ComputeBounds();
666 double* bounds = this->Bounds;
667 if (x[0] < bounds[0] || x[0] > bounds[1] || x[1] < bounds[2] || x[1] > bounds[3] ||
668 x[2] < bounds[4] || x[2] > bounds[5])
669 {
670 return 0;
671 }
672
673 // It's easiest if these computations are done in canonical space
674 this->GenerateFaces();
675
676 // This algorithm is adaptive; if there are enough faces in this
677 // polyhedron, a cell locator is built to accelerate intersections.
678 // Otherwise brute force looping over cells is used.
679 vtkIdType* faceArray = this->Faces->GetPointer(0);
680 vtkIdType nfaces = *faceArray++;
681 if (nfaces > 25)
682 {
683 this->ConstructLocator();
684 }
685
686 // We need a length to normalize the computations
687 double length = sqrt(this->Superclass::GetLength2());
688
689 // Perform in/out by shooting random rays. Multiple rays are fired
690 // to improve accuracy of the result.
691 //
692 // The variable iterNumber counts the number of rays fired and is
693 // limited by the defined variable VTK_MAX_ITER.
694 //
695 // The variable deltaVotes keeps track of the number of votes for
696 // "in" versus "out" of the surface. When deltaVotes > 0, more votes
697 // have counted for "in" than "out". When deltaVotes < 0, more votes
698 // have counted for "out" than "in". When the delta_vote exceeds or
699 // equals the defined variable VTK_VOTE_THRESHOLD, then the
700 // appropriate "in" or "out" status is returned.
701 //
702 double rayMag, ray[3], xray[3], t, pcoords[3], xint[3];
703 int i, numInts, iterNumber, deltaVotes, subId;
704 vtkIdType idx, numCells;
705 double tol = tolerance * length;
706
707 for (deltaVotes = 0, iterNumber = 1;
708 (iterNumber < VTK_MAX_ITER) && (std::abs(deltaVotes) < VTK_VOTE_THRESHOLD); iterNumber++)
709 {
710 // Define a random ray to fire.
711 do
712 {
713 for (i = 0; i < 3; i++)
714 {
715 ray[i] = vtkMath::Random(-1.0, 1.0);
716 }
717 rayMag = vtkMath::Norm(ray);
718 } while (rayMag == 0.0);
719
720 // The ray must be appropriately sized wrt the bounding box. (It has to go
721 // all the way through the bounding box.)
722 for (i = 0; i < 3; i++)
723 {
724 xray[i] = x[i] + (length / rayMag) * ray[i];
725 }
726
727 // Intersect the line with each of the candidate cells
728 numInts = 0;
729
730 if (this->LocatorConstructed)
731 {
732 // Retrieve the candidate cells from the locator
733 this->CellLocator->FindCellsAlongLine(x, xray, tol, this->CellIds);
734 numCells = this->CellIds->GetNumberOfIds();
735
736 for (idx = 0; idx < numCells; idx++)
737 {
738 this->PolyData->GetCell(this->CellIds->GetId(idx), this->Cell);
739 if (this->Cell->IntersectWithLine(x, xray, tol, t, xint, pcoords, subId))
740 {
741 // Check for vertex, edge or face intersections
742 // count the number of 0 or 1 pcoords
743 int pcount = 0;
744 for (int p = 0; p < 3; ++p)
745 {
746 if (pcoords[p] == 0.0 || pcoords[p] == 1.0)
747 {
748 pcount++;
749 }
750 }
751 // pcount = 1, exact face intersection
752 // pcount = 2, exact edge intersection
753 // pcount = 3, exact vertex intersection
754 if (pcount == 0)
755 {
756 numInts++;
757 }
758 }
759 } // for all candidate cells
760 }
761 else
762 {
763 numCells = nfaces;
764 this->ConstructPolyData();
765
766 for (idx = 0; idx < numCells; idx++)
767 {
768 this->PolyData->GetCell(idx, this->Cell);
769 if (this->Cell->IntersectWithLine(x, xray, tol, t, xint, pcoords, subId))
770 {
771 // Check for vertex, edge or face intersections
772 // count the number of 0 or 1 pcoords
773 int pcount = 0;
774 for (int p = 0; p < 3; ++p)
775 {
776 if (pcoords[p] == 0.0 || pcoords[p] == 1.0)
777 {
778 pcount++;
779 }
780 }
781 // pcount = 1, exact face intersection
782 // pcount = 2, exact edge intersection
783 // pcount = 3, exact vertex intersection
784 if (pcount == 0)
785 {
786 numInts++;
787 }
788 }
789 } // for all candidate cells
790 }
791
792 // Count the result
793 if (numInts != 0 && (numInts % 2) == 0)
794 {
795 --deltaVotes;
796 }
797 else
798 {
799 ++deltaVotes;
800 }
801 } // try another ray
802
803 // If the number of votes is positive, the point is inside
804 //
805 return (deltaVotes < 0 ? 0 : 1);
806 }
807
808 //------------------------------------------------------------------------------
809 // Determine whether or not a polyhedron is convex. This method is adapted
810 // from Devillers et al., "Checking the Convexity of Polytopes and the
811 // Planarity of Subdivisions", Computational Geometry, Volume 11, Issues 3 - 4,
812 // December 1998, Pages 187 - 208.
IsConvex()813 bool vtkPolyhedron::IsConvex()
814 {
815 double x[2][3], n[3], c[3], c0[3], c1[3], c0p[3], c1p[3], n0[3], n1[3];
816 double np[3], tmp0, tmp1;
817 vtkIdType i, w[2], edgeId, edgeFaces[2], loc, v, *face, r = 0;
818 const double eps = FLT_EPSILON;
819
820 std::vector<double> p(this->PointIds->GetNumberOfIds());
821 vtkIdVectorType d(this->PointIds->GetNumberOfIds(), 0);
822
823 // initialization
824 this->GenerateEdges();
825 this->GenerateFaces();
826 this->ConstructPolyData();
827 this->ComputeBounds();
828
829 // loop over all edges in the polyhedron
830 this->EdgeTable->InitTraversal();
831 while ((edgeId = this->EdgeTable->GetNextEdge(w[0], w[1])) >= 0)
832 {
833 // get the edge points
834 this->Points->GetPoint(w[0], x[0]);
835 this->Points->GetPoint(w[1], x[1]);
836
837 // get the local face ids
838 this->EdgeFaces->GetTypedTuple(edgeId, edgeFaces);
839
840 // get the face vertex ids for the first face
841 loc = this->FaceLocations->GetValue(edgeFaces[0]);
842 face = this->Faces->GetPointer(loc);
843
844 // compute the centroid and normal for the first face
845 vtkPolygon::ComputeCentroid(this->Points, face[0], face + 1, c0);
846 vtkPolygon::ComputeNormal(this->Points, face[0], face + 1, n0);
847
848 // get the face vertex ids for the second face
849 loc = this->FaceLocations->GetValue(edgeFaces[1]);
850 face = this->Faces->GetPointer(loc);
851
852 // compute the centroid and normal for the second face
853 vtkPolygon::ComputeCentroid(this->Points, face[0], face + 1, c1);
854 vtkPolygon::ComputeNormal(this->Points, face[0], face + 1, n1);
855
856 // check for local convexity (the average of the two centroids must be
857 // "below" both faces, as defined by their outward normals).
858 for (i = 0; i < 3; i++)
859 {
860 c[i] = (c1[i] + c0[i]) * .5;
861 c0p[i] = c[i] - c0[i];
862 c1p[i] = c[i] - c1[i];
863 }
864
865 if (vtkMath::Dot(n0, c0p) > 0. || vtkMath::Dot(n1, c1p) > 0.)
866 {
867 return false;
868 }
869
870 // check if the edge is a seam edge
871 // 1. the edge must not be vertical
872 // 2. the two faces must lie on the same side of a vertical plane
873 // 3. the upper face must not be vertical
874
875 // 1. simply check that the unit normal along the seam has x or y
876 // components
877 for (i = 0; i < 3; i++)
878 {
879 n[i] = x[1][i] - x[0][i];
880 }
881 vtkMath::Normalize(n);
882 if (std::abs(n[0]) < eps && std::abs(n[1]) < eps)
883 {
884 continue;
885 }
886
887 // 2. we need a plane through the seam and through a vector parallel to the
888 // z axis (or, more accurately, we need a vector perpendicular to this
889 // plane). This vector can be computed using the cross product between
890 // the a vector along the edge, and the vertical axis.
891 np[0] = +n[1];
892 np[1] = -n[0];
893 np[2] = 0;
894
895 for (i = 0; i < 3; i++)
896 {
897 c[i] = (x[1][i] + x[0][i]) * .5;
898 c0p[i] = c0[i] - c[i];
899 c1p[i] = c1[i] - c[i];
900 }
901
902 // if the vectors from the seam centroid to the face centroid are in the
903 // same direction relative to the plane, then condition 2 is satisfied.
904 tmp0 = vtkMath::Dot(np, c0p);
905 tmp1 = vtkMath::Dot(np, c1p);
906
907 if ((tmp0 < 0.) != (tmp1 < 0.))
908 {
909 continue;
910 }
911
912 // 3. We get the z component of the normal of the highest face
913 // If this is null, the face is in the vertical plane
914 tmp0 = c0[2] > c1[2] ? n0[2] : n1[2];
915 if (std::abs(tmp0) < eps)
916 {
917 continue;
918 }
919
920 // at this point, we know we have a seam edge. We now look at each vertex
921 // in the seam and determine whether or not it is a right-2-seam vertex. A
922 // convex polytope has exactly one right-2-seam vertex.
923 for (i = 0; i < 2; i++)
924 {
925 v = w[i];
926
927 // are there already 2 seams associated with this vertex? If so, then the
928 // projection of the polytope onto the x-y plane would have multiple seams
929 // emanating from the vertex => non-convex.
930 if (d[v] == 2)
931 {
932 return false;
933 }
934
935 // is this the first time that this vertex has been associated with a
936 // seam? If so, increment its seam count and record the x-coordinate of
937 // the adjacent edge vertex.
938 if (d[v] == 0)
939 {
940 d[v]++;
941 p[v] = x[(i + 1) % 2][0];
942 }
943 else
944 {
945 d[v]++;
946 // is v a right-2-seam vertex (i.e. is the x-value of v larger than the
947 // x-values of both u and p[v])?
948 if (x[i][0] > x[(i + 1) % 2][0] && x[i][0] > p[v])
949 {
950 // is this the first right-2-seam vertex?
951 if (r == 0)
952 {
953 r++;
954 }
955 else
956 {
957 return false;
958 }
959 }
960 }
961 }
962 }
963
964 return true;
965 }
966
967 //------------------------------------------------------------------------------
CellBoundary(int vtkNotUsed (subId),const double pcoords[3],vtkIdList * pts)968 int vtkPolyhedron::CellBoundary(int vtkNotUsed(subId), const double pcoords[3], vtkIdList* pts)
969 {
970 double x[3], n[3], o[3], v[3];
971 double dist, minDist = VTK_DOUBLE_MAX;
972 vtkIdType numFacePts = -1;
973 vtkIdType* facePts = nullptr;
974
975 // compute coordinates
976 this->ComputePositionFromParametricCoordinate(pcoords, x);
977
978 vtkPolyhedronFaceIterator faceIter(this->GetNumberOfFaces(), this->Faces->GetPointer(1));
979 while (faceIter.Id < faceIter.NumberOfPolygons)
980 {
981 if (faceIter.CurrentPolygonSize < 3)
982 {
983 vtkErrorMacro("Find a face with "
984 << faceIter.CurrentPolygonSize
985 << " vertices. Cannot return CellBoundary due to this degenerate case.");
986 break;
987 }
988
989 vtkPolygon::ComputeNormal(this->Points, faceIter.CurrentPolygonSize, faceIter.Current, n);
990 vtkMath::Normalize(n);
991 this->Points->GetPoint(faceIter.Current[0], o);
992 v[0] = x[0] - o[0];
993 v[1] = x[1] - o[1];
994 v[2] = x[2] - o[2];
995 dist = std::abs(vtkMath::Dot(v, n));
996 if (dist < minDist)
997 {
998 minDist = dist;
999 numFacePts = faceIter.CurrentPolygonSize;
1000 facePts = faceIter.Current;
1001 }
1002
1003 ++faceIter;
1004 }
1005
1006 pts->Reset();
1007 if (numFacePts > 0)
1008 {
1009 for (vtkIdType i = 0; i < numFacePts; i++)
1010 {
1011 pts->InsertNextId(this->PointIds->GetId(facePts[i]));
1012 }
1013 }
1014
1015 // determine whether point is inside of polygon
1016 if (pcoords[0] >= 0.0 && pcoords[0] <= 1.0 && pcoords[1] >= 0.0 && pcoords[1] <= 1.0 &&
1017 pcoords[2] >= 0.0 && pcoords[2] <= 1.0 &&
1018 (this->IsInside(x, std::numeric_limits<double>::infinity())))
1019 {
1020 return 1;
1021 }
1022 else
1023 {
1024 return 0;
1025 }
1026 }
1027
1028 //------------------------------------------------------------------------------
EvaluatePosition(const double x[3],double closestPoint[3],int & vtkNotUsed (subId),double pcoords[3],double & minDist2,double weights[])1029 int vtkPolyhedron::EvaluatePosition(const double x[3], double closestPoint[3],
1030 int& vtkNotUsed(subId), double pcoords[3], double& minDist2, double weights[])
1031 {
1032 // compute parametric coordinates
1033 this->ComputeParametricCoordinate(x, pcoords);
1034
1035 // construct polydata, the result is stored in this->PolyData,
1036 // the cell array is stored in this->Polys
1037 this->ConstructPolyData();
1038
1039 // Construct cell locator
1040 this->ConstructLocator();
1041
1042 // find closest point and store the squared distance
1043 vtkIdType cellId;
1044 int id;
1045 double cp[3];
1046 this->Cell->Initialize();
1047 this->CellLocator->FindClosestPoint(x, cp, this->Cell, cellId, id, minDist2);
1048
1049 if (closestPoint)
1050 {
1051 closestPoint[0] = cp[0];
1052 closestPoint[1] = cp[1];
1053 closestPoint[2] = cp[2];
1054 }
1055
1056 // get the MVC weights
1057 this->InterpolateFunctions(x, weights);
1058
1059 // set distance to be zero, if point is inside
1060 int isInside = this->IsInside(x, std::numeric_limits<double>::infinity());
1061 if (isInside)
1062 {
1063 minDist2 = 0.0;
1064 }
1065
1066 return isInside;
1067 }
1068
1069 //------------------------------------------------------------------------------
EvaluateLocation(int & vtkNotUsed (subId),const double pcoords[3],double x[3],double * weights)1070 void vtkPolyhedron::EvaluateLocation(
1071 int& vtkNotUsed(subId), const double pcoords[3], double x[3], double* weights)
1072 {
1073 this->ComputePositionFromParametricCoordinate(pcoords, x);
1074
1075 this->InterpolateFunctions(x, weights);
1076 }
1077
1078 //------------------------------------------------------------------------------
Derivatives(int vtkNotUsed (subId),const double pcoords[3],const double * values,int dim,double * derivs)1079 void vtkPolyhedron::Derivatives(
1080 int vtkNotUsed(subId), const double pcoords[3], const double* values, int dim, double* derivs)
1081 {
1082 int i, j, k, idx;
1083 for (j = 0; j < dim; j++)
1084 {
1085 for (i = 0; i < 3; i++)
1086 {
1087 derivs[j * dim + i] = 0.0;
1088 }
1089 }
1090
1091 static const double Sample_Offset_In_Parameter_Space = 0.01;
1092
1093 double x[4][3];
1094 double coord[3];
1095
1096 // compute positions of point and three offset sample points
1097 coord[0] = pcoords[0];
1098 coord[1] = pcoords[1];
1099 coord[2] = pcoords[2];
1100 this->ComputePositionFromParametricCoordinate(coord, x[0]);
1101
1102 coord[0] += Sample_Offset_In_Parameter_Space;
1103 this->ComputePositionFromParametricCoordinate(coord, x[1]);
1104 coord[0] = pcoords[0];
1105
1106 coord[1] += Sample_Offset_In_Parameter_Space;
1107 this->ComputePositionFromParametricCoordinate(coord, x[2]);
1108 coord[1] = pcoords[1];
1109
1110 coord[2] += Sample_Offset_In_Parameter_Space;
1111 this->ComputePositionFromParametricCoordinate(coord, x[3]);
1112 coord[2] = pcoords[2];
1113
1114 this->ConstructPolyData();
1115 int numVerts = this->PolyData->GetNumberOfPoints();
1116
1117 double* weights = new double[numVerts];
1118 double* sample = new double[dim * 4];
1119 // for each sample point, sample data values
1120 for (idx = 0, k = 0; k < 4; k++) // loop over three sample points
1121 {
1122 this->InterpolateFunctions(x[k], weights);
1123 for (j = 0; j < dim; j++, idx++) // over number of derivates requested
1124 {
1125 sample[idx] = 0.0;
1126 for (i = 0; i < numVerts; i++)
1127 {
1128 sample[idx] += weights[i] * values[j + i * dim];
1129 }
1130 }
1131 }
1132
1133 double v1[3], v2[3], v3[3];
1134 double l1, l2, l3;
1135 // compute differences along the two axes
1136 for (i = 0; i < 3; i++)
1137 {
1138 v1[i] = x[1][i] - x[0][i];
1139 v2[i] = x[2][i] - x[0][i];
1140 v3[i] = x[3][i] - x[0][i];
1141 }
1142 l1 = vtkMath::Normalize(v1);
1143 l2 = vtkMath::Normalize(v2);
1144 l3 = vtkMath::Normalize(v3);
1145
1146 // compute derivatives along x-y-z axes
1147 double ddx, ddy, ddz;
1148 for (j = 0; j < dim; j++)
1149 {
1150 ddx = (sample[dim + j] - sample[j]) / l1;
1151 ddy = (sample[2 * dim + j] - sample[j]) / l2;
1152 ddz = (sample[3 * dim + j] - sample[j]) / l3;
1153
1154 // project onto global x-y-z axes
1155 derivs[3 * j] = ddx * v1[0] + ddy * v2[0] + ddz * v3[0];
1156 derivs[3 * j + 1] = ddx * v1[1] + ddy * v2[1] + ddz * v3[1];
1157 derivs[3 * j + 2] = ddx * v1[2] + ddy * v2[2] + ddz * v3[2];
1158 }
1159
1160 delete[] weights;
1161 delete[] sample;
1162 }
1163
1164 //------------------------------------------------------------------------------
GetParametricCoords()1165 double* vtkPolyhedron::GetParametricCoords()
1166 {
1167 return nullptr;
1168 }
1169
1170 //------------------------------------------------------------------------------
InterpolateFunctions(const double x[3],double * sf)1171 void vtkPolyhedron::InterpolateFunctions(const double x[3], double* sf)
1172 {
1173 // construct polydata, the result is stored in this->PolyData,
1174 // the cell array is stored in this->Polys
1175 this->ConstructPolyData();
1176
1177 // compute the weights
1178 if (!this->PolyData->GetPoints())
1179 {
1180 return;
1181 }
1182 vtkMeanValueCoordinatesInterpolator::ComputeInterpolationWeights(
1183 x, this->PolyData->GetPoints(), this->Polys, sf);
1184 }
1185
1186 //------------------------------------------------------------------------------
InterpolateDerivs(const double x[3],double * derivs)1187 void vtkPolyhedron::InterpolateDerivs(const double x[3], double* derivs)
1188 {
1189 (void)x;
1190 (void)derivs;
1191 }
1192
1193 //------------------------------------------------------------------------------
Triangulate(int vtkNotUsed (index),vtkIdList * ptIds,vtkPoints * pts)1194 int vtkPolyhedron::Triangulate(int vtkNotUsed(index), vtkIdList* ptIds, vtkPoints* pts)
1195 {
1196 ptIds->Reset();
1197 pts->Reset();
1198
1199 if (!this->GetPoints() || !this->GetNumberOfPoints())
1200 {
1201 return 0;
1202 }
1203
1204 this->ComputeBounds();
1205
1206 // use ordered triangulator to triangulate the polyhedron.
1207 vtkSmartPointer<vtkOrderedTriangulator> triangulator =
1208 vtkSmartPointer<vtkOrderedTriangulator>::New();
1209
1210 triangulator->InitTriangulation(this->Bounds, this->GetNumberOfPoints());
1211 triangulator->PreSortedOff();
1212
1213 double point[3];
1214 for (vtkIdType i = 0; i < this->GetNumberOfPoints(); i++)
1215 {
1216 this->GetPoints()->GetPoint(i, point);
1217 triangulator->InsertPoint(i, point, point, 0);
1218 }
1219 triangulator->Triangulate();
1220
1221 triangulator->AddTetras(0, ptIds, pts);
1222
1223 // convert to global Ids
1224 vtkIdType* ids = ptIds->GetPointer(0);
1225 for (vtkIdType i = 0; i < ptIds->GetNumberOfIds(); i++)
1226 {
1227 ids[i] = this->PointIds->GetId(ids[i]);
1228 }
1229
1230 return 1;
1231 }
1232
1233 //------------------------------------------------------------------------------
GeneratePointToIncidentFacesAndValenceAtPoint()1234 void vtkPolyhedron::GeneratePointToIncidentFacesAndValenceAtPoint()
1235 {
1236 // Allocate memory
1237 this->PointToIncidentFaces = new vtkIdType*[this->GetNumberOfPoints()];
1238 this->ValenceAtPoint = new vtkIdType[this->GetNumberOfPoints()];
1239 // Add the faces that hold each cell local point id
1240 std::vector<std::set<vtkIdType>> setFacesOfPoint(this->GetNumberOfPoints());
1241 for (int faceIndex = 0; faceIndex < this->GetNumberOfFaces(); faceIndex++)
1242 {
1243 auto face = this->GetFace(faceIndex);
1244 // For each point of the face
1245 for (int pointIndexFace = 0; pointIndexFace < face->GetNumberOfPoints(); pointIndexFace++)
1246 {
1247 // Get the global id of the point of the face
1248 auto pointId = face->GetPointId(pointIndexFace);
1249 // Transform the global id of the point of the face to the local id of the point in the cell
1250 auto pointCellLocalId = (*this->PointIdMap)[pointId];
1251 // Insert this face in the set
1252 setFacesOfPoint[pointCellLocalId].insert(faceIndex);
1253 }
1254 }
1255 // Fill in ValenceAtPoint and PointToIncidentFaces using the set data
1256 for (int pointIndex = 0; pointIndex < this->GetNumberOfPoints(); pointIndex++)
1257 {
1258 this->ValenceAtPoint[pointIndex] = static_cast<vtkIdType>(setFacesOfPoint[pointIndex].size());
1259 this->PointToIncidentFaces[pointIndex] = new vtkIdType[ValenceAtPoint[pointIndex]];
1260 int indexInsert = 0;
1261 for (auto faceId : setFacesOfPoint[pointIndex])
1262 {
1263 this->PointToIncidentFaces[pointIndex][indexInsert++] = faceId;
1264 }
1265 }
1266 }
1267
1268 //------------------------------------------------------------------------------
GetPointToIncidentFaces(vtkIdType pointId,const vtkIdType * & faceIds)1269 vtkIdType vtkPolyhedron::GetPointToIncidentFaces(vtkIdType pointId, const vtkIdType*& faceIds)
1270 {
1271 assert(pointId < this->GetNumberOfPoints() && "pointId too large");
1272 if (this->ValenceAtPoint == nullptr)
1273 {
1274 this->GeneratePointToIncidentFacesAndValenceAtPoint();
1275 }
1276 faceIds = this->PointToIncidentFaces[pointId];
1277 return this->ValenceAtPoint[pointId];
1278 }
1279
IntersectWithContour(vtkCell * cell,vtkDataArray * pointScalars,vtkPointIdMap * pointIdMap,double value,std::function<bool (double,double)> & compare,bool & allTrue)1280 bool IntersectWithContour(vtkCell* cell, vtkDataArray* pointScalars, vtkPointIdMap* pointIdMap,
1281 double value, std::function<bool(double, double)>& compare, bool& allTrue)
1282 {
1283 allTrue = true;
1284 bool allFalse = true;
1285
1286 int nPoints = cell->GetNumberOfPoints();
1287 for (int i = 0; i < nPoints; ++i)
1288 {
1289 vtkIdType globalPid = cell->GetPointId(i);
1290 vtkIdType localPid = pointIdMap->find(globalPid)->second;
1291
1292 double pointValue = pointScalars->GetTuple1(localPid);
1293
1294 if (compare(pointValue, value))
1295 {
1296 allFalse = false;
1297 }
1298 else
1299 {
1300 allTrue = false;
1301 }
1302 }
1303
1304 return !(allTrue || allFalse);
1305 }
1306
1307 // start new contouring code //
1308
1309 /*
1310
1311 This code contains a new way of polyhedral contouring. The approach is as follows:
1312 each of the polyhedron faces is triangulated (independent on normal orientation).
1313
1314 After triangulation, the contouring will give exactly 0 or 1 lines across the
1315 (tri-)faces. This allows for a straightforward face-edge-contourpoint walking to
1316 create one or more closed contour polygons.
1317
1318 The face-edge walking starts a given contour point. Using a lookup structure,
1319 the edge of the contour point is used to find an unvisited face in the list of
1320 the two faces that border the edge. The edges of that face are then searched
1321 to find the other edge with a contour point. These two contour points then define
1322 one contour line. The walking procedure stops when the starting contour point
1323 is reached again. The collection of lines forms a closed polyhedron.
1324
1325 */
1326
CheckWatertightNonManifoldPolyhedron(vtkPolyhedron * cell,EdgeSet & originalEdges)1327 bool CheckWatertightNonManifoldPolyhedron(vtkPolyhedron* cell, EdgeSet& originalEdges)
1328 {
1329 EdgeFaceSetMap directMap;
1330 int nFaces = cell->GetNumberOfFaces();
1331 for (int i = 0; i < nFaces; ++i)
1332 {
1333 vtkCell* face = cell->GetFace(i);
1334 for (int j = 0; j < face->GetNumberOfEdges(); ++j)
1335 {
1336 Edge e(face->GetEdge(j));
1337 originalEdges.insert(e);
1338
1339 auto at = directMap.find(e);
1340 if (at == directMap.end())
1341 {
1342 std::set<vtkIdType> facesOfEdge;
1343 facesOfEdge.insert(i);
1344 directMap.insert(make_pair(e, facesOfEdge));
1345 }
1346 else
1347 {
1348 std::set<vtkIdType>& facesOfEdge = at->second;
1349 facesOfEdge.insert(i);
1350 }
1351 }
1352 }
1353
1354 size_t nEdges = cell->GetNumberOfEdges();
1355 size_t sizeMap = directMap.size();
1356 if (sizeMap != nEdges)
1357 {
1358 vtkGenericWarningMacro(
1359 << "The number of edges in the edge>face map does not match the number of edges of the cell");
1360 return false;
1361 }
1362
1363 bool ok = true;
1364
1365 for (const auto& entry : directMap)
1366 {
1367 const std::set<vtkIdType>& facesOfEdge = entry.second;
1368 if (facesOfEdge.size() != 2)
1369 {
1370 vtkGenericWarningMacro(
1371 << "The polyhedron is not watertight or non-manifold because the number of faces of edge "
1372 << entry.first.first << "-" << entry.first.second << " is not 2 but "
1373 << facesOfEdge.size());
1374 ok = false;
1375 }
1376 }
1377
1378 return ok;
1379 }
1380
1381 /*
1382
1383 When directly triangulating the polyhedron faces that are not simple triangles or quads (i.e.
1384 they're polygons), a problem can occur which gives the resulting triangulated polyhedron
1385 non-manifold triangle faces
1386
1387 For example:
1388
1389 0 ----- 1 ----- 2
1390 | | |
1391 | | |
1392 | 6 |
1393 | | |
1394 | | |
1395 3 ----- 4 ----- 5
1396
1397 this can be triangulated as (0,1,6), (0,6,3), (3,6,4) and (1,2,6), (6,2,5), (6,5,4) (that would be
1398 OK) OR triangulated as (0,1,4), (0,4,3), (1,6,4) and (1,2,5), (1,5,4), (1,6,4) (that would
1399 be NOT OK because of the duplicate (1,6,4) triangle)
1400
1401 In fact, the ear-clipping polygon triangulation can produce, depending on the geometry,
1402 the *unwanted* triangulation instead of the desired one because it prioritizes triangles with
1403 inner angles close to 60 degrees, even though it then ends with a triangle with a very large
1404 internal angle (up to 180 degrees).
1405
1406 Therefore the preferred approach is to triangulate a polygon using a fan triangulation that gives
1407 the smallest range of internal angles. This approach will always choose to triangulate starting at
1408 (6) in the example given above. If (6) is moved out-of-plane as it were (see TestPolyhedron5.cxx)
1409 then the tetrahedralization gives a face triangulation that includes the edge (1)-(4), but
1410 triangulates the face as (1-4-2)-(2-4-5). The now preferred method triangulates it as
1411 (6-2-1)-(6-2-5)-(6-5-4), thereby preserving the original shape of the polygon, even if it is
1412 slightly concave. Note that extremely concave polygons will give completely incorrect triangulations
1413 with this method, but that would also be problematic for the tetrahedralization approach.
1414
1415 */
1416 // by using an *ordered* set, the triangles are consistently ordered, independent of face normal
1417
FindLowestIndex(vtkIdType n,vtkIdType * arr)1418 int FindLowestIndex(vtkIdType n, vtkIdType* arr)
1419 {
1420 int lowest(-1);
1421 vtkIdType min(VTK_ID_MAX);
1422 for (int i = 0; i < n; ++i)
1423 {
1424 if (arr[i] < min)
1425 {
1426 lowest = i;
1427 min = arr[i];
1428 }
1429 }
1430 return lowest;
1431 }
1432
FindLowestNeighbor(vtkIdType n,vtkIdType * arr,int idx,bool & mustReverse)1433 void FindLowestNeighbor(vtkIdType n, vtkIdType* arr, int idx, bool& mustReverse)
1434 {
1435 idx += n; // add n to prevent negative remainders
1436 vtkIdType left = arr[(idx - 1) % n];
1437 vtkIdType right = arr[(idx + 1) % n];
1438 if (left < right)
1439 {
1440 mustReverse = true;
1441 }
1442 else if (left > right)
1443 {
1444 mustReverse = false;
1445 }
1446 }
1447
1448 // independent of direction of the quad, return the same triangle(s). If a quad is organized
1449 // [0,1,2,3] or [1,2,3,0] or whatever it should return the same two triangles so that two adjacent
1450 // cells that have opposite normals on a quad will have the same consistent face triangulation and
1451 // therefore the same polygonized border.
TriangulateQuad(vtkCell * quad,FaceVector & faces)1452 void TriangulateQuad(vtkCell* quad, FaceVector& faces)
1453 {
1454 std::vector<vtkIdType> consistentTri1(3), consistentTri2(2);
1455 int l = FindLowestIndex(4, quad->GetPointIds()->GetPointer(0));
1456 bool mustReverse(false);
1457 FindLowestNeighbor(4, quad->GetPointIds()->GetPointer(0), l, mustReverse);
1458
1459 if (mustReverse)
1460 {
1461 int m = l + 4; // add four to prevent negative remainders: ' (0-1)%4 => -1 '
1462 consistentTri1[0] = quad->GetPointIds()->GetId(l);
1463 consistentTri1[1] = quad->GetPointIds()->GetId((m - 1) % 4);
1464 consistentTri1[2] = quad->GetPointIds()->GetId((m - 2) % 4);
1465
1466 consistentTri2[0] = quad->GetPointIds()->GetId(l);
1467 consistentTri2[1] = quad->GetPointIds()->GetId((m - 2) % 4);
1468 consistentTri2[2] = quad->GetPointIds()->GetId((m - 3) % 4);
1469 }
1470 else
1471 {
1472 consistentTri1[0] = quad->GetPointIds()->GetId(l);
1473 consistentTri1[1] = quad->GetPointIds()->GetId((l + 1) % 4);
1474 consistentTri1[2] = quad->GetPointIds()->GetId((l + 2) % 4);
1475
1476 consistentTri2[0] = quad->GetPointIds()->GetId(l);
1477 consistentTri2[1] = quad->GetPointIds()->GetId((l + 2) % 4);
1478 consistentTri2[2] = quad->GetPointIds()->GetId((l + 3) % 4);
1479 }
1480
1481 faces.push_back(consistentTri1);
1482 faces.push_back(consistentTri2);
1483 }
1484
TriangulatePolygonAt(vtkCell * polygon,int offset,vtkIdList * triIds)1485 int TriangulatePolygonAt(vtkCell* polygon, int offset, vtkIdList* triIds)
1486 {
1487 triIds->Reset();
1488 int nPoints = polygon->GetNumberOfPoints();
1489
1490 for (int i = 0; i < nPoints - 2; ++i)
1491 {
1492 int idx0 = offset;
1493 int idx1 = (i + offset + 1) % nPoints;
1494 int idx2 = (i + offset + 2) % nPoints;
1495 triIds->InsertNextId(polygon->GetPointId(idx0));
1496 triIds->InsertNextId(polygon->GetPointId(idx1));
1497 triIds->InsertNextId(polygon->GetPointId(idx2));
1498 }
1499 return nPoints - 2;
1500 }
1501
CalculateAngles(const vtkIdType * tri,vtkPoints * phPoints,const vtkPointIdMap * pointIdMap,double & minAngle,double & maxAngle)1502 void CalculateAngles(const vtkIdType* tri, vtkPoints* phPoints, const vtkPointIdMap* pointIdMap,
1503 double& minAngle, double& maxAngle)
1504 {
1505 vtkIdType idx0 = tri[0];
1506 vtkIdType idx1 = tri[1];
1507 vtkIdType idx2 = tri[2];
1508
1509 idx0 = pointIdMap->find(idx0)->second;
1510 idx1 = pointIdMap->find(idx1)->second;
1511 idx2 = pointIdMap->find(idx2)->second;
1512
1513 double p[9];
1514 phPoints->GetPoint(idx0, p + 0);
1515 phPoints->GetPoint(idx1, p + 3);
1516 phPoints->GetPoint(idx2, p + 6);
1517
1518 minAngle = DBL_MAX;
1519 maxAngle = 0;
1520
1521 vtkVector3d left, right;
1522 for (int i = 0; i < 3; ++i)
1523 {
1524 int a = 3 * i;
1525 int b = 3 * ((i + 1) % 3);
1526 int c = 3 * ((i + 2) % 3);
1527
1528 double* p0 = p + a;
1529 double* p1 = p + b;
1530 double* p2 = p + c;
1531
1532 left.Set(p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2]);
1533 right.Set(p2[0] - p0[0], p2[1] - p0[1], p2[2] - p0[2]);
1534 left.Normalize();
1535 right.Normalize();
1536
1537 double dot = left.Dot(right);
1538 // rounding errors can occur in the vtkVector3d::Dot function,
1539 // clamp to [-1, 1] (i.e. the input range for the acos function)
1540 dot = std::min(1.0, dot);
1541 dot = std::max(-1.0, dot);
1542
1543 double angle = acos(dot) * 180.0 / vtkMath::Pi();
1544
1545 minAngle = std::min(angle, minAngle);
1546 maxAngle = std::max(angle, maxAngle);
1547 }
1548 }
1549
TriangulatePolygon(vtkCell * polygon,FaceVector & faces,vtkIdList * triIds,vtkPoints * phPoints,vtkPointIdMap * pointIdMap)1550 void TriangulatePolygon(vtkCell* polygon, FaceVector& faces, vtkIdList* triIds, vtkPoints* phPoints,
1551 vtkPointIdMap* pointIdMap)
1552 {
1553 // attempt a fan triangulation for each point on the polygon and choose the
1554 // fan triangulation with the lowest range in internal angles differing from 60 degrees
1555
1556 int nPoints = polygon->GetNumberOfPoints();
1557 std::vector<double> minAngles(nPoints, DBL_MAX), maxAngles(nPoints, 0);
1558
1559 for (int i = 0; i < nPoints; ++i)
1560 {
1561 int nTris = TriangulatePolygonAt(polygon, i, triIds);
1562 for (int j = 0; j < nTris; ++j)
1563 {
1564 double minAngle, maxAngle;
1565 CalculateAngles(triIds->GetPointer(3 * j), phPoints, pointIdMap, minAngle, maxAngle);
1566 minAngles[i] = std::min(minAngles[i], minAngle);
1567 maxAngles[i] = std::max(maxAngles[i], maxAngle);
1568 }
1569 }
1570
1571 double minRange(DBL_MAX);
1572 int choose(-1);
1573 for (int i = 0; i < nPoints; ++i)
1574 {
1575 double minDiff = std::abs(60.0 - minAngles[i]);
1576 double maxDiff = std::abs(maxAngles[i] - 60.0);
1577 double range = minDiff + maxDiff;
1578 if (range < minRange)
1579 {
1580 choose = i;
1581 minRange = range;
1582 }
1583 }
1584
1585 int nTris = TriangulatePolygonAt(polygon, choose, triIds);
1586 for (int i = 0; i < nTris; ++i)
1587 {
1588 Face tri;
1589 tri.push_back(triIds->GetId(3 * i + 0));
1590 tri.push_back(triIds->GetId(3 * i + 1));
1591 tri.push_back(triIds->GetId(3 * i + 2));
1592 faces.push_back(tri);
1593 }
1594 }
1595
TriangulateFace(vtkCell * face,FaceVector & faces,vtkIdList * triIds,vtkPoints * phPoints,vtkPointIdMap * pointIdMap)1596 void TriangulateFace(vtkCell* face, FaceVector& faces, vtkIdList* triIds, vtkPoints* phPoints,
1597 vtkPointIdMap* pointIdMap)
1598 {
1599 switch (face->GetCellType())
1600 {
1601 case VTK_TRIANGLE:
1602 {
1603 Face tri;
1604 for (int i = 0; i < 3; ++i)
1605 {
1606 tri.push_back(face->GetPointIds()->GetId(i));
1607 }
1608 faces.push_back(tri);
1609 break;
1610 }
1611 case VTK_QUAD:
1612 {
1613 TriangulateQuad(face, faces);
1614 break;
1615 }
1616 case VTK_POLYGON:
1617 {
1618 TriangulatePolygon(face, faces, triIds, phPoints, pointIdMap);
1619 break;
1620 }
1621 default:
1622 {
1623 vtkGenericWarningMacro(<< "Unable to triangulate face cell type " << face->GetCellType());
1624 }
1625 }
1626 }
1627
CheckNonManifoldTriangulation(EdgeFaceSetMap & edgeFaceMap)1628 bool CheckNonManifoldTriangulation(EdgeFaceSetMap& edgeFaceMap)
1629 {
1630 for (const auto& entry : edgeFaceMap)
1631 {
1632 if (entry.second.size() != 2)
1633 {
1634 return false;
1635 }
1636 }
1637 return true;
1638 }
1639
GetContourPoints(double value,vtkPolyhedron * cell,vtkPointIdMap * pointIdMap,FaceEdgesVector & faceEdgesVector,EdgeFaceSetMap & edgeFaceMap,EdgeSet & originalEdges,std::vector<std::vector<vtkIdType>> & oririginalFaceTriFaceMap,PointIndexEdgeMultiMap & contourPointEdgeMultiMap,EdgePointIndexMap & edgeContourPointMap,vtkIncrementalPointLocator * locator,vtkDataArray * pointScalars,vtkPointData * inPd,vtkPointData * outPd)1640 bool GetContourPoints(double value, vtkPolyhedron* cell,
1641 vtkPointIdMap* pointIdMap, // from global id to local cell id
1642 FaceEdgesVector& faceEdgesVector, EdgeFaceSetMap& edgeFaceMap, EdgeSet& originalEdges,
1643 std::vector<std::vector<vtkIdType>>& oririginalFaceTriFaceMap,
1644 PointIndexEdgeMultiMap& contourPointEdgeMultiMap, EdgePointIndexMap& edgeContourPointMap,
1645 vtkIncrementalPointLocator* locator, vtkDataArray* pointScalars, vtkPointData* inPd,
1646 vtkPointData* outPd)
1647 {
1648
1649 vtkIdType nFaces = cell->GetNumberOfFaces();
1650
1651 // this will contain the (possibly triangulated) faces
1652 // that will be contoured.
1653 FaceVector faces;
1654
1655 if (!CheckWatertightNonManifoldPolyhedron(cell, originalEdges))
1656 {
1657 return false;
1658 }
1659
1660 // temporaries for triangulation
1661 vtkNew<vtkIdList> triIds;
1662
1663 for (vtkIdType i = 0; i < nFaces; ++i)
1664 {
1665 vtkCell* face = cell->GetFace(i);
1666 if (!face)
1667 {
1668 return false;
1669 }
1670
1671 size_t nTris = faces.size();
1672 TriangulateFace(face, faces, triIds, cell->GetPoints(), pointIdMap);
1673 std::vector<vtkIdType> trisOfFace;
1674 for (size_t j = nTris; j < faces.size(); ++j)
1675 {
1676 trisOfFace.push_back((vtkIdType)j);
1677 }
1678 oririginalFaceTriFaceMap.push_back(trisOfFace);
1679 }
1680
1681 // because of the triangulation performed above,
1682 // the faces vector now contains only faces that give exactly 0 or 1 contour lines.
1683 // this enables the walking of edge-face-contourpoint tuples to give closed contour polygon(s)
1684
1685 // make the edge-face map and the face edges list
1686 nFaces = (vtkIdType)faces.size();
1687 for (int i = 0; i < nFaces; ++i)
1688 {
1689 Face& face = faces[i];
1690 size_t nFacePoints = face.size();
1691
1692 EdgeVector edges;
1693 for (size_t j = 0; j < nFacePoints; ++j)
1694 {
1695 // each edge is in global id space.
1696 Edge e(face[j], face[(j + 1) % nFacePoints]);
1697 edges.push_back(e);
1698
1699 auto at = edgeFaceMap.find(e);
1700 if (at == edgeFaceMap.end())
1701 {
1702 std::set<vtkIdType> facesOfEdge;
1703 facesOfEdge.insert(i); // this edge is connected to face i
1704 edgeFaceMap.insert(make_pair(e, facesOfEdge));
1705 }
1706 else
1707 {
1708 std::set<vtkIdType>& facesOfEdge = at->second;
1709 facesOfEdge.insert(i);
1710 }
1711 }
1712
1713 faceEdgesVector.push_back(edges);
1714 }
1715
1716 if (!CheckNonManifoldTriangulation(edgeFaceMap))
1717 {
1718 vtkGenericWarningMacro(<< "A cell with a non-manifold triangulation has been encountered. This "
1719 "cell cannot be contoured.");
1720 return false;
1721 }
1722
1723 vtkPoints* cellPoints = cell->GetPoints();
1724
1725 const double eps = 1e-6;
1726
1727 double p0[3], p1[3], cp[3]; // left, right and contour point
1728
1729 for (const auto& entry : edgeFaceMap)
1730 {
1731 const Edge& edge = entry.first;
1732
1733 // here we need to convert the global ids of the edge to
1734 // local ids to find the points and the point scalars.
1735 auto at0 = pointIdMap->find(edge.first);
1736 auto at1 = pointIdMap->find(edge.second);
1737 if (at0 == pointIdMap->end() || at1 == pointIdMap->end())
1738 {
1739 vtkGenericWarningMacro(<< "Could not find global id " << edge.first << " or " << edge.second);
1740 continue;
1741 }
1742
1743 vtkIdType id0 = at0->second;
1744 vtkIdType id1 = at1->second;
1745
1746 double v0 = pointScalars->GetTuple1(id0);
1747 double v1 = pointScalars->GetTuple1(id1);
1748
1749 // TODO: check if a face falls completely in the value being contoured.
1750 // then add face DIRECTLY.
1751
1752 // TODO: what when an edge is completely on a contour value?
1753
1754 // TODO: what when an existing point is on a contour value?
1755 // in that case the edge-face-edge walking is no longer consistent:
1756 // from the point you can walk to each of the faces that border the point,
1757 // which often is larger than two.
1758
1759 // FOR ALL ISSUES ABOVE, FOR NOW:
1760 // clamp the fraction to be in <eps, 1-eps> to
1761 // resolve any difficulties that arise from a contour lying within
1762 // machine tolerance on an existing mesh point, edge or face.
1763
1764 if ((v0 < value && v1 >= value) || (v1 < value && v0 >= value))
1765 {
1766 cellPoints->GetPoint(id0, p0);
1767 cellPoints->GetPoint(id1, p1);
1768
1769 // note that the predicate for the if-statement we're in prohibits v1 == v0 == value
1770 // that means that an edge that is exactly on the contour will never be in the contour.
1771 // instead, two points that lie just off two other edges branching off that edge will
1772 // form the contour instead. That also prevents division by zero because v1 != v0 always
1773 double f = (value - v0) / (v1 - v0);
1774
1775 f = std::max(0.0 + eps, f);
1776 f = std::min(1.0 - eps, f);
1777
1778 for (int i = 0; i < 3; ++i)
1779 {
1780 cp[i] = (1.0 - f) * p0[i] + f * p1[i];
1781 }
1782
1783 vtkIdType ptId(-1);
1784 locator->InsertUniquePoint(cp, ptId);
1785
1786 // after point addition, also add the interpolated point value
1787 outPd->InterpolateEdge(inPd, ptId, edge.first, edge.second, f);
1788
1789 // store result in the point->edge lookup structure
1790 contourPointEdgeMultiMap.insert(make_pair(ptId, edge));
1791 }
1792 }
1793
1794 // build the reverse lookup structure edge->point
1795 for (const auto& entry : contourPointEdgeMultiMap)
1796 {
1797 auto range = contourPointEdgeMultiMap.equal_range(entry.first);
1798 for (auto jt = range.first; jt != range.second; ++jt)
1799 {
1800 edgeContourPointMap.insert(make_pair(jt->second, entry.first));
1801 }
1802 }
1803
1804 return true;
1805 }
1806
CreateContours(EdgeFaceSetMap & edgeFaceMap,FaceEdgesVector & faceEdgesVector,EdgePointIndexMap & edgeContourPointMap,EdgeSet & originalEdges,std::function<void (vtkIdList *)> contourCallback)1807 int CreateContours(EdgeFaceSetMap& edgeFaceMap, FaceEdgesVector& faceEdgesVector,
1808 EdgePointIndexMap& edgeContourPointMap, EdgeSet& originalEdges,
1809 std::function<void(vtkIdList*)> contourCallback)
1810 {
1811 EdgeSet availableContourEdges;
1812 for (const auto& entry : edgeContourPointMap)
1813 {
1814 availableContourEdges.insert(entry.first);
1815 }
1816
1817 vtkNew<vtkIdList> poly;
1818 EdgeSet visited;
1819 while (!availableContourEdges.empty())
1820 {
1821 Edge start = *availableContourEdges.begin();
1822 Edge at(start);
1823 vtkIdType lastFace(-1);
1824
1825 do
1826 {
1827 vtkIdType cp = edgeContourPointMap.find(at)->second;
1828 if (originalEdges.find(at) != originalEdges.end())
1829 {
1830 poly->InsertNextId(cp);
1831 }
1832
1833 visited.insert(at);
1834
1835 const std::set<vtkIdType>& facesOfEdge = edgeFaceMap[at];
1836
1837 vtkIdType face(lastFace);
1838 for (const vtkIdType& faceOfEdge : facesOfEdge)
1839 {
1840 if (lastFace != faceOfEdge)
1841 {
1842 face = faceOfEdge;
1843 break;
1844 }
1845 }
1846
1847 if (face == lastFace)
1848 {
1849 vtkGenericWarningMacro(<< "Face navigation failed in polyhedral contouring");
1850 return EXIT_FAILURE;
1851 }
1852
1853 lastFace = face;
1854
1855 const EdgeVector& edgesOfFace = faceEdgesVector[face];
1856
1857 for (const auto& otherEdge : edgesOfFace)
1858 {
1859 if (equal_fn()(otherEdge, at))
1860 {
1861 continue;
1862 }
1863
1864 auto found = edgeContourPointMap.find(otherEdge);
1865 if (found != edgeContourPointMap.end())
1866 {
1867 at = otherEdge;
1868 break;
1869 }
1870 }
1871 } while (!equal_fn()(at, start));
1872
1873 if (poly->GetNumberOfIds() > 2)
1874 {
1875 // do something with the poly
1876 // contour: add directly to result;
1877 // clip: use poly to carve off unwanted part(s)
1878 contourCallback(poly);
1879 }
1880
1881 for (const Edge& it : visited)
1882 {
1883 availableContourEdges.erase(it);
1884 }
1885 poly->Reset();
1886 visited.clear();
1887 }
1888
1889 return EXIT_SUCCESS;
1890 }
1891
Contour(double value,vtkDataArray * pointScalars,vtkIncrementalPointLocator * locator,vtkCellArray * verts,vtkCellArray * lines,vtkCellArray * polys,vtkPointData * inPd,vtkPointData * outPd,vtkCellData * inCd,vtkIdType cellId,vtkCellData * outCd)1892 void vtkPolyhedron::Contour(double value, vtkDataArray* pointScalars,
1893 vtkIncrementalPointLocator* locator, vtkCellArray* verts, vtkCellArray* lines,
1894 vtkCellArray* polys, vtkPointData* inPd, vtkPointData* outPd, vtkCellData* inCd, vtkIdType cellId,
1895 vtkCellData* outCd)
1896 {
1897 EdgeFaceSetMap edgeFaceMap;
1898 FaceEdgesVector faceEdgesVector;
1899 PointIndexEdgeMultiMap contourPointEdgeMultiMap;
1900 EdgePointIndexMap edgeContourPointMap;
1901 EdgeSet originalEdges;
1902 std::vector<std::vector<vtkIdType>> oririginalFaceTriFaceMap;
1903
1904 if (!GetContourPoints(value, this, this->PointIdMap, faceEdgesVector, edgeFaceMap, originalEdges,
1905 oririginalFaceTriFaceMap, contourPointEdgeMultiMap, edgeContourPointMap, locator,
1906 pointScalars, inPd, outPd))
1907 {
1908 return;
1909 }
1910
1911 vtkIdType offset(0);
1912 if (verts)
1913 {
1914 offset += verts->GetNumberOfCells();
1915 }
1916 if (lines)
1917 {
1918 offset += lines->GetNumberOfCells();
1919 }
1920
1921 if (contourPointEdgeMultiMap.empty())
1922 {
1923 return; // no contours made
1924 }
1925
1926 // the callback lambda will add each polygon found polys cell array
1927 std::function<void(vtkIdList*)> cb = [=](vtkIdList* poly) {
1928 if (!poly)
1929 return;
1930
1931 vtkIdType newCellId =
1932 offset + polys->InsertNextCell(poly->GetNumberOfIds(), poly->GetPointer(0));
1933 outCd->CopyData(inCd, cellId, newCellId);
1934 };
1935
1936 CreateContours(edgeFaceMap, faceEdgesVector, edgeContourPointMap, originalEdges, cb);
1937 }
1938
1939 // start new clipping code
1940 // first some support functions, see below for the Clip(...) function
1941
PolygonAsEdges(std::vector<vtkIdType> & polygon,std::vector<Edge> & edges,std::unordered_map<Edge,int,hash_fn,equal_fn> & edgeCount)1942 void PolygonAsEdges(std::vector<vtkIdType>& polygon, std::vector<Edge>& edges,
1943 std::unordered_map<Edge, int, hash_fn, equal_fn>& edgeCount)
1944 {
1945 for (size_t i = 0; i < polygon.size(); ++i)
1946 {
1947 Edge e(polygon[i], polygon[(i + 1) % polygon.size()]);
1948 edges.push_back(e);
1949
1950 auto at = edgeCount.find(e);
1951 if (at == edgeCount.end())
1952 {
1953 edgeCount.insert(make_pair(e, 1));
1954 }
1955 else
1956 {
1957 int& counter = at->second;
1958 counter++;
1959 }
1960 }
1961 }
1962
FindNext(std::vector<Edge> & unordered,const Edge & last,std::vector<Edge>::iterator & next,Edge & nextEdge)1963 bool FindNext(
1964 std::vector<Edge>& unordered, const Edge& last, std::vector<Edge>::iterator& next, Edge& nextEdge)
1965 {
1966 for (auto it = unordered.begin(); it != unordered.end(); ++it)
1967 {
1968 if (last.second == it->first)
1969 {
1970 next = it;
1971 nextEdge = *it;
1972 return true;
1973 }
1974 else if (last.second == it->second)
1975 {
1976 nextEdge = Edge(it->second, it->first);
1977 next = it;
1978 return true;
1979 }
1980 }
1981
1982 return false;
1983 }
1984
OrderEdgePolygon(std::vector<Edge> & unordered,std::vector<std::vector<Edge>> & ordered)1985 bool OrderEdgePolygon(std::vector<Edge>& unordered, std::vector<std::vector<Edge>>& ordered)
1986 {
1987 if (unordered.empty())
1988 {
1989 return true;
1990 }
1991
1992 std::vector<Edge> edgePolygon;
1993
1994 // ! we are NOT taking a reference here on purpose because when
1995 // ! the vector 'unordered' has its first element removed, a reference would
1996 // ! point to the *NEW* first element of the vector, or be invalid if the
1997 // ! vector backing store is completely re-allocated.
1998 // ! So, don't do this: Edge& last = *unordered.begin();
1999
2000 Edge last = *unordered.begin();
2001 edgePolygon.push_back(last);
2002 unordered.erase(unordered.begin());
2003
2004 while (!unordered.empty())
2005 {
2006 std::vector<Edge>::iterator next;
2007 Edge nextEdge;
2008 if (!FindNext(unordered, last, next, nextEdge))
2009 {
2010 if (!unordered.empty())
2011 {
2012 last = *unordered.begin();
2013 }
2014 else
2015 {
2016 break;
2017 }
2018
2019 ordered.push_back(edgePolygon);
2020 edgePolygon.clear();
2021 continue;
2022 }
2023
2024 edgePolygon.push_back(nextEdge);
2025 last = nextEdge;
2026 unordered.erase(next);
2027 }
2028 ordered.push_back(edgePolygon);
2029 return true;
2030 }
2031
EdgesToPolygon(std::vector<Edge> & edges,std::vector<vtkIdType> & polygon)2032 void EdgesToPolygon(std::vector<Edge>& edges, std::vector<vtkIdType>& polygon)
2033 {
2034 for (auto it = edges.begin(); it != edges.end(); ++it)
2035 {
2036 polygon.push_back(it->first);
2037 }
2038 }
2039
EdgesToPolygons(std::vector<std::vector<Edge>> & edgePolygons,std::vector<std::vector<vtkIdType>> & polygons)2040 void EdgesToPolygons(
2041 std::vector<std::vector<Edge>>& edgePolygons, std::vector<std::vector<vtkIdType>>& polygons)
2042 {
2043 for (auto it = edgePolygons.begin(); it != edgePolygons.end(); ++it)
2044 {
2045 std::vector<Edge>& edgePolygon = *it;
2046 std::vector<vtkIdType> polygon;
2047 EdgesToPolygon(edgePolygon, polygon);
2048 polygons.push_back(polygon);
2049 }
2050 }
2051
PruneContourPoints(std::vector<std::vector<vtkIdType>> & merged,EdgeSet & originalEdges,PointIndexEdgeMultiMap & contourPointEdgeMultiMap)2052 void PruneContourPoints(std::vector<std::vector<vtkIdType>>& merged, EdgeSet& originalEdges,
2053 PointIndexEdgeMultiMap& contourPointEdgeMultiMap)
2054 {
2055 for (auto it = merged.begin(); it != merged.end(); ++it)
2056 {
2057 std::vector<vtkIdType>& polygon = *it;
2058 // don't use size_t because the index i will get to -1 in the loop below
2059 // and size_t is *UNSIGNED*
2060 int i = (int)polygon.size() - 1;
2061 for (; i >= 0; --i)
2062 {
2063 auto at = contourPointEdgeMultiMap.find(polygon[i]);
2064 if (at != contourPointEdgeMultiMap.end())
2065 {
2066 bool doErase(true);
2067 auto eq = contourPointEdgeMultiMap.equal_range(polygon[i]);
2068 for (auto jt = eq.first; jt != eq.second; ++jt)
2069 {
2070 const Edge& edgeOfContourPoint = jt->second;
2071 if (originalEdges.find(edgeOfContourPoint) != originalEdges.end())
2072 {
2073 doErase = false;
2074 break;
2075 }
2076 }
2077
2078 if (doErase)
2079 {
2080 // the contour point is on a non-original edge: remove it from the polygon.
2081 polygon.erase(polygon.begin() + i);
2082 }
2083 }
2084 }
2085 }
2086 }
2087
MergeTriFacePolygons(std::vector<std::vector<vtkIdType>> & toMerge,std::vector<std::vector<vtkIdType>> & merged,EdgeSet & originalEdges,PointIndexEdgeMultiMap & contourPointEdgeMultiMap)2088 void MergeTriFacePolygons(std::vector<std::vector<vtkIdType>>& toMerge,
2089 std::vector<std::vector<vtkIdType>>& merged, EdgeSet& originalEdges,
2090 PointIndexEdgeMultiMap& contourPointEdgeMultiMap)
2091 {
2092 // this is a five-step procedure:
2093
2094 // 1) convert from vector<vtkIdType> to vector<Edge>
2095 // 2) remove duplicate edges;
2096 // 3) order the remaining edges head-to-tail;
2097 // 4) convert back from vector<Edge> to vector<vtkIdType>
2098 // 5) prune contour points that are not on original edges.
2099
2100 // step 1: convert from std::vector<vtkIdType> to std::vector<Edge>
2101 std::vector<std::vector<Edge>> polygonsAsEdges;
2102 std::unordered_map<Edge, int, hash_fn, equal_fn> edgeCount;
2103 for (auto it = toMerge.begin(); it != toMerge.end(); ++it)
2104 {
2105 std::vector<Edge> edgesPolygon;
2106 PolygonAsEdges(*it, edgesPolygon, edgeCount);
2107 polygonsAsEdges.push_back(edgesPolygon);
2108 }
2109
2110 // step 2: remove duplicate edges.
2111 for (auto it = polygonsAsEdges.begin(); it != polygonsAsEdges.end(); ++it)
2112 {
2113 std::vector<Edge>& edgesPolygon = *it;
2114 // don't use size_t because the index i will get to -1 in the loop below
2115 // and size_t is *UNSIGNED* => overflow
2116 int i = (int)edgesPolygon.size() - 1;
2117 for (; i >= 0; --i)
2118 {
2119 int ec = edgeCount.find(edgesPolygon[i])->second;
2120 if (ec == 2)
2121 {
2122 edgesPolygon.erase(edgesPolygon.begin() + i);
2123 }
2124 }
2125 }
2126
2127 // step 3: throw remaining edges together
2128 std::vector<Edge> withoutDuplicates;
2129 for (auto it = polygonsAsEdges.begin(); it != polygonsAsEdges.end(); ++it)
2130 {
2131 std::vector<Edge>& edgesPolygon = *it;
2132 for (auto jt = edgesPolygon.begin(); jt != edgesPolygon.end(); ++jt)
2133 {
2134 withoutDuplicates.push_back(*jt);
2135 }
2136 }
2137
2138 // step 3: and merge them
2139 std::vector<std::vector<Edge>> result;
2140 OrderEdgePolygon(withoutDuplicates, result);
2141
2142 // step 4: convert back to std::vector<vtkIdType> polygons
2143 EdgesToPolygons(result, merged);
2144
2145 // step 5: prune contour points that are not on original edges.
2146 PruneContourPoints(merged, originalEdges, contourPointEdgeMultiMap);
2147 }
2148
MergeTriFacePolygons(vtkPolyhedron * cell,std::unordered_map<vtkIdType,std::vector<vtkIdType>> & triFacePolygonMap,std::vector<std::vector<vtkIdType>> & oririginalFaceTriFaceMap,PointIndexEdgeMultiMap & contourPointEdgeMultiMap,EdgeSet & originalEdges,std::vector<std::vector<vtkIdType>> & polygons)2149 void MergeTriFacePolygons(vtkPolyhedron* cell,
2150 std::unordered_map<vtkIdType, std::vector<vtkIdType>>& triFacePolygonMap,
2151 std::vector<std::vector<vtkIdType>>& oririginalFaceTriFaceMap,
2152 PointIndexEdgeMultiMap& contourPointEdgeMultiMap, EdgeSet& originalEdges,
2153 std::vector<std::vector<vtkIdType>>& polygons)
2154 {
2155 // for each *original* face, find the list of triangulated faces
2156 // and use these to get the list of polygons on the original face
2157 int nFaces = cell->GetNumberOfFaces();
2158 for (int i = 0; i < nFaces; ++i)
2159 {
2160 const std::vector<vtkIdType>& triFacesOfOriginalFace = oririginalFaceTriFaceMap[i];
2161
2162 std::vector<std::vector<vtkIdType>> facePolygons;
2163 for (auto it = triFacesOfOriginalFace.begin(); it != triFacesOfOriginalFace.end(); ++it)
2164 {
2165 vtkIdType triFace = *it;
2166 auto at = triFacePolygonMap.find(triFace);
2167 if (at != triFacePolygonMap.end())
2168 facePolygons.push_back(at->second);
2169 }
2170
2171 if (!facePolygons.empty())
2172 {
2173 std::vector<std::vector<vtkIdType>> mergedPolygons;
2174 MergeTriFacePolygons(facePolygons, mergedPolygons, originalEdges, contourPointEdgeMultiMap);
2175 for (auto it = mergedPolygons.begin(); it != mergedPolygons.end(); ++it)
2176 {
2177 polygons.push_back(*it);
2178 }
2179 }
2180 }
2181 }
2182
Clip(double value,vtkDataArray * pointScalars,vtkIncrementalPointLocator * locator,vtkCellArray * connectivity,vtkPointData * inPd,vtkPointData * outPd,vtkCellData * inCd,vtkIdType cellId,vtkCellData * outCd,int insideOut)2183 void vtkPolyhedron::Clip(double value, vtkDataArray* pointScalars,
2184 vtkIncrementalPointLocator* locator, vtkCellArray* connectivity, vtkPointData* inPd,
2185 vtkPointData* outPd, vtkCellData* inCd, vtkIdType cellId, vtkCellData* outCd, int insideOut)
2186 {
2187 // set the compare function
2188 std::function<bool(double, double)> c = [insideOut](double a, double b) {
2189 if (insideOut)
2190 return std::less_equal<double>()(a, b);
2191
2192 return std::greater_equal<double>()(a, b);
2193 };
2194
2195 bool all(true);
2196
2197 // check if polyhedron is all in
2198 bool intersect = IntersectWithContour(this, pointScalars, this->PointIdMap, value, c, all);
2199 if (!intersect && all)
2200 {
2201 double x[3];
2202
2203 vtkNew<vtkIdList> faceStream;
2204 int nFaces = this->GetNumberOfFaces();
2205 faceStream->InsertNextId(nFaces);
2206 for (int i = 0; i < nFaces; ++i)
2207 {
2208 vtkCell* face = this->GetFace(i);
2209 int nFacePoints = (int)face->GetNumberOfPoints();
2210 faceStream->InsertNextId(nFacePoints);
2211 for (int j = 0; j < nFacePoints; ++j)
2212 {
2213 face->GetPoints()->GetPoint(j, x);
2214
2215 vtkIdType id(-1);
2216 locator->InsertUniquePoint(x, id);
2217 faceStream->InsertNextId(id);
2218 outPd->CopyData(inPd, face->GetPointId(j), id);
2219 }
2220 }
2221 if (nFaces > 0)
2222 {
2223 vtkIdType newCellId = connectivity->InsertNextCell(faceStream);
2224 outCd->CopyData(inCd, cellId, newCellId);
2225 }
2226 return;
2227 }
2228
2229 EdgeFaceSetMap edgeFaceMap;
2230 FaceEdgesVector faceEdgesVector;
2231 PointIndexEdgeMultiMap contourPointEdgeMultiMap;
2232 EdgePointIndexMap edgeContourPointMap;
2233 EdgeSet originalEdges;
2234 std::vector<std::vector<vtkIdType>> oririginalFaceTriFaceMap;
2235
2236 if (!GetContourPoints(value, this, this->PointIdMap, faceEdgesVector, edgeFaceMap, originalEdges,
2237 oririginalFaceTriFaceMap, contourPointEdgeMultiMap, edgeContourPointMap, locator,
2238 pointScalars, inPd, outPd))
2239 {
2240 return;
2241 }
2242
2243 if (contourPointEdgeMultiMap.empty())
2244 {
2245 return;
2246 }
2247
2248 std::unordered_map<vtkIdType, std::vector<vtkIdType>> triFacePolygonMap;
2249
2250 vtkPoints* cellPoints = this->GetPoints();
2251
2252 // for all (triangulated) faces, walk the edges and insert (+) points and contour points
2253 // note: the edges are oriented head-to-tail and neighbor-to-neighbor, i.e. [0-1][1-2][2-0]
2254 for (size_t i = 0; i < faceEdgesVector.size(); ++i)
2255 {
2256 const EdgeVector& edges = faceEdgesVector[i];
2257
2258 std::vector<vtkIdType> polygon;
2259 for (auto edgeIt = edges.begin(); edgeIt != edges.end(); ++edgeIt)
2260 {
2261 const Edge& edge = *edgeIt;
2262 vtkIdType v0 = edge.first;
2263 auto localIdIt = this->PointIdMap->find(v0);
2264 if (localIdIt == this->PointIdMap->end())
2265 {
2266 vtkGenericWarningMacro(<< "Could not find global id " << v0);
2267 continue;
2268 }
2269 vtkIdType localId = localIdIt->second;
2270
2271 double val0 = pointScalars->GetTuple1(localId);
2272 if (c(val0, value))
2273 {
2274 vtkIdType id(-1);
2275 locator->InsertUniquePoint(cellPoints->GetPoint(localId), id);
2276 // we have added a point, so add point data to the output too
2277 // that has to be done in global id space
2278 outPd->CopyData(inPd, v0, id);
2279 polygon.push_back(id);
2280 }
2281
2282 // if the current edge contains a contour point, add that as well
2283 // note: due to the edge ordering this works.
2284 auto at = edgeContourPointMap.find(edge);
2285 if (at != edgeContourPointMap.end())
2286 {
2287 polygon.push_back(at->second);
2288 }
2289 }
2290
2291 // if a polygon was identified (if all face points are all + or all -, there is no polygon)
2292 if (!polygon.empty())
2293 {
2294 triFacePolygonMap.insert(make_pair(static_cast<vtkIdType>(i), polygon));
2295 }
2296 }
2297
2298 std::vector<std::vector<vtkIdType>> polygons;
2299 MergeTriFacePolygons(this, triFacePolygonMap, oririginalFaceTriFaceMap, contourPointEdgeMultiMap,
2300 originalEdges, polygons);
2301
2302 // next, get the contour polygons.
2303
2304 // inside the callback lambda function defined below, we can only use pointers to capture
2305 // variables
2306 std::vector<std::vector<vtkIdType>>* pPolygons = &polygons;
2307
2308 std::function<void(vtkIdList*)> cb = [=](vtkIdList* poly) {
2309 vtkIdType nIds = poly->GetNumberOfIds();
2310 std::vector<vtkIdType> polygon;
2311 polygon.reserve(nIds);
2312 for (int i = 0; i < nIds; ++i)
2313 {
2314 polygon.push_back(poly->GetId(i));
2315 }
2316 if (!polygon.empty())
2317 pPolygons->push_back(polygon);
2318 };
2319
2320 CreateContours(edgeFaceMap, faceEdgesVector, edgeContourPointMap, originalEdges, cb);
2321
2322 // this next bit finds closed polyhedra by looking at disjoint sets of point ids
2323 // that hold the polyhedra. Note that if two closed polyhedra share one point
2324 // that they are identified as one closed polyhedron with two closed parts.
2325 while (!polygons.empty())
2326 {
2327 // the set of point ids that form a closed polyhedron
2328 std::unordered_set<vtkIdType> polyhedralIdSet;
2329
2330 // this list holds the polygons by moving references
2331 // in the polygons list of polyhedral faces that
2332 // belong to the polyhedron being built.
2333 std::vector<std::vector<vtkIdType>> polyhedralFaceList;
2334
2335 // while one face is added, keep looping all faces that
2336 // were not yet added. The face last added can make faces that were
2337 // skipped earlier be valid candidates now. At a certain point, no
2338 // faces can be added anymore, and the polyhedron is finished.
2339 bool add = true;
2340 while (add)
2341 {
2342 add = false;
2343 auto it = polygons.begin();
2344 while (it != polygons.end())
2345 {
2346 // If there are empty polygons, we erase them
2347 while (it != polygons.end() && it->empty())
2348 {
2349 it = polygons.erase(it);
2350 }
2351 if (it == polygons.end())
2352 {
2353 // All polygons were empty
2354 break;
2355 }
2356 if (polyhedralIdSet.empty())
2357 {
2358 // Insert seed polygon in the polyhedron
2359 polyhedralIdSet.insert(it->begin(), it->end());
2360 continue;
2361 }
2362
2363 const std::vector<vtkIdType>& nextPolygon = *it;
2364 auto polygon_it = nextPolygon.begin();
2365 bool insertedNextPolygon = false;
2366 for (; polygon_it != nextPolygon.end(); ++polygon_it)
2367 {
2368 // Check if the next polygon has any common point with the seed polygon
2369 if (polyhedralIdSet.find(*polygon_it) != polyhedralIdSet.end())
2370 {
2371 polyhedralIdSet.insert(nextPolygon.begin(), nextPolygon.end());
2372 polyhedralFaceList.emplace_back(std::move(*it));
2373 it = polygons.erase(it);
2374 // We might have missed a polygon earlier because
2375 // polyhedralIdSet has new ids now
2376 // this flag allows to scan again the list polygons
2377 add = true;
2378 insertedNextPolygon = true;
2379 // We found a polygon, we can look for another one now
2380 break;
2381 }
2382 }
2383 if (it == polygons.end())
2384 {
2385 break;
2386 }
2387 if (!insertedNextPolygon)
2388 {
2389 ++it;
2390 }
2391 }
2392 }
2393 if (!polyhedralFaceList.empty())
2394 {
2395 // next, build the face stream for the polyhedron.
2396 vtkNew<vtkIdList> polyhedron;
2397 // first entry: # of faces:
2398 polyhedron->InsertNextId(static_cast<vtkIdType>(polyhedralFaceList.size()));
2399 for (const auto& polyFace : polyhedralFaceList)
2400 {
2401 // each face entry starts with # points in that face
2402 polyhedron->InsertNextId(static_cast<vtkIdType>(polyFace.size()));
2403 for (const auto& id : polyFace)
2404 {
2405 // then all global face point ids
2406 polyhedron->InsertNextId(id);
2407 }
2408 }
2409
2410 vtkIdType newCellId = connectivity->InsertNextCell(polyhedron);
2411 // we've added a cell, so add cell data too
2412 outCd->CopyData(inCd, cellId, newCellId);
2413 }
2414 }
2415 }
2416
2417 //------------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)2418 void vtkPolyhedron::PrintSelf(ostream& os, vtkIndent indent)
2419 {
2420 this->Superclass::PrintSelf(os, indent);
2421
2422 os << indent << "Triangle:\n";
2423 this->Triangle->PrintSelf(os, indent.GetNextIndent());
2424
2425 os << indent << "Polygon:\n";
2426 this->Polygon->PrintSelf(os, indent.GetNextIndent());
2427
2428 os << indent << "Tetra:\n";
2429 this->Tetra->PrintSelf(os, indent.GetNextIndent());
2430
2431 os << indent << "Faces:\n";
2432 this->GlobalFaces->PrintSelf(os, indent.GetNextIndent());
2433 }
2434