vtkIntersectionPolyDataFilter.cxx 86.8 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkIntersectionPolyDataFilter.cxx

  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notice for more information.

=========================================================================*/
15
16
17
18
19
20
/** @file vtkIntersectionPolyDataFilter.cxx
 *  @brief This is a filter that performs the intersection of surfaces
 *  @author Adam Updegrove
 *  @author updega2@gmail.com
 */

21
22
23
24
#include "vtkIntersectionPolyDataFilter.h"

#include "vtkCellArray.h"
#include "vtkCellData.h"
25
#include "vtkCleanPolyData.h"
26
27
#include "vtkDelaunay2D.h"
#include "vtkDoubleArray.h"
28
#include "vtkIdList.h"
29
30
31
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkLine.h"
32
#include "vtkLongArray.h"
33
#include "vtkOBBTree.h"
34
#include "vtkObjectFactory.h"
35
36
37
#include "vtkPlane.h"
#include "vtkPointData.h"
#include "vtkPointLocator.h"
38
39
40
#include "vtkPoints.h"
#include "vtkPolyDataNormals.h"
#include "vtkPolygon.h"
41
42
43
#include "vtkSmartPointer.h"
#include "vtkSortDataArray.h"
#include "vtkTriangle.h"
44
45
46
47
#include "vtkTriangleFilter.h"
#include "vtkTransform.h"
#include "vtkTransformPolyDataFilter.h"
#include "vtkUnstructuredGrid.h"
48
49

#include <map>
50
51
52
#include <sstream>
#include <iostream>
#include <list>
53
54

//----------------------------------------------------------------------------
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
// Helper typedefs and data structures.
namespace {

struct simPoint
{
  vtkIdType id;
  double pt[3];
  int concave;
};

struct simPolygon
{
  std::list<simPoint> points;
  int orientation;
};

}

73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
typedef std::multimap< vtkIdType, vtkIdType >    IntersectionMapType;
typedef IntersectionMapType::iterator            IntersectionMapIteratorType;

//typedef std::pair< vtkIdType, vtkIdType >            CellEdgePairType;
typedef struct _CellEdgeLine {
  vtkIdType CellId;
  vtkIdType EdgeId;
  vtkIdType LineId;
} CellEdgeLineType;

typedef std::multimap< vtkIdType, CellEdgeLineType > PointEdgeMapType;
typedef PointEdgeMapType::iterator                   PointEdgeMapIteratorType;


//----------------------------------------------------------------------------
// Private implementation to hide STL.
//----------------------------------------------------------------------------
class vtkIntersectionPolyDataFilter::Impl
{
public:
  Impl();
  virtual ~Impl();

96
  //Finds all triangle triangle intersections between two input OOBTrees
97
98
99
  static int FindTriangleIntersections(vtkOBBNode *node0, vtkOBBNode *node1,
                                       vtkMatrix4x4 *transform, void *arg);

100
  //Runs the split mesh for the designated input surface
101
102
103
104
105
  int SplitMesh(int inputIndex, vtkPolyData *output,
                vtkPolyData *intersectionLines);

protected:

106
  //Split cells into polygons created by intersection lines
107
  vtkCellArray* SplitCell(vtkPolyData *input, vtkIdType cellId,
108
109
110
111
                                     vtkIdType *cellPts,
                                     IntersectionMapType *map,
                                     vtkPolyData *interLines, int inputIndex,
                                     int numCurrCells);
112

113
114
  //Function to add point to check edge list for remeshing step
  int AddToPointEdgeMap(int index, vtkIdType ptId, double x[3],
115
116
117
118
                         vtkPolyData *mesh, vtkIdType cellId,
                         vtkIdType edgeId, vtkIdType lineId,
                         vtkIdType triPts[3]);

119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
  //Function to add information about the new cell data
  void AddToNewCellMap(int inputIndex, int interPtCount, int interPts[3],
                       vtkPolyData *interLines,int numCurrCells);

  //Function inside SplitCell to get the smaller triangle loops
  int GetLoops(vtkPolyData *pd, std::vector<simPolygon> *loops);

  //Get individual polygon loop of splitting cell
  int GetSingleLoop(vtkPolyData *pd,simPolygon *loop, vtkIdType nextCell,
                    bool *interPtBool, bool *lineBool);

  //Follow a loop orienation to iterate around a split polygon
  int FollowLoopOrientation(vtkPolyData *pd, simPolygon *loop,
                            vtkIdType *nextCell,
                            vtkIdType nextPt, vtkIdType prevPt,
                            vtkIdList *pointCells);

  //Set the loop orientation based on CW CCW geometric test
  void SetLoopOrientation(vtkPolyData *pd, simPolygon *loop,
                          vtkIdType *nextCell, vtkIdType nextPt,
                          vtkIdType prevPt, vtkIdList *pointCells);

  //Get the loop orienation is already given
  int GetLoopOrientation(vtkPolyData *pd, vtkIdType cell, vtkIdType ptId1,
      vtkIdType ptId2);

  //Orient the triangle based on the transform for remeshing
  void Orient(vtkPolyData *pd, vtkTransform *transform, vtkPolyData *boundary,
                  vtkPolygon *boundarypoly);

  //Checks to make sure multiple lines are not added to the same triangle
  //that needs to re-triangulated
  int CheckLine(vtkPolyData *pd, vtkIdType ptId1, vtkIdType ptId2);

