vtkIntersectionPolyDataFilter.cxx

/*=========================================================================

  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.

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

#include "vtkIntersectionPolyDataFilter.h"

#include "vtkCellArray.h"
#include "vtkCellData.h"
#include "vtkCleanPolyData.h"
#include "vtkDelaunay2D.h"
#include "vtkDoubleArray.h"
#include "vtkIdList.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkLine.h"
#include "vtkLongArray.h"
#include "vtkOBBTree.h"
#include "vtkObjectFactory.h"
#include "vtkPlane.h"
#include "vtkPointData.h"
#include "vtkPointLocator.h"
#include "vtkPoints.h"
#include "vtkPolyDataNormals.h"
#include "vtkPolygon.h"
#include "vtkSmartPointer.h"
#include "vtkSortDataArray.h"
#include "vtkTriangle.h"
#include "vtkTriangleFilter.h"
#include "vtkTransform.h"
#include "vtkTransformPolyDataFilter.h"
#include "vtkUnstructuredGrid.h"

#include <map>
#include <sstream>
#include <iostream>
#include <list>

//----------------------------------------------------------------------------
// Helper typedefs and data structures.
namespace {

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

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

}

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();

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

  //Runs the split mesh for the designated input surface
  int SplitMesh(int inputIndex, vtkPolyData *output,
                vtkPolyData *intersectionLines);

protected:

  //Split cells into polygons created by intersection lines
  vtkCellArray* SplitCell(vtkPolyData *input, vtkIdType cellId,
                                     vtkIdType *cellPts,
                                     IntersectionMapType *map,
                                     vtkPolyData *interLines, int inputIndex,
                                     int numCurrCells);

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

  //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);

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

  // Stores the intersection lines.
  vtkCellArray        *IntersectionLines;

  vtkIdTypeArray      *SurfaceId;
  vtkIdTypeArray      *NewCellIds[2];

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

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

  // 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];
  vtkIntArray         *BoundaryPoints[2];

  // 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];
  IntersectionMapType *IntersectionPtsMap[2];
  IntersectionMapType *PointMapper;

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

  //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;

protected:
  Impl(const Impl&); // purposely not implemented
  void operator=(const Impl&); // purposely not implemented

};

//----------------------------------------------------------------------------
vtkIntersectionPolyDataFilter::Impl::Impl() :
  OBBTree1(0), IntersectionLines(0), SurfaceId(0), PointMerger(0)
{
  for (int i = 0; i < 2; i++)
    {
    this->Mesh[i]                 = NULL;
    this->CellIds[i]              = NULL;
    this->IntersectionMap[i]      = new IntersectionMapType();
    this->IntersectionPtsMap[i]   = new IntersectionMapType();
    this->PointEdgeMap[i]         = new PointEdgeMapType();
    }
  this->PointMapper             = new IntersectionMapType();
  this->SplittingPD             = vtkPolyData::New();
  this->TransformSign = 0;
  this->Tol = 1e-6;
}

//----------------------------------------------------------------------------
vtkIntersectionPolyDataFilter::Impl::~Impl()
{
  for (int i = 0; i < 2; i++)
    {
    delete this->IntersectionMap[i];
    delete this->IntersectionPtsMap[i];
    delete this->PointEdgeMap[i];
    }
  delete this->PointMapper;
  this->SplittingPD->Delete();
}


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

  //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
  int numCells0 = node0->Cells->GetNumberOfIds();

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

    //Make sure the cell is a triangle
    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];
            double surfaceid[2];
            int intersects =
              vtkIntersectionPolyDataFilter::TriangleTriangleIntersection
              (triPts0[0], triPts0[1], triPts0[2],
               triPts1[0], triPts1[1], triPts1[2],
               coplanar, outpt0, outpt1, surfaceid, tolerance);

            if (coplanar)
              {
              // Coplanar triangle intersection is not handled.
              // This intersection will not be included in the output. TODO
              //vtkDebugMacro(<<"Coplanar");
              intersects = 0;
              continue;
              }

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

              vtkIdType ptId0, ptId1;
              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])
                {
                intersectionSurfaceId->InsertValue(ptId0, 3);
                }
              else
                {
                if (unique[0])
                {
                intersectionSurfaceId->InsertValue(ptId0, surfaceid[0]);
                }
              else
                {
                if (intersectionSurfaceId->GetValue(ptId0) != 3)
                  {
                  intersectionSurfaceId->InsertValue(ptId0, surfaceid[0]);
                  }

                }
              if (unique[1])
                {
                intersectionSurfaceId->InsertValue(ptId1, surfaceid[1]);
                }
              else
                {
                if (intersectionSurfaceId->GetValue(ptId1) != 3)
                  {
                  intersectionSurfaceId->InsertValue(ptId1, surfaceid[1]);
                  }
                }
              }

            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();
              int newLine = info->CheckLine(lineTest, ptId0, ptId1);
              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);
              intersectionLines->InsertCellPoint(ptId0);
              intersectionLines->InsertCellPoint(ptId1);

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

              info->PointCellIds[0]->InsertValue(ptId0, cellId0);
              info->PointCellIds[0]->InsertValue(ptId1, cellId0);
              info->PointCellIds[1]->InsertValue(ptId0, cellId1);
              info->PointCellIds[1]->InsertValue(ptId1, cellId1);

              info->IntersectionMap[0]->
                insert(std::make_pair(cellId0, lineId));
              info->IntersectionMap[1]->
                insert(std::make_pair(cellId1, lineId));

              // Check which edges of cellId0 and cellId1 outpt0 and
              // outpt1 are on, if any.
              int isOnEdge=0;
              int m0p0=0, m0p1=0, m1p0=0, m1p1=0;
              for (vtkIdType edgeId = 0; edgeId < 3; edgeId++)
                {
                isOnEdge = info->AddToPointEdgeMap(0, ptId0, outpt0,
                    mesh0, cellId0, edgeId, lineId, triPtIds0);
                if (isOnEdge != -1)
                  {
                  m0p0++;
                  }
                isOnEdge = info->AddToPointEdgeMap(0, ptId1, outpt1,
                    mesh0, cellId0, edgeId, lineId, triPtIds0);
                if (isOnEdge != -1)
                  {
                  m0p1++;
                  }
                isOnEdge = info->AddToPointEdgeMap(1, ptId0, outpt0,
                    mesh1, cellId1, edgeId, lineId, triPtIds1);
                if (isOnEdge != -1)
                  {
                  m1p0++;
                  }
                isOnEdge = info->AddToPointEdgeMap(1, ptId1, outpt1,
                    mesh1, cellId1, edgeId, lineId, triPtIds1);
                if (isOnEdge != -1)
                  {
                  m1p1++;
                  }
                }
              //Special cases caught by tolerance and not from the Point
              //Merger
              if (m0p0 > 0 && m1p0 > 0)
                {
                intersectionSurfaceId->InsertValue(ptId0, 3);
                }
              if (m0p1 > 0 && m1p1 > 0)
                {
                intersectionSurfaceId->InsertValue(ptId1, 3);
                }
              }
            //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));
              }
            }
            }
          }
        }
      }
    }
  return 1;
}


//----------------------------------------------------------------------------
int vtkIntersectionPolyDataFilter::Impl
::SplitMesh(int inputIndex, vtkPolyData *output, vtkPolyData *intersectionLines)
{
  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();
  outPD->CopyAllocate(inPD, input->GetNumberOfPoints());

  // 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);
    this->BoundaryPoints[inputIndex]->InsertValue(ptId, 0);
    }

  // 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
    vtkIdType cellId = this->PointCellIds[inputIndex]->GetValue(id);

    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);
    this->BoundaryPoints[inputIndex]->InsertValue(newPtId, 0);
    }

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

  if (input->GetPolys()->GetNumberOfCells() > 0)
    {
    vtkCellArray *cells = input->GetPolys();
    vtkIdType newId = output->GetNumberOfCells();

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

    newPolys->EstimateSize(cells->GetNumberOfCells(), 3);
    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++)
      {
      if (nptsX != 3)
        {
        vtkGenericWarningMacro(<< "vtkIntersectionPolyDataFilter only works"
                                << " with triangle meshes.");
        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.
      bool needsSplit = intersectionMap->find(cellIdX)
        != intersectionMap->end();
      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);

          if (intersectionMap->find(nbrId) != intersectionMap->end())
            {
            needsSplit = true;
            }
          } // for (vtkIdType nbr = 0; ...
        } // for (vtkIdType pt = 0; ...

      // Splitting occurs here
      if (!needsSplit)
        {
        // Just insert the cell and copy the cell data
        newId = newPolys->InsertNextCell(3, pts);
        outCD->CopyData(inCD, cellIdX, newId);
        }
      else
        {
        //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();
        vtkCellArray *splitCells = this->SplitCell
          (input, cellIdX, pts, intersectionMap, splitLines,
           inputIndex,numCurrCells);
        if (splitCells == NULL)
          {
          vtkDebugWithObjectMacro(this->ParentFilter, <<"Error in splitting cell!");
          return 0;
          }


        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);

        vtkIdType npts, *ptIds, subCellId;
        splitCells->InitTraversal();
        for (subCellId = 0;splitCells->GetNextCell(npts, ptIds);subCellId++)
          {
          // 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);

          if (vtkMath::Dot(normal, subCellNormal) > 0)
            {
            newId = newPolys->InsertNextCell(npts, ptIds);
            }
          else
            {
            newId = newPolys->InsertNextCell(npts);
            for (int i = 0; i < npts; i++)
              {
              newPolys->InsertCellPoint(ptIds[ npts-i-1 ]);
              }
            }

          outCD->CopyData(inCD, cellIdX, newId); // Duplicate cell data
          }
        splitCells->Delete();
        }
      } // for (cells->InitTraversal(); ...
    } //if inputGetPolys()->GetNumberOfCells() > 1 ...

  return 1;
}

vtkCellArray* vtkIntersectionPolyDataFilter::Impl
::SplitCell(vtkPolyData *input, vtkIdType cellId, vtkIdType *cellPts,
            IntersectionMapType *map,
            vtkPolyData *interLines, int inputIndex,
            int numCurrCells)
{
  // 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};

  // Gather points from the cell
  vtkSmartPointer< vtkPoints > points = vtkSmartPointer< vtkPoints >::New();
  vtkSmartPointer< vtkPointLocator > merger =
    vtkSmartPointer< vtkPointLocator >::New();
  merger->SetTolerance(this->Tol);
  merger->InitPointInsertion(points, input->GetBounds());

  double xyz[3];
  for (int i = 0; i < 3; i++)
    {
    if (cellPts[i] >= input->GetNumberOfPoints())
      {
      vtkGenericWarningMacro(<< "invalid point read 1");
      }
    input->GetPoint(cellPts[i], xyz);
    merger->InsertNextPoint(xyz);
    cellBoundaryPt->InsertNextValue(1);
    }

  // Set up line cells and array to track the just the intersecting lines
  // on the cell.
  vtkSmartPointer< vtkCellArray > lines =
    vtkSmartPointer< vtkCellArray >::New();
  vtkSmartPointer< vtkCellArray > interceptlines =
    vtkSmartPointer< vtkCellArray >::New();

  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;

  IntersectionMapIteratorType iterLower = map->lower_bound(cellId);
  IntersectionMapIteratorType iterUpper = map->upper_bound(cellId);
  //Get all the lines associated with the original cell
  while (iterLower != iterUpper)
    {
    vtkIdType lineId = iterLower->second;
    vtkIdType nLinePts, *linePtIds;
    interLines->GetLines()->GetCell(3*lineId, nLinePts, linePtIds);

    interceptlines->InsertNextCell(2);
    lines->InsertNextCell(2);
    //Loop through the points of each line
    for (vtkIdType i = 0; i < nLinePts; i++)
      {
      std::map< vtkIdType, vtkIdType >::iterator location =
        ptIdMap.find(linePtIds[i]);
      //If point already isn't in list
      if (location == ptIdMap.end())
        {
        interLines->GetPoint(linePtIds[i], xyz);
        if (linePtIds[i] >= interLines->GetNumberOfPoints())
          {
          vtkGenericWarningMacro(<< "invalid point read 2");
          }
        //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)
            {
            cellBoundaryPt->InsertValue(ptIdMap[linePtIds[i]], 1);
            }
          else
            {
            cellBoundaryPt->InsertValue(ptIdMap[linePtIds[i]], 0);
            }
          }
        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] ]);
        }
      //Point is already in list, so run through checks with its value
      else
        {
        interceptlines->InsertCellPoint(location->second);
        lines->InsertCellPoint(location->second);
        if (location->second < 3)
          {
          CellPointOnInterLine[location->second] = 1;
          }
        }
      }
    ++iterLower;
    }

  // Now check the neighbors of the cell
  IntersectionMapIteratorType ptIterLower;
  IntersectionMapIteratorType ptIterUpper;
  IntersectionMapIteratorType cellIterLower;
  IntersectionMapIteratorType cellIterUpper;
  vtkSmartPointer< vtkIdList > nbrCellIds =
    vtkSmartPointer< vtkIdList >::New();
  for (vtkIdType i = 0; i < 3; i++)
    {
    //Get Points belonging to each edge of this cell
    vtkIdType edgePtId0 = cellPts[i];
    vtkIdType edgePtId1 = cellPts[(i+1) % 3];

    double edgePt0[3], edgePt1[3];
    if (edgePtId0 >= input->GetNumberOfPoints())
      {
      vtkGenericWarningMacro(<< "invalid point read 3");
      }
    if (edgePtId1 >= input->GetNumberOfPoints())
      {
      vtkGenericWarningMacro(<< "invalid point read 4");
      }
    input->GetPoint(edgePtId0, edgePt0);
    input->GetPoint(edgePtId1, edgePt1);

    nbrCellIds->Reset();
    input->GetCellEdgeNeighbors(cellId, edgePtId0, edgePtId1, nbrCellIds);
    //Loop through attached neighbor cells and check for split edges
    for (vtkIdType j = 0; j < nbrCellIds->GetNumberOfIds(); j++)
      {
      vtkIdType nbrCellId = nbrCellIds->GetId(j);
      iterLower = map->lower_bound(nbrCellId);
      iterUpper = map->upper_bound(nbrCellId);
      while (iterLower != iterUpper)
        {
        vtkIdType lineId = iterLower->second;
        vtkIdType nLinePts, *linePtIds;
        interLines->GetLines()->GetCell(3*lineId, nLinePts, linePtIds);
        for (vtkIdType k = 0; k < nLinePts; k++)
          {
          if (linePtIds[k] >= interLines->GetNumberOfPoints())
            {
            vtkGenericWarningMacro(<< "invalid point read 5");
            }
          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)
            {
            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)
              {
              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)
                    {
                    cellBoundaryPt->InsertValue(ptIdMap[linePtIds[k]], 1);
                    }

                  else
                    {
                    cellBoundaryPt->InsertValue(ptIdMap[linePtIds[k]], 0);
                    }

                  }
                else
                  {
                  if (ptIdMap[ linePtIds[k] ] < 3)
                    {
                    CellPointOnInterLine[ptIdMap[ linePtIds[k] ]] = 1;
                    }
                  }
                }
              ++cellIterLower;
              }
            ++ptIterLower;
            }
          }
        ++iterLower;
        }
      }
    }

  // Set up reverse ID map
  std::map< vtkIdType, vtkIdType > reverseIdMap;
  std::map< vtkIdType, vtkIdType > reverseLineIdMap;
  std::map< vtkIdType, vtkIdType >::iterator iter = ptIdMap.begin();
  while (iter != ptIdMap.end())
    {
    // 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.
    if (reverseIdMap.find(iter->second) == reverseIdMap.end())
      {
      reverseIdMap[ iter->second ] = iter->first + input->GetNumberOfPoints();
      }
    if (reverseLineIdMap.find(iter->second) == reverseLineIdMap.end())
      {
      reverseLineIdMap[ iter->second ] = iter->first;
      }
    ++iter;
    }

  double v0[3], v1[3], n[3], c[3];
  vtkTriangle::TriangleCenter(p0, p1, p2, c);
  vtkTriangle::ComputeNormal(p0, p1, p2, n);
  vtkMath::Perpendiculars(n, v0, v1, 0.0);

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

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

    interPtBool[ptId] = false;
    if (cellBoundaryPt->GetValue(ptId))
      {
      // Point is on line. Add its id to id list and add its angle to
      // angle list.
      edgePtIdList->InsertNextValue(ptId);
      double d[3];
      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)
      {
      this->BoundaryPoints[inputIndex]->InsertValue(reverseIdMap[ptId], 1);
      }
    else if (CellPointOnInterLine[ptId])
      {
      this->BoundaryPoints[inputIndex]->InsertValue(cellPts[ptId], 1);
      }
    else
      {
      this->BoundaryPoints[inputIndex]->InsertValue(cellPts[ptId], 0);
      }

    }