1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkLoopBooleanPolyDataFilter.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 /** @file vtkLoopBooleanPolyDataFilter.cxx
16 * @brief This is the filter to perform boolean operations
17 * @author Adam Updegrove
18 * @author updega2@gmail.com
19 */
20
21 #include "vtkLoopBooleanPolyDataFilter.h"
22
23 #include "vtkAppendPolyData.h"
24 #include "vtkCellArray.h"
25 #include "vtkCellData.h"
26 #include "vtkCleanPolyData.h"
27 #include "vtkDoubleArray.h"
28 #include "vtkFloatArray.h"
29 #include "vtkGenericCell.h"
30 #include "vtkInformation.h"
31 #include "vtkInformationVector.h"
32 #include "vtkIntersectionPolyDataFilter.h"
33 #include "vtkMath.h"
34 #include "vtkMergeCells.h"
35 #include "vtkObjectFactory.h"
36 #include "vtkPointData.h"
37 #include "vtkPointLocator.h"
38 #include "vtkPolyDataNormals.h"
39 #include "vtkSmartPointer.h"
40 #include "vtkTransform.h"
41 #include "vtkTransformPolyDataFilter.h"
42 #include "vtkTriangle.h"
43 #include "vtkUnstructuredGrid.h"
44
45 #include <iostream>
46 #include <list>
47 #include <sstream>
48 #include <string>
49
50 //------------------------------------------------------------------------------
51 // Helper typedefs and data structures.
52 namespace
53 {
54
55 struct simLine
56 {
57 vtkIdType id;
58 vtkIdType pt1;
59 vtkIdType pt2;
60 };
61
62 struct simLoop
63 {
64 std::list<simLine> cells;
65 vtkIdType startPt;
66 vtkIdType endPt;
67 int loopType;
68 };
69
70 }
71
72 // Implementation function
73 class vtkLoopBooleanPolyDataFilter::Impl
74 {
75 public:
76 Impl();
77 virtual ~Impl();
78
79 void Initialize();
80 void SetCheckArrays();
81 void SetBoundaryArrays();
82 void ResetCheckArrays();
83 void GetBooleanRegions(int inputIndex, std::vector<simLoop>* loops);
84 void DetermineIntersection(std::vector<simLoop>* loops);
85 void PerformBoolean(vtkPolyData* output, int booleanOperation);
86 void ThresholdRegions(vtkPolyData** surfaces);
87
88 protected:
89 int RunLoopFind(vtkIdType interPt, vtkIdType nextCell, bool* usedPt, simLoop* loop);
90
91 int RunLoopTest(vtkIdType interPt, vtkIdType nextCell, simLoop* loop, bool* usedPt);
92
93 int GetCellOrientation(vtkPolyData* pd, vtkIdType cellId, vtkIdType p0, vtkIdType p1, int index);
94
95 int FindRegion(int inputIndex, int fillnumber, int start, int fill);
96
97 int FindRegionTipToe(int inputIndex, int fillnumber, int fill);
98
99 public:
100 int IntersectionCase;
101
102 vtkPolyData* Mesh[2];
103 vtkPolyData* IntersectionLines;
104
105 vtkIntArray* BoundaryPointArray[2];
106 vtkIntArray* BoundaryCellArray[2];
107 vtkIntArray* BooleanArray[2];
108 vtkIntArray* NewCellIds[2];
109
110 vtkIdType* Checked[2];
111 vtkIdType* CheckedCarefully[2];
112 vtkIdType* PointMapper[2];
113 vtkIdType* ReversePointMapper[2];
114
115 vtkIdList* CheckCells;
116 vtkIdList* CheckCells2;
117 vtkIdList* CheckCellsCareful;
118 vtkIdList* CheckCellsCareful2;
119
120 // Pointer to overarching filter
121 vtkLoopBooleanPolyDataFilter* ParentFilter;
122 };
123
Impl()124 vtkLoopBooleanPolyDataFilter::Impl::Impl()
125 : CheckCells(nullptr)
126 , CheckCells2(nullptr)
127 , CheckCellsCareful(nullptr)
128 , CheckCellsCareful2(nullptr)
129 {
130 for (int i = 0; i < 2; i++)
131 {
132 this->Mesh[i] = vtkPolyData::New();
133
134 this->BooleanArray[i] = vtkIntArray::New();
135 this->BoundaryPointArray[i] = vtkIntArray::New();
136 this->BoundaryCellArray[i] = vtkIntArray::New();
137 this->NewCellIds[i] = vtkIntArray::New();
138
139 this->Checked[i] = nullptr;
140 this->CheckedCarefully[i] = nullptr;
141 this->PointMapper[i] = nullptr;
142 this->ReversePointMapper[i] = nullptr;
143 }
144 this->IntersectionLines = vtkPolyData::New();
145 this->CheckCells = vtkIdList::New();
146 this->CheckCells2 = vtkIdList::New();
147 this->CheckCellsCareful = vtkIdList::New();
148 this->CheckCellsCareful2 = vtkIdList::New();
149
150 // Intersection Case:
151 // 0 -> Only hard closed intersection loops
152 // 1 -> At least one soft closed intersection loop
153 // 2 -> At least one open intersection loop
154 this->IntersectionCase = 0;
155 }
156
~Impl()157 vtkLoopBooleanPolyDataFilter::Impl::~Impl()
158 {
159 for (int i = 0; i < 2; i++)
160 {
161 this->Mesh[i]->Delete();
162 this->BooleanArray[i]->Delete();
163 this->BoundaryPointArray[i]->Delete();
164 this->BoundaryCellArray[i]->Delete();
165 this->NewCellIds[i]->Delete();
166
167 delete[] this->Checked[i];
168 delete[] this->CheckedCarefully[i];
169 delete[] this->PointMapper[i];
170 delete[] this->ReversePointMapper[i];
171 }
172 this->IntersectionLines->Delete();
173 this->CheckCells->Delete();
174 this->CheckCells2->Delete();
175 this->CheckCellsCareful->Delete();
176 this->CheckCellsCareful2->Delete();
177 }
178
179 // Flood fill algorithm to find region of mesh separated by intersection lines
FindRegion(int inputIndex,int fillnumber,int start,int fill)180 int vtkLoopBooleanPolyDataFilter::Impl::FindRegion(
181 int inputIndex, int fillnumber, int start, int fill)
182 {
183 vtkDebugWithObjectMacro(this->ParentFilter, << "Finding region with fill " << fillnumber
184 << " of mesh " << inputIndex << " with cellID "
185 << this->CheckCells->GetId(0));
186
187 // Id List to store neighbor cells for each set of nodes and a cell
188 vtkSmartPointer<vtkIdList> neighbors = vtkSmartPointer<vtkIdList>::New();
189 vtkSmartPointer<vtkIdList> tmp = vtkSmartPointer<vtkIdList>::New();
190
191 vtkIdType numCheckCells;
192 // Get neighboring cell for each pair of points in current cell
193 // While there are still cells to be checked, find neighbor cells
194 while ((numCheckCells = this->CheckCells->GetNumberOfIds()) > 0)
195 {
196 for (int c = 0; c < numCheckCells; c++)
197 {
198 vtkIdType cellId = this->CheckCells->GetId(c);
199 // Get the three points of the cell
200 const vtkIdType* pts = nullptr;
201 vtkIdType npts = 0;
202 this->Mesh[inputIndex]->GetCellPoints(cellId, npts, pts);
203 if (this->Checked[inputIndex][cellId] == 0)
204 {
205 // Mark cell as checked and insert the fillnumber value to cell
206 if (fill)
207 {
208 this->BooleanArray[inputIndex]->InsertValue(cellId, fillnumber);
209 }
210 this->Checked[inputIndex][cellId] = 1;
211 for (int i = 0; i < npts; i++)
212 {
213 vtkIdType p1 = pts[i];
214 // Get the cells attached to each point
215 this->Mesh[inputIndex]->GetPointCells(p1, neighbors);
216 vtkIdType numNeighbors = neighbors->GetNumberOfIds();
217
218 // For each neighboring cell
219 for (int j = 0; j < numNeighbors; j++)
220 {
221 // If this cell is close to a boundary
222 if (this->BoundaryCellArray[inputIndex]->GetValue(neighbors->GetId(j)))
223 {
224 // If this cell hasn't been checked already
225 if (this->CheckedCarefully[inputIndex][neighbors->GetId(j)] == 0)
226 {
227 // Add this cell to the careful check cells list and run
228 // the region finding tip toe code
229 this->CheckCellsCareful->InsertNextId(neighbors->GetId(j));
230 if (fill)
231 {
232 this->FindRegionTipToe(inputIndex, fillnumber, 1);
233 }
234 else
235 {
236 this->FindRegionTipToe(inputIndex, fillnumber, 0);
237 }
238 this->CheckCellsCareful->Reset();
239 this->CheckCellsCareful2->Reset();
240 }
241 }
242 // Cell needs to be added to check list
243 else
244 {
245 this->CheckCells2->InsertNextId(neighbors->GetId(j));
246 }
247 }
248 }
249 }
250 // This statement is for if the start cell is a boundary cell
251 else if (this->CheckedCarefully[inputIndex][cellId] == 0 && start)
252 {
253 start = 0;
254 this->CheckCells->Reset();
255 this->CheckCellsCareful->InsertNextId(cellId);
256 if (fill)
257 {
258 this->FindRegionTipToe(inputIndex, fillnumber, 1);
259 }
260 else
261 {
262 this->FindRegionTipToe(inputIndex, fillnumber, 0);
263 }
264 }
265 }
266
267 // Swap the current check list to the full check list and continue
268 tmp = this->CheckCells;
269 this->CheckCells = this->CheckCells2;
270 this->CheckCells2 = tmp;
271 tmp->Reset();
272 }
273 return 1;
274 }
275
276 // This is the slow version of the flood fill algorithm that is initiated
277 // when we get close to a boundary to ensure we don't cross the line
FindRegionTipToe(int inputIndex,int fillnumber,int fill)278 int vtkLoopBooleanPolyDataFilter::Impl::FindRegionTipToe(int inputIndex, int fillnumber, int fill)
279 {
280 // Id List to store neighbor cells for each set of nodes and a cell
281 vtkSmartPointer<vtkIdList> tmp = vtkSmartPointer<vtkIdList>::New();
282 vtkSmartPointer<vtkIdList> neighborIds = vtkSmartPointer<vtkIdList>::New();
283
284 vtkIdType numCheckCells;
285 // Get neighboring cell for each pair of points in current cell
286 // While there are still cells to be checked
287 while ((numCheckCells = this->CheckCellsCareful->GetNumberOfIds()) > 0)
288 {
289 for (int c = 0; c < numCheckCells; c++)
290 {
291 neighborIds->Reset();
292 vtkIdType cellId = this->CheckCellsCareful->GetId(c);
293 // Get the three points of the cell
294 const vtkIdType* pts = nullptr;
295 vtkIdType npts = 0;
296 this->Mesh[inputIndex]->GetCellPoints(cellId, npts, pts);
297 // Update this cell to have been checked carefully and assign it
298 // with the fillnumber scalar
299 if (this->CheckedCarefully[inputIndex][cellId] == 0)
300 {
301 if (fill)
302 {
303 this->BooleanArray[inputIndex]->InsertValue(cellId, fillnumber);
304 }
305 this->CheckedCarefully[inputIndex][cellId] = 1;
306 // For each edge of the cell
307 vtkDebugWithObjectMacro(this->ParentFilter, << "Checking edges of cell " << cellId);
308 for (int i = 0; i < npts; i++)
309 {
310 vtkIdType p1 = pts[i];
311 vtkIdType p2 = pts[(i + 1) % (npts)];
312
313 vtkSmartPointer<vtkIdList> neighbors = vtkSmartPointer<vtkIdList>::New();
314 // Initial check to make sure the cell is in fact a face cell
315 this->Mesh[inputIndex]->GetCellEdgeNeighbors(cellId, p1, p2, neighbors);
316 vtkIdType numNeighbors = neighbors->GetNumberOfIds();
317
318 // Check to make sure it is an oustide surface cell,
319 // i.e. one neighbor
320 if (numNeighbors == 1)
321 {
322 int count = 0;
323 // Check to see if cell is on the boundary,
324 // if it is get adjacent lines
325 if (this->BoundaryPointArray[inputIndex]->GetValue(p1) == 1)
326 {
327 count++;
328 }
329
330 if (this->BoundaryPointArray[inputIndex]->GetValue(p2) == 1)
331 {
332 count++;
333 }
334
335 vtkIdType neighbor = neighbors->GetId(0);
336 // if cell is not on the boundary, add new cell to check list
337 if (count < 2)
338 {
339 neighborIds->InsertNextId(neighbor);
340 }
341 // if cell is on boundary, check to make sure it isn't
342 // false positive; don't add to check list. This is done by
343 // getting the boundary lines attached to each point, then
344 // intersecting the two lists. If the result is zero, then this
345 // is a false positive
346 else
347 {
348 // Variables for the boundary cells adjacent to the boundary point
349 vtkSmartPointer<vtkIdList> bLinesOne = vtkSmartPointer<vtkIdList>::New();
350 vtkSmartPointer<vtkIdList> bLinesTwo = vtkSmartPointer<vtkIdList>::New();
351
352 vtkIdType bPt1 = PointMapper[inputIndex][p1];
353 this->IntersectionLines->GetPointCells(bPt1, bLinesOne);
354
355 vtkIdType bPt2 = PointMapper[inputIndex][p2];
356 this->IntersectionLines->GetPointCells(bPt2, bLinesTwo);
357
358 bLinesOne->IntersectWith(bLinesTwo);
359 // Cell is false positive. Add to check list.
360 if (bLinesOne->GetNumberOfIds() == 0)
361 {
362 vtkDebugWithObjectMacro(this->ParentFilter, << "False positive! " << neighbor);
363 neighborIds->InsertNextId(neighbor);
364 }
365 else
366 {
367 vtkDebugWithObjectMacro(
368 this->ParentFilter, << "I have not been added because false");
369 }
370 }
371 }
372 else
373 {
374 vtkDebugWithObjectMacro(this->ParentFilter, << "NumNei is not 1");
375 vtkDebugWithObjectMacro(this->ParentFilter, << "Number of Neighbors " << numNeighbors);
376 vtkDebugWithObjectMacro(this->ParentFilter, << "Cell is " << cellId);
377 for (int k = 0; k < numNeighbors; k++)
378 {
379 vtkDebugWithObjectMacro(this->ParentFilter, << "Id!!! " << neighbors->GetId(k));
380 }
381 }
382 }
383
384 vtkIdType numIds = neighborIds->GetNumberOfIds();
385 if (numIds > 0)
386 {
387 // Add all Ids in current list to global list of Ids
388 for (int k = 0; k < numIds; k++)
389 {
390 vtkIdType neighborId = neighborIds->GetId(k);
391 if (this->CheckedCarefully[inputIndex][neighborId] == 0)
392 {
393 this->CheckCellsCareful2->InsertNextId(neighborId);
394 }
395 else if (this->Checked[inputIndex][neighborId] == 0)
396 {
397 this->CheckCells2->InsertNextId(neighborId);
398 }
399 }
400 }
401 }
402 }
403
404 // Add current list of checked cells to the full list and continue
405 tmp = this->CheckCellsCareful;
406 this->CheckCellsCareful = this->CheckCellsCareful2;
407 this->CheckCellsCareful2 = tmp;
408 tmp->Reset();
409 }
410 return 1;
411 }
412
Initialize()413 void vtkLoopBooleanPolyDataFilter::Impl::Initialize()
414 {
415 for (int i = 0; i < 2; i++)
416 {
417 if (this->Mesh[i]->GetNumberOfPoints() == 0 || this->Mesh[i]->GetNumberOfCells() == 0)
418 {
419 vtkGenericWarningMacro(<< "Mesh has zero points or cells and "
420 << "cannot run filter");
421 return;
422 }
423
424 // Get the number of Polys for scalar allocation
425 int numPolys = this->Mesh[i]->GetNumberOfPolys();
426 int numPts = this->Mesh[i]->GetNumberOfPoints();
427 int numLinePts = this->IntersectionLines->GetNumberOfPoints();
428
429 // Allocate space for each Boundary Array and the fill array
430 this->BoundaryPointArray[i]->SetNumberOfTuples(numPts);
431 this->BoundaryCellArray[i]->SetNumberOfTuples(numPolys);
432 this->BooleanArray[i]->SetNumberOfTuples(numPolys);
433 this->Checked[i] = new vtkIdType[numPolys];
434 this->CheckedCarefully[i] = new vtkIdType[numPolys];
435 this->PointMapper[i] = new vtkIdType[numPolys];
436 this->ReversePointMapper[i] = new vtkIdType[numLinePts];
437
438 for (int j = 0; j < numPts; j++)
439 {
440 this->BoundaryPointArray[i]->InsertValue(j, 0);
441 }
442 for (int j = 0; j < numPolys; j++)
443 {
444 this->BoundaryCellArray[i]->InsertValue(j, 0);
445 this->BooleanArray[i]->InsertValue(j, 0);
446 this->Checked[i][j] = 0;
447 this->CheckedCarefully[i][j] = 0;
448 this->PointMapper[i][j] = -1;
449 }
450 for (int j = 0; j < numLinePts; j++)
451 {
452 this->ReversePointMapper[i][j] = -1;
453 }
454 }
455 this->NewCellIds[0]->DeepCopy(this->IntersectionLines->GetCellData()->GetArray("NewCell0ID"));
456 this->NewCellIds[1]->DeepCopy(this->IntersectionLines->GetCellData()->GetArray("NewCell1ID"));
457
458 this->BooleanArray[0]->SetName("BooleanRegion");
459 this->BooleanArray[1]->SetName("BooleanRegion");
460 this->Mesh[0]->GetCellData()->AddArray(this->BooleanArray[0]);
461 this->Mesh[0]->GetCellData()->SetActiveScalars("BooleanRegion");
462 this->Mesh[1]->GetCellData()->AddArray(this->BooleanArray[1]);
463 this->Mesh[1]->GetCellData()->SetActiveScalars("BooleanRegion");
464
465 this->BoundaryCellArray[0]->SetName("BoundaryCells");
466 this->BoundaryCellArray[1]->SetName("BoundaryCells");
467 this->Mesh[0]->GetCellData()->AddArray(this->BoundaryCellArray[0]);
468 this->Mesh[0]->GetCellData()->SetActiveScalars("BoundaryCells");
469 this->Mesh[1]->GetCellData()->AddArray(this->BoundaryCellArray[1]);
470 this->Mesh[1]->GetCellData()->SetActiveScalars("BoundaryCells");
471
472 this->BoundaryPointArray[0]->SetName("BoundaryPoints");
473 this->BoundaryPointArray[1]->SetName("BoundaryPoints");
474 this->Mesh[0]->GetPointData()->AddArray(this->BoundaryPointArray[0]);
475 this->Mesh[0]->GetPointData()->SetActiveScalars("BoundaryPoints");
476 this->Mesh[1]->GetPointData()->AddArray(this->BoundaryPointArray[1]);
477 this->Mesh[1]->GetPointData()->SetActiveScalars("BoundaryPoints");
478 }
479
480 // Function to find the regions on each input separated by the intersection
481 // lines
GetBooleanRegions(int inputIndex,std::vector<simLoop> * loops)482 void vtkLoopBooleanPolyDataFilter::Impl::GetBooleanRegions(
483 int inputIndex, std::vector<simLoop>* loops)
484 {
485 vtkSmartPointer<vtkPolyData> tmpPolyData = vtkSmartPointer<vtkPolyData>::New();
486 tmpPolyData->DeepCopy(this->Mesh[inputIndex]);
487 tmpPolyData->BuildLinks();
488
489 std::vector<simLoop>::iterator loopit;
490 std::list<simLine>::iterator cellit;
491
492 // For each intersection loop
493 for (loopit = loops->begin(); loopit != loops->end(); ++loopit)
494 {
495 std::list<simLine> loopcells = (loopit)->cells;
496 // Go through each cell in the loop
497 for (cellit = loopcells.begin(); cellit != loopcells.end(); ++cellit)
498 {
499 simLine nextLine;
500 nextLine = *cellit;
501 vtkIdType nextCell = nextLine.id;
502 vtkIdType p1 = nextLine.pt1;
503 vtkIdType p2 = nextLine.pt2;
504 vtkIdType outputCellId0 = this->NewCellIds[inputIndex]->GetComponent(nextCell, 0);
505 vtkIdType outputCellId1 = this->NewCellIds[inputIndex]->GetComponent(nextCell, 1);
506 // If the cell has not been given an orientation from the flood fill
507 // algorithm and it has an id from vtkIntersectionPolyDataFilter
508 if (outputCellId0 != -1)
509 {
510 // If the cell hasn't been touched
511 if (this->CheckedCarefully[inputIndex][outputCellId0] == 0)
512 {
513 int sign1 = this->GetCellOrientation(tmpPolyData, outputCellId0, p1, p2, inputIndex);
514 // If cell orientation is found
515 if (sign1 != 0)
516 {
517 this->CheckCells->InsertNextId(outputCellId0);
518 this->FindRegion(inputIndex, sign1, 1, 1);
519 this->CheckCells->Reset();
520 this->CheckCells2->Reset();
521 this->CheckCellsCareful->Reset();
522 this->CheckCellsCareful2->Reset();
523 }
524 }
525 }
526 // Check cell on other side of intersection line
527 if (outputCellId1 != -1)
528 {
529 // If the cell hasn't been touched
530 if (this->CheckedCarefully[inputIndex][outputCellId1] == 0)
531 {
532 int sign2 = this->GetCellOrientation(tmpPolyData, outputCellId1, p1, p2, inputIndex);
533 // If cell orientation is found
534 if (sign2 != 0)
535 {
536 this->CheckCells->InsertNextId(outputCellId1);
537 this->FindRegion(inputIndex, sign2, 1, 1);
538 this->CheckCells->Reset();
539 this->CheckCells2->Reset();
540 this->CheckCellsCareful->Reset();
541 this->CheckCellsCareful2->Reset();
542 }
543 }
544 }
545 }
546 }
547 }
548
549 // Get Cell orientation so we know which value to flood fill a region with
GetCellOrientation(vtkPolyData * pd,vtkIdType cellId,vtkIdType p0,vtkIdType p1,int index)550 int vtkLoopBooleanPolyDataFilter::Impl::GetCellOrientation(
551 vtkPolyData* pd, vtkIdType cellId, vtkIdType p0, vtkIdType p1, int index)
552 {
553 vtkDebugWithObjectMacro(this->ParentFilter, << "CellId: " << cellId);
554 vtkIdType npts;
555 const vtkIdType* pts;
556 pd->BuildLinks();
557 pd->GetCellPoints(cellId, npts, pts);
558 // pt0Id and pt1Id are from intersectionLines PolyData and I am trying
559 // to compare these to the point ids in pd.
560 vtkIdType cellPtId0 = this->ReversePointMapper[index][p0];
561 vtkIdType cellPtId1 = this->ReversePointMapper[index][p1];
562 double points[3][3];
563 vtkIdType cellPtId2 = 0;
564 for (int j = 0; j < npts; j++)
565 {
566 pd->GetPoint(pts[j], points[j]);
567 if (cellPtId0 != pts[j] && cellPtId1 != pts[j])
568 {
569 cellPtId2 = pts[j];
570 }
571 }
572 vtkSmartPointer<vtkPoints> cellPts = vtkSmartPointer<vtkPoints>::New();
573 cellPts->InsertNextPoint(pd->GetPoint(cellPtId0));
574 cellPts->InsertNextPoint(pd->GetPoint(cellPtId1));
575 cellPts->InsertNextPoint(pd->GetPoint(cellPtId2));
576
577 vtkSmartPointer<vtkPolyData> cellPD = vtkSmartPointer<vtkPolyData>::New();
578 cellPD->SetPoints(cellPts);
579
580 vtkSmartPointer<vtkCellArray> cellLines = vtkSmartPointer<vtkCellArray>::New();
581 for (int j = 0; j < npts; j++)
582 {
583 int spot1 = j;
584 int spot2 = (j + 1) % 3;
585 cellLines->InsertNextCell(2);
586 cellLines->InsertCellPoint(spot1);
587 cellLines->InsertCellPoint(spot2);
588 }
589 cellPD->SetLines(cellLines);
590
591 // Set up a transform that will rotate the points to the
592 // XY-plane (normal aligned with z-axis).
593 vtkSmartPointer<vtkTransform> transform = vtkSmartPointer<vtkTransform>::New();
594 double zaxis[3] = { 0.0, 0.0, 1.0 };
595 double rotationAxis[3], normal[3], center[3], rotationAngle;
596
597 vtkTriangle::ComputeNormal(points[0], points[1], points[2], normal);
598
599 double dotZAxis = vtkMath::Dot(normal, zaxis);
600 if (fabs(1.0 - dotZAxis) < 1e-6)
601 {
602 // Aligned with z-axis
603 rotationAxis[0] = 1.0;
604 rotationAxis[1] = 0.0;
605 rotationAxis[2] = 0.0;
606 rotationAngle = 0.0;
607 }
608 else if (fabs(1.0 + dotZAxis) < 1e-6)
609 {
610 // Co-linear with z-axis, but reversed sense.
611 // Aligned with z-axis
612 rotationAxis[0] = 1.0;
613 rotationAxis[1] = 0.0;
614 rotationAxis[2] = 0.0;
615 rotationAngle = 180.0;
616 }
617 else
618 {
619 // The general case
620 vtkMath::Cross(normal, zaxis, rotationAxis);
621 vtkMath::Normalize(rotationAxis);
622 rotationAngle = vtkMath::DegreesFromRadians(acos(vtkMath::Dot(zaxis, normal)));
623 }
624
625 transform->PreMultiply();
626 transform->Identity();
627
628 vtkDebugWithObjectMacro(this->ParentFilter, << "ROTATION ANGLE " << rotationAngle);
629 transform->RotateWXYZ(rotationAngle, rotationAxis[0], rotationAxis[1], rotationAxis[2]);
630
631 vtkTriangle::TriangleCenter(points[0], points[1], points[2], center);
632 transform->Translate(-center[0], -center[1], -center[2]);
633
634 vtkSmartPointer<vtkTransformPolyDataFilter> transformer =
635 vtkSmartPointer<vtkTransformPolyDataFilter>::New();
636 transformer->SetInputData(cellPD);
637 transformer->SetTransform(transform);
638 transformer->Update();
639
640 vtkSmartPointer<vtkPolyData> transPD = vtkSmartPointer<vtkPolyData>::New();
641 transPD = transformer->GetOutput();
642 transPD->BuildLinks();
643
644 double area = 0;
645 double tedgept1[3];
646 double tedgept2[3];
647 vtkIdType newpt;
648 for (newpt = 0; newpt < transPD->GetNumberOfPoints() - 1; newpt++)
649 {
650 transPD->GetPoint(newpt, tedgept1);
651 transPD->GetPoint(newpt + 1, tedgept2);
652 area = area + (tedgept1[0] * tedgept2[1]) - (tedgept2[0] * tedgept1[1]);
653 }
654 transPD->GetPoint(newpt, tedgept1);
655 transPD->GetPoint(0, tedgept2);
656 area = area + (tedgept1[0] * tedgept2[1]) - (tedgept2[0] * tedgept1[1]);
657
658 int value = 0;
659 double tolerance = 1e-6;
660 if (area < 0 && fabs(area) > tolerance)
661 {
662 value = -1;
663 }
664 else if (area > 0 && fabs(area) > tolerance)
665 {
666 value = 1;
667 }
668 else
669 {
670 vtkDebugWithObjectMacro(this->ParentFilter, << "Line pts are " << p0 << " and " << p1);
671 vtkDebugWithObjectMacro(
672 this->ParentFilter, << "PD pts are " << cellPtId0 << " and " << cellPtId1);
673 value = 0;
674 }
675
676 return value;
677 }
678
679 // Reset the find region arrays to test another region
ResetCheckArrays()680 void vtkLoopBooleanPolyDataFilter::Impl::ResetCheckArrays()
681 {
682 for (int i = 0; i < 2; i++)
683 {
684 int numPolys = this->Mesh[i]->GetNumberOfCells();
685 for (vtkIdType cellId = 0; cellId < numPolys; cellId++)
686 {
687 if (this->BoundaryCellArray[i]->GetValue(cellId) == 1)
688 {
689 this->Checked[i][cellId] = 1;
690 this->CheckedCarefully[i][cellId] = 0;
691 }
692 else
693 {
694 this->Checked[i][cellId] = 0;
695 this->CheckedCarefully[i][cellId] = 1;
696 }
697 }
698 }
699 }
700
701 // Set the boundary arrays on the mesh
SetBoundaryArrays()702 void vtkLoopBooleanPolyDataFilter::Impl::SetBoundaryArrays()
703 {
704 // Variables used in the function
705 // Point locator to find points on mesh that are the points on the boundary
706 // lines
707 vtkSmartPointer<vtkPointLocator> pointLocator1 = vtkSmartPointer<vtkPointLocator>::New();
708 vtkSmartPointer<vtkPointLocator> pointLocator2 = vtkSmartPointer<vtkPointLocator>::New();
709 pointLocator1->SetDataSet(this->Mesh[0]);
710 pointLocator1->BuildLocator();
711 pointLocator2->SetDataSet(this->Mesh[1]);
712 pointLocator2->BuildLocator();
713
714 int numPoints = this->IntersectionLines->GetNumberOfPoints();
715
716 for (vtkIdType pointId = 0; pointId < numPoints; pointId++)
717 {
718 double pt[3];
719 this->IntersectionLines->GetPoint(pointId, pt);
720 // Find point on mesh
721 vtkIdType bp1 = pointLocator1->FindClosestPoint(pt);
722 this->PointMapper[0][bp1] = pointId;
723 this->ReversePointMapper[0][pointId] = bp1;
724 this->BoundaryPointArray[0]->InsertValue(bp1, 1);
725 vtkSmartPointer<vtkIdList> bpCellIds1 = vtkSmartPointer<vtkIdList>::New();
726 this->Mesh[0]->GetPointCells(bp1, bpCellIds1);
727 // Set the point mapping array
728 // Assign each cell attached to this point as a boundary cell
729 for (int i = 0; i < bpCellIds1->GetNumberOfIds(); i++)
730 {
731 this->BoundaryCellArray[0]->InsertValue(bpCellIds1->GetId(i), 1);
732 this->Checked[0][bpCellIds1->GetId(i)] = 1;
733 }
734
735 vtkIdType bp2 = pointLocator2->FindClosestPoint(pt);
736 this->PointMapper[1][bp2] = pointId;
737 this->ReversePointMapper[1][pointId] = bp2;
738 this->BoundaryPointArray[1]->InsertValue(bp2, 1);
739 vtkSmartPointer<vtkIdList> bpCellIds2 = vtkSmartPointer<vtkIdList>::New();
740 this->Mesh[1]->GetPointCells(bp2, bpCellIds2);
741 // Set the point mapping array
742 // Assign each cell attached to this point as a boundary cell
743 for (int i = 0; i < bpCellIds2->GetNumberOfIds(); i++)
744 {
745 this->BoundaryCellArray[1]->InsertValue(bpCellIds2->GetId(i), 1);
746 this->Checked[1][bpCellIds2->GetId(i)] = 1;
747 }
748 }
749 }
750
751 // Set the original finding region check arrays
SetCheckArrays()752 void vtkLoopBooleanPolyDataFilter::Impl::SetCheckArrays()
753 {
754 for (int i = 0; i < 2; i++)
755 {
756 // Get the number of Polys for scalar allocation
757 int numPolys = this->Mesh[i]->GetNumberOfPolys();
758
759 for (int j = 0; j < numPolys; j++)
760 {
761 if (this->Checked[i][j] == 0)
762 {
763 this->CheckedCarefully[i][j] = 1;
764 }
765 else
766 {
767 this->CheckedCarefully[i][j] = 0;
768 }
769 }
770 }
771 }
772 //------------------------------------------------------------------------------
773
774 vtkStandardNewMacro(vtkLoopBooleanPolyDataFilter);
775
776 //------------------------------------------------------------------------------
vtkLoopBooleanPolyDataFilter()777 vtkLoopBooleanPolyDataFilter::vtkLoopBooleanPolyDataFilter()
778 {
779 this->Operation = VTK_UNION;
780
781 this->SetNumberOfInputPorts(2);
782 this->SetNumberOfOutputPorts(2);
783 this->NoIntersectionOutput = 1;
784
785 this->NumberOfIntersectionPoints = 0;
786 this->NumberOfIntersectionLines = 0;
787
788 this->Status = 1;
789 this->Tolerance = 1e-6;
790 }
791
792 //------------------------------------------------------------------------------
793 vtkLoopBooleanPolyDataFilter::~vtkLoopBooleanPolyDataFilter() = default;
794
795 //------------------------------------------------------------------------------
RequestData(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)796 int vtkLoopBooleanPolyDataFilter::RequestData(vtkInformation* vtkNotUsed(request),
797 vtkInformationVector** inputVector, vtkInformationVector* outputVector)
798 {
799 vtkInformation* inInfo0 = inputVector[0]->GetInformationObject(0);
800 vtkInformation* inInfo1 = inputVector[1]->GetInformationObject(0);
801 vtkInformation* outInfo0 = outputVector->GetInformationObject(0);
802 vtkInformation* outInfo1 = outputVector->GetInformationObject(1);
803
804 if (!inInfo0 || !inInfo1 || !outInfo0 || !outInfo1)
805 {
806 this->Status = 0;
807 return 0;
808 }
809
810 vtkPolyData* input0 = vtkPolyData::SafeDownCast(inInfo0->Get(vtkDataObject::DATA_OBJECT()));
811 vtkPolyData* input1 = vtkPolyData::SafeDownCast(inInfo1->Get(vtkDataObject::DATA_OBJECT()));
812
813 vtkPolyData* outputSurface =
814 vtkPolyData::SafeDownCast(outInfo0->Get(vtkDataObject::DATA_OBJECT()));
815 vtkPolyData* outputIntersection =
816 vtkPolyData::SafeDownCast(outInfo1->Get(vtkDataObject::DATA_OBJECT()));
817
818 if (!input0 || !input1 || !outputSurface || !outputIntersection)
819 {
820 this->Status = 0;
821 return 0;
822 }
823
824 // Get intersected versions
825 vtkSmartPointer<vtkIntersectionPolyDataFilter> polydataIntersection =
826 vtkSmartPointer<vtkIntersectionPolyDataFilter>::New();
827 polydataIntersection->SetInputConnection(0, this->GetInputConnection(0, 0));
828 polydataIntersection->SetInputConnection(1, this->GetInputConnection(1, 0));
829 polydataIntersection->SplitFirstOutputOn();
830 polydataIntersection->SplitSecondOutputOn();
831 polydataIntersection->SetTolerance(this->Tolerance);
832 polydataIntersection->Update();
833 if (polydataIntersection->GetStatus() != 1)
834 {
835 this->Status = 0;
836 return 0;
837 }
838
839 this->NumberOfIntersectionPoints = polydataIntersection->GetNumberOfIntersectionPoints();
840 this->NumberOfIntersectionLines = polydataIntersection->GetNumberOfIntersectionLines();
841
842 vtkDebugMacro(<< "Intersection is Done!!!");
843
844 vtkLoopBooleanPolyDataFilter::Impl* impl = new vtkLoopBooleanPolyDataFilter::Impl();
845 impl->ParentFilter = this;
846 impl->Mesh[0]->DeepCopy(polydataIntersection->GetOutput(1));
847 impl->Mesh[0]->BuildLinks();
848 impl->Mesh[1]->DeepCopy(polydataIntersection->GetOutput(2));
849 impl->Mesh[1]->BuildLinks();
850 impl->IntersectionLines->ShallowCopy(polydataIntersection->GetOutput(0));
851
852 if (this->NumberOfIntersectionPoints == 0 || this->NumberOfIntersectionLines == 0)
853 {
854 vtkWarningMacro(<< "No intersections!");
855 if (this->NoIntersectionOutput == 0)
856 {
857 delete impl;
858 return 1;
859 }
860 else
861 {
862 for (int i = 0; i < 2; i++)
863 {
864 // Get the number of Polys for scalar allocation
865 int numPolys = impl->Mesh[i]->GetNumberOfPolys();
866 int numPts = impl->Mesh[i]->GetNumberOfPoints();
867 for (int j = 0; j < numPts; j++)
868 {
869 impl->BoundaryPointArray[i]->InsertValue(j, 0);
870 }
871
872 for (int j = 0; j < numPolys; j++)
873 {
874 impl->BoundaryCellArray[i]->InsertValue(j, 0);
875 }
876
877 impl->BoundaryCellArray[i]->SetName("BoundaryCells");
878 impl->Mesh[i]->GetCellData()->AddArray(impl->BoundaryCellArray[i]);
879 impl->Mesh[i]->GetCellData()->SetActiveScalars("BoundaryCells");
880
881 impl->BoundaryPointArray[i]->SetName("BoundaryPoints");
882 impl->Mesh[i]->GetPointData()->AddArray(impl->BoundaryPointArray[i]);
883 impl->Mesh[i]->GetPointData()->SetActiveScalars("BoundaryPoints");
884 }
885 if (this->NoIntersectionOutput == 1)
886 {
887 vtkDebugMacro(<< "Only returning first surface");
888 outputSurface->DeepCopy(impl->Mesh[0]);
889 }
890 else if (this->NoIntersectionOutput == 2)
891 {
892 vtkDebugMacro(<< "Only returning second surface");
893 outputSurface->DeepCopy(impl->Mesh[1]);
894 }
895 else
896 {
897 vtkDebugMacro(<< "Keeping both");
898
899 vtkSmartPointer<vtkAppendPolyData> appender = vtkSmartPointer<vtkAppendPolyData>::New();
900 appender->AddInputData(impl->Mesh[0]);
901 appender->AddInputData(impl->Mesh[1]);
902 appender->Update();
903 outputSurface->DeepCopy(appender->GetOutput());
904 }
905 delete impl;
906 return 1;
907 }
908 }
909
910 double badtri1[2], badtri2[2];
911 double freeedge1[2], freeedge2[2];
912 impl->Mesh[0]->GetCellData()->GetArray("BadTriangle")->GetRange(badtri1, 0);
913 impl->Mesh[0]->GetCellData()->GetArray("FreeEdge")->GetRange(freeedge1, 0);
914
915 impl->Mesh[1]->GetCellData()->GetArray("BadTriangle")->GetRange(badtri2, 0);
916 impl->Mesh[1]->GetCellData()->GetArray("FreeEdge")->GetRange(freeedge2, 0);
917
918 // Set the check and boundary arrays for region finding
919 vtkDebugMacro(<< "Initializing");
920 impl->Initialize();
921 vtkDebugMacro(<< "Setting Bound Arrays");
922 impl->SetBoundaryArrays();
923 vtkDebugMacro(<< "Setting Check Arrays");
924 impl->SetCheckArrays();
925
926 // Determine the intersection type and obtain the intersection loops
927 // to give to Boolean Region finding
928 vtkDebugMacro(<< "Determining Intersection Type");
929 std::vector<simLoop> loops;
930 impl->DetermineIntersection(&loops);
931
932 // Get the regions bounded by the intersection lines and give correct
933 // orientation
934 impl->GetBooleanRegions(0, &loops);
935 vtkDebugMacro(<< "DONE WITH 1");
936 impl->GetBooleanRegions(1, &loops);
937 vtkDebugMacro(<< "DONE WITH 2");
938
939 // Combine certain orientations based on the operation desired
940 impl->PerformBoolean(outputSurface, this->Operation);
941
942 // Number of bad triangles and free edges (Should be zero for watertight,
943 // manifold surfaces!
944 vtkDebugMacro(<< "SURFACE 1 BAD TRI MIN: " << badtri1[0] << " MAX: " << badtri1[1]);
945 vtkDebugMacro(<< "SURFACE 1 FREE EDGE MIN: " << freeedge1[0] << " MAX: " << freeedge1[1]);
946 vtkDebugMacro(<< "SURFACE 2 BAD TRI MIN: " << badtri2[0] << " MAX: " << badtri2[1]);
947 vtkDebugMacro(<< "SURFACE 2 FREE EDGE MIN: " << freeedge2[0] << " MAX: " << freeedge2[1]);
948
949 double fullbadtri[2], fullfreeedge[2], dummy[2];
950 vtkIntersectionPolyDataFilter::CleanAndCheckSurface(outputSurface, dummy, this->Tolerance);
951 outputSurface->GetCellData()->GetArray("BadTriangle")->GetRange(fullbadtri, 0);
952 outputSurface->GetCellData()->GetArray("FreeEdge")->GetRange(fullfreeedge, 0);
953
954 // Add Normals
955 vtkSmartPointer<vtkPolyDataNormals> normaler = vtkSmartPointer<vtkPolyDataNormals>::New();
956 normaler->SetInputData(outputSurface);
957 normaler->AutoOrientNormalsOn();
958 normaler->Update();
959 outputSurface->DeepCopy(normaler->GetOutput());
960
961 vtkDebugMacro(<< "FULL SURFACE BAD TRI MIN: " << fullbadtri[0]);
962 vtkDebugMacro(<< " MAX: " << fullbadtri[1]);
963 vtkDebugMacro(<< "FULL SURFACE FREE EDGE MIN: " << fullfreeedge[0]);
964 vtkDebugMacro(<< " MAX: " << fullfreeedge[1]);
965
966 delete impl;
967 return 1;
968 }
969
970 //------------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)971 void vtkLoopBooleanPolyDataFilter::PrintSelf(ostream& os, vtkIndent indent)
972 {
973 this->Superclass::PrintSelf(os, indent);
974
975 os << indent << "Operation: ";
976 switch (this->Operation)
977 {
978 case VTK_UNION:
979 os << "UNION";
980 break;
981
982 case VTK_INTERSECTION:
983 os << "INTERSECTION";
984 break;
985
986 case VTK_DIFFERENCE:
987 os << "DIFFERENCE";
988 break;
989 }
990 os << "\n";
991 os << indent << "No Intersection Output: " << this->NoIntersectionOutput << "\n";
992 os << indent << "Tolerance: " << this->Tolerance << "\n";
993 os << indent << "NumberOfIntersectionPoints: " << this->NumberOfIntersectionPoints << "\n";
994 os << indent << "NumberOfIntersectionLines: " << this->NumberOfIntersectionLines << "\n";
995 }
996
997 //------------------------------------------------------------------------------
FillInputPortInformation(int port,vtkInformation * info)998 int vtkLoopBooleanPolyDataFilter::FillInputPortInformation(int port, vtkInformation* info)
999 {
1000 if (!this->Superclass::FillInputPortInformation(port, info))
1001 {
1002 return 0;
1003 }
1004 if (port == 0)
1005 {
1006 info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkPolyData");
1007 }
1008 else if (port == 1)
1009 {
1010 info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkPolyData");
1011 info->Set(vtkAlgorithm::INPUT_IS_OPTIONAL(), 0);
1012 }
1013 return 1;
1014 }
1015
1016 //------------------------------------------------------------------------------
1017
1018 // Determine type of intersection
DetermineIntersection(std::vector<simLoop> * loops)1019 void vtkLoopBooleanPolyDataFilter::Impl::DetermineIntersection(std::vector<simLoop>* loops)
1020 {
1021 int numInterPts = this->IntersectionLines->GetNumberOfPoints();
1022 bool* usedPt;
1023 usedPt = new bool[numInterPts];
1024 for (vtkIdType interPt = 0; interPt < numInterPts; interPt++)
1025 {
1026 usedPt[interPt] = false;
1027 }
1028
1029 for (vtkIdType interPt = 0; interPt < numInterPts; interPt++)
1030 {
1031 if (usedPt[interPt] == false)
1032 {
1033 simLoop newloop;
1034 vtkSmartPointer<vtkIdList> cellIds = vtkSmartPointer<vtkIdList>::New();
1035 this->IntersectionLines->GetPointCells(interPt, cellIds);
1036 if (cellIds->GetNumberOfIds() > 2)
1037 {
1038 vtkDebugWithObjectMacro(
1039 this->ParentFilter, << "Number Of Cells is greater than 2 for first point " << interPt);
1040 }
1041 else if (cellIds->GetNumberOfIds() < 2)
1042 {
1043 vtkDebugWithObjectMacro(
1044 this->ParentFilter, << "Number Of Cells is less than 2 for point " << interPt);
1045 }
1046
1047 vtkIdType nextCell = cellIds->GetId(0);
1048
1049 // Run through intersection lines to get loops!
1050 newloop.startPt = interPt;
1051 int caseId = 0;
1052 caseId = this->RunLoopFind(interPt, nextCell, usedPt, &newloop);
1053 if (caseId != -1)
1054 {
1055 // If the intersection loop is open
1056 if (this->IntersectionCase == 2)
1057 {
1058 vtkIdType nextPt = caseId;
1059 vtkDebugWithObjectMacro(this->ParentFilter, << "End point of open loop is " << nextPt);
1060 newloop.endPt = nextPt;
1061 newloop.loopType = 2;
1062 nextCell = cellIds->GetId(1);
1063 vtkIdType newId = this->RunLoopFind(interPt, nextCell, usedPt, &newloop);
1064 newloop.startPt = newId;
1065 // Save start and end point in custom data structure for loop
1066 }
1067 else
1068 {
1069 newloop.loopType = 1;
1070 }
1071 }
1072 usedPt[interPt] = true;
1073 loops->push_back(newloop);
1074 }
1075 }
1076 vtkDebugWithObjectMacro(this->ParentFilter, << "Number Of Loops: " << loops->size());
1077
1078 delete[] usedPt;
1079 }
1080
1081 // Function for if the intersection is soft closed or open
RunLoopFind(vtkIdType interPt,vtkIdType nextCell,bool * usedPt,simLoop * loop)1082 int vtkLoopBooleanPolyDataFilter::Impl::RunLoopFind(
1083 vtkIdType interPt, vtkIdType nextCell, bool* usedPt, simLoop* loop)
1084 {
1085 vtkIdType prevPt = interPt;
1086 vtkIdType nextPt = interPt;
1087 vtkSmartPointer<vtkIdList> pointIds = vtkSmartPointer<vtkIdList>::New();
1088 vtkSmartPointer<vtkIdList> cellIds = vtkSmartPointer<vtkIdList>::New();
1089
1090 IntersectionLines->GetCellPoints(nextCell, pointIds);
1091 if (pointIds->GetNumberOfIds() > 2)
1092 {
1093 vtkDebugWithObjectMacro(
1094 this->ParentFilter, << "Number Of Points is greater than 2 for first cell " << nextCell);
1095 }
1096 else if (pointIds->GetNumberOfIds() < 2)
1097 {
1098 vtkDebugWithObjectMacro(
1099 this->ParentFilter, << "Number Of Points is less than 2 for first cell " << nextCell);
1100 }
1101
1102 if (pointIds->GetId(0) == nextPt)
1103 {
1104 nextPt = pointIds->GetId(1);
1105 }
1106 else
1107 {
1108 nextPt = pointIds->GetId(0);
1109 }
1110 simLine newline;
1111 newline.pt1 = prevPt;
1112 newline.pt2 = nextPt;
1113 newline.id = nextCell;
1114 loop->cells.push_back(newline);
1115
1116 usedPt[nextPt] = true;
1117 while (nextPt != loop->cells.front().pt1)
1118 {
1119 IntersectionLines->GetPointCells(nextPt, cellIds);
1120 if (cellIds->GetNumberOfIds() > 2)
1121 {
1122 IntersectionCase = 1;
1123 vtkDebugWithObjectMacro(
1124 this->ParentFilter, << "Number Of Cells is greater than 2 for point " << nextPt);
1125 usedPt[nextPt] = false;
1126 nextCell = this->RunLoopTest(nextPt, nextCell, loop, usedPt);
1127 if (nextCell == -1)
1128 {
1129 break;
1130 }
1131 vtkDebugWithObjectMacro(this->ParentFilter, << "Next cell is " << nextCell);
1132 }
1133 else if (cellIds->GetNumberOfIds() < 2)
1134 {
1135 vtkDebugWithObjectMacro(
1136 this->ParentFilter, << "Number Of Cells is less than 2 for point " << nextPt);
1137 IntersectionCase = 2;
1138 return nextPt;
1139 }
1140 else
1141 {
1142 if (cellIds->GetId(0) == nextCell)
1143 {
1144 nextCell = cellIds->GetId(1);
1145 }
1146 else
1147 {
1148 nextCell = cellIds->GetId(0);
1149 }
1150 }
1151
1152 IntersectionLines->GetCellPoints(nextCell, pointIds);
1153 if (pointIds->GetNumberOfIds() > 2)
1154 {
1155 vtkDebugWithObjectMacro(
1156 this->ParentFilter, << "Number Of Points is greater than 2 for cell " << nextCell);
1157 }
1158 else if (pointIds->GetNumberOfIds() < 2)
1159 {
1160 vtkDebugWithObjectMacro(
1161 this->ParentFilter, << "Number Of Points is less than 2 for first cell " << nextCell);
1162 }
1163 prevPt = nextPt;
1164 if (pointIds->GetId(0) == nextPt)
1165 {
1166 nextPt = pointIds->GetId(1);
1167 }
1168 else
1169 {
1170 nextPt = pointIds->GetId(0);
1171 }
1172 usedPt[nextPt] = true;
1173
1174 simLine newestline;
1175 newestline.pt1 = prevPt;
1176 newestline.pt2 = nextPt;
1177 newestline.id = nextCell;
1178 loop->cells.push_back(newestline);
1179 }
1180 loop->endPt = nextPt;
1181 loop->loopType = 0;
1182 vtkDebugWithObjectMacro(this->ParentFilter, << "Start and End Point are " << nextPt);
1183
1184 return -1;
1185 }
1186
1187 // Tests an orientation in a specified region
RunLoopTest(vtkIdType interPt,vtkIdType nextCell,simLoop * loop,bool * usedPt)1188 int vtkLoopBooleanPolyDataFilter::Impl::RunLoopTest(
1189 vtkIdType interPt, vtkIdType nextCell, simLoop* loop, bool* usedPt)
1190 {
1191 // This test is only if the intersection has soft closed loops
1192 vtkDebugWithObjectMacro(this->ParentFilter, << "Running Loop Test to find right loop");
1193 vtkIdType stopCell = nextCell;
1194 vtkIdType prevPt = interPt;
1195 vtkIdType nextPt = interPt;
1196 vtkSmartPointer<vtkIdList> pointIds = vtkSmartPointer<vtkIdList>::New();
1197 vtkSmartPointer<vtkPolyData> tmpPolyData = vtkSmartPointer<vtkPolyData>::New();
1198 int input = 0;
1199 tmpPolyData->DeepCopy(this->Mesh[input]);
1200 tmpPolyData->BuildLinks();
1201 std::list<simLine>::iterator cellit;
1202
1203 vtkSmartPointer<vtkIdList> cellIds = vtkSmartPointer<vtkIdList>::New();
1204 IntersectionLines->GetPointCells(nextPt, cellIds);
1205 vtkDebugWithObjectMacro(this->ParentFilter,
1206 << "Number of cells should be more than two!! " << cellIds->GetNumberOfIds());
1207 for (vtkIdType i = 0; i < cellIds->GetNumberOfIds(); i++)
1208 {
1209 int numRegionsFound = 0;
1210 vtkIdType cellId = cellIds->GetId(i);
1211 vtkDebugWithObjectMacro(this->ParentFilter, << "Testing cell " << cellId);
1212 IntersectionLines->GetCellPoints(cellId, pointIds);
1213 if (pointIds->GetNumberOfIds() > 2)
1214 {
1215 vtkDebugWithObjectMacro(
1216 this->ParentFilter, << "Number Of Points is greater than 2 for first cell " << nextCell);
1217 }
1218 else if (pointIds->GetNumberOfIds() < 2)
1219 {
1220 vtkDebugWithObjectMacro(
1221 this->ParentFilter, << "Number Of Points is less than 2 for first cell " << nextCell);
1222 }
1223
1224 if (pointIds->GetId(0) == interPt)
1225 {
1226 nextPt = pointIds->GetId(1);
1227 }
1228 else
1229 {
1230 nextPt = pointIds->GetId(0);
1231 }
1232
1233 if (usedPt[nextPt] == true)
1234 {
1235 vtkDebugWithObjectMacro(this->ParentFilter, << "Bad One");
1236 }
1237 if (cellId != stopCell && usedPt[nextPt] != true)
1238 {
1239 simLine newline;
1240 newline.id = cellId;
1241 newline.pt1 = prevPt;
1242 newline.pt2 = nextPt;
1243 loop->cells.push_back(newline);
1244 vtkDebugWithObjectMacro(this->ParentFilter, << "Cell id is: " << cellId);
1245 for (cellit = loop->cells.begin(); cellit != loop->cells.end(); ++cellit)
1246 {
1247 simLine nextLine;
1248 nextLine = *cellit;
1249 nextCell = nextLine.id;
1250 vtkIdType p1 = nextLine.pt1;
1251 vtkIdType p2 = nextLine.pt2;
1252 vtkDebugWithObjectMacro(this->ParentFilter, << "Line cell is " << nextCell);
1253 vtkIdType outputCellId0 = this->NewCellIds[input]->GetComponent(nextCell, 0);
1254 vtkIdType outputCellId1 = this->NewCellIds[input]->GetComponent(nextCell, 1);
1255 if (outputCellId0 != -1)
1256 {
1257 if (this->CheckedCarefully[input][outputCellId0] == 0)
1258 {
1259 int sign1 = this->GetCellOrientation(tmpPolyData, outputCellId0, p1, p2, input);
1260 if (sign1 == -1)
1261 {
1262 numRegionsFound++;
1263 this->CheckCells->InsertNextId(outputCellId0);
1264 this->FindRegion(input, sign1, 1, 0);
1265 this->CheckCells->Reset();
1266 this->CheckCells2->Reset();
1267 this->CheckCellsCareful->Reset();
1268 this->CheckCellsCareful2->Reset();
1269 }
1270 }
1271 }
1272 if (outputCellId1 != -1)
1273 {
1274 if (this->CheckedCarefully[input][outputCellId1] == 0)
1275 {
1276 int sign2 = this->GetCellOrientation(tmpPolyData, outputCellId1, p1, p2, input);
1277 if (sign2 == -1)
1278 {
1279 numRegionsFound++;
1280 this->CheckCells->InsertNextId(outputCellId1);
1281 this->FindRegion(input, sign2, 1, 0);
1282 this->CheckCells->Reset();
1283 this->CheckCells2->Reset();
1284 this->CheckCellsCareful->Reset();
1285 this->CheckCellsCareful2->Reset();
1286 }
1287 }
1288 }
1289 }
1290 loop->cells.pop_back();
1291 this->ResetCheckArrays();
1292 vtkDebugWithObjectMacro(
1293 this->ParentFilter, << "Number of Regions Found: " << numRegionsFound);
1294 if (numRegionsFound == 1)
1295 {
1296 vtkDebugWithObjectMacro(this->ParentFilter, << "Legitimate Loop found");
1297 return cellId;
1298 }
1299 }
1300 }
1301 vtkDebugWithObjectMacro(this->ParentFilter, << "Start and End Point are " << nextPt);
1302
1303 return -1;
1304 }
1305
1306 // Combine the correct regions for output boolean
PerformBoolean(vtkPolyData * output,int booleanOperation)1307 void vtkLoopBooleanPolyDataFilter::Impl::PerformBoolean(vtkPolyData* output, int booleanOperation)
1308 {
1309 // vtkSmartPointer<vtkThreshold> thresholder =
1310 // vtkSmartPointer<vtkThreshold>::New();
1311 // vtkSmartPointer<vtkDataSetSurfaceFilter> surfacer =
1312 // vtkSmartPointer<vtkDataSetSurfaceFilter>::New();
1313
1314 vtkPolyData* surfaces[4];
1315 for (int i = 0; i < 4; i++)
1316 {
1317 surfaces[i] = vtkPolyData::New();
1318 }
1319
1320 this->ThresholdRegions(surfaces);
1321 // thresholder->SetInputData(this->Mesh[0]);
1322 // thresholder->SetInputArrayToProcess(0, 0, 0, 1, "BooleanRegion");
1323 // thresholder->SetThresholdFunction(vtkThreshold::THRESHOLD_BETWEEN);
1324 // thresholder->SetLowerThreshold(-1.0);
1325 // thresholder->SetUpperThreshold(-1.0);
1326 // thresholder->Update();
1327 // surfacer->SetInputData(thresholder->GetOutput());
1328 // surfacer->Update();
1329 // surfaces[0]->DeepCopy(surfacer->GetOutput());
1330
1331 // thresholder->SetInputData(this->Mesh[0]);
1332 // thresholder->SetInputArrayToProcess(0, 0, 0, 1, "BooleanRegion");
1333 // thresholder->SetThresholdFunction(vtkThreshold::THRESHOLD_BETWEEN);
1334 // thresholder->SetLowerThreshold(1.0);
1335 // thresholder->SetUpperThreshold(1.0);
1336 // thresholder->Update();
1337 // surfacer->SetInputData(thresholder->GetOutput());
1338 // surfacer->Update();
1339 // surfaces[1]->DeepCopy(surfacer->GetOutput());
1340
1341 // thresholder->SetInputData(this->Mesh[1]);
1342 // thresholder->SetInputArrayToProcess(0, 0, 0, 1, "BooleanRegion");
1343 // thresholder->SetThresholdFunction(vtkThreshold::THRESHOLD_BETWEEN);
1344 // thresholder->SetLowerThreshold(1.0);
1345 // thresholder->SetUpperThreshold(1.0);
1346 // thresholder->Update();
1347 // surfacer->SetInputData(thresholder->GetOutput());
1348 // surfacer->Update();
1349 // surfaces[2]->DeepCopy(surfacer->GetOutput());
1350
1351 // thresholder->SetInputData(this->Mesh[1]);
1352 // thresholder->SetInputArrayToProcess(0, 0, 0, 1, "BooleanRegion");
1353 // thresholder->SetThresholdFunction(vtkThreshold::THRESHOLD_BETWEEN);
1354 // thresholder->SetLowerThreshold(-1.0);
1355 // thresholder->SetUpperThreshold(-1.0);
1356 // thresholder->Update();
1357 // surfacer->SetInputData(thresholder->GetOutput());
1358 // surfacer->Update();
1359 // surfaces[3]->DeepCopy(surfacer->GetOutput());
1360
1361 vtkSmartPointer<vtkAppendPolyData> appender = vtkSmartPointer<vtkAppendPolyData>::New();
1362
1363 // If open intersection case, make sure correct region is taken
1364 if (this->IntersectionCase == 2)
1365 {
1366 vtkSmartPointer<vtkPolyData> tmp = vtkSmartPointer<vtkPolyData>::New();
1367 int numCells[4];
1368 std::list<vtkIdType> nocellregion;
1369 for (int i = 0; i < 4; i++)
1370 {
1371 numCells[i] = surfaces[i]->GetNumberOfCells();
1372 if (numCells[i] == 0)
1373 {
1374 nocellregion.push_back(i);
1375 }
1376 }
1377 if (!nocellregion.empty())
1378 {
1379 if (nocellregion.front() == 0)
1380 {
1381 tmp->DeepCopy(surfaces[1]);
1382 surfaces[1]->DeepCopy(surfaces[0]);
1383 surfaces[0]->DeepCopy(tmp);
1384 }
1385 if (nocellregion.back() == 2)
1386 {
1387 tmp->DeepCopy(surfaces[3]);
1388 surfaces[3]->DeepCopy(surfaces[2]);
1389 surfaces[2]->DeepCopy(tmp);
1390 }
1391 }
1392 }
1393 if (booleanOperation == 0)
1394 {
1395 appender->AddInputData(surfaces[0]);
1396 appender->AddInputData(surfaces[2]);
1397 }
1398 if (booleanOperation == 1)
1399 {
1400 appender->AddInputData(surfaces[1]);
1401 appender->AddInputData(surfaces[3]);
1402 }
1403 if (booleanOperation == 2)
1404 {
1405 appender->AddInputData(surfaces[0]);
1406 appender->AddInputData(surfaces[3]);
1407 }
1408 appender->Update();
1409
1410 output->DeepCopy(appender->GetOutput());
1411
1412 for (int i = 0; i < 4; i++)
1413 {
1414 surfaces[i]->Delete();
1415 }
1416 }
1417
ThresholdRegions(vtkPolyData ** surfaces)1418 void vtkLoopBooleanPolyDataFilter::Impl::ThresholdRegions(vtkPolyData** surfaces)
1419 {
1420 vtkPoints* points[4];
1421 vtkCellArray* cells[4];
1422 vtkIntArray* boundaryPoints[4];
1423 vtkIntArray* boundaryCells[4];
1424 vtkIntArray* booleanCells[4];
1425
1426 for (int i = 0; i < 4; i++)
1427 {
1428 points[i] = vtkPoints::New();
1429 cells[i] = vtkCellArray::New();
1430 boundaryPoints[i] = vtkIntArray::New();
1431 boundaryCells[i] = vtkIntArray::New();
1432 booleanCells[i] = vtkIntArray::New();
1433 }
1434
1435 for (int i = 0; i < 2; i++)
1436 {
1437 int numCells = this->Mesh[i]->GetNumberOfCells();
1438 for (int j = 0; j < numCells; j++)
1439 {
1440 int value = this->BooleanArray[i]->GetValue(j);
1441 vtkIdType npts;
1442 const vtkIdType* pts;
1443 this->Mesh[i]->GetCellPoints(j, npts, pts);
1444 if (value < 0)
1445 {
1446 vtkSmartPointer<vtkIdList> newPointIds = vtkSmartPointer<vtkIdList>::New();
1447 newPointIds->SetNumberOfIds(3);
1448 for (int k = 0; k < npts; k++)
1449 {
1450 double pt[3];
1451 this->Mesh[i]->GetPoint(pts[k], pt);
1452 vtkIdType newId = points[3 * i]->InsertNextPoint(pt);
1453 newPointIds->SetId(k, newId);
1454 boundaryPoints[3 * i]->InsertValue(newId, this->BoundaryPointArray[i]->GetValue(pts[k]));
1455 }
1456 vtkIdType cellId = cells[3 * i]->InsertNextCell(newPointIds);
1457 boundaryCells[3 * i]->InsertValue(cellId, this->BoundaryCellArray[i]->GetValue(j));
1458 booleanCells[3 * i]->InsertValue(cellId, this->BooleanArray[i]->GetValue(j));
1459 }
1460 if (value > 0)
1461 {
1462 vtkSmartPointer<vtkIdList> newPointIds = vtkSmartPointer<vtkIdList>::New();
1463 newPointIds->SetNumberOfIds(3);
1464 for (int k = 0; k < npts; k++)
1465 {
1466 double pt[3];
1467 this->Mesh[i]->GetPoint(pts[k], pt);
1468 vtkIdType newId = points[i + 1]->InsertNextPoint(pt);
1469 newPointIds->SetId(k, newId);
1470 boundaryPoints[i + 1]->InsertValue(newId, this->BoundaryPointArray[i]->GetValue(pts[k]));
1471 }
1472 vtkIdType cellId = cells[i + 1]->InsertNextCell(newPointIds);
1473 boundaryCells[i + 1]->InsertValue(cellId, this->BoundaryCellArray[i]->GetValue(j));
1474 booleanCells[i + 1]->InsertValue(cellId, this->BooleanArray[i]->GetValue(j));
1475 }
1476 }
1477 }
1478 for (int i = 0; i < 4; i++)
1479 {
1480 surfaces[i]->SetPoints(points[i]);
1481 surfaces[i]->SetPolys(cells[i]);
1482 surfaces[i]->BuildLinks();
1483 boundaryPoints[i]->SetName("BoundaryPoints");
1484 surfaces[i]->GetPointData()->AddArray(boundaryPoints[i]);
1485 boundaryCells[i]->SetName("BoundaryCells");
1486 surfaces[i]->GetCellData()->AddArray(boundaryCells[i]);
1487 booleanCells[i]->SetName("BooleanRegion");
1488 surfaces[i]->GetCellData()->AddArray(booleanCells[i]);
1489
1490 points[i]->Delete();
1491 cells[i]->Delete();
1492 boundaryPoints[i]->Delete();
1493 boundaryCells[i]->Delete();
1494 booleanCells[i]->Delete();
1495 }
1496 }
1497