  //Gets a transform to the XY plane for three points comprising a triangle
  int GetTransform(vtkTransform *transform, vtkPoints *points);
155
156
157
158
159
160
161
162

public:
  vtkPolyData         *Mesh[2];
  vtkOBBTree          *OBBTree1;

  // Stores the intersection lines.
  vtkCellArray        *IntersectionLines;

163
164
165
  vtkIdTypeArray      *SurfaceId;
  vtkIdTypeArray      *NewCellIds[2];

166
167
168
169
  // Cell data that indicates in which cell each intersection
  // lies. One array for each output surface.
  vtkIdTypeArray      *CellIds[2];

170
171
  // Cell data that indicates on which surface the intersection point lies.

172
173
174
175
176
  // Map from points to the cells that contain them. Used for point
  // data interpolation. For points on the edge between two cells, it
  // does not matter which cell is recorded bcause the interpolation
  // will be the same.  One array for each output surface.
  vtkIdTypeArray      *PointCellIds[2];
177
  vtkIntArray         *BoundaryPoints[2];
178
179
180
181
182
183
184

  // Merging filter used to convert intersection lines from "line
  // soup" to connected polylines.
  vtkPointLocator     *PointMerger;

  // Map from cell ID to intersection line.
  IntersectionMapType *IntersectionMap[2];
185
186
  IntersectionMapType *IntersectionPtsMap[2];
  IntersectionMapType *PointMapper;
187
188
189
190
191

  // Map from point to an edge on which it resides, the ID of the
  // cell, and the ID of the line.
  PointEdgeMapType    *PointEdgeMap[2];

192
193
194
195
196
197
198
199
200
  //vtkPolyData to hold current splitting cell. Used to double check area
  //of small area cells
  vtkPolyData *SplittingPD;
  int         TransformSign;
  double      Tol;

  // Pointer to overarching filter
  vtkIntersectionPolyDataFilter *ParentFilter;

201
202
203
204
205
206
207
208
protected:
  Impl(const Impl&); // purposely not implemented
  void operator=(const Impl&); // purposely not implemented

};

//----------------------------------------------------------------------------
vtkIntersectionPolyDataFilter::Impl::Impl() :
209
  OBBTree1(0), IntersectionLines(0), SurfaceId(0), PointMerger(0)
210
211
{
  for (int i = 0; i < 2; i++)
212
    {
213
214
215
216
217
    this->Mesh[i]                 = NULL;
    this->CellIds[i]              = NULL;
    this->IntersectionMap[i]      = new IntersectionMapType();
    this->IntersectionPtsMap[i]   = new IntersectionMapType();
    this->PointEdgeMap[i]         = new PointEdgeMapType();
218
    }
219
220
221
222
  this->PointMapper             = new IntersectionMapType();
  this->SplittingPD             = vtkPolyData::New();
  this->TransformSign = 0;
  this->Tol = 1e-6;
223
224
225
226
227
228
}

//----------------------------------------------------------------------------
vtkIntersectionPolyDataFilter::Impl::~Impl()
{
  for (int i = 0; i < 2; i++)
229
    {
230
    delete this->IntersectionMap[i];
231
    delete this->IntersectionPtsMap[i];
232
    delete this->PointEdgeMap[i];
233
    }
234
235
  delete this->PointMapper;
  this->SplittingPD->Delete();
236
237
238
239
240
241
242
243
244
245
246
}


