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 #include "vtkPolyhedron.h"
16
17 #include "vtkCellArray.h"
18 #include "vtkIdTypeArray.h"
19 #include "vtkDoubleArray.h"
20 #include "vtkMath.h"
21 #include "vtkObjectFactory.h"
22 #include "vtkOrderedTriangulator.h"
23 #include "vtkPointData.h"
24 #include "vtkPoints.h"
25 #include "vtkTetra.h"
26 #include "vtkTriangle.h"
27 #include "vtkQuad.h"
28 #include "vtkPolygon.h"
29 #include "vtkLine.h"
30 #include "vtkEdgeTable.h"
31 #include "vtkPolyData.h"
32 #include "vtkCellLocator.h"
33 #include "vtkGenericCell.h"
34 #include "vtkPointLocator.h"
35 #include "vtkMeanValueCoordinatesInterpolator.h"
36 #include "vtkSmartPointer.h"
37 #include "vtkMergePoints.h"
38 #include "vtkCellData.h"
39 #include "vtkDataArray.h"
40 #include "vtkType.h"
41
42 #include <map>
43 #include <vector>
44 #include <set>
45 #include <list>
46 #include <limits>
47
48 vtkStandardNewMacro(vtkPolyhedron);
49
50 // Special typedef
51 typedef std::vector<vtkIdType> vtkIdVectorType;
52 class vtkPointIdMap : public std::map<vtkIdType,vtkIdType>{};
53 class vtkIdToIdMapType : public std::map<vtkIdType, vtkIdType>{};
54 class vtkIdToIdVectorMapType : public std::map<vtkIdType, vtkIdVectorType>{};
55 typedef std::map<vtkIdType,vtkIdType*>::iterator PointIdMapIterator;
56 typedef vtkIdToIdVectorMapType::iterator vtkIdToIdVectorMapIteratorType;
57 typedef std::pair<vtkIdType, vtkIdVectorType> vtkIdToIdVectorPairType;
58 typedef std::pair<vtkIdType, vtkIdType> vtkIdToIdPairType;
59 typedef std::set<vtkIdType> vtkIdSetType;
60
61 // Special class for iterating through polyhedron faces
62 //----------------------------------------------------------------------------
63 class vtkPolyhedronFaceIterator
64 {
65 public:
66 vtkIdType CurrentPolygonSize;
67 vtkIdType *Polygon;
68 vtkIdType *Current;
69 vtkIdType NumberOfPolygons;
70 vtkIdType Id;
71
vtkPolyhedronFaceIterator(vtkIdType numFaces,vtkIdType * t)72 vtkPolyhedronFaceIterator(vtkIdType numFaces, vtkIdType *t)
73 {
74 this->CurrentPolygonSize = t[0];
75 this->Polygon = t;
76 this->Current = t+1;
77 this->NumberOfPolygons = numFaces;
78 this->Id = 0;
79 }
operator ++()80 vtkIdType* operator++()
81 {
82 this->Current += this->CurrentPolygonSize + 1;
83 this->Polygon = this->Current - 1;
84 this->Id++;
85 if (this->Id < this->NumberOfPolygons)
86 {
87 this->CurrentPolygonSize = this->Polygon[0];
88 }
89 else
90 {
91 this->CurrentPolygonSize = VTK_ID_MAX;
92 }
93 return this->Current;
94 }
95 };
96
97
98 // Special class for iterating through vertices on a polygon face
99 //----------------------------------------------------------------------------
100 class vtkPolygonVertexIterator
101 {
102 public:
103 vtkIdType *Current;
104 vtkIdType NumberOfVertices;
105 vtkIdType Id;
106 // 1 or 0 for iterating along its original direction or reverse
107 vtkIdType IterDirection;
108
vtkPolygonVertexIterator(vtkIdType numVertices,vtkIdType startVertex,vtkIdType * startVertexPointer,vtkIdType nextVertex)109 vtkPolygonVertexIterator(vtkIdType numVertices, vtkIdType startVertex,
110 vtkIdType *startVertexPointer, vtkIdType nextVertex)
111 {
112 this->Current = startVertexPointer;
113 this->NumberOfVertices = numVertices;
114 this->Id = startVertex;
115 this->IterDirection = 1;
116 vtkIdType nextId = this->Id + 1;
117 vtkIdType *next = this->Current + 1;
118 if (nextId == this->NumberOfVertices)
119 {
120 next -= this->NumberOfVertices;
121 }
122 if (*next != nextVertex)
123 {
124 this->IterDirection = 0;
125 }
126 }
127
operator ++()128 vtkIdType* operator++()
129 {
130 if (this->IterDirection)
131 {
132 this->Id++;
133 this->Current++;
134 if (this->Id == this->NumberOfVertices)
135 {
136 this->Id = 0;
137 this->Current -= this->NumberOfVertices;
138 }
139 }
140 else
141 {
142 this->Id--;
143 this->Current--;
144 if (this->Id == -1)
145 {
146 this->Id = this->NumberOfVertices - 1;
147 this->Current += this->NumberOfVertices;
148 }
149 }
150 return this->Current;
151 }
152 };
153
154 //----------------------------------------------------------------------------
155 class vtkPolyhedron::vtkInternal
156 {
157 public:
158 vtkIdTypeArray * FacesBackup;
159 vtkEdgeTable * EdgeTableBackup;
160
vtkInternal()161 vtkInternal()
162 {
163 this->FacesBackup = NULL;
164 this->EdgeTableBackup = NULL;
165 }
166
~vtkInternal()167 ~vtkInternal()
168 {
169 this->FacesBackup = NULL;
170 this->EdgeTableBackup = NULL;
171 }
172
173 //----------------------------------------------------------------------------
174 // Here we use a point merger to try to prevent the problem of duplicated
175 // points in the input.
RemoveDuplicatedPointsFromFaceArrayAndEdgeTable(vtkPoints * points,vtkIdTypeArray * & faces,vtkEdgeTable * & edgeTable,double * bounds)176 void RemoveDuplicatedPointsFromFaceArrayAndEdgeTable(vtkPoints * points,
177 vtkIdTypeArray * & faces,
178 vtkEdgeTable * & edgeTable,
179 double *bounds)
180 {
181 const double eps = 0.000001;
182 vtkSmartPointer<vtkPoints> newPoints = vtkSmartPointer<vtkPoints>::New();
183 vtkSmartPointer<vtkPointLocator> merge = vtkSmartPointer<vtkPointLocator>::New();
184 merge->SetTolerance(eps);
185 merge->InitPointInsertion(newPoints, bounds);
186 bool foundDupPoint = false;
187 vtkIdType pid = -1;
188 vtkIdToIdMapType pidMap0;
189 for (vtkIdType i = 0; i < points->GetNumberOfPoints(); i++)
190 {
191 if (!merge->InsertUniquePoint(points->GetPoint(i), pid))
192 {
193 foundDupPoint = true;
194 }
195 if (pidMap0.find(pid) == pidMap0.end())
196 {
197 pidMap0.insert(vtkIdToIdPairType(pid,i));
198 }
199 }
200
201 // update face array and edge table if necessary.
202 if (foundDupPoint)
203 {
204 vtkIdToIdMapType pidMap;
205 for (vtkIdType i = 0; i < points->GetNumberOfPoints(); i++)
206 {
207 pid = merge->IsInsertedPoint(points->GetPoint(i));
208 pidMap.insert(vtkIdToIdPairType(i, pidMap0.find(pid)->second));
209 }
210
211 this->FacesBackup = faces;
212 this->EdgeTableBackup = edgeTable;
213
214 vtkIdType nfaces = 0;
215 vtkIdType insertId = 0;
216
217 faces = vtkIdTypeArray::New();
218 faces->SetNumberOfTuples(points->GetNumberOfPoints()*10);
219 faces->InsertComponent(insertId++, 0, 0); // allocate space for nfaces
220 edgeTable = vtkEdgeTable::New();
221 edgeTable->InitEdgeInsertion(points->GetNumberOfPoints());
222
223 vtkPolyhedronFaceIterator
224 faceIter(this->FacesBackup->GetValue(0), this->FacesBackup->GetPointer(1));
225 while (faceIter.Id < faceIter.NumberOfPolygons)
226 {
227 vtkIdVectorType vVector;
228 for (vtkIdType i = 0; i < faceIter.CurrentPolygonSize; i++)
229 {
230 pid = pidMap.find(faceIter.Current[i])->second;
231 vVector.push_back(pid);
232 }
233 bool dupPointRemoved = true;
234 while (dupPointRemoved && vVector.size() > 2)
235 {
236 dupPointRemoved = false;
237 if (vVector[0] == vVector[vVector.size()-1])
238 {
239 vVector.erase(vVector.begin()+vVector.size()-1);
240 dupPointRemoved = true;
241 }
242 for (size_t i = 1; i < vVector.size(); i++)
243 {
244 if (vVector[i] == vVector[i-1])
245 {
246 vVector.erase(vVector.begin()+i);
247 dupPointRemoved = true;
248 }
249 }
250 }
251 if (vVector.size() < 3)
252 {
253 ++faceIter;
254 continue;
255 }
256
257 nfaces++;
258
259 faces->InsertComponent(insertId++, 0, vVector.size());
260 for (size_t i = 0; i < vVector.size(); i++)
261 {
262 faces->InsertComponent(insertId++, 0, vVector[i]);
263 }
264 if (edgeTable->IsEdge(vVector[0],vVector[vVector.size()-1]) == (-1))
265 {
266 edgeTable->InsertEdge(vVector[0],vVector[vVector.size()-1]);
267 }
268 for (size_t i = 1; i < vVector.size(); i++)
269 {
270 if (edgeTable->IsEdge(vVector[i],vVector[i-1]) == (-1))
271 {
272 edgeTable->InsertEdge(vVector[i],vVector[i-1]);
273 }
274 }
275
276 ++faceIter;
277 }
278
279 faces->SetComponent(0,0,nfaces);
280 }
281 else
282 {
283 this->FacesBackup = NULL;
284 this->EdgeTableBackup = NULL;
285 }
286 }
287
288 //----------------------------------------------------------------------------
289 // Here we use a point merger to try to prevent the problem of duplicated
290 // points in the input.
RestoreFaceArrayAndEdgeTable(vtkIdTypeArray * & faces,vtkEdgeTable * & edgeTable)291 void RestoreFaceArrayAndEdgeTable(vtkIdTypeArray * & faces,
292 vtkEdgeTable * & edgeTable)
293 {
294 if (this->FacesBackup)
295 {
296 faces->Delete();
297 faces = this->FacesBackup;
298 }
299 if (this->EdgeTableBackup)
300 {
301 edgeTable->Delete();
302 edgeTable = this->EdgeTableBackup;
303 }
304 }
305
306 //----------------------------------------------------------------------------
307 // insert new id element in between two existing adjacent id elements.
308 // this is a convenient function. no check whether the input elements
309 // exist in the vector. no check for element adjacency.
InsertNewIdToIdVector(vtkIdVectorType & idVector,vtkIdType id,vtkIdType id0,vtkIdType id1)310 int InsertNewIdToIdVector(vtkIdVectorType & idVector, vtkIdType id,
311 vtkIdType id0, vtkIdType id1)
312 {
313 if (idVector.size() < 2)
314 {
315 return 0;
316 }
317
318 size_t num = idVector.size();
319 if ((idVector[0] == id0 && idVector[num-1] == id1)
320 ||(idVector[0] == id1 && idVector[num-1] == id0))
321 {
322 idVector.push_back(id);
323 return 1;
324 }
325
326 vtkIdVectorType::iterator iter = idVector.begin();
327 for (; iter != idVector.end(); ++iter)
328 {
329 if (*iter == id0 || *iter == id1)
330 {
331 ++iter;
332 idVector.insert(iter, id);
333 return 1;
334 }
335 }
336
337 return 0;
338 };
339
340 // Convinient function used by clip. The id is the vector index of the positive
341 // point, id0 is the vector index of the start point, and id1 is the vector index
342 // of the end point.
343 //----------------------------------------------------------------------------
EraseSegmentFromIdVector(vtkIdVectorType & idVector,vtkIdType id,vtkIdType id0,vtkIdType id1)344 int EraseSegmentFromIdVector(vtkIdVectorType & idVector, vtkIdType id,
345 vtkIdType id0, vtkIdType id1)
346 {
347 // three possible cases
348 // first case: 0 -- id0 -- id -- id1 -- size-1
349 if (id0 < id && id < id1)
350 {
351 idVector.erase(idVector.begin() + id0 + 1, idVector.begin() + id1);
352 }
353 // second case: 0 -- id1 -- id0 -- id -- size-1
354 // third case: 0 -- id -- id1 -- id0 -- size-1
355 else if (id1 < id0 && (id0 < id || id < id1))
356 {
357 idVector.erase(idVector.begin() + id0 + 1, idVector.end());
358 idVector.erase(idVector.begin(), idVector.begin() + id1);
359 }
360 else
361 {
362 // we should never get here.
363 return 0;
364 }
365 return 1;
366 };
367
368 // convert the point ids from map.first to map.second
369 //----------------------------------------------------------------------------
ConvertPointIds(vtkIdType npts,vtkIdType * pts,vtkIdToIdMapType & map,vtkIdType reverse=0)370 int ConvertPointIds(vtkIdType npts, vtkIdType * pts,
371 vtkIdToIdMapType & map, vtkIdType reverse = 0)
372 {
373 for (vtkIdType i = 0; i < npts; i++)
374 {
375 vtkIdType id = reverse ? npts-1-i : i;
376 vtkIdToIdMapType::iterator iter = map.find(pts[id]);
377 if (iter == map.end())
378 {
379 return 0;
380 }
381 pts[id] = iter->second;
382 }
383 return 1;
384 };
385
386 //----------------------------------------------------------------------------
387 // The connected contour points are found by (1) locating the current
388 // contour point in the face loop, (2) looping through face point:
389 // meet a positive point, keep going.
390 // meet a contour point, store it and stop marching in this direction.
391 // meet a negative point, stop marching in this direction.
392 // meet the same point from both directions, stop.
393 // This loop may find zero, one or two connected contour points.
FindConnectedContourPointsOnFace(vtkIdVectorType & facePtsVector,vtkIdVectorType & faceContourPtsVec,vtkIdType currContourPoint,vtkIdVectorType & pointLabelVec,vtkIdSetType & connectedContourPtsSet,vtkIdSetType & unConnectedContourPtsSet)394 void FindConnectedContourPointsOnFace(vtkIdVectorType & facePtsVector,
395 vtkIdVectorType & faceContourPtsVec,
396 vtkIdType currContourPoint,
397 vtkIdVectorType & pointLabelVec,
398 vtkIdSetType & connectedContourPtsSet,
399 vtkIdSetType & unConnectedContourPtsSet)
400 {
401 vtkIdType numFacePoints = static_cast<vtkIdType>(facePtsVector.size());
402 if (numFacePoints < 3)
403 {
404 return;
405 }
406 if (faceContourPtsVec.size() < 2)
407 {
408 return;
409 }
410 // locate the id of the startContourPt inside the face loop
411 vtkIdType startPt = -1;
412 for (vtkIdType i = 0; i < numFacePoints; i++)
413 {
414 if (currContourPoint == facePtsVector[i])
415 {
416 startPt = i;
417 break;
418 }
419 }
420
421 if (startPt < 0 || startPt >= numFacePoints)
422 {
423 return;
424 }
425
426 vtkIdType leftEndPt = -1; // face loop index
427 vtkIdType rightEndPt = -1; // face loop index
428 vtkIdType leftEndPoint = -1; // point id
429 vtkIdType rightEndPoint = -1; // point id
430 vtkIdType leftEndPassPositivePoint = 0;
431 vtkIdType rightEndPassPositivePoint = 0;
432 // search in one direction.
433 vtkIdType endPt = startPt - 1;
434 for (; endPt != startPt; endPt--)
435 {
436 if (endPt < 0)
437 {
438 endPt = numFacePoints - 1;
439 if (endPt == startPt)
440 {
441 break;
442 }
443 }
444 if (pointLabelVec[facePtsVector[endPt]] == -1)//negative point reached. stop
445 {
446 break;
447 }
448 else if (pointLabelVec[facePtsVector[endPt]] == 0)//contour pt reached. stop
449 {
450 leftEndPt = endPt;
451 leftEndPoint = facePtsVector[endPt];
452 break;
453 }
454 else
455 {
456 leftEndPassPositivePoint = 1;
457 }
458 // positive pt reached. continue.
459 }
460
461 // check if already loop through the entire face
462 if (endPt != startPt)
463 {
464 vtkIdType prevEndPt = endPt;
465
466 // search in the other direction
467 for (endPt = startPt + 1; endPt != prevEndPt; endPt++)
468 {
469 if (endPt > numFacePoints - 1)
470 {
471 endPt = 0;
472 if (endPt == prevEndPt)
473 {
474 break;
475 }
476 if (endPt == startPt)
477 {
478 break;
479 }
480 }
481 if (pointLabelVec[facePtsVector[endPt]] == -1)//negative point reached. stop
482 {
483 break;
484 }
485 else if (pointLabelVec[facePtsVector[endPt]] == 0)//contour pt reached. stop
486 {
487 rightEndPt = endPt;
488 rightEndPoint = facePtsVector[endPt];
489 break;
490 }
491 else
492 {
493 rightEndPassPositivePoint = 1;
494 }
495 }
496 }
497
498 // need to check a special case where startPt, leftEndPoint and rightEndPoint
499 // are directly connected or connected by a series of other contour points,
500 // and startPt is at one end of the contour strip. We can check this situation
501 // using leftEndPassPositivePoint and leftEndPassPositivePoint. If both are
502 // 1, then the three points are not on a contour strip. If both are 0, then
503 // startPt is not at one end of the contour strip.
504 if (leftEndPoint >= 0 && rightEndPoint >=0 && leftEndPoint != rightEndPoint)
505 {
506 if (leftEndPassPositivePoint != rightEndPassPositivePoint)
507 {
508 bool foundNonContourPoint = false;
509 for (endPt = leftEndPt - 1; endPt != rightEndPt; endPt--)
510 {
511 if (endPt < 0)
512 {
513 endPt = numFacePoints - 1;
514 if (endPt == rightEndPt)
515 {
516 break;
517 }
518 }
519 if (pointLabelVec[facePtsVector[endPt]] != 0)
520 {
521 foundNonContourPoint = true;
522 break;
523 }
524 }
525 if (!foundNonContourPoint)// startPt on one end of the contour strip
526 {
527 if (leftEndPassPositivePoint)
528 {
529 leftEndPoint = -1;
530 }
531 else
532 {
533 rightEndPoint = -1;
534 }
535 }
536 }
537 }
538
539 if (leftEndPoint >= 0)
540 {
541 connectedContourPtsSet.insert(leftEndPoint);
542 }
543 if (rightEndPoint >= 0)
544 {
545 connectedContourPtsSet.insert(rightEndPoint);
546 }
547 for (size_t i = 0; i < faceContourPtsVec.size(); i++)
548 {
549 if (faceContourPtsVec[i] != leftEndPoint &&
550 faceContourPtsVec[i] != rightEndPoint &&
551 faceContourPtsVec[i] != currContourPoint)
552 {
553 unConnectedContourPtsSet.insert(faceContourPtsVec[i]);
554 }
555 }
556 };
557
558 //----------------------------------------------------------------------------
RemoveIdFromIdToIdVectorMap(vtkIdToIdVectorMapType & map,vtkIdType id)559 void RemoveIdFromIdToIdVectorMap(vtkIdToIdVectorMapType & map, vtkIdType id)
560 {
561 vtkIdToIdVectorMapIteratorType mit = map.begin();
562 for (; mit != map.end(); ++mit)
563 {
564 vtkIdVectorType::iterator vit = mit->second.begin();
565 for (; vit != mit->second.end(); ++vit)
566 {
567 if ((*vit) == id)
568 {
569 mit->second.erase(vit);
570 break;
571 }
572 }
573 }
574 };
575
576 //----------------------------------------------------------------------------
577 // For each contour point, extract its adjacent faces, then extract other
578 // contour points on the same face that can be connected to the current
579 // points.
580 // The connected contour points are found by (1) locating the current
581 // contour point in the face loop, (2) looping through face point:
582 // meet a positive point, keep going.
583 // meet a contour point, store it and stop marching in this direction.
584 // meet a negative point, stop marching in this direction.
585 // meet the same point from both directions, stop.
586 // This loop may find zero, one or two connected contour points.
ExtractContourConnectivities(vtkIdToIdVectorMapType & ceMap,vtkIdSetType & cpSet,vtkIdVectorType & pointLabelVector,vtkIdToIdVectorMapType & pointToFacesMap,vtkIdToIdVectorMapType & faceToPointsMap,vtkIdToIdVectorMapType & faceToContourPointsMap)587 int ExtractContourConnectivities(
588 vtkIdToIdVectorMapType & ceMap,
589 vtkIdSetType & cpSet,
590 vtkIdVectorType & pointLabelVector,
591 vtkIdToIdVectorMapType & pointToFacesMap,
592 vtkIdToIdVectorMapType & faceToPointsMap,
593 vtkIdToIdVectorMapType & faceToContourPointsMap)
594 {
595 int maxConnectivity = 0;
596 if (cpSet.empty())
597 {
598 return 0;
599 }
600
601 vtkIdSetType contourBranchesSet;
602 vtkIdSetType nonContourBranchesSet;
603 vtkIdVectorType contourBranchesVector;
604 vtkIdSetType::iterator cpSetIt;
605 vtkIdToIdVectorMapType::iterator fcpMapIt, fvMapIt, ceMapIt, ceMapIt1;
606 for (cpSetIt = cpSet.begin(); cpSetIt != cpSet.end(); /*manual increment*/)
607 {
608 contourBranchesSet.clear();
609 nonContourBranchesSet.clear();
610 contourBranchesVector.clear();
611 vtkIdType pid = *cpSetIt;
612 vtkIdVectorType fVector = pointToFacesMap.find(pid)->second;
613 for (size_t i = 0; i < fVector.size(); i++)
614 {
615 // find adjacent faces that contain contour points
616 fcpMapIt = faceToContourPointsMap.find(fVector[i]);
617 if (fcpMapIt == faceToContourPointsMap.end())
618 {
619 continue;
620 }
621 fvMapIt = faceToPointsMap.find(fVector[i]);
622 if (fvMapIt == faceToPointsMap.end())
623 {
624 cout << "Cannot find point ids of a face. We should never get "
625 "here. Contouring aborted." << endl;
626 return 0;
627 }
628
629
630 // find connected contour points and store them in the set. Notice that
631 // some weird topology will classify a point as a connected contour point
632 // in one face and a non-connected contour point in some other face. we
633 // will extract the union.
634 FindConnectedContourPointsOnFace(
635 fvMapIt->second, fcpMapIt->second, pid,
636 pointLabelVector, contourBranchesSet, nonContourBranchesSet);
637 }
638
639 if (!contourBranchesSet.empty())
640 {
641 vtkIdSetType::iterator ccpSetIt = contourBranchesSet.begin();
642 for (; ccpSetIt != contourBranchesSet.end(); ++ccpSetIt)
643 {
644 if (nonContourBranchesSet.find(*ccpSetIt) == nonContourBranchesSet.end())
645 {
646 contourBranchesVector.push_back(*ccpSetIt);
647 }
648 }
649 }
650
651 if (contourBranchesVector.size() >= 2)
652 {
653 ceMap.insert(
654 vtkIdToIdVectorPairType(pid, contourBranchesVector));
655 ++cpSetIt;
656 }
657 else // throw away point contour or edge contour.
658 {
659 if (cpSetIt != cpSet.begin())
660 {
661 vtkIdSetType::iterator tempIt = cpSetIt;
662 --cpSetIt;
663 cpSet.erase(tempIt);
664 ++cpSetIt;
665 }
666 else
667 {
668 cpSet.erase(cpSetIt);
669 cpSetIt = cpSet.begin();
670 }
671 }
672 }
673
674 // sanity check, all edges should be listed twice
675 for (ceMapIt = ceMap.begin(); ceMapIt != ceMap.end(); ++ceMapIt)
676 {
677 vtkIdVectorType edges = ceMapIt->second;
678 for (size_t i = 0; i < edges.size(); i++)
679 {
680 bool foundMatch = false;
681 ceMapIt1 = ceMap.find(edges[i]);
682 if (ceMapIt1 != ceMap.end())
683 {
684 for (size_t j = 0; j < ceMapIt1->second.size(); j++)
685 {
686 if (ceMapIt->first == ceMapIt1->second[j])
687 {
688 foundMatch = true;
689 break;
690 }
691 }
692 }
693 if (!foundMatch)
694 {
695 edges.erase(edges.begin()+i);
696 i--;
697 }
698 }
699 ceMapIt->second = edges;
700 }
701
702 // clean 0 or 1-connected contour from ceMap
703 for (ceMapIt = ceMap.begin(); ceMapIt != ceMap.end(); /*manual increment*/)
704 {
705 if (ceMapIt->second.size() >= 2)
706 {
707 ++ceMapIt;
708 continue;
709 }
710
711 cpSetIt = cpSet.find(ceMapIt->first);
712 if (cpSetIt != cpSet.end())
713 {
714 cpSet.erase(cpSetIt);
715 }
716
717 if (ceMapIt != ceMap.begin())
718 {
719 vtkIdToIdVectorMapType::iterator tempIt = ceMapIt;
720 --ceMapIt;
721 ceMap.erase(tempIt);
722 ++ceMapIt;
723 }
724 else
725 {
726 ceMap.erase(ceMapIt);
727 ceMapIt = ceMap.begin();
728 }
729 }
730
731 // set maxConnectivity.
732 for (ceMapIt = ceMap.begin(); ceMapIt != ceMap.end(); ++ceMapIt)
733 {
734 if (static_cast<int>(ceMapIt->second.size()) > maxConnectivity)
735 {
736 maxConnectivity = static_cast<int>(ceMapIt->second.size());
737 }
738 }
739
740 return maxConnectivity;
741 };
742
743 //----------------------------------------------------------------------------
744 // Use eigenvalues to determine the dimension of the input contour points.
745 // This chunk of code is mostly copied from vtkOBBTree::ComputeOBB()
746 // Function returns 0 if input is a single point, 1 if co-linear,
747 // 2 if co-planar, 3 if 3D. It also returns the center as well as the normal
748 // (the eigenvector with the smallest eigenvalue) of the input contour pointset.
CheckContourDimensions(vtkPoints * points,vtkIdType npts,const vtkIdType * ptIds,double * normal,double * center)749 static int CheckContourDimensions(vtkPoints* points, vtkIdType npts,
750 const vtkIdType * ptIds,
751 double * normal, double * center)
752 {
753 static const double eigenvalueRatioThresh = 0.001;
754
755 if (npts < 3)
756 {
757 // Defensively return zeros here for normal and center.
758 if (normal)
759 {
760 normal[0] = 0.0;
761 normal[1] = 0.0;
762 normal[2] = 0.0;
763 }
764 if (center)
765 {
766 center[0] = 0.0;
767 center[1] = 0.0;
768 center[2] = 0.0;
769 }
770 return npts - 1;
771 }
772
773 vtkIdType i, j;
774 double x[3], mean[3], xp[3], *v[3], v0[3], v1[3], v2[3];
775 double *a[3], a0[3], a1[3], a2[3], eigValue[3];
776
777 // Compute mean
778 mean[0] = mean[1] = mean[2] = 0.0;
779 for (i=0; i < npts; i++ )
780 {
781 points->GetPoint(ptIds[i], x);
782 mean[0] += x[0];
783 mean[1] += x[1];
784 mean[2] += x[2];
785 }
786 for (i=0; i < 3; i++)
787 {
788 mean[i] /= npts;
789 }
790
791 // Compute covariance matrix
792 a[0] = a0; a[1] = a1; a[2] = a2;
793 for (i=0; i < 3; i++)
794 {
795 a0[i] = a1[i] = a2[i] = 0.0;
796 }
797
798 for (j = 0; j < npts; j++ )
799 {
800 points->GetPoint(ptIds[j], x);
801 xp[0] = x[0] - mean[0]; xp[1] = x[1] - mean[1]; xp[2] = x[2] - mean[2];
802 for (i = 0; i < 3; i++)
803 {
804 a0[i] += xp[0] * xp[i];
805 a1[i] += xp[1] * xp[i];
806 a2[i] += xp[2] * xp[i];
807 }
808 }//for all points
809
810 for (i=0; i < 3; i++)
811 {
812 a0[i] /= npts;
813 a1[i] /= npts;
814 a2[i] /= npts;
815 }
816
817 // Extract axes (i.e., eigenvectors) from covariance matrix.
818 v[0] = v0; v[1] = v1; v[2] = v2;
819 vtkMath::Jacobi(a,eigValue,v);
820
821 int ret = 3;
822
823 if ((eigValue[2] / eigValue[0]) < eigenvalueRatioThresh)
824 {
825 ret--;
826 }
827 if ((eigValue[1] / eigValue[0]) < eigenvalueRatioThresh)
828 {
829 ret--;
830 }
831
832 if (normal)
833 {
834 for (i =0; i < 3; i++)
835 {
836 double norm = vtkMath::Norm(a[i], 3);
837 if (norm > 0.000001)
838 {
839 break;
840 }
841 }
842 if (i < 3)
843 {
844 normal[0] = v2[0];
845 normal[1] = v2[1];
846 normal[2] = v2[2];
847 }
848 else
849 {
850 points->GetPoint(ptIds[0], v0);
851 points->GetPoint(ptIds[1], v1);
852 v0[0] = v0[0] - mean[0];
853 v0[1] = v0[1] - mean[1];
854 v0[2] = v0[2] - mean[2];
855 v1[0] = v1[0] - mean[0];
856 v1[1] = v1[1] - mean[1];
857 v1[2] = v1[2] - mean[2];
858 vtkMath::Normalize(v0);
859 vtkMath::Normalize(v1);
860 vtkMath::Cross(v0, v1, normal);
861 vtkMath::Normalize(normal);
862 }
863 }
864 if (center)
865 {
866 center[0] = mean[0];
867 center[1] = mean[1];
868 center[2] = mean[2];
869 }
870
871 return ret;
872 };
873
874 //----------------------------------------------------------------------------
875 // For each contour point, compute the normal (pointing to the positive side),
876 // then sort the other contour points connected to it, such that the connecting
877 // edges are ordered contour-clockwise when viewed from the normal direction.
878
879 // Input ceMap shows that a contour point (map->first) is connected to a number
880 // of other contour points (map->second). It does not distinguish boundary
881 // edges from internal edges. The following function also update ceMap such that
882 // a boundary edge a-->b (assuming traversing from the counter-clockwise
883 // direction) is only stored once ({a, [b, ...]}). an internal edge a<-->b is
884 // stored twice ({a, [b, ...] and {b, [a, ...]}}.
885
886 // Current implementation of this function assumes planar contours, we only
887 // compute normal once and reuse it for all other contour points.
888 // TODO: for non-planar cut, need to compute normal for each contour point. We
889 // then project edges onto a tangent plane and sort them.
OrderMultiConnectedContourPoints(vtkIdToIdVectorMapType & cpMap,vtkIdToIdVectorMapType & cpBackupMap,vtkIdSetType & cpSet,vtkPoints * points)890 static void OrderMultiConnectedContourPoints(vtkIdToIdVectorMapType & cpMap,
891 vtkIdToIdVectorMapType & cpBackupMap,
892 vtkIdSetType & cpSet,
893 vtkPoints * points)
894 {
895 double p[3], x0[3], x1[3], e0[3], e1[3], nn[3];
896 vtkIdSetType::iterator setIt;
897 vtkIdVectorType pids;
898 for (setIt = cpSet.begin(); setIt != cpSet.end(); ++setIt)
899 {
900 pids.push_back(*setIt);
901 }
902
903 // return if the input contour points are 1D. Note: the function also
904 // computes normal n and center o.
905 double o[3], n[3];
906 if (CheckContourDimensions(
907 points, static_cast<vtkIdType>(pids.size()), &(pids[0]), n, o) < 2)
908 {
909 return;
910 }
911 vtkMath::Normalize(n);
912
913 // locate an extreme point in a direction normal to the normal. this
914 // extreme point is a convex vertex.
915 vtkIdToIdVectorMapType::iterator mapIt = cpMap.begin();
916 points->GetPoint(mapIt->first, p);
917 e0[0] = p[0] - o[0];
918 e0[1] = p[1] - o[1];
919 e0[2] = p[2] - o[2];
920 vtkMath::Normalize(e0);
921 vtkMath::Cross(e0, n, nn);
922 vtkMath::Normalize(nn);
923
924 double maxDistance = VTK_DOUBLE_MIN;
925 vtkIdType maxPid = -1;
926 for (; mapIt != cpMap.end(); ++mapIt)
927 {
928 points->GetPoint(mapIt->first, p);
929 e0[0] = p[0] - o[0];
930 e0[1] = p[1] - o[1];
931 e0[2] = p[2] - o[2];
932 double distance = vtkMath::Dot(nn, e0);
933 if (distance > maxDistance)
934 {
935 maxDistance = distance;
936 maxPid = mapIt->first;
937 }
938 }
939
940 // Order edges of the contour point contour-clockwise. Note that a boundary
941 // point has two boundary edges. We will remove the incoming boundary edge
942 // and store the outgoing boundary edge at the end (after all internal edges).
943 // incoming and outgoing boudnary edges are defined when they are traversed
944 // counter-clockwisely.
945 std::vector<double> extremePointAngles; // record the angles of extreme point
946 vtkIdVectorType edges;
947 size_t edgesSize = 0;
948 const double eps = 0.0000001;
949 for (mapIt = cpMap.begin(); mapIt != cpMap.end(); ++mapIt)
950 {
951 edges = mapIt->second;
952 edgesSize = edges.size();
953
954 // If the contour point is 2-connected we don't need to order them.
955 if (edgesSize >=3 || mapIt->first == maxPid)
956 {
957 // get the current first edge
958 points->GetPoint(mapIt->first, p);
959 points->GetPoint(edges[0], x0);
960 e0[0] = x0[0] - p[0];
961 e0[1] = x0[1] - p[1];
962 e0[2] = x0[2] - p[2];
963 vtkMath::Normalize(e0);
964 vtkMath::Cross(e0, n, x0);
965 vtkMath::Cross(n, x0, e0);
966 vtkMath::Normalize(e0);
967
968 // compute the angles from other edges to the first edge
969 std::vector<double> angles;
970 angles.push_back(0);
971 const double maxDotProduct = 0.95;
972 for (size_t i = 1; i < edgesSize; i++)
973 {
974 points->GetPoint(edges[i], x1);
975 e1[0] = x1[0] - p[0];
976 e1[1] = x1[1] - p[1];
977 e1[2] = x1[2] - p[2];
978 vtkMath::Normalize(e1);
979 vtkMath::Cross(e1, n, x1);
980 vtkMath::Cross(n, x1, e1);
981 vtkMath::Normalize(e1);
982 double dotproduct = vtkMath::Dot(e0, e1);
983 double angle = acos(dotproduct);
984 if (dotproduct < maxDotProduct && dotproduct > -maxDotProduct)
985 {
986 vtkMath::Cross(e0, e1, nn);
987 if (vtkMath::Dot(n, nn) < 0)
988 {
989 angle = 2.0*vtkMath::Pi() - angle;
990 }
991 }
992 else if (dotproduct > maxDotProduct)
993 {
994 vtkMath::Cross(e0, n, nn);
995 angle = acos(vtkMath::Dot(nn, e1)) - vtkMath::Pi()/2.0;
996 }
997 else if (dotproduct < -maxDotProduct)
998 {
999 vtkMath::Cross(n, e0, nn);
1000 angle = acos(vtkMath::Dot(nn, e1)) + vtkMath::Pi()/2.0;
1001 }
1002 if (angle < -eps)
1003 {
1004 angle += 2.0*vtkMath::Pi();
1005 }
1006 if (angle > 2.0*vtkMath::Pi()+eps)
1007 {
1008 angle -= 2.0*vtkMath::Pi();
1009 }
1010 angles.push_back(angle);
1011 }
1012
1013 // sort edges
1014 for (size_t i = 1; i < edgesSize-1; i++)
1015 {
1016 for (size_t j = i+1; j < edgesSize; j++)
1017 {
1018 if (angles[i] > angles[j])
1019 {
1020 vtkIdType temp = edges[i];
1021 edges[i] = edges[j];
1022 edges[j] = temp;
1023 double angle = angles[i];
1024 angles[i] = angles[j];
1025 angles[j] = angle;
1026 }
1027 }
1028 }
1029
1030 mapIt->second = edges;
1031
1032 if (mapIt->first == maxPid)
1033 {
1034 extremePointAngles = angles;
1035 }
1036 }
1037 }
1038
1039 // store the sorted map.
1040 cpBackupMap = cpMap;
1041
1042 // find the incoming and outgoing boundary edges of the extreme point. we use
1043 // the observation: if the outgoing boundary edge is chosen as the reference
1044 // edge. the angle between all other edges and the outgoing boundary edges
1045 // will be in [0, pi]. the incoming boundary edge will be the one that is
1046 // previous to the outgoing boundary edge.
1047 mapIt = cpMap.find(maxPid);
1048 edges = mapIt->second;
1049 edgesSize = edges.size();
1050 if (extremePointAngles.size() != edgesSize)
1051 {
1052 cout << "The size of the edge array does not match the size of the "
1053 "angle array. We should never get here." << endl;
1054 return;
1055 }
1056 vtkIdType outBoundary = -1;
1057 vtkIdType inBoundary = -1;
1058 for (size_t i = 0; i < edgesSize; i++)
1059 {
1060 double angle0 = extremePointAngles[i];
1061 size_t j = 0;
1062 for (; j < edgesSize; j++)
1063 {
1064 double angle = extremePointAngles[j] - angle0;
1065 if (angle < 0)
1066 {
1067 angle = angle + 2.0*vtkMath::Pi();
1068 }
1069 if (angle > vtkMath::Pi())
1070 {
1071 break;
1072 }
1073 }
1074 if (j == edgesSize)
1075 {
1076 outBoundary = static_cast<vtkIdType>(i);
1077 inBoundary = outBoundary - 1 < 0 ?
1078 static_cast<vtkIdType>(edgesSize) - 1 : outBoundary - 1;
1079 break;
1080 }
1081 }
1082
1083 vtkIdType prevPid = maxPid;
1084 vtkIdType currPid = edges[outBoundary];
1085
1086 // remove incoming boundary edge.
1087 edges.erase(edges.begin() + inBoundary);
1088 cpMap.find(maxPid)->second = edges;
1089
1090 // traverse the contour graph to remove all incoming boundary edges.
1091 while (currPid != maxPid)
1092 {
1093 edges = cpMap.find(currPid)->second;
1094 edgesSize = edges.size();
1095 size_t i;
1096 bool foundPrevPid = false;
1097 for (i = 0; i < edgesSize; i++)
1098 {
1099 if (edges[i] == prevPid)
1100 {
1101 inBoundary = static_cast<vtkIdType>(i);
1102 outBoundary = inBoundary + 1 >= static_cast<vtkIdType>(edgesSize) ?
1103 0 : inBoundary + 1;
1104 foundPrevPid = true;
1105 break;
1106 }
1107 }
1108 if (!foundPrevPid) // traversing failed.
1109 {
1110 return;
1111 }
1112 prevPid = currPid;
1113 currPid = edges[outBoundary];
1114 edges.erase(edges.begin() + inBoundary);
1115 cpMap.find(prevPid)->second = edges;
1116 }
1117 };
1118
1119 //-----------------------------------------------------------------------------
OrderTwoConnectedContourPoints(vtkIdToIdVectorMapType & cpMap,vtkIdToIdVectorMapType & cpBackupMap)1120 void OrderTwoConnectedContourPoints(vtkIdToIdVectorMapType & cpMap,
1121 vtkIdToIdVectorMapType & cpBackupMap)
1122 {
1123 // backup the map.
1124 cpBackupMap = cpMap;
1125
1126 // traverse edges
1127 vtkIdToIdVectorMapType::iterator mapIt = cpMap.begin();
1128 vtkIdVectorType edges = mapIt->second;
1129 vtkIdType startPid = mapIt->first;
1130
1131 // choose one as incoming edge and one as outgoing edge
1132 vtkIdType outBoundary = 0;
1133 vtkIdType inBoundary = 1;
1134
1135 // find next point
1136 vtkIdType prevPid = mapIt->first;
1137 vtkIdType currPid = edges[outBoundary];
1138
1139 // remove incoming boundary edge.
1140 edges.erase(edges.begin() + inBoundary);
1141 cpMap.find(startPid)->second = edges;
1142
1143 // traverse the edge graph to remove all incoming boundary edges.
1144 while (currPid != startPid)
1145 {
1146 mapIt = cpMap.find(currPid);
1147 if (mapIt == cpMap.end())
1148 {
1149 cout << "Find an unexpected case. The input polyhedron cell may not be a "
1150 << "water tight or the polygonal faces may not be planar. Contouring "
1151 << "will continue, but this cell may not be processed correctly." << endl;
1152 break;
1153 }
1154 edges = mapIt->second;
1155 if (edges[0] == prevPid)
1156 {
1157 inBoundary = 0;
1158 outBoundary = 1;
1159 }
1160 else
1161 {
1162 inBoundary = 1;
1163 outBoundary = 0;
1164 }
1165 prevPid = currPid;
1166 currPid = edges[outBoundary];
1167 edges.erase(edges.begin() + inBoundary);
1168 cpMap.find(prevPid)->second = edges;
1169 }
1170 };
1171
1172 //----------------------------------------------------------------------------
1173 // This function is called when InternalContour() finds an unexpected case
1174 // (typically caused by a non-watertight cell). In this case, we will ignore
1175 // the existing edges between contours. Instead, simply order them as a polygon
1176 // around the center point.
OrderDisconnectedContourPoints(vtkIdSetType & cpSet,vtkPoints * points,vtkIdVectorType & pointLabelVector,vtkIdVectorType & polygon)1177 static int OrderDisconnectedContourPoints(vtkIdSetType & cpSet,
1178 vtkPoints * points,
1179 vtkIdVectorType & pointLabelVector,
1180 vtkIdVectorType & polygon)
1181 {
1182 polygon.clear();
1183 if (cpSet.empty())
1184 {
1185 return 0;
1186 }
1187
1188 double x[3], e0[3], e[3], nn[3];
1189 vtkIdSetType::iterator setIt;
1190 for (setIt = cpSet.begin(); setIt != cpSet.end(); ++setIt)
1191 {
1192 polygon.push_back(*setIt);
1193 }
1194
1195 // return if the input contour points are 1D. Note: the function also
1196 // computes normal n and center o.
1197 double o[3], n[3];
1198 if (CheckContourDimensions(
1199 points, static_cast<vtkIdType>(polygon.size()), &(polygon[0]), n, o) < 2)
1200 {
1201 return 0;
1202 }
1203
1204 // make sure normal n points to the positive side
1205 vtkIdType numPoints = static_cast<vtkIdType>(pointLabelVector.size());
1206 for (vtkIdType i = 0; i < numPoints; i++)
1207 {
1208 if (pointLabelVector[i] == 1)
1209 {
1210 points->GetPoint(i, x);
1211 e[0] = x[0] - o[0];
1212 e[1] = x[1] - o[1];
1213 e[2] = x[2] - o[2];
1214 if (vtkMath::Dot(e, n) < 0)
1215 {
1216 n[0] = -n[0];
1217 n[1] = -n[1];
1218 n[2] = -n[2];
1219 }
1220 break;
1221 }
1222 else if (pointLabelVector[i] == -1)
1223 {
1224 points->GetPoint(i, x);
1225 e[0] = x[0] - o[0];
1226 e[1] = x[1] - o[1];
1227 e[2] = x[2] - o[2];
1228 if (vtkMath::Dot(e, n) > 0)
1229 {
1230 n[0] = -n[0];
1231 n[1] = -n[1];
1232 n[2] = -n[2];
1233 }
1234 break;
1235 }
1236 }
1237
1238 // now loop over contour points to order them.
1239 std::vector<double> angles;
1240 angles.push_back(0.0);
1241
1242 // choose to start from the first point
1243 points->GetPoint(polygon[0], x);
1244 e0[0] = x[0] - o[0];
1245 e0[1] = x[1] - o[1];
1246 e0[2] = x[2] - o[2];
1247 vtkMath::Cross(e0, n, nn);
1248 vtkMath::Cross(n, nn, e0);
1249 vtkMath::Normalize(e0);
1250
1251 // compute the angles from other edges to the first edge
1252 for (size_t i = 1; i < polygon.size(); i++)
1253 {
1254 points->GetPoint(polygon[i], x);
1255 e[0] = x[0] - o[0];
1256 e[1] = x[1] - o[1];
1257 e[2] = x[2] - o[2];
1258 vtkMath::Cross(e, n, nn);
1259 vtkMath::Cross(n, nn, e);
1260 vtkMath::Normalize(e);
1261
1262 const double maxDotProduct = 0.95;
1263 double dotproduct = vtkMath::Dot(e0, e);
1264 double angle = acos(dotproduct);
1265 if (dotproduct < maxDotProduct && dotproduct > -maxDotProduct)
1266 {
1267 vtkMath::Cross(e0, e, nn);
1268 if (vtkMath::Dot(n, nn) < 0)
1269 {
1270 angle += vtkMath::Pi();
1271 }
1272 }
1273 else if (dotproduct > maxDotProduct)
1274 {
1275 vtkMath::Cross(e0, n, nn);
1276 angle = acos(vtkMath::Dot(nn, e)) - vtkMath::Pi()/2.0;
1277 }
1278 else
1279 {
1280 vtkMath::Cross(n, e0, nn);
1281 angle = acos(vtkMath::Dot(nn, e)) + vtkMath::Pi()/2.0;
1282 }
1283 angles.push_back(angle);
1284 }
1285
1286 // sort contour points
1287 for (size_t i = 1; i < polygon.size(); i++)
1288 {
1289 for (size_t j = i+1; j < polygon.size(); j++)
1290 {
1291 if (angles[i] > angles[j])
1292 {
1293 vtkIdType temp = polygon[i];
1294 polygon[i] = polygon[j];
1295 polygon[j] = temp;
1296 double angle = angles[i];
1297 angles[i] = angles[j];
1298 angles[j] = angle;
1299 }
1300 }
1301 }
1302
1303 return 1;
1304 }
1305
1306 //----------------------------------------------------------------------------
1307 // Note: the triangulation results are inserted into the input cellArray, which
1308 // does not need to be empty.
Triangulate3DContour(vtkIdType npts,vtkIdType * pts,vtkCellArray * cellArray)1309 void Triangulate3DContour(vtkIdType npts, vtkIdType * pts,
1310 vtkCellArray *cellArray)
1311 {
1312 vtkIdType start = 0;
1313 vtkIdType end = npts-1;
1314 vtkIdType ids[3];
1315
1316 while (start < end)
1317 {
1318 ids[0] = pts[start++];
1319 ids[1] = pts[start];
1320 ids[2] = pts[end];
1321 cellArray->InsertNextCell(3, ids);
1322
1323 if (start >= end - 1)
1324 {
1325 return;
1326 }
1327
1328 ids[0] = pts[end];
1329 ids[1] = pts[start];
1330 ids[2] = pts[--end];
1331 cellArray->InsertNextCell(3, ids);
1332 }
1333 };
1334
1335 }; //end vtkInternal class
1336
1337
1338 //----------------------------------------------------------------------------
1339 // Construct the hexahedron with eight points.
vtkPolyhedron()1340 vtkPolyhedron::vtkPolyhedron()
1341 {
1342 this->Line = vtkLine::New();
1343 this->Triangle = vtkTriangle::New();
1344 this->Quad = vtkQuad::New();
1345 this->Polygon = vtkPolygon::New();
1346 this->Tetra = vtkTetra::New();
1347 this->GlobalFaces = vtkIdTypeArray::New();
1348 this->FaceLocations = vtkIdTypeArray::New();
1349 this->PointIdMap = new vtkPointIdMap;
1350
1351 this->EdgesGenerated = 0;
1352 this->EdgeTable = vtkEdgeTable::New();
1353 this->Edges = vtkIdTypeArray::New();
1354 this->Edges->SetNumberOfComponents(2);
1355
1356 this->FacesGenerated = 0;
1357 this->Faces = vtkIdTypeArray::New();
1358
1359 this->BoundsComputed = 0;
1360
1361 this->PolyDataConstructed = 0;
1362 this->PolyData = vtkPolyData::New();
1363 this->Polys = vtkCellArray::New();
1364 //this->Polys->Register(this);
1365 //this->Polys->Delete();
1366 this->PolyConnectivity = vtkIdTypeArray::New();
1367 this->LocatorConstructed = 0;
1368 this->CellLocator = vtkCellLocator::New();
1369 this->CellIds = vtkIdList::New();
1370 this->Cell = vtkGenericCell::New();
1371
1372 this->Internal = new vtkInternal();
1373 }
1374
1375 //----------------------------------------------------------------------------
~vtkPolyhedron()1376 vtkPolyhedron::~vtkPolyhedron()
1377 {
1378 this->Line->Delete();
1379 this->Triangle->Delete();
1380 this->Quad->Delete();
1381 this->Polygon->Delete();
1382 this->Tetra->Delete();
1383 this->GlobalFaces->Delete();
1384 this->FaceLocations->Delete();
1385 delete this->PointIdMap;
1386 this->EdgeTable->Delete();
1387 this->Edges->Delete();
1388 this->Faces->Delete();
1389 this->PolyData->Delete();
1390 this->Polys->Delete();
1391 this->PolyConnectivity->Delete();
1392 this->CellLocator->Delete();
1393 this->CellIds->Delete();
1394 this->Cell->Delete();
1395 delete this->Internal;
1396 }
1397
1398 //----------------------------------------------------------------------------
ComputeBounds()1399 void vtkPolyhedron::ComputeBounds()
1400 {
1401 if ( this->BoundsComputed )
1402 {
1403 return;
1404 }
1405
1406 this->Superclass::GetBounds(); //stored in this->Bounds
1407 this->BoundsComputed = 1;
1408 }
1409
1410 //----------------------------------------------------------------------------
ConstructPolyData()1411 void vtkPolyhedron::ConstructPolyData()
1412 {
1413 if (this->PolyDataConstructed)
1414 {
1415 return;
1416 }
1417
1418 // Here's a trick, we're going to use the Faces array as the connectivity
1419 // array. Note that the Faces have an added nfaces value at the beginning
1420 // of the array. Other than that,it's a vtkCellArray. So we play games
1421 // with the pointers.
1422 this->GenerateFaces();
1423
1424 if (this->Faces->GetNumberOfTuples() == 0)
1425 {
1426 return;
1427 }
1428
1429 this->PolyConnectivity->SetNumberOfTuples(this->Faces->GetMaxId()-1);
1430 this->PolyConnectivity->
1431 SetArray(this->Faces->GetPointer(1), this->Faces->GetMaxId()-1, 1);
1432 this->Polys->SetNumberOfCells(*(this->Faces->GetPointer(0)));
1433 this->Polys->
1434 SetCells(*(this->Faces->GetPointer(0)), this->PolyConnectivity);
1435
1436 // Standard setup
1437 this->PolyData->Initialize();
1438 this->PolyData->SetPoints(this->Points);
1439 this->PolyData->SetPolys(this->Polys);
1440
1441 this->PolyDataConstructed = 1;
1442 }
1443
GetPolyData()1444 vtkPolyData* vtkPolyhedron::GetPolyData()
1445 {
1446 if (!this->PolyDataConstructed)
1447 {
1448 this->ConstructPolyData();
1449 }
1450
1451 return this->PolyData;
1452 }
1453 //----------------------------------------------------------------------------
ConstructLocator()1454 void vtkPolyhedron::ConstructLocator()
1455 {
1456 if (this->LocatorConstructed)
1457 {
1458 return;
1459 }
1460
1461 this->ConstructPolyData();
1462
1463 // With the polydata set up, we can assign it to the locator
1464 this->CellLocator->Initialize();
1465 this->CellLocator->SetDataSet(this->PolyData);
1466 this->CellLocator->BuildLocator();
1467
1468 this->LocatorConstructed = 1;
1469 }
1470
1471
1472 //----------------------------------------------------------------------------
ComputeParametricCoordinate(double x[3],double pc[3])1473 void vtkPolyhedron::ComputeParametricCoordinate(double x[3], double pc[3])
1474 {
1475 this->ComputeBounds();
1476 double *bounds = this->Bounds;
1477
1478 pc[0] = (x[0] - bounds[0]) / (bounds[1] - bounds[0]);
1479 pc[1] = (x[1] - bounds[2]) / (bounds[3] - bounds[2]);
1480 pc[2] = (x[2] - bounds[4]) / (bounds[5] - bounds[4]);
1481 }
1482
1483 //----------------------------------------------------------------------------
1484 void vtkPolyhedron::
ComputePositionFromParametricCoordinate(double pc[3],double x[3])1485 ComputePositionFromParametricCoordinate(double pc[3], double x[3])
1486 {
1487 this->ComputeBounds();
1488 double *bounds = this->Bounds;
1489 x[0] = ( 1 - pc[0] )* bounds[0] + pc[0] * bounds[1];
1490 x[1] = ( 1 - pc[1] )* bounds[2] + pc[1] * bounds[3];
1491 x[2] = ( 1 - pc[2] )* bounds[4] + pc[2] * bounds[5];
1492 }
1493
1494 //----------------------------------------------------------------------------
1495 // Should be called by GetCell() prior to any other method invocation and after the
1496 // points, point ids, and faces have been loaded.
Initialize()1497 void vtkPolyhedron::Initialize()
1498 {
1499 // Clear out any remaining memory.
1500 this->PointIdMap->clear();
1501
1502 // We need to create a reverse map from the point ids to their canonical cell
1503 // ids. This is a fancy way of saying that we have to be able to rapidly go
1504 // from a PointId[i] to the location i in the cell.
1505 vtkIdType i, id, numPointIds = this->PointIds->GetNumberOfIds();
1506 for (i=0; i < numPointIds; ++i)
1507 {
1508 id = this->PointIds->GetId(i);
1509 (*this->PointIdMap)[id] = i;
1510 }
1511
1512 // Edges have to be reset
1513 this->EdgesGenerated = 0;
1514 this->EdgeTable->Reset();
1515 this->Edges->Reset();
1516 this->Faces->Reset();
1517
1518 // Polys have to be reset
1519 this->Polys->Reset();
1520 this->PolyConnectivity->Reset();
1521
1522 // Faces may need renumbering later. This means converting the face ids from
1523 // global ids to local, canonical ids.
1524 this->FacesGenerated = 0;
1525
1526 // No bounds have been computed as of yet.
1527 this->BoundsComputed = 0;
1528
1529 // No supplemental geometric stuff created
1530 this->PolyDataConstructed = 0;
1531 this->LocatorConstructed = 0;
1532 }
1533
1534 //----------------------------------------------------------------------------
GetNumberOfEdges()1535 int vtkPolyhedron::GetNumberOfEdges()
1536 {
1537 // Make sure edges have been generated.
1538 if ( ! this->EdgesGenerated )
1539 {
1540 this->GenerateEdges();
1541 }
1542
1543 return static_cast<int>(this->Edges->GetNumberOfTuples());
1544 }
1545
1546 //----------------------------------------------------------------------------
1547 // This method requires that GenerateEdges() is invoked beforehand.
GetEdge(int edgeId)1548 vtkCell *vtkPolyhedron::GetEdge(int edgeId)
1549 {
1550 // Make sure edges have been generated.
1551 if ( ! this->EdgesGenerated )
1552 {
1553 this->GenerateEdges();
1554 }
1555
1556 // Make sure requested edge is within range
1557 vtkIdType numEdges = this->Edges->GetNumberOfTuples();
1558
1559 if ( edgeId < 0 || edgeId >= numEdges )
1560 {
1561 return NULL;
1562 }
1563
1564 // Return the requested edge
1565 vtkIdType edge[2];
1566 this->Edges->GetTupleValue(edgeId,edge);
1567
1568 // Recall that edge tuples are stored in canonical numbering
1569 for (int i=0; i<2; i++)
1570 {
1571 this->Line->PointIds->SetId(i,this->PointIds->GetId(edge[i]));
1572 this->Line->Points->SetPoint(i,this->Points->GetPoint(edge[i]));
1573 }
1574
1575 return this->Line;
1576 }
1577
1578 //----------------------------------------------------------------------------
GenerateEdges()1579 int vtkPolyhedron::GenerateEdges()
1580 {
1581 if ( this->EdgesGenerated )
1582 {
1583 return this->Edges->GetNumberOfTuples();
1584 }
1585
1586 //check the number of faces and return if there aren't any
1587 if ( this->GlobalFaces->GetNumberOfTuples() == 0 ||
1588 this->GlobalFaces->GetValue(0) <= 0 )
1589 {
1590 return 0;
1591 }
1592
1593 // Loop over all faces, inserting edges into the table
1594 vtkIdType *faces = this->GlobalFaces->GetPointer(0);
1595 vtkIdType nfaces = faces[0];
1596 vtkIdType *face = faces + 1;
1597 vtkIdType fid, i, edge[2], npts;
1598
1599 this->EdgeTable->InitEdgeInsertion(this->Points->GetNumberOfPoints());
1600 for (fid=0; fid < nfaces; ++fid)
1601 {
1602 npts = face[0];
1603 for (i=1; i <= npts; ++i)
1604 {
1605 edge[0] = (*this->PointIdMap)[face[i]];
1606 edge[1] = (*this->PointIdMap)[(i != npts ? face[i+1] : face[1])];
1607 if ( this->EdgeTable->IsEdge(edge[0],edge[1]) == (-1) )
1608 {
1609 this->EdgeTable->InsertEdge(edge[0],edge[1]);
1610 this->Edges->InsertNextTupleValue(edge);
1611 }
1612 }
1613 face += face[0] + 1;
1614 } //for all faces
1615
1616 // Okay all done
1617 this->EdgesGenerated = 1;
1618 return this->Edges->GetNumberOfTuples();
1619 }
1620
1621 //----------------------------------------------------------------------------
GetNumberOfFaces()1622 int vtkPolyhedron::GetNumberOfFaces()
1623 {
1624 // Make sure faces have been generated.
1625 if ( ! this->FacesGenerated )
1626 {
1627 this->GenerateFaces();
1628 }
1629
1630 if (this->GlobalFaces->GetNumberOfTuples() == 0)
1631 {
1632 return 0;
1633 }
1634
1635 return static_cast<int>(this->GlobalFaces->GetValue(0));
1636 }
1637
1638 //----------------------------------------------------------------------------
GenerateFaces()1639 void vtkPolyhedron::GenerateFaces()
1640 {
1641 if ( this->FacesGenerated )
1642 {
1643 return;
1644 }
1645
1646 if (this->GlobalFaces->GetNumberOfTuples() == 0)
1647 {
1648 return;
1649 }
1650
1651 // Basically we just run through the faces and change the global ids to the
1652 // canonical ids using the PointIdMap.
1653 this->Faces->SetNumberOfTuples(this->GlobalFaces->GetNumberOfTuples());
1654 vtkIdType *gFaces = this->GlobalFaces->GetPointer(0);
1655 vtkIdType *faces = this->Faces->GetPointer(0);
1656 vtkIdType nfaces = gFaces[0]; faces[0] = nfaces;
1657 vtkIdType *gFace = gFaces + 1;
1658 vtkIdType *face = faces + 1;
1659 vtkIdType fid, i, id, npts;
1660
1661 for (fid=0; fid < nfaces; ++fid)
1662 {
1663 npts = gFace[0];
1664 face[0] = npts;
1665 for (i=1; i <= npts; ++i)
1666 {
1667 id = (*this->PointIdMap)[gFace[i]];
1668 face[i] = id;
1669 }
1670 gFace += gFace[0] + 1;
1671 face += face[0] + 1;
1672 } //for all faces
1673
1674
1675 // Okay we've done the deed
1676 this->FacesGenerated = 1;
1677 }
1678
1679 //----------------------------------------------------------------------------
GetFace(int faceId)1680 vtkCell *vtkPolyhedron::GetFace(int faceId)
1681 {
1682 if ( faceId < 0 || faceId >= this->GlobalFaces->GetValue(0) )
1683 {
1684 return NULL;
1685 }
1686
1687 this->GenerateFaces();
1688
1689 // Okay load up the polygon
1690 vtkIdType i, p, loc = this->FaceLocations->GetValue(faceId);
1691 vtkIdType *face = this->GlobalFaces->GetPointer(loc);
1692
1693 this->Polygon->PointIds->SetNumberOfIds(face[0]);
1694 this->Polygon->Points->SetNumberOfPoints(face[0]);
1695
1696 // grab faces in global id space
1697 for (i=0; i < face[0]; ++i)
1698 {
1699 this->Polygon->PointIds->SetId(i,face[i+1]);
1700 p = (*this->PointIdMap)[face[i+1]];
1701 this->Polygon->Points->SetPoint(i,this->Points->GetPoint(p));
1702 }
1703
1704 return this->Polygon;
1705 }
1706
1707 //----------------------------------------------------------------------------
1708 // Specify the faces for this cell.
SetFaces(vtkIdType * faces)1709 void vtkPolyhedron::SetFaces(vtkIdType *faces)
1710 {
1711 // Set up face structure
1712 this->GlobalFaces->Reset();
1713 this->FaceLocations->Reset();
1714
1715 if (!faces)
1716 {
1717 return;
1718 }
1719
1720 vtkIdType nfaces = faces[0];
1721 this->FaceLocations->SetNumberOfValues(nfaces);
1722
1723 this->GlobalFaces->InsertNextValue(nfaces);
1724 vtkIdType *face = faces + 1;
1725 vtkIdType faceLoc = 1;
1726 vtkIdType i, fid, npts;
1727
1728 for (fid=0; fid < nfaces; ++fid)
1729 {
1730 npts = face[0];
1731 this->GlobalFaces->InsertNextValue(npts);
1732 for (i=1; i<=npts; ++i)
1733 {
1734 this->GlobalFaces->InsertNextValue(face[i]);
1735 }
1736 this->FaceLocations->SetValue(fid,faceLoc);
1737
1738 faceLoc += face[0] + 1;
1739 face = faces + faceLoc;
1740 } //for all faces
1741 }
1742
1743 //----------------------------------------------------------------------------
1744 // Return the list of faces for this cell.
GetFaces()1745 vtkIdType *vtkPolyhedron::GetFaces()
1746 {
1747 if (!this->GlobalFaces->GetNumberOfTuples())
1748 {
1749 return NULL;
1750 }
1751
1752 return this->GlobalFaces->GetPointer(0);
1753 }
1754
1755 //----------------------------------------------------------------------------
IntersectWithLine(double p1[3],double p2[3],double tol,double & tMin,double xMin[3],double pc[3],int & subId)1756 int vtkPolyhedron::IntersectWithLine(double p1[3], double p2[3], double tol,
1757 double& tMin, double xMin[3],
1758 double pc[3], int& subId)
1759 {
1760 // It's easiest if this is done in canonical space
1761 this->GenerateFaces();
1762
1763 // Loop over all the faces, intersecting them in turn.
1764 vtkIdType *face = this->Faces->GetPointer(0);
1765 vtkIdType nfaces = *face++;
1766 vtkIdType npts, i, fid, numHits=0;
1767 double t=VTK_FLOAT_MAX;
1768 double x[3];
1769
1770 tMin=VTK_FLOAT_MAX;
1771 for (fid=0; fid < nfaces; ++fid)
1772 {
1773 npts = face[0];
1774 vtkIdType hit = 0;
1775 switch (npts)
1776 {
1777 case 3: //triangle
1778 for (i=0; i<3; i++)
1779 {
1780 this->Triangle->Points->SetPoint(i,this->Points->GetPoint(face[i+1]));
1781 this->Triangle->PointIds->SetId(i,face[i+1]);
1782 }
1783 hit = this->Triangle->IntersectWithLine(p1,p2,tol,t,x,pc,subId);
1784 break;
1785 case 4: //quad
1786 for (i=0; i<4; i++)
1787 {
1788 this->Quad->Points->SetPoint(i,this->Points->GetPoint(face[i+1]));
1789 this->Quad->PointIds->SetId(i,face[i+1]);
1790 }
1791 hit = this->Quad->IntersectWithLine(p1,p2,tol,t,x,pc,subId);
1792 break;
1793 default: //general polygon
1794 this->Polygon->GetPoints()->SetNumberOfPoints(npts);
1795 this->Polygon->GetPointIds()->SetNumberOfIds(npts);
1796 for (i=0; i<npts; i++)
1797 {
1798 this->Polygon->Points->SetPoint(i,this->Points->GetPoint(face[i+1]));
1799 this->Polygon->PointIds->SetId(i,face[i+1]);
1800 }
1801 hit = this->Polygon->IntersectWithLine(p1,p2,tol,t,x,pc,subId);
1802 break;
1803 }
1804
1805 // Update minimum hit
1806 if ( hit )
1807 {
1808 numHits++;
1809 if ( t < tMin )
1810 {
1811 tMin = t;
1812 xMin[0] = x[0]; xMin[1] = x[1]; xMin[2] = x[2];
1813 }
1814 }
1815
1816 face += face[0] + 1;
1817 }//for all faces
1818
1819 // Compute parametric coordinates
1820 this->ComputeParametricCoordinate(xMin,pc);
1821
1822 return numHits;
1823 }
1824
1825 #define VTK_MAX_ITER 10 //Maximum iterations for ray-firing
1826 #define VTK_VOTE_THRESHOLD 3
1827
1828 //----------------------------------------------------------------------------
1829 // Shoot random rays and count the number of intersections
IsInside(double x[3],double tolerance)1830 int vtkPolyhedron::IsInside(double x[3], double tolerance)
1831 {
1832 // do a quick bounds check
1833 this->ComputeBounds();
1834 double *bounds = this->Bounds;
1835 if ( x[0] < bounds[0] || x[0] > bounds[1] ||
1836 x[1] < bounds[2] || x[1] > bounds[3] ||
1837 x[2] < bounds[4] || x[2] > bounds[5])
1838 {
1839 return 0;
1840 }
1841
1842 // It's easiest if these computations are done in canonical space
1843 this->GenerateFaces();
1844
1845 // This algorithm is adaptive; if there are enough faces in this
1846 // polyhedron, a cell locator is built to accelerate intersections.
1847 // Otherwise brute force looping over cells is used.
1848 vtkIdType *faceArray = this->Faces->GetPointer(0);
1849 vtkIdType nfaces = *faceArray++;
1850 if ( nfaces > 25 )
1851 {
1852 this->ConstructLocator();
1853 }
1854
1855 // We need a length to normalize the computations
1856 double length = sqrt( this->Superclass::GetLength2() );
1857
1858 // Perform in/out by shooting random rays. Multiple rays are fired
1859 // to improve accuracy of the result.
1860 //
1861 // The variable iterNumber counts the number of rays fired and is
1862 // limited by the defined variable VTK_MAX_ITER.
1863 //
1864 // The variable deltaVotes keeps track of the number of votes for
1865 // "in" versus "out" of the surface. When deltaVotes > 0, more votes
1866 // have counted for "in" than "out". When deltaVotes < 0, more votes
1867 // have counted for "out" than "in". When the delta_vote exceeds or
1868 // equals the defined variable VTK_VOTE_THRESHOLD, then the
1869 // appropriate "in" or "out" status is returned.
1870 //
1871 double rayMag, ray[3], xray[3], t, pcoords[3], xint[3];
1872 int i, numInts, iterNumber, deltaVotes, subId;
1873 vtkIdType idx, numCells;
1874 double tol = tolerance * length;
1875
1876 for (deltaVotes = 0, iterNumber = 1;
1877 (iterNumber < VTK_MAX_ITER) && (abs(deltaVotes) < VTK_VOTE_THRESHOLD);
1878 iterNumber++)
1879 {
1880 // Define a random ray to fire.
1881 do
1882 {
1883 for (i=0; i<3; i++)
1884 {
1885 ray[i] = vtkMath::Random(-1.0,1.0);
1886 }
1887 rayMag = vtkMath::Norm(ray);
1888 }
1889 while (rayMag == 0.0);
1890
1891 // The ray must be appropriately sized wrt the bounding box. (It has to go
1892 // all the way through the bounding box.)
1893 for (i=0; i<3; i++)
1894 {
1895 xray[i] = x[i] + (length/rayMag)*ray[i];
1896 }
1897
1898 // Intersect the line with each of the candidate cells
1899 numInts = 0;
1900
1901 if ( this->LocatorConstructed )
1902 {
1903 // Retrieve the candidate cells from the locator
1904 this->CellLocator->FindCellsAlongLine(x,xray,tol,this->CellIds);
1905 numCells = this->CellIds->GetNumberOfIds();
1906
1907 for ( idx=0; idx < numCells; idx++ )
1908 {
1909 this->PolyData->GetCell(this->CellIds->GetId(idx), this->Cell);
1910 if ( this->Cell->IntersectWithLine(x, xray, tol, t, xint, pcoords, subId) )
1911 {
1912 numInts++;
1913 }
1914 } //for all candidate cells
1915 }
1916 else
1917 {
1918 numCells = nfaces;
1919 this->ConstructPolyData();
1920
1921 for ( idx=0; idx < numCells; idx++ )
1922 {
1923 this->PolyData->GetCell(idx, this->Cell);
1924 if ( this->Cell->IntersectWithLine(x, xray, tol, t, xint, pcoords, subId) )
1925 {
1926 numInts++;
1927 }
1928 } //for all candidate cells
1929 }
1930
1931 // Count the result
1932 if ( (numInts % 2) == 0)
1933 {
1934 --deltaVotes;
1935 }
1936 else
1937 {
1938 ++deltaVotes;
1939 }
1940 } //try another ray
1941
1942 // If the number of votes is positive, the point is inside
1943 //
1944 return ( deltaVotes < 0 ? 0 : 1 );
1945 }
1946
1947 #undef VTK_MAX_ITER
1948 #undef VTK_VOTE_THRESHOLD
1949
1950
1951 //----------------------------------------------------------------------------
CellBoundary(int vtkNotUsed (subId),double pcoords[3],vtkIdList * pts)1952 int vtkPolyhedron::CellBoundary(int vtkNotUsed(subId), double pcoords[3],
1953 vtkIdList *pts)
1954 {
1955 double x[3], n[3], o[3], v[3];
1956 double dist, minDist = VTK_DOUBLE_MAX;
1957 vtkIdType numFacePts = -1;
1958 vtkIdType * facePts = 0;
1959
1960 // compute coordinates
1961 this->ComputePositionFromParametricCoordinate(pcoords, x);
1962
1963 vtkPolyhedronFaceIterator
1964 faceIter(this->GetNumberOfFaces(), this->Faces->GetPointer(1));
1965 while (faceIter.Id < faceIter.NumberOfPolygons)
1966 {
1967 if (faceIter.CurrentPolygonSize < 3)
1968 {
1969 continue;
1970 }
1971
1972 vtkPolygon::ComputeNormal(this->Points, faceIter.CurrentPolygonSize,
1973 faceIter.Current, n);
1974 vtkMath::Normalize(n);
1975 this->Points->GetPoint(faceIter.Current[0], o);
1976 v[0] = x[0] - o[0];
1977 v[1] = x[1] - o[1];
1978 v[2] = x[2] - o[2];
1979 dist = fabs(vtkMath::Dot(v, n));
1980 if (dist < minDist)
1981 {
1982 minDist = dist;
1983 numFacePts = faceIter.CurrentPolygonSize;
1984 facePts = faceIter.Current;
1985 }
1986
1987 ++faceIter;
1988 }
1989
1990 pts->Reset();
1991 if (numFacePts > 0)
1992 {
1993 for (vtkIdType i = 0; i < numFacePts; i++)
1994 {
1995 pts->InsertNextId(this->PointIds->GetId(facePts[i]));
1996 }
1997 }
1998
1999 // determine whether point is inside of polygon
2000 if ( pcoords[0] >= 0.0 && pcoords[0] <= 1.0 &&
2001 pcoords[1] >= 0.0 && pcoords[1] <= 1.0 &&
2002 pcoords[2] >= 0.0 && pcoords[2] <= 1.0 &&
2003 (this->IsInside(x, std::numeric_limits<double>::infinity())) )
2004 {
2005 return 1;
2006 }
2007 else
2008 {
2009 return 0;
2010 }
2011 }
2012
2013 //----------------------------------------------------------------------------
EvaluatePosition(double x[3],double * closestPoint,int & vtkNotUsed (subId),double pcoords[3],double & minDist2,double * weights)2014 int vtkPolyhedron::EvaluatePosition( double x[3], double * closestPoint,
2015 int & vtkNotUsed(subId), double pcoords[3],
2016 double & minDist2, double * weights )
2017 {
2018 // compute parametric coordinates
2019 this->ComputeParametricCoordinate(x, pcoords);
2020
2021 // construct polydata, the result is stored in this->PolyData,
2022 // the cell array is stored in this->Polys
2023 this->ConstructPolyData();
2024
2025 // Construct cell locator
2026 this->ConstructLocator();
2027
2028 // find closest point and store the squared distance
2029 vtkIdType cellId;
2030 int id;
2031 double cp[3];
2032 this->Cell->Initialize();
2033 this->CellLocator->FindClosestPoint(
2034 x, cp, this->Cell, cellId, id, minDist2 );
2035
2036 if (closestPoint)
2037 {
2038 closestPoint[0] = cp[0];
2039 closestPoint[1] = cp[1];
2040 closestPoint[2] = cp[2];
2041 }
2042
2043 // get the MVC weights
2044 this->InterpolateFunctions(x, weights);
2045
2046 // set distance to be zero, if point is inside
2047 int isInside = this->IsInside(x, std::numeric_limits<double>::infinity());
2048 if (isInside)
2049 {
2050 minDist2 = 0.0;
2051 }
2052
2053 return isInside;
2054 }
2055
2056 //----------------------------------------------------------------------------
EvaluateLocation(int & vtkNotUsed (subId),double pcoords[3],double x[3],double * weights)2057 void vtkPolyhedron::EvaluateLocation( int & vtkNotUsed(subId), double pcoords[3],
2058 double x[3], double * weights )
2059 {
2060 this->ComputePositionFromParametricCoordinate(pcoords, x);
2061
2062 this->InterpolateFunctions(x, weights);
2063 }
2064
2065 //----------------------------------------------------------------------------
Derivatives(int vtkNotUsed (subId),double pcoords[3],double * values,int dim,double * derivs)2066 void vtkPolyhedron::Derivatives(int vtkNotUsed(subId), double pcoords[3],
2067 double *values, int dim, double *derivs)
2068 {
2069 int i, j, k, idx;
2070 for ( j = 0; j < dim; j++ )
2071 {
2072 for ( i = 0; i < 3; i++ )
2073 {
2074 derivs[j*dim + i] = 0.0;
2075 }
2076 }
2077
2078 static const double Sample_Offset_In_Parameter_Space = 0.01;
2079
2080 double x[4][3];
2081 double coord[3];
2082
2083 //compute positions of point and three offset sample points
2084 coord[0] = pcoords[0];
2085 coord[1] = pcoords[1];
2086 coord[2] = pcoords[2];
2087 this->ComputePositionFromParametricCoordinate(coord, x[0]);
2088
2089 coord[0] += Sample_Offset_In_Parameter_Space;
2090 this->ComputePositionFromParametricCoordinate(coord, x[1]);
2091 coord[0] = pcoords[0];
2092
2093 coord[1] += Sample_Offset_In_Parameter_Space;
2094 this->ComputePositionFromParametricCoordinate(coord, x[2]);
2095 coord[1] = pcoords[1];
2096
2097 coord[2] += Sample_Offset_In_Parameter_Space;
2098 this->ComputePositionFromParametricCoordinate(coord, x[3]);
2099 coord[2] = pcoords[2];
2100
2101 this->ConstructPolyData();
2102 int numVerts = this->PolyData->GetNumberOfPoints();
2103
2104 double *weights = new double[numVerts];
2105 double *sample = new double[dim*4];
2106 //for each sample point, sample data values
2107 for ( idx = 0, k = 0; k < 4; k++ ) //loop over three sample points
2108 {
2109 this->InterpolateFunctions(x[k],weights);
2110 for ( j = 0; j < dim; j++, idx++) //over number of derivates requested
2111 {
2112 sample[idx] = 0.0;
2113 for ( i = 0; i < numVerts; i++ )
2114 {
2115 sample[idx] += weights[i] * values[j + i*dim];
2116 }
2117 }
2118 }
2119
2120 double v1[3], v2[3], v3[3];
2121 double l1, l2, l3;
2122 //compute differences along the two axes
2123 for ( i = 0; i < 3; i++ )
2124 {
2125 v1[i] = x[1][i] - x[0][i];
2126 v2[i] = x[2][i] - x[0][i];
2127 v3[i] = x[3][i] - x[0][i];
2128 }
2129 l1 = vtkMath::Normalize(v1);
2130 l2 = vtkMath::Normalize(v2);
2131 l3 = vtkMath::Normalize(v3);
2132
2133
2134 //compute derivatives along x-y-z axes
2135 double ddx, ddy, ddz;
2136 for ( j = 0; j < dim; j++ )
2137 {
2138 ddx = (sample[ dim+j] - sample[j]) / l1;
2139 ddy = (sample[2*dim+j] - sample[j]) / l2;
2140 ddz = (sample[3*dim+j] - sample[j]) / l3;
2141
2142 //project onto global x-y-z axes
2143 derivs[3*j] = ddx*v1[0] + ddy*v2[0] + ddz*v3[0];
2144 derivs[3*j + 1] = ddx*v1[1] + ddy*v2[1] + ddz*v3[1];
2145 derivs[3*j + 2] = ddx*v1[2] + ddy*v2[2] + ddz*v3[2];
2146 }
2147
2148 delete [] weights;
2149 delete [] sample;
2150 }
2151
2152 //----------------------------------------------------------------------------
GetParametricCoords()2153 double *vtkPolyhedron::GetParametricCoords()
2154 {
2155 return NULL;
2156 }
2157
2158 //----------------------------------------------------------------------------
InterpolateFunctions(double x[3],double * sf)2159 void vtkPolyhedron::InterpolateFunctions(double x[3], double *sf)
2160 {
2161 // construct polydata, the result is stored in this->PolyData,
2162 // the cell array is stored in this->Polys
2163 this->ConstructPolyData();
2164
2165 // compute the weights
2166 if (!this->PolyData->GetPoints())
2167 {
2168 return;
2169 }
2170 vtkMeanValueCoordinatesInterpolator::ComputeInterpolationWeights(
2171 x, this->PolyData->GetPoints(), this->Polys, sf);
2172 }
2173
2174 //----------------------------------------------------------------------------
InterpolateDerivs(double x[3],double * derivs)2175 void vtkPolyhedron::InterpolateDerivs(double x[3], double *derivs)
2176 {
2177 (void)x;
2178 (void)derivs;
2179 }
2180
2181 //----------------------------------------------------------------------------
Triangulate(int vtkNotUsed (index),vtkIdList * ptIds,vtkPoints * pts)2182 int vtkPolyhedron::Triangulate(int vtkNotUsed(index), vtkIdList *ptIds,
2183 vtkPoints *pts)
2184 {
2185 ptIds->Reset();
2186 pts->Reset();
2187
2188 if (!this->GetPoints() || !this->GetNumberOfPoints())
2189 {
2190 return 0;
2191 }
2192
2193 this->ComputeBounds();
2194
2195 // use ordered triangulator to triangulate the polyhedron.
2196 vtkSmartPointer<vtkOrderedTriangulator> triangulator =
2197 vtkSmartPointer<vtkOrderedTriangulator>::New();
2198
2199 triangulator->InitTriangulation(this->Bounds, this->GetNumberOfPoints());
2200 triangulator->PreSortedOff();
2201
2202 double point[3];
2203 for (vtkIdType i = 0; i < this->GetNumberOfPoints(); i++)
2204 {
2205 this->GetPoints()->GetPoint(i, point);
2206 triangulator->InsertPoint(i, point, point, 0);
2207 }
2208 triangulator->Triangulate();
2209
2210 triangulator->AddTetras(0, ptIds, pts);
2211
2212 // convert to global Ids
2213 vtkIdType* ids = ptIds->GetPointer(0);
2214 for (vtkIdType i = 0; i < ptIds->GetNumberOfIds(); i++)
2215 {
2216 ids[i] = this->PointIds->GetId(ids[i]);
2217 }
2218
2219 return 1;
2220 }
2221
2222 //----------------------------------------------------------------------------
IntersectWithContour(double value,int insideOut,vtkDataArray * inScalars)2223 int vtkPolyhedron::IntersectWithContour(double value,
2224 int insideOut,
2225 vtkDataArray *inScalars)
2226 {
2227 const double eps = 0.000001;
2228 bool allPositive = true;
2229 bool allNegative = true;
2230 for (vtkIdType pid = 0; pid < this->Points->GetNumberOfPoints(); pid++)
2231 {
2232 double v = inScalars->GetComponent(pid,0);
2233 if (v < value + eps)
2234 {
2235 allPositive = false;
2236 }
2237 else if (v > value - eps)
2238 {
2239 allNegative = false;
2240 }
2241 }
2242
2243 if ((allPositive && insideOut) || (allNegative && !insideOut))
2244 {
2245 return 2;
2246 }
2247
2248 if (allPositive || allNegative)
2249 {
2250 return 1;
2251 }
2252
2253 return 0;
2254 }
2255
2256 //----------------------------------------------------------------------------
2257 // Internal implementation of contouring algorithm
2258 // NOTE: inScalars are in canonoical id space. while inPd are in global id space.
InternalContour(double value,int insideOut,vtkIncrementalPointLocator * locator,vtkDataArray * inScalars,vtkDataArray * outScalars,vtkPointData * inPd,vtkPointData * outPd,vtkCellArray * contourPolys,vtkIdToIdVectorMapType & faceToPointsMap,vtkIdToIdVectorMapType & pointToFacesMap,vtkIdToIdMapType & pointIdMap)2259 int vtkPolyhedron::InternalContour(double value,
2260 int insideOut,
2261 vtkIncrementalPointLocator *locator,
2262 vtkDataArray *inScalars,
2263 vtkDataArray *outScalars,
2264 vtkPointData *inPd,
2265 vtkPointData *outPd,
2266 vtkCellArray *contourPolys,
2267 vtkIdToIdVectorMapType & faceToPointsMap,
2268 vtkIdToIdVectorMapType & pointToFacesMap,
2269 vtkIdToIdMapType & pointIdMap)
2270 {
2271 const double eps = 0.000001;
2272 double x0[3], x1[3], x[3];
2273 double v0, v1, v, t;
2274
2275 vtkIdType p0, p1, pid, fid, outPid, globalP0, globalP1;
2276 void * ptr = NULL;
2277
2278 pointToFacesMap.clear();
2279 faceToPointsMap.clear();
2280 pointIdMap.clear();
2281
2282 vtkIdVectorType pointLabelVector;
2283 for (pid = 0; pid < this->Points->GetNumberOfPoints(); pid++)
2284 {
2285 v = inScalars->GetComponent(pid,0);
2286 if (v < value + eps)
2287 {
2288 if (v > value - eps)
2289 {
2290 pointLabelVector.push_back(0);
2291 }
2292 else
2293 {
2294 pointLabelVector.push_back(-1);
2295 }
2296 }
2297 else if (v > value - eps)
2298 {
2299 if (v < value + eps)
2300 {
2301 pointLabelVector.push_back(0);
2302 }
2303 else
2304 {
2305 pointLabelVector.push_back(1);
2306 }
2307 }
2308 }
2309
2310 vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
2311 points->DeepCopy(this->Points);
2312
2313 if (outScalars)
2314 {
2315 for (vtkIdType i = 0; i < inScalars->GetNumberOfTuples(); i++)
2316 {
2317 outScalars->InsertNextTuple1(inScalars->GetTuple1(i));
2318 }
2319 }
2320
2321 // construct a face to contour points map
2322 vtkIdToIdVectorMapType faceToContourPointsMap;
2323
2324 vtkIdToIdVectorMapIteratorType vfMapIt, vfMapIt0, vfMapIt1;
2325 vtkIdToIdVectorMapIteratorType fvMapIt, fcpMapIt, fcpMapItTemp;
2326
2327 // loop through all faces to construct PointToFacesMap and FaceToPointsMap
2328 vtkPolyhedronFaceIterator
2329 faceIter(this->Faces->GetValue(0), this->Faces->GetPointer(1));
2330 while (faceIter.Id < faceIter.NumberOfPolygons)
2331 {
2332 // the rest code of this function assumes that a face contains at least
2333 // three vertices. return if find a single-vertex or double-vertex face.
2334 if (faceIter.CurrentPolygonSize < 3)
2335 {
2336 vtkErrorMacro("Find a face with " << faceIter.CurrentPolygonSize <<
2337 " vertices. Contouring aborted due to this degenrate case.");
2338 return -1;
2339 }
2340
2341 fid = faceIter.Id;
2342 vtkIdVectorType vVector;
2343 for (vtkIdType i = 0; i < faceIter.CurrentPolygonSize; i++)
2344 {
2345 pid = faceIter.Current[i];
2346 vfMapIt = pointToFacesMap.find(pid);
2347 if (vfMapIt != pointToFacesMap.end())
2348 {
2349 vfMapIt->second.push_back(fid);
2350 }
2351 else
2352 {
2353 vtkIdVectorType fVector;
2354 fVector.push_back(fid);
2355 pointToFacesMap.insert(vtkIdToIdVectorPairType(pid, fVector));
2356 }
2357 vVector.push_back(pid);
2358 }
2359
2360 faceToPointsMap.insert(vtkIdToIdVectorPairType(fid, vVector));
2361 ++faceIter;
2362 }
2363
2364 // loop through all edges to find contour points and store them in the point
2365 // locator. if the contour points are new (not overlap with any of original
2366 // vertex), update PointToFacesMap, FaceToPointsMap and FaceToContourPointsMap.
2367 vtkIdSetType cpSet; // contour point set
2368 this->EdgeTable->InitTraversal();
2369 while (this->EdgeTable->GetNextEdge(p0, p1, ptr))
2370 {
2371 // If both vertices are positive or negative, we do nothing and continue;
2372 if ((pointLabelVector[p0] == 1 && pointLabelVector[p1] == 1) ||
2373 (pointLabelVector[p0] == -1 && pointLabelVector[p1] == -1))
2374 {
2375 continue;
2376 }
2377
2378 globalP0 = this->PointIds->GetId(p0);
2379 globalP1 = this->PointIds->GetId(p1);
2380
2381 v0 = inScalars->GetComponent(p0,0);
2382 v1 = inScalars->GetComponent(p1,0);
2383
2384 points->GetPoint(p0, x0);
2385 points->GetPoint(p1, x1);
2386
2387 // If one or two of the vertices are contour points, we maintain the face
2388 // to contour point map then continue
2389 if (!pointLabelVector[p0] || !pointLabelVector[p1])
2390 {
2391 vtkIdType contourVertexIds[2];
2392 contourVertexIds[0] = -1;
2393 contourVertexIds[1] = -1;
2394 if (pointLabelVector[p0] == 0)
2395 {
2396 if (cpSet.insert(p0).second) // check if the point already exist in set
2397 {
2398 if (locator->InsertUniquePoint(x0, outPid))
2399 {
2400 outPd->CopyData(inPd, globalP0, outPid);
2401 }
2402 pointIdMap.insert(vtkIdToIdPairType(p0, outPid));
2403 contourVertexIds[0] = p0;
2404 }
2405 }
2406 if (pointLabelVector[p1] == 0)
2407 {
2408 if (cpSet.insert(p1).second) // check if the point already exist in set
2409 {
2410 if (locator->InsertUniquePoint(x1, outPid))
2411 {
2412 outPd->CopyData(inPd, globalP1, outPid);
2413 }
2414 pointIdMap.insert(vtkIdToIdPairType(p1, outPid));
2415 contourVertexIds[1] = p1;
2416 }
2417 }
2418
2419 for (int i = 0; i < 2; i++)
2420 {
2421 if (contourVertexIds[i] < 0)
2422 {
2423 continue;
2424 }
2425
2426 vfMapIt = pointToFacesMap.find(contourVertexIds[i]);
2427 if (vfMapIt == pointToFacesMap.end())
2428 {
2429 vtkErrorMacro("Cannot locate adjacent faces of a vertex. We should "
2430 "never get here. Contouring continue but result maybe wrong.");
2431 continue;
2432 }
2433
2434 for (size_t k = 0; k < vfMapIt->second.size(); k++)
2435 {
2436 vtkIdType contourFaceId = vfMapIt->second[k];
2437 fcpMapIt = faceToContourPointsMap.find(contourFaceId);
2438 if (fcpMapIt != faceToContourPointsMap.end())
2439 {
2440 fcpMapIt->second.push_back(contourVertexIds[i]);
2441 }
2442 else
2443 {
2444 vtkIdVectorType contourPointVector;
2445 contourPointVector.push_back(contourVertexIds[i]);
2446 faceToContourPointsMap.insert(
2447 vtkIdToIdVectorPairType(contourFaceId, contourPointVector));
2448 }
2449 }
2450 }
2451
2452 continue;
2453 }
2454
2455 // If two edge vertices are one positive and one negative. We need to
2456 // insert new contour points on this edge.
2457
2458 t = (value - v0)/(v1 - v0);
2459 x[0] = (1 - t) * x0[0] + t * x1[0];
2460 x[1] = (1 - t) * x0[1] + t * x1[1];
2461 x[2] = (1 - t) * x0[2] + t * x1[2];
2462
2463 pid = points->InsertNextPoint(x);
2464 // update pointLabelVector: we know the pid will be the number of existing
2465 // point (original verices plus previously inserted contour points)
2466 pointLabelVector.push_back(0);
2467
2468 // update PointToFacesMap: there should be two and only two faces adjacent
2469 // to the newly inserted contour point.
2470 vfMapIt0 = pointToFacesMap.find(p0);
2471 vfMapIt1 = pointToFacesMap.find(p1);
2472 vtkIdVectorType fVector;
2473 vtkIdVectorType fVector0 = vfMapIt0->second;
2474 vtkIdVectorType fVector1 = vfMapIt1->second;
2475 for (size_t i = 0; i < fVector0.size(); i++)
2476 {
2477 for (size_t j = 0; j < fVector1.size(); j++)
2478 {
2479 if (fVector0[i] == fVector1[j])
2480 {
2481 fVector.push_back(fVector0[i]);
2482 }
2483 }
2484 }
2485 if (fVector.size() != 2)
2486 {
2487 continue;
2488 }
2489 pointToFacesMap.insert(vtkIdToIdVectorPairType(pid, fVector));
2490
2491 // update FaceToPointsMap: insert the new point to the adjacent faces,
2492 // but still need to keep the order
2493 for (int k = 0; k < 2; k++)
2494 {
2495 fvMapIt = faceToPointsMap.find(fVector[k]);
2496 this->Internal->InsertNewIdToIdVector(fvMapIt->second, pid, p0, p1);
2497 }
2498
2499 // update FaceToContourPointsMap: insert the new point to the adjacent faces
2500 for (int k = 0; k < 2; k++)
2501 {
2502 fcpMapIt = faceToContourPointsMap.find(fVector[k]);
2503 if (fcpMapIt != faceToContourPointsMap.end())
2504 {
2505 fcpMapIt->second.push_back(pid);
2506 }
2507 else
2508 {
2509 vtkIdVectorType contourPointVector;
2510 contourPointVector.push_back(pid);
2511 faceToContourPointsMap.insert(
2512 vtkIdToIdVectorPairType(fVector[k], contourPointVector));
2513 }
2514 }
2515
2516 // Maintain point data. only add to locator when it has never been added
2517 // as contour point of previous processed cells.
2518 if (locator->InsertUniquePoint(x, outPid) && outPd)
2519 {
2520 outPd->InterpolateEdge(inPd,outPid,globalP0,globalP1,t);
2521 }
2522
2523 // A point unique to merge may not be unique to locator, since it may have
2524 // been inserted to locator as contour point of previous processed cells.
2525 if (outScalars)
2526 {
2527 outScalars->InsertTuple1(pid, value);
2528 }
2529
2530 pointIdMap.insert(vtkIdToIdPairType(pid, outPid));
2531
2532 cpSet.insert(pid);
2533 }
2534
2535 // Extract valid edges between contour points. We store edge information in a
2536 // edge map ceMap. The key (first field) of ceMap is contour point Pd. The
2537 // second field of ceMap is a vector of Ids of connected contour points. This
2538 // process may remove point from cpSet if that point is only connected to one
2539 // other contour point and therefore form a edge face.
2540 vtkIdToIdVectorMapType ceMap; // edge map
2541 int maxConnectivity = this->Internal->ExtractContourConnectivities(
2542 ceMap, cpSet, pointLabelVector, pointToFacesMap,
2543 faceToPointsMap, faceToContourPointsMap);
2544
2545 // special handling of point or line cases.
2546 if (cpSet.size() < 3 || ceMap.size() < 3)
2547 {
2548 for (size_t i = 0; i < pointLabelVector.size(); i++)
2549 {
2550 if (pointLabelVector[i] == 1)
2551 {
2552 return insideOut ? 2 : 1;
2553 }
2554 if (pointLabelVector[i] == -1)
2555 {
2556 return insideOut ? 1 : 2;
2557 }
2558 }
2559 return -1;
2560 }
2561
2562 // The following process needs to know whether a contour point is boundary
2563 // point and therefore contain two boundary edges.
2564 // This information is important. As we traverse the edges to extract polygons
2565 // contour edges only need to be traversed once, while internal edges need to
2566 // be traversed twice.
2567 // Note that in the simple case, where all contour points are 2-connected and
2568 // all contour edges are boundary edges. The result contour only contains one
2569 // single polygon. Otherwise, there are both boundary points (connected to two
2570 // boundary edges and zero or one or multiple internal edges) and internal
2571 // points (only connected to internal edges). The result contour contains
2572 // multiple polygons. In the latter case, we will need to distinguish boundary
2573 // contour points and interior contour points.
2574
2575 // ceMap only shows that a contour point (map->first) is connected to a number
2576 // of other contour points (map->second). the following function computes
2577 // the normal of the contour point (map->first) and then sorts the connected
2578 // contour points (map->second) such that the connected edges are ordered
2579 // counter-clockwise. the sorted edge graph is stored in ceBackupMap.
2580 // the following function also distinguishes boundary edges from internal ones
2581 // a boundary edge a-->b (assuming traversing from the counter-clockwise
2582 // direction) is only stored once ({a, [b, ...]}). an internal edge a<-->b is
2583 // stored twice ({a, [b, ...] and {b, [a, ...]}}. this graph is stored in
2584 // the updated ceMap.
2585 vtkIdToIdVectorMapType ceBackupMap;
2586 if (maxConnectivity > 2)
2587 {
2588 this->Internal->OrderMultiConnectedContourPoints(ceMap, ceBackupMap,
2589 cpSet, points);
2590 }
2591 else
2592 {
2593 this->Internal->OrderTwoConnectedContourPoints(ceMap, ceBackupMap);
2594 }
2595
2596 // cpSet and ceMap defines the contour graph. We now need to travel through
2597 // the graph to extract non-overlapping polygons. The polygons can share
2598 // edges but none of them is a subset of another one.
2599 // Here we use the order of the edges. Specifically, when a contour point
2600 // is visited, we will choose the outgoing edge to be the edge previous to the
2601 // incoming edge in the ceBackupMap.
2602 std::vector<vtkIdVectorType> polygonVector;
2603 vtkIdToIdVectorMapType::iterator ceMapIt, ceBackupMapIt;
2604 vtkIdSetType::iterator cpSetIt = cpSet.end();
2605 vtkIdVectorType::iterator cpVectorIt;
2606
2607 // backup ceMap. During graph travasal, we will remove edges from contour point
2608 // which can mess up the ordering.
2609 vtkIdSetType cpBackupSet = cpSet;
2610 bool unexpectedCell = false;
2611 while (!cpSet.empty())
2612 {
2613 vtkIdType startPid = *(cpSet.begin());
2614
2615 // check if the point still have untravelled outgoing edges.
2616 ceMapIt = ceMap.find(startPid);
2617 if (ceMapIt == ceMap.end())
2618 {
2619 cpSet.erase(cpSetIt);
2620 continue;
2621 }
2622
2623 vtkIdType currPid = startPid;
2624 vtkIdType prevPid = -1;
2625 vtkIdType nextPid = -1;
2626
2627 // vector to record points on a contour polygon
2628 vtkIdVectorType cpLoop;
2629
2630 // continue to find the next contour point.
2631 while (!cpLoop.empty() || prevPid == -1)
2632 {
2633
2634 // when back to the start point, break the loop.
2635 if (!cpLoop.empty() && currPid == startPid)
2636 {
2637 break;
2638 }
2639
2640 cpSetIt = cpSet.find(currPid);
2641 ceMapIt = ceMap.find(currPid);
2642
2643 // we should never arrive to a deadend.
2644 if (ceMapIt == ceMap.end() || cpSetIt == cpSet.end())
2645 {
2646 unexpectedCell = true;
2647 break;
2648 }
2649
2650 // add current point to the polygon loop
2651 cpLoop.push_back(currPid);
2652
2653 // get the current available outgoing edges
2654 vtkIdVectorType edges = ceMapIt->second;
2655
2656 // choose the next point to travel. the outgoing edge is chosen to be the
2657 // one previous to the incoming edge.
2658 if (prevPid == -1)
2659 {
2660 nextPid = edges[0];
2661 }
2662 else
2663 {
2664 if (edges.size() == 1)
2665 {
2666 nextPid = edges[0];
2667 }
2668 else if (edges.size() == 2)
2669 {
2670 nextPid = edges[0] == prevPid ? edges[1] : edges[0];
2671 }
2672 else
2673 {
2674 vtkIdVectorType backupEdges = ceBackupMap.find(currPid)->second;
2675 for (size_t i = 0; i < backupEdges.size(); i++)
2676 {
2677 if (backupEdges[i] == prevPid)
2678 {
2679 if (i == 0)
2680 {
2681 nextPid = backupEdges[backupEdges.size() - 1];
2682 }
2683 else
2684 {
2685 nextPid = backupEdges[i-1];
2686 }
2687 break;
2688 }
2689 }
2690 }
2691 }
2692
2693 // remove the outgoing edge
2694 bool foundEdge = false;
2695 for (size_t i = 0; i < edges.size(); i++)
2696 {
2697 if (edges[i] == nextPid)
2698 {
2699 foundEdge = true;
2700 edges.erase(edges.begin()+i);
2701 }
2702 }
2703
2704 // the next edge shouldn't have been travelled and thus been removed
2705 if (!foundEdge)
2706 {
2707 unexpectedCell = true;
2708 break;
2709 }
2710
2711 // removing point from ceMap and cpSet if all its edges have been visited.
2712 if (edges.empty())
2713 {
2714 ceMap.erase(ceMapIt);
2715 cpSet.erase(cpSetIt);
2716 }
2717 else
2718 {
2719 ceMapIt->second = edges;
2720 }
2721
2722 // move on
2723 prevPid = currPid;
2724 currPid = nextPid;
2725 nextPid = -1;
2726
2727 }// end_inner_while_loop
2728
2729 if (unexpectedCell)
2730 {
2731 //vtkWarningMacro("Find an unexpected case. The input polyhedron cell may "
2732 //"not be a water tight cell. Or the contouring function is non-planar and "
2733 //"intersects more than two edges and/or vertices on one face of the input "
2734 //"polyhedron cell. Contouring will continue, but this cell will be not be "
2735 //"processed.");
2736
2737 polygonVector.clear();
2738 vtkIdVectorType polygon;
2739 if (this->Internal->OrderDisconnectedContourPoints(cpBackupSet,
2740 points, pointLabelVector, polygon))
2741 {
2742 polygonVector.push_back(polygon);
2743 }
2744 break;
2745 }
2746
2747 if (!cpLoop.empty())
2748 {
2749 // record polygon loop.
2750 polygonVector.push_back(cpLoop);
2751 }
2752 } // end_outer_while_loop
2753
2754 //
2755 // Finally, add contour polygons to the output
2756 for (size_t i = 0; i < polygonVector.size(); i++)
2757 {
2758 vtkIdVectorType polygon = polygonVector[i];
2759
2760 vtkIdType npts = static_cast<vtkIdType>(polygon.size());
2761 vtkIdType *pts = &(polygon[0]);
2762
2763 if (npts < 3) // skip point or line contour
2764 {
2765 continue;
2766 }
2767
2768 // check the dimensionality of the contour
2769 int ret = this->Internal->
2770 CheckContourDimensions(points, npts, pts, NULL, NULL);
2771
2772 if (ret <= 1) // skip single point or co-linear points
2773 {
2774 }
2775 else if (ret == 2) // planar polygon, add directly
2776 {
2777 contourPolys->InsertNextCell(npts, pts);
2778 }
2779 else // 3D points, need to triangulate the original polygon
2780 {
2781 this->Internal->Triangulate3DContour(npts, pts, contourPolys);
2782 }
2783 }
2784
2785 return 0;
2786 }
2787
2788 //----------------------------------------------------------------------------
Contour(double value,vtkDataArray * pointScalars,vtkIncrementalPointLocator * locator,vtkCellArray * verts,vtkCellArray * lines,vtkCellArray * polys,vtkPointData * inPd,vtkPointData * outPd,vtkCellData * inCd,vtkIdType cellId,vtkCellData * outCd)2789 void vtkPolyhedron::Contour(double value,
2790 vtkDataArray *pointScalars,
2791 vtkIncrementalPointLocator *locator,
2792 vtkCellArray *verts,
2793 vtkCellArray *lines,
2794 vtkCellArray *polys,
2795 vtkPointData *inPd, vtkPointData *outPd,
2796 vtkCellData *inCd, vtkIdType cellId,
2797 vtkCellData *outCd)
2798 {
2799 vtkIdToIdVectorMapType faceToPointsMap;
2800 vtkIdToIdVectorMapType pointToFacesMap;
2801 vtkIdToIdMapType pointIdMap; //local one, not this->PointIdMap
2802 vtkIdType offset = 0;
2803 if (verts)
2804 {
2805 offset += verts->GetNumberOfCells();
2806 }
2807 if (lines)
2808 {
2809 offset += lines->GetNumberOfCells();
2810 }
2811
2812 // initialization
2813 this->GenerateEdges();
2814 this->GenerateFaces();
2815 this->ConstructPolyData();
2816 this->ComputeBounds();
2817
2818 vtkIdVectorType pointLabelVector;
2819 if (this->IntersectWithContour(value, 0, pointScalars))
2820 {
2821 return;
2822 }
2823
2824 this->Internal->RemoveDuplicatedPointsFromFaceArrayAndEdgeTable(
2825 this->Points, this->Faces, this->EdgeTable, this->Bounds);
2826
2827 vtkSmartPointer<vtkCellArray> contourPolys =
2828 vtkSmartPointer<vtkCellArray>::New();
2829
2830 int ret = this->InternalContour(value, 0, locator, pointScalars,
2831 NULL, inPd, outPd, contourPolys,
2832 faceToPointsMap, pointToFacesMap, pointIdMap);
2833 if (ret != 0)
2834 {
2835 this->Internal->RestoreFaceArrayAndEdgeTable(this->Faces, this->EdgeTable);
2836 return;
2837 }
2838
2839 vtkIdType npts = 0;
2840 vtkIdType *pts = 0;
2841 contourPolys->InitTraversal();
2842 while (contourPolys->GetNextCell(npts, pts))
2843 {
2844 if (!this->Internal->ConvertPointIds(npts, pts, pointIdMap))
2845 {
2846 vtkErrorMacro("Cannot find the id of an output point. We should never "
2847 "get here. Contouring aborted.");
2848 this->Internal->RestoreFaceArrayAndEdgeTable(this->Faces, this->EdgeTable);
2849 return;
2850 }
2851
2852 vtkIdType newCellId = offset + polys->InsertNextCell(npts, pts);
2853 outCd->CopyData(inCd, cellId, newCellId);
2854 }
2855
2856 this->Internal->RestoreFaceArrayAndEdgeTable(this->Faces, this->EdgeTable);
2857 }
2858
2859 //----------------------------------------------------------------------------
Clip(double value,vtkDataArray * pointScalars,vtkIncrementalPointLocator * locator,vtkCellArray * connectivity,vtkPointData * inPd,vtkPointData * outPd,vtkCellData * inCd,vtkIdType cellId,vtkCellData * outCd,int insideOut)2860 void vtkPolyhedron::Clip(double value,
2861 vtkDataArray *pointScalars,
2862 vtkIncrementalPointLocator *locator,
2863 vtkCellArray *connectivity,
2864 vtkPointData *inPd, vtkPointData *outPd,
2865 vtkCellData *inCd, vtkIdType cellId,
2866 vtkCellData *outCd, int insideOut)
2867 {
2868 vtkIdToIdVectorMapType faceToPointsMap;
2869 vtkIdToIdVectorMapType pointToFacesMap;
2870 vtkIdToIdMapType pointIdMap; //local one, not this->PointIdMap
2871 vtkIdType newPid, newCellId;
2872
2873 vtkIdType npts = 0;
2874 vtkIdType *pts = 0;
2875
2876 // initialization
2877 this->GenerateEdges();
2878 this->GenerateFaces();
2879 this->ConstructPolyData();
2880 this->ComputeBounds();
2881
2882 // vector to store cell connectivity
2883 vtkIdVectorType cellVector;
2884
2885 // vector to store which side of the clip function the polyhedron vertices are
2886 vtkIdVectorType pointLabelVector;
2887
2888 // check if polyhedron is all in
2889 if (this->IntersectWithContour(value, insideOut, pointScalars) == 1)
2890 {
2891 cellVector.push_back(this->Faces->GetValue(0));
2892
2893 // loop through all faces to add them into cellVector
2894 vtkPolyhedronFaceIterator
2895 faceIter(this->Faces->GetValue(0), this->Faces->GetPointer(1));
2896 while (faceIter.Id < faceIter.NumberOfPolygons)
2897 {
2898 vtkIdVectorType pids;
2899 for (vtkIdType i = 0; i < faceIter.CurrentPolygonSize; i++)
2900 {
2901 vtkIdType pid = faceIter.Current[i];
2902 if (locator->InsertUniquePoint(this->Points->GetPoint(pid), newPid))
2903 {
2904 vtkIdType globalPid = this->PointIds->GetId(pid);
2905 outPd->CopyData(inPd, globalPid, newPid);
2906 }
2907 pids.push_back(pid);
2908 pointIdMap.insert(vtkIdToIdPairType(pid, newPid));
2909 }
2910
2911 npts = static_cast<vtkIdType>(pids.size());
2912 if (npts == 0)
2913 {
2914 ++faceIter;
2915 continue;
2916 }
2917 pts = &(pids[0]);
2918 if (!this->Internal->ConvertPointIds(npts, pts, pointIdMap))
2919 {
2920 vtkErrorMacro("Cannot find the id of an output point. We should never "
2921 "get here. Clipping aborted.");
2922 this->Internal->RestoreFaceArrayAndEdgeTable(this->Faces, this->EdgeTable);
2923 return;
2924 }
2925 cellVector.push_back(npts);
2926 cellVector.insert(cellVector.end(), pts, pts+npts);
2927
2928 ++faceIter;
2929 }
2930 if (!cellVector.empty())
2931 {
2932 newCellId = connectivity->InsertNextCell(
2933 static_cast<vtkIdType>(cellVector.size()), &(cellVector[0]));
2934 outCd->CopyData(inCd, cellId, newCellId);
2935 }
2936 return;
2937 }
2938
2939 this->Internal->RemoveDuplicatedPointsFromFaceArrayAndEdgeTable(
2940 this->Points, this->Faces, this->EdgeTable, this->Bounds);
2941
2942 vtkSmartPointer<vtkDoubleArray> contourScalars =
2943 vtkSmartPointer<vtkDoubleArray>::New();
2944 vtkSmartPointer<vtkCellArray> contourPolys =
2945 vtkSmartPointer<vtkCellArray>::New();
2946
2947 int ret = this->InternalContour(value, insideOut, locator, pointScalars,
2948 contourScalars, inPd, outPd, contourPolys,
2949 faceToPointsMap, pointToFacesMap, pointIdMap);
2950
2951 // error occurs
2952 if (ret == -1)
2953 {
2954 this->Internal->RestoreFaceArrayAndEdgeTable(this->Faces, this->EdgeTable);
2955 return;
2956 }
2957
2958 // polyhedron is all outside
2959 if (ret == 2)
2960 {
2961 this->Internal->RestoreFaceArrayAndEdgeTable(this->Faces, this->EdgeTable);
2962 return;
2963 }
2964
2965 // polyhedron is all inside
2966 // FIXME: Documentation needed:
2967 // 1. How this can ever happen given the IntersectWithContour call above?
2968 // 2. If it can happen, how+why is it different than the code above that
2969 // copies the cell to the output?
2970 if (ret == 1)
2971 {
2972 cellVector.push_back(this->Faces->GetValue(0));
2973
2974 // loop through all faces to add them into cellVector
2975 vtkPolyhedronFaceIterator
2976 faceIter(this->Faces->GetValue(0), this->Faces->GetPointer(1));
2977 while (faceIter.Id < faceIter.NumberOfPolygons)
2978 {
2979 vtkIdVectorType pids;
2980 for (vtkIdType i = 0; i < faceIter.CurrentPolygonSize; i++)
2981 {
2982 vtkIdType pid = faceIter.Current[i];
2983 if (locator->InsertUniquePoint(this->Points->GetPoint(pid), newPid))
2984 {
2985 vtkIdType globalPid = this->PointIds->GetId(pid);
2986 outPd->CopyData(inPd, globalPid, newPid);
2987 }
2988 pids.push_back(pid);
2989 pointIdMap.insert(vtkIdToIdPairType(pid, newPid));
2990 }
2991
2992 npts = static_cast<vtkIdType>(pids.size());
2993 if (npts == 0)
2994 {
2995 if (faceIter.Id < faceIter.NumberOfPolygons - 1)
2996 {
2997 ++faceIter;
2998 continue;
2999 }
3000 else
3001 {
3002 break;
3003 }
3004 }
3005 pts = &(pids[0]);
3006 if (!this->Internal->ConvertPointIds(npts, pts, pointIdMap))
3007 {
3008 vtkErrorMacro("Cannot find the id of an output point. We should never "
3009 "get here. Clipping aborted.");
3010 this->Internal->RestoreFaceArrayAndEdgeTable(this->Faces, this->EdgeTable);
3011 return;
3012 }
3013 cellVector.push_back(npts);
3014 cellVector.insert(cellVector.end(), pts, pts+npts);
3015
3016 ++faceIter;
3017 }
3018 if (!cellVector.empty())
3019 {
3020 newCellId = connectivity->InsertNextCell(
3021 static_cast<vtkIdType>(cellVector.size()), &(cellVector[0]));
3022 outCd->CopyData(inCd, cellId, newCellId);
3023 }
3024 this->Internal->RestoreFaceArrayAndEdgeTable(this->Faces, this->EdgeTable);
3025 return;
3026 }
3027
3028 // prepare visited array for all faces
3029 bool* visited = new bool [this->Faces->GetValue(0)];
3030 for (int i = 0; i < this->Faces->GetValue(0); i++)
3031 {
3032 visited[i] = false;
3033 }
3034
3035 const double eps = 0.000001;
3036
3037 // Main algorithm: go through all positive points (points on the right side
3038 // of the contour). These do not include contour points.
3039 // For each point on the right side, find all of its adjacent faces. There
3040 // maybe two types of faces, (1) faces with all positive points, or
3041 // (2) faces with positive negative and contour points. For case (1), we will
3042 // keep the original face and add it into the result polyhedron. For case (2),
3043 // we will subdivide the original face, and add the subface that includes
3044 // positive points into the result polyhedron.
3045 std::vector<vtkIdVectorType> faces;
3046 vtkIdToIdVectorMapIteratorType pfMapIt, fpMapIt;
3047 for (vtkIdType pid = 0; pid < this->Points->GetNumberOfPoints(); pid++)
3048 {
3049 // find if a point is a positive point
3050 double v = contourScalars->GetComponent(pid,0);
3051 if ( (insideOut && (v > value-eps)) || ((!insideOut) && (v < value+eps)) )
3052 {
3053 continue;
3054 }
3055
3056 // find adjacent faces of the positive point
3057 pfMapIt = pointToFacesMap.find(pid);
3058 if (pfMapIt == pointToFacesMap.end())
3059 {
3060 continue;
3061 }
3062 vtkIdVectorType fids = pfMapIt->second;
3063
3064 // for each adjacent face
3065 for (size_t i = 0; i < fids.size(); i++)
3066 {
3067 vtkIdType fid = fids[i];
3068 if (visited[fid])
3069 {
3070 continue;
3071 }
3072
3073 fpMapIt = faceToPointsMap.find(fid);
3074 if (fpMapIt == faceToPointsMap.end())
3075 {
3076 vtkErrorMacro("Cannot locate points on a face. We should "
3077 "never get here. Clipping continues but may generate wrong result.");
3078 continue;
3079 }
3080 vtkIdVectorType pids = fpMapIt->second;
3081 vtkIdType numFacePoints = static_cast<vtkIdType>(pids.size());
3082
3083 // locate the positive point inside the id vector.
3084 vtkIdType positivePt = -1;
3085 for (vtkIdType j = 0; j < numFacePoints; j++)
3086 {
3087 if (pid == pids[j])
3088 {
3089 positivePt = j;
3090 break;
3091 }
3092 }
3093
3094 // positive point not found: this can happen when the current face
3095 // has been partially visited before, and some points have been removed from
3096 // its point vector.
3097 if (positivePt < 0 || positivePt >= numFacePoints)
3098 {
3099 continue;
3100 }
3101
3102 // a new id vector to hold ids of points on new surface patch
3103 vtkIdVectorType newpids;
3104 newpids.push_back(pid);
3105
3106 // step through the ajacent points on both sides of the positive point.
3107 // stop when a contour point or a negative point is hit.
3108 bool startFound = false;
3109 bool endFound = false;
3110
3111 vtkIdType startPt = positivePt - 1;
3112 vtkIdType endPt = positivePt + 1;
3113 for (vtkIdType k = 0; k < numFacePoints; k++)
3114 {
3115 if (startFound && endFound)
3116 {
3117 break;
3118 }
3119
3120 if (!startFound)
3121 {
3122 if (startPt < 0)
3123 {
3124 startPt = numFacePoints - 1;
3125 }
3126
3127 newpids.insert(newpids.begin(), pids[startPt]);
3128 v = contourScalars->GetComponent(pids[startPt],0);
3129 if ((insideOut && (v > value-eps)) || ((!insideOut) && (v < value+eps)))
3130 {
3131 startFound = true;
3132 if ((insideOut && (v > value+eps)) || ((!insideOut) && (v < value-eps)))
3133 {
3134 vtkWarningMacro("A positive point is directly connected to a "
3135 "negative point with no contour point in between. We should "
3136 "never get here.");
3137 if (startPt == numFacePoints-1)
3138 {
3139 startPt = 0;
3140 }
3141 else
3142 {
3143 startPt++;
3144 }
3145 newpids.erase(newpids.begin());
3146 }
3147 }
3148 else
3149 {
3150 startPt--;
3151 }
3152 }
3153
3154 if (!endFound)
3155 {
3156 if (endPt > numFacePoints - 1)
3157 {
3158 endPt = 0;
3159 }
3160
3161 newpids.push_back(pids[endPt]);
3162 v = contourScalars->GetComponent(pids[endPt],0);
3163 if ((insideOut && (v > value-eps)) || ((!insideOut) && (v < value+eps)))
3164 {
3165 endFound = true;
3166 if ((insideOut && (v > value+eps)) || ((!insideOut) && (v < value-eps)))
3167 {
3168 vtkWarningMacro("A positive point is directly connected to a "
3169 "negative point with no contour point in between. We should "
3170 "never get here.");
3171 if (endPt == 0)
3172 {
3173 endPt = numFacePoints-1;
3174 }
3175 else
3176 {
3177 endPt--;
3178 }
3179 newpids.pop_back();
3180 }
3181 }
3182 else
3183 {
3184 endPt++;
3185 }
3186 }
3187 }// end inner for loop for finding start and end points
3188
3189 // if face are entirely positive, add it directly into the face list
3190 if (!startFound && !endFound)
3191 {
3192 visited[fid] = true;
3193 faces.push_back(pids);
3194 }
3195
3196 // if face contain contour points
3197 else if (startFound && endFound)
3198 {
3199 // a point or a line
3200 if (newpids.size() < 3)
3201 {
3202 visited[fid] = true;
3203 }
3204 // if face only contains one contour point, this is a special case that
3205 // may only happen when one of the original vertex is a contour point.
3206 // we will add this face to the result polyhedron.
3207 else if (startPt == endPt)
3208 {
3209 visited[fid] = true;
3210 faces.push_back(pids);
3211 }
3212 // Face contain at least two contour points. In this case, we will create
3213 // a new face patch whose close boundary is start point -->contour point
3214 // --> end point --> start point. Notice that the face may contain other
3215 // positive points and contour points. So we will not label the face as
3216 // visited. Instead, we will erase the chunk from start point to end
3217 // point from the point id vector of the face. So that the other part
3218 // can still be visited in the future.
3219 else
3220 {
3221 if (!this->Internal->EraseSegmentFromIdVector(
3222 pids, positivePt, startPt, endPt))
3223 {
3224 vtkErrorMacro("Erase segment from Id vector failed. We should "
3225 "never get here.");
3226 visited[fid] = true;
3227 continue;
3228 }
3229 if (pids.size()<=2) // all but two contour points are left
3230 {
3231 pids.clear();
3232 visited[fid] = true;
3233 }
3234 fpMapIt->second = pids;
3235 faces.push_back(newpids);
3236 }
3237 }
3238
3239 // only find start or only find end. this should never happen
3240 else
3241 {
3242 visited[fid] = true;
3243 vtkErrorMacro("We should never get here. Locating contour points failed. "
3244 "Clipping continues but may generate wrong result.");
3245 }
3246 } // end for each face
3247
3248 } // end for_pid
3249
3250 delete [] visited;
3251
3252 // not a valid output when the clip plane passes through the cell boundary
3253 // faces.
3254 if (faces.empty())
3255 {
3256 this->Internal->RestoreFaceArrayAndEdgeTable(this->Faces, this->EdgeTable);
3257 return;
3258 }
3259
3260 vtkIdType numAllFaces = contourPolys->GetNumberOfCells() +
3261 static_cast<vtkIdType>(faces.size());
3262 cellVector.push_back(numAllFaces);
3263
3264 // add contour faces
3265 contourPolys->InitTraversal();
3266 while (contourPolys->GetNextCell(npts, pts))
3267 {
3268 if (!this->Internal->ConvertPointIds(npts, pts, pointIdMap, insideOut))
3269 {
3270 vtkErrorMacro("Cannot find the id of an output point. We should never "
3271 "get here. Clipping aborted.");
3272 this->Internal->RestoreFaceArrayAndEdgeTable(this->Faces, this->EdgeTable);
3273 return;
3274 }
3275 cellVector.push_back(npts);
3276 cellVector.insert(cellVector.end(), pts, pts+npts);
3277 }
3278
3279 // add other faces
3280 for (size_t i = 0; i < faces.size(); i++)
3281 {
3282 vtkIdVectorType pids = faces[i];
3283 for (size_t j = 0; j < pids.size(); j++)
3284 {
3285 vtkIdType pid = pids[j];
3286 vtkIdToIdMapType::iterator iter = pointIdMap.find(pid);
3287 if (iter == pointIdMap.end()) // must be original points
3288 {
3289 if (locator->InsertUniquePoint(this->Points->GetPoint(pid), newPid))
3290 {
3291 vtkIdType globalPid = this->PointIds->GetId(pid);
3292 outPd->CopyData(inPd, globalPid, newPid);
3293 }
3294 pointIdMap.insert(vtkIdToIdPairType(pid, newPid));
3295 }
3296 }
3297
3298 npts = static_cast<vtkIdType>(pids.size());
3299 pts = &(pids[0]);
3300 if (!this->Internal->ConvertPointIds(npts, pts, pointIdMap))
3301 {
3302 vtkErrorMacro("Cannot find the id of an output point. We should never "
3303 "get here. Clipping aborted.");
3304 this->Internal->RestoreFaceArrayAndEdgeTable(this->Faces, this->EdgeTable);
3305 return;
3306 }
3307 cellVector.push_back(npts);
3308 cellVector.insert(cellVector.end(), pts, pts+npts);
3309 }
3310
3311 newCellId = connectivity->InsertNextCell(
3312 static_cast<vtkIdType>(cellVector.size()), &(cellVector[0]));
3313 outCd->CopyData(inCd, cellId, newCellId);
3314
3315 this->Internal->RestoreFaceArrayAndEdgeTable(this->Faces, this->EdgeTable);
3316 }
3317
3318 //----------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)3319 void vtkPolyhedron::PrintSelf(ostream& os, vtkIndent indent)
3320 {
3321 this->Superclass::PrintSelf(os,indent);
3322
3323 os << indent << "Triangle:\n";
3324 this->Triangle->PrintSelf(os,indent.GetNextIndent());
3325
3326 os << indent << "Polygon:\n";
3327 this->Polygon->PrintSelf(os,indent.GetNextIndent());
3328
3329 os << indent << "Tetra:\n";
3330 this->Tetra->PrintSelf(os,indent.GetNextIndent());
3331
3332 os << indent << "Faces:\n";
3333 this->GlobalFaces->PrintSelf(os,indent.GetNextIndent());
3334
3335 }
3336