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