//----------------------------------------------------------------------------
int vtkIntersectionPolyDataFilter::Impl
::FindTriangleIntersections(vtkOBBNode *node0, vtkOBBNode *node1,
                            vtkMatrix4x4 *transform, void *arg)
{
  vtkIntersectionPolyDataFilter::Impl *info =
    reinterpret_cast<vtkIntersectionPolyDataFilter::Impl*>(arg);

247
248
249
250
251
252
253
254
255
256
257
258
  //Set up local structures to hold Impl array information
  vtkPolyData     *mesh0                 = info->Mesh[0];
  vtkPolyData     *mesh1                 = info->Mesh[1];
  vtkOBBTree      *obbTree1              = info->OBBTree1;
  vtkCellArray    *intersectionLines     = info->IntersectionLines;
  vtkIdTypeArray  *intersectionSurfaceId = info->SurfaceId;
  vtkIdTypeArray  *intersectionCellIds0  = info->CellIds[0];
  vtkIdTypeArray  *intersectionCellIds1  = info->CellIds[1];
  vtkPointLocator *pointMerger           = info->PointMerger;
  double tolerance                       = info->Tol;

  //The number of cells in OBBTree
259
260
261
262
263
264
265
  int numCells0 = node0->Cells->GetNumberOfIds();

  for (vtkIdType id0 = 0; id0 < numCells0; id0++)
    {
    vtkIdType cellId0 = node0->Cells->GetId(id0);
    int type0 = mesh0->GetCellType(cellId0);

266
    //Make sure the cell is a triangle
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
    if (type0 == VTK_TRIANGLE)
      {
      vtkIdType npts0, *triPtIds0;
      mesh0->GetCellPoints(cellId0, npts0, triPtIds0);
      double triPts0[3][3];
      for (vtkIdType id = 0; id < npts0; id++)
        {
        mesh0->GetPoint(triPtIds0[id], triPts0[id]);
        }

      if (obbTree1->TriangleIntersectsNode
          (node1, triPts0[0], triPts0[1], triPts0[2], transform))
        {
        int numCells1 = node1->Cells->GetNumberOfIds();
        for (vtkIdType id1 = 0; id1 < numCells1; id1++)
          {
          vtkIdType cellId1 = node1->Cells->GetId(id1);
          int type1 = mesh1->GetCellType(cellId1);
          if (type1 == VTK_TRIANGLE)
            {
            // See if the two cells actually intersect. If they do,
            // add an entry into the intersection maps and add an
            // intersection line.
            vtkIdType npts1, *triPtIds1;
            mesh1->GetCellPoints(cellId1, npts1, triPtIds1);

            double triPts1[3][3];
            for (vtkIdType id = 0; id < npts1; id++)
              {
              mesh1->GetPoint(triPtIds1[id], triPts1[id]);
              }

            int coplanar = 0;
            double outpt0[3], outpt1[3];
301
            double surfaceid[2];
302
303
304
305
            int intersects =
              vtkIntersectionPolyDataFilter::TriangleTriangleIntersection
              (triPts0[0], triPts0[1], triPts0[2],
               triPts1[0], triPts1[1], triPts1[2],
306
               coplanar, outpt0, outpt1, surfaceid, tolerance);
307

308
            if (coplanar)
309
310
              {
              // Coplanar triangle intersection is not handled.
311
312
313
              // This intersection will not be included in the output. TODO
              //vtkDebugMacro(<<"Coplanar");
              intersects = 0;
314
315
316
              continue;
              }

317
318
319
            //If actual intersection, add point and cell to edge, line,
            //and surface maps!
            if (intersects)
320
321
322
323
              {
              vtkIdType lineId = intersectionLines->GetNumberOfCells();

              vtkIdType ptId0, ptId1;
324
325
326
327
328
329
330
331
332
333
334
335
              int unique[2];
              unique[0] = pointMerger->InsertUniquePoint(outpt0, ptId0);
              unique[1] = pointMerger->InsertUniquePoint(outpt1, ptId1);

              int addline = 1;
              if (ptId0 == ptId1)
                {
                addline = 0;
                }

              if (ptId0 == ptId1 && surfaceid[0] != surfaceid[1])
                {
Adam Updegrove's avatar
Adam Updegrove committed
336
                intersectionSurfaceId->InsertValue(ptId0, 3);
337
338
339
340
341
                }
              else
                {
                if (unique[0])
                {
Adam Updegrove's avatar
Adam Updegrove committed
342
                intersectionSurfaceId->InsertValue(ptId0, surfaceid[0]);
343
344
345
346
347
                }
              else
                {
                if (intersectionSurfaceId->GetValue(ptId0) != 3)
                  {
Adam Updegrove's avatar
Adam Updegrove committed
348
                  intersectionSurfaceId->InsertValue(ptId0, surfaceid[0]);
349
350
351
352
353
                  }

                }
              if (unique[1])
                {
Adam Updegrove's avatar
Adam Updegrove committed
354
                intersectionSurfaceId->InsertValue(ptId1, surfaceid[1]);
355
356
357
358
359
                }
              else
                {
                if (intersectionSurfaceId->GetValue(ptId1) != 3)
                  {
Adam Updegrove's avatar
Adam Updegrove committed
360
                  intersectionSurfaceId->InsertValue(ptId1, surfaceid[1]);
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
                  }
                }
              }

            info->IntersectionPtsMap[0]->
              insert(std::make_pair(ptId0, cellId0));
            info->IntersectionPtsMap[1]->
              insert(std::make_pair(ptId0, cellId1));
            info->IntersectionPtsMap[0]->
              insert(std::make_pair(ptId1, cellId0));
            info->IntersectionPtsMap[1]->
              insert(std::make_pair(ptId1, cellId1));

            //Check to see if duplicate line. Line can only be a duplicate
            //line if both points are not unique and they don't
            //equal eachother
            if (!unique[0] && !unique[1] && ptId0 != ptId1)
              {
              vtkSmartPointer<vtkPolyData> lineTest =
                vtkSmartPointer<vtkPolyData>::New();
              lineTest->SetPoints(pointMerger->GetPoints());
              lineTest->SetLines(intersectionLines);
              lineTest->BuildLinks();
Adam Updegrove's avatar
Adam Updegrove committed
384
              int newLine = info->CheckLine(lineTest, ptId0, ptId1);
385
386
387
388
389
390
391
392
393
394
395
              if (newLine == 0)
                {
                addline = 0;
                }
              }
            if (addline)
              {
              //If the line is new and does not consist of two identical
              //points, add the line to the intersection and update
              //mapping information
              intersectionLines->InsertNextCell(2);
396
397
398
399
400
401
              intersectionLines->InsertCellPoint(ptId0);
              intersectionLines->InsertCellPoint(ptId1);

              intersectionCellIds0->InsertNextValue(cellId0);
              intersectionCellIds1->InsertNextValue(cellId1);

402
403
404
405
              info->PointCellIds[0]->InsertValue(ptId0, cellId0);
              info->PointCellIds[0]->InsertValue(ptId1, cellId0);
              info->PointCellIds[1]->InsertValue(ptId0, cellId1);
              info->PointCellIds[1]->InsertValue(ptId1, cellId1);
406

407
408
409
410
              info->IntersectionMap[0]->
                insert(std::make_pair(cellId0, lineId));
              info->IntersectionMap[1]->
                insert(std::make_pair(cellId1, lineId));
411
412
413

              // Check which edges of cellId0 and cellId1 outpt0 and
              // outpt1 are on, if any.
414
              int isOnEdge=0;
Adam Updegrove's avatar
Adam Updegrove committed
415
              int m0p0=0, m0p1=0, m1p0=0, m1p1=0;
416
417
              for (vtkIdType edgeId = 0; edgeId < 3; edgeId++)
                {
418
                isOnEdge = info->AddToPointEdgeMap(0, ptId0, outpt0,
Adam Updegrove's avatar
Adam Updegrove committed
419
                    mesh0, cellId0, edgeId, lineId, triPtIds0);
420
421
422
423
424
                if (isOnEdge != -1)
                  {
                  m0p0++;
                  }
                isOnEdge = info->AddToPointEdgeMap(0, ptId1, outpt1,
Adam Updegrove's avatar
Adam Updegrove committed
425
                    mesh0, cellId0, edgeId, lineId, triPtIds0);
426
427
428
429
430
                if (isOnEdge != -1)
                  {
                  m0p1++;
                  }
                isOnEdge = info->AddToPointEdgeMap(1, ptId0, outpt0,
Adam Updegrove's avatar
Adam Updegrove committed
431
                    mesh1, cellId1, edgeId, lineId, triPtIds1);
432
433
434
435
436
                if (isOnEdge != -1)
                  {
                  m1p0++;
                  }
                isOnEdge = info->AddToPointEdgeMap(1, ptId1, outpt1,
Adam Updegrove's avatar
Adam Updegrove committed
437
                    mesh1, cellId1, edgeId, lineId, triPtIds1);
438
439
440
441
                if (isOnEdge != -1)
                  {
                  m1p1++;
                  }
442
                }
443
444
445
446
              //Special cases caught by tolerance and not from the Point
              //Merger
              if (m0p0 > 0 && m1p0 > 0)
                {
Adam Updegrove's avatar
Adam Updegrove committed
447
                intersectionSurfaceId->InsertValue(ptId0, 3);
448
449
450
                }
              if (m0p1 > 0 && m1p1 > 0)
                {
Adam Updegrove's avatar
Adam Updegrove committed
451
                intersectionSurfaceId->InsertValue(ptId1, 3);
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
                }
              }
            //Add information about origin surface to std::maps for
            //checks later
            if (intersectionSurfaceId->GetValue(ptId0) == 1)
              {
              info->IntersectionPtsMap[0]->
                insert(std::make_pair(ptId0, cellId0));
              }
            else if (intersectionSurfaceId->GetValue(ptId0) == 2)
              {
              info->IntersectionPtsMap[1]->
                insert(std::make_pair(ptId0, cellId1));
              }
            else
              {
              info->IntersectionPtsMap[0]->
                insert(std::make_pair(ptId0, cellId0));
              info->IntersectionPtsMap[1]->
                insert(std::make_pair(ptId0, cellId1));
              }
            if (intersectionSurfaceId->GetValue(ptId1) == 1)
              {
              info->IntersectionPtsMap[0]->
                insert(std::make_pair(ptId1, cellId0));
              }
            else if (intersectionSurfaceId->GetValue(ptId1) == 2)
              {
              info->IntersectionPtsMap[1]->
                insert(std::make_pair(ptId1, cellId1));
              }
            else
              {
              info->IntersectionPtsMap[0]->
                insert(std::make_pair(ptId1, cellId0));
              info->IntersectionPtsMap[1]->
                insert(std::make_pair(ptId1, cellId1));
489
490
              }
            }
491
            }
492
493
494
495
          }
        }
      }
    }
496
  return 1;
497
498
499
500
501
}


//----------------------------------------------------------------------------
int vtkIntersectionPolyDataFilter::Impl
Adam Updegrove's avatar
Adam Updegrove committed
502
::SplitMesh(int inputIndex, vtkPolyData *output, vtkPolyData *intersectionLines)
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
{
  vtkPolyData *input = this->Mesh[inputIndex];
  IntersectionMapType *intersectionMap = this->IntersectionMap[inputIndex];
  vtkCellData *inCD  = input->GetCellData();
  vtkCellData *outCD = output->GetCellData();
  vtkIdType numCells = input->GetNumberOfCells();
  vtkIdType cellIdX = 0;

  //
  // Process points
  //
  vtkIdType inputNumPoints = input->GetPoints()->GetNumberOfPoints();
  vtkSmartPointer< vtkPoints > points = vtkSmartPointer< vtkPoints >::New();
  points->Allocate(100);
  output->SetPoints(points);

  //
  // Split intersection lines. The lines structure is constructed
  // using a vtkPointLocator. However, some lines may have an endpoint
  // on a cell edge that has no neighbor. We need to duplicate a line
  // point in such a case and update the point ID in the line cell.
  //
  vtkSmartPointer< vtkPolyData > splitLines =
    vtkSmartPointer <vtkPolyData >::New();
  splitLines->DeepCopy(intersectionLines);

  vtkPointData *inPD  = input->GetPointData();
  vtkPointData *outPD = output->GetPointData();
531
  outPD->CopyAllocate(inPD, input->GetNumberOfPoints());
532
533
534
535
536
537
538
539

  // Copy over the point data from the input
  for (vtkIdType ptId = 0; ptId < inputNumPoints; ptId++)
    {
    double pt[3];
    input->GetPoints()->GetPoint(ptId, pt);
    output->GetPoints()->InsertNextPoint(pt);
    outPD->CopyData(inPD, ptId, ptId);
Adam Updegrove's avatar
Adam Updegrove committed
540
    this->BoundaryPoints[inputIndex]->InsertValue(ptId, 0);
541
542
543
544
545
546
547
548
549
550
551
    }

  // Copy the points from splitLines to the output, interpolating the
  // data as we go.
  for (vtkIdType id = 0; id < splitLines->GetNumberOfPoints(); id++)
    {
    double pt[3];
    splitLines->GetPoint(id, pt);
    vtkIdType newPtId = output->GetPoints()->InsertNextPoint(pt);

    // Retrieve the cell ID from splitLines
552
    vtkIdType cellId = this->PointCellIds[inputIndex]->GetValue(id);
553
554
555
556
557
558
559

    double closestPt[3], pcoords[3], dist2, weights[3];
    int subId;
    vtkCell *cell = input->GetCell(cellId);
    cell->EvaluatePosition(pt, closestPt, subId, pcoords, dist2, weights);
    outPD->InterpolatePoint(input->GetPointData(), newPtId, cell->PointIds,
                            weights);
Adam Updegrove's avatar
Adam Updegrove committed
560
    this->BoundaryPoints[inputIndex]->InsertValue(newPtId, 0);
561
562
563
564
565
566
567
    }

  //
  // Process cells
  //
  outCD->CopyAllocate(inCD, numCells);

568
  if (input->GetPolys()->GetNumberOfCells() > 0)
569
570
    {
    vtkCellArray *cells = input->GetPolys();
571
    vtkIdType newId = output->GetNumberOfCells();
572
573
574
575

    vtkSmartPointer< vtkCellArray > newPolys =
      vtkSmartPointer< vtkCellArray >::New();

Adam Updegrove's avatar
Adam Updegrove committed
576
    newPolys->EstimateSize(cells->GetNumberOfCells(), 3);
577
578
579
580
581
582
583
584
585
586
    output->SetPolys(newPolys);

    vtkSmartPointer< vtkIdList > edgeNeighbors =
      vtkSmartPointer< vtkIdList >::New();
    vtkIdType nptsX = 0;
    vtkIdType *pts = 0;
    vtkSmartPointer< vtkIdList > cellsToCheck =
      vtkSmartPointer< vtkIdList >::New();
    for (cells->InitTraversal(); cells->GetNextCell(nptsX, pts); cellIdX++)
      {
587
      if (nptsX != 3)
588
        {
589
590
        vtkGenericWarningMacro(<< "vtkIntersectionPolyDataFilter only works"
                                << " with triangle meshes.");
591
592
593
594
595
596
597
598
599
600
601
602
        continue;
        }

      cellsToCheck->Reset();
      cellsToCheck->Allocate(nptsX+1);
      cellsToCheck->InsertNextId(cellIdX);

      // Collect the cells relevant for splitting this cell.  If the
      // cell is in the intersection map, split. If not, one of its
      // edges may be split by an intersection line that splits a
      // neighbor cell. Mark the cell as needing a split if this is
      // the case.
603
604
      bool needsSplit = intersectionMap->find(cellIdX)
        != intersectionMap->end();
605
606
607
608
609
610
611
612
613
614
615
      for (vtkIdType ptId = 0; ptId < nptsX; ptId++)
        {
        vtkIdType pt0Id = pts[ptId];
        vtkIdType pt1Id = pts[(ptId+1) % nptsX];
        edgeNeighbors->Reset();
        input->GetCellEdgeNeighbors(cellIdX, pt0Id, pt1Id, edgeNeighbors);
        for (vtkIdType nbr = 0; nbr < edgeNeighbors->GetNumberOfIds(); nbr++)
          {
          vtkIdType nbrId = edgeNeighbors->GetId(nbr);
          cellsToCheck->InsertNextId(nbrId);

616
          if (intersectionMap->find(nbrId) != intersectionMap->end())
617
618
619
620
621
622
623
            {
            needsSplit = true;
            }
          } // for (vtkIdType nbr = 0; ...
        } // for (vtkIdType pt = 0; ...

      // Splitting occurs here
624
      if (!needsSplit)
625
626
627
628
629
630
631
        {
        // Just insert the cell and copy the cell data
        newId = newPolys->InsertNextCell(3, pts);
        outCD->CopyData(inCD, cellIdX, newId);
        }
      else
        {
632
633
634
        //Total number of cells so that we know the id numbers of the new
        //cells added and we can add it to the new cell id mapping
        int numCurrCells = newPolys->GetNumberOfCells();
635
        vtkCellArray *splitCells = this->SplitCell
636
637
638
639
          (input, cellIdX, pts, intersectionMap, splitLines,
           inputIndex,numCurrCells);
        if (splitCells == NULL)
          {
Adam Updegrove's avatar
Adam Updegrove committed
640
          vtkDebugWithObjectMacro(this->ParentFilter, <<"Error in splitting cell!");
641
642
643
          return 0;
          }

644
645
646
647
648
649
650
651

        double pt0[3], pt1[3], pt2[3], normal[3];
        points->GetPoint(pts[0], pt0);
        points->GetPoint(pts[1], pt1);
        points->GetPoint(pts[2], pt2);
        vtkTriangle::ComputeNormal(pt0, pt1, pt2, normal);
        vtkMath::Normalize(normal);

652
653
654
        vtkIdType npts, *ptIds, subCellId;
        splitCells->InitTraversal();
        for (subCellId = 0;splitCells->GetNextCell(npts, ptIds);subCellId++)
655
656
657
658
659
660
661
662
663
664
          {
          // Check for reversed cells. I'm not sure why, but in some
          // cases, cells are reversed.
          double subCellNormal[3];
          points->GetPoint(ptIds[0], pt0);
          points->GetPoint(ptIds[1], pt1);
          points->GetPoint(ptIds[2], pt2);
          vtkTriangle::ComputeNormal(pt0, pt1, pt2, subCellNormal);
          vtkMath::Normalize(subCellNormal);

665
          if (vtkMath::Dot(normal, subCellNormal) > 0)
666
667
668
669
670
671
            {
            newId = newPolys->InsertNextCell(npts, ptIds);
            }
          else
            {
            newId = newPolys->InsertNextCell(npts);
672
            for (int i = 0; i < npts; i++)
673
              {
674
              newPolys->InsertCellPoint(ptIds[ npts-i-1 ]);
675
676
677
678
679
680
681
682
              }
            }

          outCD->CopyData(inCD, cellIdX, newId); // Duplicate cell data
          }
        splitCells->Delete();
        }
      } // for (cells->InitTraversal(); ...
683
    } //if inputGetPolys()->GetNumberOfCells() > 1 ...
684
685
686
687
688
689

  return 1;
}

vtkCellArray* vtkIntersectionPolyDataFilter::Impl
::SplitCell(vtkPolyData *input, vtkIdType cellId, vtkIdType *cellPts,
690
691
692
            IntersectionMapType *map,
            vtkPolyData *interLines, int inputIndex,
            int numCurrCells)
693
{
694
695
696
697
698
699
700
701
702
703
704
705
706
  // Copy down the SurfaceID array that tells which surface the point belongs
  // to
  vtkIdTypeArray *surfaceMapper;
  surfaceMapper = vtkIdTypeArray::SafeDownCast(
    interLines->GetPointData()->GetArray("SurfaceID"));

  //Array to keep track of which points are on the boundary of the cell
  vtkSmartPointer<vtkIdTypeArray> cellBoundaryPt =
    vtkSmartPointer<vtkIdTypeArray>::New();
  //Array to tell wheter the original cell points lie on the intersecting
  //line
  int CellPointOnInterLine[3] = {0,0,0};

707
708
709
710
  // Gather points from the cell
  vtkSmartPointer< vtkPoints > points = vtkSmartPointer< vtkPoints >::New();
  vtkSmartPointer< vtkPointLocator > merger =
    vtkSmartPointer< vtkPointLocator >::New();
711
712
  merger->SetTolerance(this->Tol);
  merger->InitPointInsertion(points, input->GetBounds());
713
714

  double xyz[3];
715
  for (int i = 0; i < 3; i++)
716
    {
717
    if (cellPts[i] >= input->GetNumberOfPoints())
718
719
720
721
      {
      vtkGenericWarningMacro(<< "invalid point read 1");
      }
    input->GetPoint(cellPts[i], xyz);
722
723
    merger->InsertNextPoint(xyz);
    cellBoundaryPt->InsertNextValue(1);
724
725
    }

726
727
  // Set up line cells and array to track the just the intersecting lines
  // on the cell.
728
729
  vtkSmartPointer< vtkCellArray > lines =
    vtkSmartPointer< vtkCellArray >::New();
730
731
  vtkSmartPointer< vtkCellArray > interceptlines =
    vtkSmartPointer< vtkCellArray >::New();
732
733
734
735
736
737
738
739
740
741
742

  double p0[3], p1[3], p2[3];
  input->GetPoint(cellPts[0], p0);
  input->GetPoint(cellPts[1], p1);
  input->GetPoint(cellPts[2], p2);

  // This maps the point IDs for the vtkPolyData passed to
  // vtkDelaunay2D back to the original IDs in interLines. NOTE: The
  // point IDs from the cell are not stored here.
  std::map< vtkIdType, vtkIdType > ptIdMap;

743
744
745
746
  IntersectionMapIteratorType iterLower = map->lower_bound(cellId);
  IntersectionMapIteratorType iterUpper = map->upper_bound(cellId);
  //Get all the lines associated with the original cell
  while (iterLower != iterUpper)
747
748
749
    {
    vtkIdType lineId = iterLower->second;
    vtkIdType nLinePts, *linePtIds;
750
751
752
    interLines->GetLines()->GetCell(3*lineId, nLinePts, linePtIds);

    interceptlines->InsertNextCell(2);
753
    lines->InsertNextCell(2);
754
    //Loop through the points of each line
755
756
757
    for (vtkIdType i = 0; i < nLinePts; i++)
      {
      std::map< vtkIdType, vtkIdType >::iterator location =
758
759
760
        ptIdMap.find(linePtIds[i]);
      //If point already isn't in list
      if (location == ptIdMap.end())
761
        {
762
763
        interLines->GetPoint(linePtIds[i], xyz);
        if (linePtIds[i] >= interLines->GetNumberOfPoints())
764
765
766
          {
          vtkGenericWarningMacro(<< "invalid point read 2");
          }
767
768
769
770
771
772
773
774
775
776
        //Check to see if point is unique
        int unique=merger->InsertUniquePoint(xyz, ptIdMap[ linePtIds[i] ]);
        if (unique)
          {
          //If point is unique, check to see if it is actually a
          //point originating from this input surface or on both surfaces
          //Don't mark as boundary point if it originates from other surface
          if (surfaceMapper->GetValue(linePtIds[i]) == inputIndex + 1 ||
              surfaceMapper->GetValue(linePtIds[i]) == 3)
            {
Adam Updegrove's avatar
Adam Updegrove committed
777
            cellBoundaryPt->InsertValue(ptIdMap[linePtIds[i]], 1);
778
779
780
            }
          else
            {
Adam Updegrove's avatar
Adam Updegrove committed
781
            cellBoundaryPt->InsertValue(ptIdMap[linePtIds[i]], 0);
782
783
784
785
786
787
788
789
790
791
792
793
794
            }
          }
        else
          {
          //Obviously if the pointid is less than three, it is one of the
          //original cell points and can be added to the inter cell point arr
          if (ptIdMap[linePtIds[i] ] < 3)
            {
            CellPointOnInterLine[ptIdMap[linePtIds[i]]] = 1;
            }
          }
        interceptlines->InsertCellPoint(ptIdMap[ linePtIds[i] ]);
        lines->InsertCellPoint(ptIdMap[ linePtIds[i] ]);
795
        }
796
      //Point is already in list, so run through checks with its value
797
798
      else
        {
799
800
801
802
803
804
        interceptlines->InsertCellPoint(location->second);
        lines->InsertCellPoint(location->second);
        if (location->second < 3)
          {
          CellPointOnInterLine[location->second] = 1;
          }
805
806
807
808
809
810
        }
      }
    ++iterLower;
    }

  // Now check the neighbors of the cell
811
812
813
814
  IntersectionMapIteratorType ptIterLower;
  IntersectionMapIteratorType ptIterUpper;
  IntersectionMapIteratorType cellIterLower;
  IntersectionMapIteratorType cellIterUpper;
815
816
817
818
  vtkSmartPointer< vtkIdList > nbrCellIds =
    vtkSmartPointer< vtkIdList >::New();
  for (vtkIdType i = 0; i < 3; i++)
    {
819
    //Get Points belonging to each edge of this cell
820
821
822
823
    vtkIdType edgePtId0 = cellPts[i];
    vtkIdType edgePtId1 = cellPts[(i+1) % 3];

    double edgePt0[3], edgePt1[3];
824
    if (edgePtId0 >= input->GetNumberOfPoints())
825
      {
826
      vtkGenericWarningMacro(<< "invalid point read 3");
827
      }
828
    if (edgePtId1 >= input->GetNumberOfPoints())
829
      {
830
      vtkGenericWarningMacro(<< "invalid point read 4");
831
832
833
834
835
836
      }
    input->GetPoint(edgePtId0, edgePt0);
    input->GetPoint(edgePtId1, edgePt1);

    nbrCellIds->Reset();
    input->GetCellEdgeNeighbors(cellId, edgePtId0, edgePtId1, nbrCellIds);
837
    //Loop through attached neighbor cells and check for split edges
838
839
    for (vtkIdType j = 0; j < nbrCellIds->GetNumberOfIds(); j++)
      {
840
841
842
843
      vtkIdType nbrCellId = nbrCellIds->GetId(j);
      iterLower = map->lower_bound(nbrCellId);
      iterUpper = map->upper_bound(nbrCellId);
      while (iterLower != iterUpper)
844
845
846
        {
        vtkIdType lineId = iterLower->second;
        vtkIdType nLinePts, *linePtIds;
847
        interLines->GetLines()->GetCell(3*lineId, nLinePts, linePtIds);
848
849
        for (vtkIdType k = 0; k < nLinePts; k++)
          {
850
          if (linePtIds[k] >= interLines->GetNumberOfPoints())
851
            {
852
            vtkGenericWarningMacro(<< "invalid point read 5");
853
            }
854
855
856
857
858
          interLines->GetPoint(linePtIds[k], xyz);
          ptIterLower = this->PointMapper->lower_bound(linePtIds[k]);
          ptIterUpper = this->PointMapper->upper_bound(linePtIds[k]);
          //Find all points withing this neighbor cell
          while (ptIterLower != ptIterUpper)
859
            {
860
861
862
863
864
865
866
            vtkIdType mappedPtId = ptIterLower->second;
            cellIterLower = this->IntersectionPtsMap[inputIndex]->
              lower_bound(mappedPtId);
            cellIterUpper = this->IntersectionPtsMap[inputIndex]->
              upper_bound(mappedPtId);
            //Check all cell values associated with this point
            while (cellIterLower != cellIterUpper)
867
              {
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
              vtkIdType checkCellId = cellIterLower->second;

              //If this cell id is the same as the current cell id, this
              //means the point is a split edge, need to add to list!!
              if (checkCellId == cellId)
                {
                int unique=0;
                if (ptIdMap.find(linePtIds[k]) == ptIdMap.end())
                  {
                  unique = merger->InsertUniquePoint(xyz,
                      ptIdMap[ linePtIds[k] ]);
                  }

                else
                  {
                  //Point is less than 3, original cell point
                  if (ptIdMap[ linePtIds[k] ] < 3)
                    {
                    CellPointOnInterLine[ptIdMap[ linePtIds[k] ]] = 1;
                    }

                  }

                if (unique)
                  {
                  //Check to see what surface point originates from. Don't
                  //mark if point is from other surface
                  if (surfaceMapper->GetValue(linePtIds[k]) == inputIndex + 1
                      || surfaceMapper->GetValue(linePtIds[k]) == 3)
                    {
Adam Updegrove's avatar
Adam Updegrove committed
898
                    cellBoundaryPt->InsertValue(ptIdMap[linePtIds[k]], 1);
899
900
901
902
                    }

                  else
                    {
Adam Updegrove's avatar
Adam Updegrove committed
903
                    cellBoundaryPt->InsertValue(ptIdMap[linePtIds[k]], 0);
904
905
906
907
908
909
910
911
912
913
914
915
                    }

                  }
                else
                  {
                  if (ptIdMap[ linePtIds[k] ] < 3)
                    {
                    CellPointOnInterLine[ptIdMap[ linePtIds[k] ]] = 1;
                    }
                  }
                }
              ++cellIterLower;
916
              }
917
            ++ptIterLower;
918
919
920
921
922
923
924
925
926
            }
          }
        ++iterLower;
        }
      }
    }

  // Set up reverse ID map
  std::map< vtkIdType, vtkIdType > reverseIdMap;
927
  std::map< vtkIdType, vtkIdType > reverseLineIdMap;
928
  std::map< vtkIdType, vtkIdType >::iterator iter = ptIdMap.begin();
929
  while (iter != ptIdMap.end())
930
931
932
933
934
935
    {
    // If we have more than one point mapping back to the same point
    // in the input mesh, just use the first one. This will give a
    // preference for using cell points when an intersection line shares
    // a point with a a cell and prevent introducing accidental holes
    // in the mesh.
936
    if (reverseIdMap.find(iter->second) == reverseIdMap.end())
937
938
939
      {
      reverseIdMap[ iter->second ] = iter->first + input->GetNumberOfPoints();
      }
940
941
942
943
    if (reverseLineIdMap.find(iter->second) == reverseLineIdMap.end())
      {
      reverseLineIdMap[ iter->second ] = iter->first;
      }
944
945
946
947
    ++iter;
    }

  double v0[3], v1[3], n[3], c[3];
948
949
950
  vtkTriangle::TriangleCenter(p0, p1, p2, c);
  vtkTriangle::ComputeNormal(p0, p1, p2, n);
  vtkMath::Perpendiculars(n, v0, v1, 0.0);
951
952
953
954

  // For each point on an edge, compute it's relative angle about n.
  vtkSmartPointer< vtkIdTypeArray > edgePtIdList =
    vtkSmartPointer< vtkIdTypeArray >::New();
955
956
957
  vtkSmartPointer< vtkIdTypeArray > interPtIdList =
    vtkSmartPointer< vtkIdTypeArray >::New();
  edgePtIdList->Allocate(points->GetNumberOfPoints());
958
959
  vtkSmartPointer< vtkDoubleArray > angleList =
    vtkSmartPointer< vtkDoubleArray >::New();
960
961
962
  angleList->Allocate(points->GetNumberOfPoints());
  bool *interPtBool;
  interPtBool = new bool[points->GetNumberOfPoints()];
963

964
  for (vtkIdType ptId = 0; ptId < points->GetNumberOfPoints(); ptId++)
965
    {
966
967
    double x[3];
    points->GetPoint(ptId, x);
968

969
970
    interPtBool[ptId] = false;
    if (cellBoundaryPt->GetValue(ptId))
971
972
973
      {
      // Point is on line. Add its id to id list and add its angle to
      // angle list.
974
      edgePtIdList->InsertNextValue(ptId);
975
      double d[3];
976
977
978
979
980
981
982
983
984
985
986
987
988
      vtkMath::Subtract(x, c, d);
      angleList->InsertNextValue(atan2(vtkMath::Dot(d, v0),
                                         vtkMath::Dot(d, v1)));
        if (ptId > 2)
          {
          //Intersection Point!
          interPtIdList->InsertNextValue(ptId);
          interPtBool[ptId] = true;
          }
      }
    //Setting the boundary points
    if (ptId > 2)
      {
Adam Updegrove's avatar
Adam Updegrove committed
989
      this->BoundaryPoints[inputIndex]->InsertValue(reverseIdMap[ptId], 1);
990
991
992
      }
    else if (CellPointOnInterLine[ptId])
      {
Adam Updegrove's avatar
Adam Updegrove committed
993
      this->BoundaryPoints[inputIndex]->InsertValue(cellPts[ptId], 1);
994
995
996
      }
    else
      {
Adam Updegrove's avatar
Adam Updegrove committed
997
      this->BoundaryPoints[inputIndex]->InsertValue(cellPts[ptId], 0);
998
999
      }

1000
    }