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

  Program:   Visualization Toolkit
  Module:    vtkCookieCutter.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.

=========================================================================*/
#include "vtkCookieCutter.h"

#include "vtkCellData.h"
#include "vtkExecutive.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkLine.h"
#include "vtkMergePoints.h"
#include "vtkPolygon.h"
#include "vtkStreamingDemandDrivenPipeline.h"

#include <set>
#include <vector>

vtkStandardNewMacro(vtkCookieCutter);
vtkCxxSetObjectMacro(vtkCookieCutter, Locator, vtkIncrementalPointLocator);

// Helper functions------------------------------------------------------------

// Precision in the parametric coordinate system. Note that intersection
// points, and/or parallel segments routinely occur during processing. This
// tolerance is used to merge nearly coincident points.
#define VTK_DEGENERATE_TOL 0.001

namespace
{

// Note on the definition of parametric coordinates: Given a sequence of
// line segments (vi,vi+1) that form a primitive (e.g., polyline or
// polygon), the parametric coordinate t along the primitive is
// [i,i+1). Any point on a segment (like an intersection point on the
// segment) has parametric coordinate i+t, where 0 <= t < 1.

// Infrastructure for cropping----------------------------------------------
// Note that the classification Class data member is has multiple meanings,
// being applied to points as well as line segments. Initially the vertices
// are classified. Then the segments [i,(i+1)) are classified by combining
// the vertex classifications (and possibly other information like in/out
// tests). For example, the INSIDE or OUTSIDE classification refers to a
// segment that lies inside or outside the boundary of the loop
// (respectively), or to a point that may be inside or outside. The
// INTERSECTION classification refers to a point generated by the
// intersection of two edges. MULT_INTS refers to multiple >1 merged
// intersection points. The ON classification refers to a point (maybe
// merged) that has at least one ON point. The ON classification for a
// segment also holds if both of the end points are ON and the OnId of both
// are the same.
struct SortPoint
{
  // vertex and segment classification
  enum ClassEnum
  {
    VERTEX = 0,
    OUTSIDE = 1,
    INSIDE = 2,
    INTERSECTION = 4,
    MULT_INTS = 8,
    ON = 16
  };
  double T;       // parametric coordinate along primitive (line or polygon)
  int Class;      // classification of vertex or segment
  vtkIdType Id;   // Id assigned to INTERSECTION or ON point, used later to build loops
  vtkIdType OnId; // Id of edge pair if classification of point is ON
  double X[3];    // position of point
  SortPoint(double t, int c, vtkIdType id, vtkIdType onId, double x[3])
    : T(t)
    , Class(c)
    , Id(id)
    , OnId(onId)
  {
    this->X[0] = x[0];
    this->X[1] = x[1];
    this->X[2] = x[2];
  }
};

// Special sort operation on primitive parametric coordinate T-------------
bool PointSorter(SortPoint const& lhs, SortPoint const& rhs)
{
  return lhs.T < rhs.T;
}

// Vectors are used to hold points, eventually sorted
typedef std::vector<SortPoint> SortedPointsType;

// Used to analyze and merge nearly coincident points----------------------
struct MergeRange
{ // points to merge [StartIdx,EndIdx) possibly modulo around loop
  int StartIdx;
  int EndIdx;
  MergeRange(int start, int end)
    : StartIdx(start)
    , EndIdx(end)
  {
  }
};

// Vectors are used to hold merge regions
typedef std::vector<MergeRange> MergeRangeType;

// Convenience function classifies a segment of a loop or polyline--------
int ClassifySegment(SortedPointsType& sortedPoints, int i, int j, vtkIdType npts, double* p,
  double bds[6], double n[3])
{
  double midSegment[3];
  midSegment[0] = 0.5 * (sortedPoints[i].X[0] + sortedPoints[j].X[0]);
  midSegment[1] = 0.5 * (sortedPoints[i].X[1] + sortedPoints[j].X[1]);
  midSegment[2] = 0.5 * (sortedPoints[i].X[2] + sortedPoints[j].X[2]);
  if (vtkPolygon::PointInPolygon(midSegment, npts, p, bds, n) == 1)
  {
    return SortPoint::INSIDE;
  }
  else
  {
    return SortPoint::OUTSIDE;
  }
}

// Merge the list of coincident intersections along a polyline-------------
// The points are assumed sorted in parametric coordinates, not closed.
int CleanSortedPolyline(SortedPointsType& sortedPoints)
{
  int num, i, ip, j;
  double t, tp;
  bool needsAnalysis = false;
  vtkIdType maxId;

  // First do a quick check. If there are no degenerate situations just
  // return.
  num = static_cast<int>(sortedPoints.size());
  maxId = sortedPoints[0].Id;
  for (i = 0; i < (num - 1); ++i)
  {
    ip = i + 1;
    t = sortedPoints[i].T;
    tp = sortedPoints[ip].T;

    maxId = (sortedPoints[ip].Id > maxId ? sortedPoints[ip].Id : maxId);

    if (fabs(tp - t) <= VTK_DEGENERATE_TOL)
    {
      needsAnalysis = true;
    }
  }

  if (!needsAnalysis)
  {
    return maxId + 1;
  }

  // If here a deeper analysis is required. We need to segment groups of
  // points and/or vertices; deleting some and updating others.
  int start = 0, end = 0;
  MergeRangeType mergeRange;

  // segment nearly coincident points into groups [start,end)
  int imax = num - 1;
  for (i = 0; i < num;)
  {
    t = sortedPoints[i].T;
    if (end == (num - 1))
    { // last point may require special treatment
      mergeRange.push_back(MergeRange(end, end + 1));
      break;
    }
    ip = i + 1;
    tp = sortedPoints[ip].T;

    // special treatment for first point in the loop
    while (fabs(tp - t) <= VTK_DEGENERATE_TOL && ip < imax)
    {
      ip++;
      tp = sortedPoints[ip].T;
    }
    end = ip;

    // Add new segment
    mergeRange.push_back(MergeRange(start, end));

    // Move to next segment
    i = start = end;
  }

  // Now perform the analysis and update the sorted points
  SortedPointsType newSortedPoints;
  vtkIdType minId, onId, minX = 0;
  int spc, nints, sze, numMR = static_cast<int>(mergeRange.size());
  double minT;
  for (i = 0; i < numMR; ++i)
  {
    start = mergeRange[i].StartIdx;
    end = mergeRange[i].EndIdx;
    sze = (end > start ? end - start : (end + num) - start);
    if (sze == 1) // just copy vertex
    {
      newSortedPoints.push_back(sortedPoints[start]);
      continue;
    }

    minId = onId = VTK_ID_MAX;
    minT = sortedPoints[start].T;
    minX = start;
    spc = SortPoint::VERTEX;
    for (nints = 0, j = 0; j < sze; ++j)
    {
      ip = start + j;
      if (sortedPoints[ip].Id >= 0)
      {
        nints++;
        minId =
          (sortedPoints[ip].Id >= 0 && sortedPoints[ip].Id < minId ? sortedPoints[ip].Id : minId);
        onId = (sortedPoints[ip].OnId >= 0 && sortedPoints[ip].OnId < onId ? sortedPoints[ip].OnId
                                                                           : onId);
        if (sortedPoints[ip].T < minT)
        {
          minT = sortedPoints[ip].T;
          minX = ip;
        }
        spc |= sortedPoints[ip].Class;
      }
    }

    if (spc == SortPoint::INTERSECTION && nints > 1)
    {
      spc = SortPoint::MULT_INTS;
    }
    newSortedPoints.push_back(SortPoint(minT, spc, minId, onId, sortedPoints[minX].X));
  } // across all merge ranges

  // Update the sorted points array, now clean!
  sortedPoints = newSortedPoints;

  // Max id is used for allocating some space shortly
  num = static_cast<int>(sortedPoints.size());
  for (maxId = 0, i = 0; i < num; ++i)
  {
    maxId = (sortedPoints[i].Id > maxId ? sortedPoints[i].Id : maxId);
  }

  return maxId + 1;
}

// Clean up the list of intersections around the polygon-------------------
// The points are assumed sorted in parametric coordinates, closed
// loop. Nearly coincident points are merged.
//
void CleanSortedPolygon(vtkIdType npts, SortedPointsType& sortedPoints)
{
  int num, i, im, ip, j;
  double t, tm, tp = 0, moduloOffset = static_cast<double>(npts);
  bool needsAnalysis = false;

  // First do a quick check. If there are no degenerate situations just
  // return.
  num = static_cast<int>(sortedPoints.size());
  for (i = 0; i < num; ++i)
  {
    ip = (i + 1) % num;
    t = sortedPoints[i].T;
    if ((tp = sortedPoints[ip].T) < t)
    { // in case of modulo wrap
      tp += moduloOffset;
    }

    if (fabs(tp - t) <= VTK_DEGENERATE_TOL)
    {
      needsAnalysis = true;
    }
  }

  if (!needsAnalysis)
  {
    return;
  }

  // If here a deeper analysis is required. We need to segment groups of
  // points and/or vertices; eventually deleting some and updating the
  // classification of the remaining point.
  int start = 0, end = 0;
  MergeRangeType mergeRange;

  // Segment nearly coincident points into groups [start,end)
  int imax = num;
  for (i = 0; i < imax;)
  {
    if (i == (imax - 1))
    {
      mergeRange.push_back(MergeRange(i, imax));
      break;
    }

    t = sortedPoints[i].T;
    ip = (i + 1) % num;
    tp = sortedPoints[ip].T;

    // Special treatment for first point in the loop (due to modulo wrap)
    if (i == 0)
    {
      im = num - 1;
      tm = npts - sortedPoints[im].T;
      while (fabs(t - tm) <= VTK_DEGENERATE_TOL)
      {
        im--;
        tm = npts - sortedPoints[im].T;
      }
      start = (im + 1) % num;
      imax = (start == 0 ? num : start);
    }

    // Now proceed in forward direction
    while (fabs(tp - t) <= VTK_DEGENERATE_TOL && ip < imax)
    {
      ip++;
      tp = sortedPoints[ip % num].T;
    }
    end = ip;

    // Add new segment
    mergeRange.push_back(MergeRange(start, end));

    // Move to next segment
    i = start = end;
  }

  // Now perform the analysis and update the sorted points
  SortedPointsType newSortedPoints;
  vtkIdType minId, onId, minX = 0;
  int spc, nints, sze, numMR = static_cast<int>(mergeRange.size());
  double minT;
  for (i = 0; i < numMR; ++i)
  {
    start = mergeRange[i].StartIdx;
    end = mergeRange[i].EndIdx;
    sze = (end > start ? end - start : (end + num) - start);
    if (sze == 1) // just copy vertex
    {
      newSortedPoints.push_back(sortedPoints[start]);
      continue;
    }

    // Aggregate the final classification, point id, and other information
    // for merged point.
    minId = onId = VTK_ID_MAX;
    minT = sortedPoints[start].T;
    minX = start;
    spc = SortPoint::VERTEX;
    for (nints = 0, j = 0; j < sze; ++j)
    {
      ip = (start + j) % num;
      if (sortedPoints[ip].Class != SortPoint::VERTEX)
      {
        nints++;
        minId =
          (sortedPoints[ip].Id >= 0 && sortedPoints[ip].Id < minId ? sortedPoints[ip].Id : minId);
        onId = (sortedPoints[ip].OnId >= 0 && sortedPoints[ip].OnId < onId ? sortedPoints[ip].OnId
                                                                           : onId);
        spc |= sortedPoints[ip].Class;
        if (sortedPoints[ip].T < minT)
        {
          minT = sortedPoints[ip].T;
          minX = ip;
        }
      }
    }

    if (spc == SortPoint::INTERSECTION && nints > 1)
    {
      spc = SortPoint::MULT_INTS;
    }
    newSortedPoints.push_back(SortPoint(minT, spc, minId, onId, sortedPoints[minX].X));
  } // across all merge ranges

  // Update the sorted points array, now clean!
  sortedPoints = newSortedPoints;
}

// Classify polyline segments----------------------------------------------
// Entering this function the points are classified as either VERTEX,
// INTERSECTION, MULT_INTS, or ON. Combine this information to classify each
// segment of the polyline. Note that whenever possible a previous segment
// classification is propagated to the next segment to avoid extra in/out tests.
int ClassifyPolyline(
  SortedPointsType& sortedPoints, vtkIdType npts, double* p, double bds[6], double n[3])
{
  // Check for Degenerate case. Can happen when all points are merged to one.
  int num = static_cast<int>(sortedPoints.size());
  if (num < 2)
  {
    sortedPoints[0].Class = SortPoint::OUTSIDE;
    return 1;
  }

  // Classify the initial segment.
  int prevSegClass;
  if (sortedPoints[0].Class >= SortPoint::ON && sortedPoints[1].Class >= SortPoint::ON &&
    sortedPoints[0].OnId == sortedPoints[1].OnId)
  {
    prevSegClass = SortPoint::ON;
  }
  else // perform in/out test
  {
    prevSegClass = ClassifySegment(sortedPoints, 0, 1, npts, p, bds, n);
    sortedPoints[0].Class = prevSegClass;
  }

  // Okay now work along the polyline, propagating classification when possible
  // and otherwise determining the classification of each segment.
  int i, ip;
  int currClass, nextClass;
  for (i = 1; i < (num - 1); ++i)
  {
    ip = i + 1;
    currClass = sortedPoints[i].Class;
    nextClass = sortedPoints[ip].Class;

    if (currClass == SortPoint::VERTEX)
    { // propagate forward
      ;
    }
    else if (currClass == SortPoint::INTERSECTION)
    { // Flip classification
      prevSegClass = (prevSegClass == SortPoint::INSIDE ? SortPoint::OUTSIDE : SortPoint::INSIDE);
    }
    else if (currClass >= SortPoint::ON && nextClass >= SortPoint::ON &&
      sortedPoints[i].OnId == sortedPoints[ip].OnId)
    { // Segment is on loop segment
      prevSegClass = SortPoint::ON;
    }
    else
    { // complicated, do an in/out check
      prevSegClass = ClassifySegment(sortedPoints, i, ip, npts, p, bds, n);
    }
    sortedPoints[i].Class = prevSegClass;
  }

  return 1;
}

// Classify polygon segments----------------------------------------------
// Entering this function the points are classified as either VERTEX,
// INTERSECTION, MULT_INTS, or ON. Combine this information to classify each
// segment of the polygon. Note that whenever possible a previous segment
// classification is propagated to the next segment to avoid extra in/out tests.
// This deals with closed, modulo processing. Return a non-zero value if the
// loop is complex, i.e., has a "ON" segment classification.
int ClassifyPolygon(
  SortedPointsType& sortedPoints, vtkIdType npts, double* p, double bds[6], double n[3])
{
  // Check for Degenerate case. Can happen when all points are merged to one or two.
  int num = static_cast<int>(sortedPoints.size());
  if (num < 3)
  {
    for (int i = 0; i < num; ++i)
    {
      sortedPoints[i].Class = SortPoint::OUTSIDE;
    }
    return 0;
  }

  // Classify the initial segment.
  int segClass, startClass = sortedPoints[0].Class, hasOnClassification = 0;
  if (sortedPoints[0].Class >= SortPoint::ON && sortedPoints[1].Class >= SortPoint::ON &&
    sortedPoints[0].OnId == sortedPoints[1].OnId)
  {
    segClass = SortPoint::ON;
    hasOnClassification = 1;
  }
  else // perform in/out test
  {
    segClass = ClassifySegment(sortedPoints, 0, 1, npts, p, bds, n);
  }
  sortedPoints[0].Class = segClass;

  // Okay now work around the polygon, propagating segment classification
  // when possible and otherwise determining the classification of each
  // segment.
  int i, ip;
  int currClass, nextClass;
  for (i = 1; i < num; ++i)
  {
    currClass = sortedPoints[i].Class;
    ip = (i + 1) % num;
    nextClass = (ip == 0 ? startClass : sortedPoints[ip].Class); // remember start

    if (currClass == SortPoint::VERTEX)
    { // Propagate the classification forward
      ;
    }
    else if (currClass == SortPoint::INTERSECTION)
    { // Flip the classification
      segClass = (segClass == SortPoint::INSIDE ? SortPoint::OUTSIDE : SortPoint::INSIDE);
    }
    else if (currClass >= SortPoint::ON && nextClass >= SortPoint::ON &&
      sortedPoints[i].OnId == sortedPoints[ip].OnId)
    { // Segment is on loop segment
      segClass = SortPoint::ON;
      hasOnClassification = 1;
    }
    else
    { // complicated, do an in/out check
      segClass = ClassifySegment(sortedPoints, i, ip, npts, p, bds, n);
    }
    sortedPoints[i].Class = segClass;
  }

  return hasOnClassification;
}

// A helper class used to crop vtkCookieCutter input polys and lines.
class vtkCookieCutterHelper
{
public:
  vtkCookieCutterHelper(vtkIncrementalPointLocator* locator, vtkPoints* inPts,
    vtkCellData* inCellData, vtkPoints* outPts, vtkCellArray* outLines, vtkCellArray* outPolys,
    vtkCellData* outCellData);

  virtual ~vtkCookieCutterHelper();

  // Process a polyline
  void CropLine(vtkIdType cellId, vtkIdType cellOffset, vtkIdType npts, const vtkIdType* pts,
    vtkPolygon* loop, double* l, double loopBds[6], double n[3]);

  // Process a polygon
  void CropPoly(vtkIdType cellId, vtkIdType cellOffset, vtkPolygon* poly, vtkIdType npts,
    const vtkIdType* pts, vtkPolygon* loop, double* l, double loopBds[6], double n[3]);

protected:
  // Convenience method
  inline void InsertPoint(double x[3], vtkCellArray* ca);

  // Check and clean up potentially non manifold situations
  // Used to clean up complex sets of lines assumed to form one or more loops
  // (i.e., polygons). At this point, these lines are classified "ON" or are
  // classified "INSIDE". Make sure they are manifold; trim out non-manifold
  // loops.
  int ResolveTopology(vtkPolyData* pd);

public:
  vtkIncrementalPointLocator* Locator;
  vtkPoints* InPoints;
  vtkCellData* InCellData;
  vtkPoints* OutPoints;
  vtkCellArray* OutLines;
  vtkCellArray* OutPolys;
  vtkCellData* OutCellData;
};

vtkCookieCutterHelper::vtkCookieCutterHelper(vtkIncrementalPointLocator* locator, vtkPoints* inPts,
  vtkCellData* inCellData, vtkPoints* outPts, vtkCellArray* outLines, vtkCellArray* outPolys,
  vtkCellData* outCellData)
  : Locator(locator)
  , InPoints(inPts)
  , InCellData(inCellData)
  , OutPoints(outPts)
  , OutLines(outLines)
  , OutPolys(outPolys)
  , OutCellData(outCellData)
{
}

vtkCookieCutterHelper::~vtkCookieCutterHelper() {}

void vtkCookieCutterHelper::InsertPoint(double x[3], vtkCellArray* ca)
{
  vtkIdType ptId;
  this->Locator->InsertUniquePoint(x, ptId);
  ca->InsertCellPoint(ptId);
}

void vtkCookieCutterHelper::CropLine(vtkIdType cellId, vtkIdType cellOffset, vtkIdType npts,
  const vtkIdType* pts, vtkPolygon* loop, double* l, double loopBds[6], double n[3])
{
  // Make sure that this is a valid line
  if (npts < 2)
  {
    return;
  }

  // Create a vector of points with parametric coordinates, etc. which will
  // be sorted later. The tricky part is getting the classification of each
  // point correctly.
  // First insert all of the vertices defining the polyline
  vtkIdType i, j, newCellId;
  double t, u, v, x[3], x0[3], x1[3], y0[3], y1[3];
  int result;
  SortedPointsType sortedPoints;
  for (i = 0; i < npts; ++i)
  {
    t = static_cast<double>(i);
    this->InPoints->GetPoint(pts[i], x);
    sortedPoints.push_back(SortPoint(t, SortPoint::VERTEX, -1, -1, x));
  }

  // Now insert any intersection points
  vtkIdType numInts = 0, numLoopPts = loop->Points->GetNumberOfPoints();
  for (numInts = 0, i = 0; i < (npts - 1); ++i)
  {
    this->InPoints->GetPoint(pts[i], x0);
    this->InPoints->GetPoint(pts[i + 1], x1);

    // Traverse polygon loop intersecting each polygon segment
    for (j = 0; j < numLoopPts; ++j)
    {
      loop->Points->GetPoint(j, y0);
      loop->Points->GetPoint((j + 1) % numLoopPts, y1);
      if ((result = vtkLine::Intersection(x0, x1, y0, y1, u, v)) == 2)
      {
        x[0] = x0[0] + u * (x1[0] - x0[0]);
        x[1] = x0[1] + u * (x1[1] - x0[1]);
        x[2] = x0[2] + u * (x1[2] - x0[2]);
        u += static_cast<double>(i);
        v += static_cast<double>(j);
        sortedPoints.push_back(SortPoint(u, SortPoint::INTERSECTION, numInts, -1, x));
        numInts++;
      }
      else if (result == 3) // parallel lines
      {
        double x10[3], c[3], c2[3];
        vtkMath::Subtract(x1, x0, x10);
        double tol2 = 1.0e-08 * sqrt(vtkMath::Dot(x10, x10));
        if (vtkLine::DistanceBetweenLines(x0, x1, y0, y1, c, c2, u, v) > tol2)
        {
          continue; // no intersection just move ahead
        }
        else
        {
          vtkIdType onId = i * numLoopPts + j;
          vtkLine::DistanceToLine(x0, y0, y1, u, c);
          if (-VTK_DEGENERATE_TOL <= u && u <= (1.0 + VTK_DEGENERATE_TOL))
          {
            sortedPoints.push_back(
              SortPoint(static_cast<double>(i), SortPoint::ON, numInts, onId, c));
            numInts++;
          }
          vtkLine::DistanceToLine(x1, y0, y1, u, c);
          if (-VTK_DEGENERATE_TOL <= u && u <= (1.0 + VTK_DEGENERATE_TOL))
          {
            sortedPoints.push_back(
              SortPoint(static_cast<double>(i) + 1.0, SortPoint::ON, numInts, onId, c));
            numInts++;
          }
          vtkLine::DistanceToLine(y0, x0, x1, u, c);
          if (-VTK_DEGENERATE_TOL <= u && u <= (1.0 + VTK_DEGENERATE_TOL))
          {
            sortedPoints.push_back(
              SortPoint(static_cast<double>(i) + u, SortPoint::ON, numInts, onId, c));
            numInts++;
          }
          vtkLine::DistanceToLine(y1, x0, x1, u, c);
          if (-VTK_DEGENERATE_TOL <= u && u <= (1.0 + VTK_DEGENERATE_TOL))
          {
            sortedPoints.push_back(
              SortPoint(static_cast<double>(i) + u, SortPoint::ON, numInts, onId, c));
            numInts++;
          }
        } // within tolerance of other line
      }   // parallel line
    }     // intersect all line segments that form the loop
  }       // for all line segments that make up this polyline

  // Sort in parametric space
  std::sort(sortedPoints.begin(), sortedPoints.end(), &PointSorter);

  // Clean up coincident points
  CleanSortedPolyline(sortedPoints);

  // Classify the segments of the polyline
  ClassifyPolyline(sortedPoints, numLoopPts, l, loopBds, n);

  // If here, then pieces of the intersected polyline are output. Note that
  // a "ON" classification for a polyline should be output.
  vtkIdType num = static_cast<vtkIdType>(sortedPoints.size());
  vtkIdType startIdx = 0, endIdx = 0, numInsertedPts;
  while (startIdx < (num - 1))
  {
    // move to the beginning of the next interior strand.
    while (startIdx < (num - 1) && sortedPoints[startIdx].Class == SortPoint::OUTSIDE)
    {
      startIdx++;
    }
    if (startIdx >= (num - 1))
    {
      continue; // all done
    }

    // Find the end of the run, i.e., link together interior strands.
    endIdx = startIdx + 1;
    while (endIdx < (num - 1) &&
      (sortedPoints[endIdx].Class == SortPoint::INSIDE ||
        sortedPoints[endIdx].Class == SortPoint::ON))
    {
      endIdx++;
    }

    // Output line segments
    numInsertedPts = endIdx - startIdx + 1;
    newCellId = this->OutLines->InsertNextCell(numInsertedPts) + cellOffset;
    this->OutCellData->CopyData(this->InCellData, cellId, newCellId);
    for (i = startIdx; i <= endIdx; ++i)
    {
      InsertPoint(sortedPoints[i].X, this->OutLines);
    }
    startIdx = endIdx;
  } // over all sorted points
} // CropLine

int vtkCookieCutterHelper::ResolveTopology(vtkPolyData* pd)
{
  vtkIdType numLines = pd->GetNumberOfLines();
  if (numLines < 3) // non-manifold, on-edge lines don't form a loop
  {
    return 0;
  }

  // Make a pass and make sure all points are manifold, or not used at all
  vtkPoints* pts = pd->GetPoints();
  vtkIdType pid, *cells, numPts = pts->GetNumberOfPoints();
  int numBoundary = 0, numNonManifold = 0;
  vtkIdType ncells;
  for (pid = 0; pid < numPts; ++pid)
  {
    pd->GetPointCells(pid, ncells, cells);
    if (ncells == 0 || ncells == 2)
    {
      ; // skip is ok
    }
    else if (ncells == 1)
    {
      numBoundary++;
    }
    else
    {
      numNonManifold++;
    }
  } // over all points

  // Return if everything is cool
  if (numBoundary == 0 && numNonManifold == 0)
  {
    return 1;
  }

  // Special case: non-manifolds must be balanced with boundary.
  // Otherwise things are really messed up.
  else if (numBoundary != numNonManifold)
  {
    return 0;
  }

  // At this point, we have to erode away the boundary segments and
  // leave just manifolds. It's rare but does happen.
  bool resolved = false;
  while (!resolved)
  {
    resolved = true;
    for (pid = 0; pid < numPts; ++pid)
    {
      pd->GetPointCells(pid, ncells, cells);
      if (ncells == 1)
      {
        pd->RemoveCellReference(cells[0]);
        resolved = false;
      }
    } // over all points
  }   // while not resolved

  // Do a last sanity check
  for (pid = 0; pid < numPts; ++pid)
  {
    pd->GetPointCells(pid, ncells, cells);
    if (ncells == 1 || ncells > 2)
    {
      return 0;
    }
  } // over all points
  return 1;
}

void vtkCookieCutterHelper::CropPoly(vtkIdType cellId, vtkIdType cellOffset, vtkPolygon* poly,
  vtkIdType npts, const vtkIdType* pts, vtkPolygon* loop, double* l, double loopBds[6], double n[3])
{
  // Make sure that this is a valid polygon
  if (npts < 3)
  {
    return;
  }

  // Make sure that the polygons actually overlap in the x-y plane
  double polyBds[6];
  double* p = static_cast<double*>(poly->Points->GetVoidPointer(0));
  poly->GetBounds(polyBds);
  if (loopBds[0] > polyBds[1] || loopBds[1] < polyBds[0] || loopBds[2] > polyBds[3] ||
    loopBds[3] < polyBds[2])
  {
    return;
  }

  // Run around the polygon inserting intersection points into two
  // (eventually sorted) arrays. These sorted arrays will alternate inside
  // paths with outside paths, sort of like a "braid". To construct polygon
  // loops, the braid is woven together to form the loops.
  vtkIdType i, j, newCellId, numPts = 0, numInts = 0;
  double t, u, v, x[3], x0[3], x1[3], y0[3], y1[3];
  int result;
  SortedPointsType loopPoints, polyPoints;
  vtkIdType numLoopPts = loop->Points->GetNumberOfPoints();
  for (i = 0; i < npts; ++i)
  {
    t = static_cast<double>(i);
    this->InPoints->GetPoint(pts[i], x);
    polyPoints.push_back(SortPoint(t, SortPoint::VERTEX, numPts, -1, x));
    numPts++;
  }
  for (i = 0; i < numLoopPts; ++i)
  {
    t = static_cast<double>(i);
    loop->Points->GetPoint(i, x);
    loopPoints.push_back(SortPoint(t, SortPoint::VERTEX, numPts, -1, x));
    numPts++;
  }

  // Now insert intersection points
  for (i = 0; i < npts; ++i)
  {
    this->InPoints->GetPoint(pts[i], x0);
    this->InPoints->GetPoint(pts[(i + 1) % npts], x1);

    // Traverse polygon loop intersecting each polygon segment
    for (j = 0; j < numLoopPts; ++j)
    {
      loop->Points->GetPoint(j, y0);
      loop->Points->GetPoint((j + 1) % numLoopPts, y1);
      if ((result = vtkLine::Intersection(x0, x1, y0, y1, u, v)) == 2)
      {
        x[0] = x0[0] + u * (x1[0] - x0[0]);
        x[1] = x0[1] + u * (x1[1] - x0[1]);
        x[2] = x0[2] + u * (x1[2] - x0[2]);
        u += static_cast<double>(i);
        v += static_cast<double>(j);
        polyPoints.push_back(SortPoint(u, SortPoint::INTERSECTION, numPts, -1, x));
        loopPoints.push_back(SortPoint(v, SortPoint::INTERSECTION, numPts, -1, x));
        numInts++;
        numPts++;
      }
      else if (result == 3) // parallel lines
      {
        double x10[3], c[3], c2[3];
        vtkMath::Subtract(x1, x0, x10);
        double tol2 = 1.0e-08 * sqrt(vtkMath::Dot(x10, x10));
        if (vtkLine::DistanceBetweenLines(x0, x1, y0, y1, c, c2, u, v) > tol2)
        {
          continue; // no intersection just move ahead
        }
        else
        {
          vtkIdType onId = i * numLoopPts + j;
          vtkLine::DistanceToLine(x0, y0, y1, u, c);
          if (-VTK_DEGENERATE_TOL <= u && u <= (1.0 + VTK_DEGENERATE_TOL))
          {
            polyPoints.push_back(SortPoint(static_cast<double>(i), SortPoint::ON, numPts, onId, c));
            loopPoints.push_back(
              SortPoint(static_cast<double>(j) + u, SortPoint::ON, numPts, onId, c));
            numInts++;
            numPts++;
          }
          vtkLine::DistanceToLine(x1, y0, y1, u, c);
          if (-VTK_DEGENERATE_TOL <= u && u <= (1.0 + VTK_DEGENERATE_TOL))
          {
            polyPoints.push_back(
              SortPoint(static_cast<double>(i) + 1.0, SortPoint::ON, numPts, onId, c));
            loopPoints.push_back(
              SortPoint(static_cast<double>(j) + u, SortPoint::ON, numPts, onId, c));
            numInts++;
            numPts++;
          }
          vtkLine::DistanceToLine(y0, x0, x1, u, c);
          if (-VTK_DEGENERATE_TOL <= u && u <= (1.0 + VTK_DEGENERATE_TOL))
          {
            polyPoints.push_back(
              SortPoint(static_cast<double>(i) + u, SortPoint::ON, numPts, onId, c));
            loopPoints.push_back(SortPoint(static_cast<double>(j), SortPoint::ON, numPts, onId, c));
            numInts++;
            numPts++;
          }
          vtkLine::DistanceToLine(y1, x0, x1, u, c);
          if (-VTK_DEGENERATE_TOL <= u && u <= (1.0 + VTK_DEGENERATE_TOL))
          {
            polyPoints.push_back(
              SortPoint(static_cast<double>(i) + u, SortPoint::ON, numPts, onId, c));
            loopPoints.push_back(
              SortPoint(static_cast<double>(j) + 1.0, SortPoint::ON, numPts, onId, c));
            numInts++;
            numPts++;
          }
        }
      }
    } // intersect all line segments that form the loop
  }   // for all line segments that make up this polygon

  // Sort in parametric coordinates around the intersected polygon and
  // loop.
  std::sort(polyPoints.begin(), polyPoints.end(), &PointSorter);
  std::sort(loopPoints.begin(), loopPoints.end(), &PointSorter);

  // Clean loops of nearly coincident points
  CleanSortedPolygon(npts, polyPoints);
  CleanSortedPolygon(numLoopPts, loopPoints);

  // If here, we identify the strands from the polygon and the loop. A
  // strand is the "inside" region between two intersection points. Strands
  // will be braided later to form loops. Note that potentially non-manifold
  // conditions (classified "ON") require a more complex loop building process.
  ClassifyPolygon(polyPoints, numLoopPts, l, loopBds, n);
  ClassifyPolygon(loopPoints, npts, p, polyBds, n);

  // Insert "INSIDE" or "ON" edge segments into polydata (polylines) and
  // build adjacency information. We are using vtkPolyData because it does
  // everything we want, although there is a lot of allocation / deallocation
  // going on which is a potential area of speed improvement.
  vtkNew<vtkPoints> pDataPts;
  pDataPts->SetNumberOfPoints(numPts);
  vtkNew<vtkCellArray> pDataLines;
  vtkNew<vtkPolyData> pData;
  pData->SetPoints(pDataPts);
  pData->SetLines(pDataLines);

  SortedPointsType* loops[2];
  loops[0] = &polyPoints;
  loops[1] = &loopPoints;
  SortedPointsType* curLoop;
  int loopIdx;
  size_t sze, ii, ip;
  std::set<vtkIdType> edgeExist;
  vtkIdType eId, id0, id1;
  for (loopIdx = 0; loopIdx < 2; ++loopIdx)
  {
    curLoop = loops[loopIdx];
    sze = curLoop->size();
    for (ii = 0; ii < sze; ++ii)
    {
      ip = (ii + 1) % sze;
      id0 = (*curLoop)[ii].Id;
      id1 = (*curLoop)[ip].Id;
      // Since the number of edges is small, create an unique edge id from
      // the two vertex ids (id0,id1) where id0 < id1.
      eId = (id0 < id1 ? id1 * numPts + id0 : id0 * numPts + id1);
      if (((*curLoop)[ii].Class == SortPoint::INSIDE || (*curLoop)[ii].Class == SortPoint::ON) &&
        edgeExist.find(eId) == edgeExist.end())
      {
        edgeExist.insert(eId);
        pDataPts->SetPoint(id0, (*curLoop)[ii].X);
        pDataPts->SetPoint(id1, (*curLoop)[ip].X);
        pDataLines->InsertNextCell(2);
        pDataLines->InsertCellPoint(id0);
        pDataLines->InsertCellPoint(id1);
      }
    }
  }
  pData->BuildLinks();

  // Check the topology of the edges and ensure that it is valid.  If there
  // are "ON" classifications, may need to remove potential non-manifold
  // edges. Bail out if can't fix any problems.
  if (!this->ResolveTopology(pData))
  {
    return;
  }

  // Build loops (i.e., output polygons)
  vtkIdType numInsertedPts, *cells, thisCell, nextId, startId;
  vtkIdType ncells;
  vtkIdType thisNPts;
  const vtkIdType* thisPts;
  std::vector<char> visited(numPts, 0);
  // Each unvisited, connected point generates a loop
  for (i = 0; i < numPts; ++i)
  {
    if (!visited[i])
    {
      startId = nextId = i;
      pData->GetPointCells(nextId, ncells, cells);
      if (ncells != 2) // unused or merged point
      {
        continue;
      }
      numInsertedPts = 0;
      thisCell = cells[0];
      newCellId = this->OutPolys->InsertNextCell(numInsertedPts) + cellOffset;
      this->OutCellData->CopyData(this->InCellData, cellId, newCellId);

      do
      {
        visited[nextId] = 1;
        numInsertedPts++;
        InsertPoint(pDataPts->GetPoint(nextId), this->OutPolys);
        pData->GetCellPoints(thisCell, thisNPts, thisPts);
        nextId = (thisPts[0] != nextId ? thisPts[0] : thisPts[1]);
        pData->GetPointCells(nextId, ncells, cells);
        thisCell = (cells[0] != thisCell ? cells[0] : cells[1]);
      } while (nextId != startId);
      this->OutPolys->UpdateCellCount(numInsertedPts);
    }
  }
} // That was easy! CropPoly

} // anonymous namespace

//----------------------------------------------------------------------------
// Instantiate object with empty loop.
vtkCookieCutter::vtkCookieCutter()
{
  this->Locator = nullptr;
  this->SetNumberOfInputPorts(2);
}

//----------------------------------------------------------------------------
vtkCookieCutter::~vtkCookieCutter()
{
  this->SetLocator(nullptr);
}

//----------------------------------------------------------------------------
void vtkCookieCutter::SetLoopsConnection(vtkAlgorithmOutput* algOutput)
{
  this->SetInputConnection(1, algOutput);
}

//----------------------------------------------------------------------------
vtkAlgorithmOutput* vtkCookieCutter::GetLoopsConnection()
{
  return this->GetInputConnection(1, 0);
}

//----------------------------------------------------------------------------
void vtkCookieCutter::SetLoopsData(vtkDataObject* input)
{
  this->SetInputData(1, input);
}

//----------------------------------------------------------------------------
vtkDataObject* vtkCookieCutter::GetLoops()
{
  if (this->GetNumberOfInputConnections(1) < 1)
  {
    return nullptr;
  }

  return this->GetExecutive()->GetInputData(1, 0);
}

//----------------------------------------------------------------------------
int vtkCookieCutter::RequestData(vtkInformation* vtkNotUsed(request),
  vtkInformationVector** inputVector, vtkInformationVector* outputVector)
{
  // get the info objects
  vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation* loopInfo = inputVector[1]->GetInformationObject(0);
  vtkInformation* outInfo = outputVector->GetInformationObject(0);

  // get the input and output
  vtkPolyData* input = vtkPolyData::SafeDownCast(inInfo->Get(vtkDataObject::DATA_OBJECT()));
  vtkPolyData* loops = vtkPolyData::SafeDownCast(loopInfo->Get(vtkDataObject::DATA_OBJECT()));
  vtkPolyData* output = vtkPolyData::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));

  // Initialize and check data
  vtkDebugMacro(<< "Cookie cutting...");

  if (input->GetNumberOfPoints() < 1)
  {
    vtkErrorMacro("Input contains no points");
    return 1;
  }

  if (!loops || loops->GetNumberOfCells() < 1)
  {
    vtkErrorMacro("Please define a polygonal loop with at least three points");
    return 1;
  }

  // Create output data objects and prepare for processing
  vtkPoints* inPts = input->GetPoints();
  vtkPoints* loopPts = loops->GetPoints();
  vtkPoints* outPts = inPts->NewInstance();

  vtkCellData* inCellData = input->GetCellData();
  vtkCellData* outCellData = output->GetCellData();
  outCellData->CopyAllocate(inCellData);

  vtkCellArray* inVerts = input->GetVerts();
  vtkCellArray* outVerts = vtkCellArray::New();
  outVerts->AllocateCopy(inVerts);

  vtkCellArray* inLines = input->GetLines();
  vtkCellArray* outLines = vtkCellArray::New();
  outLines->AllocateCopy(inLines);

  vtkCellArray* inPolys = input->GetPolys();
  vtkCellArray* inStrips = input->GetStrips();
  vtkCellArray* outPolys = vtkCellArray::New();
  outPolys->AllocateCopy(inPolys);

  // Locator used to merge potentially duplicate points
  if (this->Locator == nullptr)
  {
    this->CreateDefaultLocator();
  }
  this->Locator->InitPointInsertion(outPts, loops->GetBounds());

  // Setup the lines and polys cropping helper
  vtkCookieCutterHelper* helper = new vtkCookieCutterHelper(
    this->Locator, inPts, inCellData, outPts, outLines, outPolys, outCellData);

  // Initialize and create polygon representing the loop
  double n[3], bds[6];
  vtkPolygon* loop = vtkPolygon::New();
  vtkPolygon* poly = vtkPolygon::New();

  // Process loops from second input. Note that the cell id in vtkPolyData
  // starts with verts, then lines, then polys, then strips.
  vtkIdType npts;
  const vtkIdType* pts;
  vtkIdType cellId = 0;
  vtkCellArray* loopPolys = loops->GetPolys();
  for (loopPolys->InitTraversal(); loopPolys->GetNextCell(npts, pts);)
  {
    vtkIdType numLoopPts = npts;
    if (numLoopPts < 3)
    {
      continue; // need a valid polygon loop, skip this one
    }
    loop->Initialize(npts, pts, loopPts);
    loop->GetBounds(bds);
    vtkPolygon::ComputeNormal(loop->Points, n);
    double* l = static_cast<double*>(loop->Points->GetVoidPointer(0));

    // Start by processing the verts. A simple in/out check.
    double x[3];
    if (inVerts->GetNumberOfCells() > 0)
    {
      for (inVerts->InitTraversal(); inVerts->GetNextCell(npts, pts); ++cellId)
      {
        for (vtkIdType i = 0; i < npts; ++i)
        {
          inPts->GetPoint(pts[i], x);
          if (vtkPolygon::PointInPolygon(x, numLoopPts, l, bds, n) == 1)
          {
            vtkIdType ptId;
            this->Locator->InsertUniquePoint(x, ptId);
            vtkIdType newCellId = outVerts->InsertNextCell(1, &ptId);
            outCellData->CopyData(inCellData, cellId, newCellId);
          }
        }
      } // for all verts
    }   // if vert cells

    // Now process lines
    if (inLines->GetNumberOfCells() > 0)
    {
      vtkIdType cellOffset = outVerts->GetNumberOfCells();
      for (inLines->InitTraversal(); inLines->GetNextCell(npts, pts); ++cellId)
      {
        helper->CropLine(cellId, cellOffset, npts, pts, loop, l, bds, n);
      }
    } // if line cells

    // Now process polygons
    if (inPolys->GetNumberOfCells() > 0)
    {
      vtkIdType cellOffset = outVerts->GetNumberOfCells() + outLines->GetNumberOfCells();
      for (inPolys->InitTraversal(); inPolys->GetNextCell(npts, pts); ++cellId)
      {
        poly->Initialize(npts, pts, inPts);
        helper->CropPoly(cellId, cellOffset, poly, npts, pts, loop, l, bds, n);
      }
    } // if polygonal cells

    // Now process triangle strips
    if (inStrips->GetNumberOfCells() > 0)
    {
      vtkIdType cellOffset = outVerts->GetNumberOfCells() + outLines->GetNumberOfCells();
      for (inStrips->InitTraversal(); inStrips->GetNextCell(npts, pts); ++cellId)
      {
        vtkIdType numTriPts = 3, triPts[3];
        for (vtkIdType i = 0; i < (npts - 2); ++i)
        {
          triPts[0] = pts[i];
          triPts[1] = pts[i + 1];
          triPts[2] = pts[i + 2];
          poly->Initialize(3, triPts, inPts);
          helper->CropPoly(cellId, cellOffset, poly, numTriPts, triPts, loop, l, bds, n);
        }
      }
    } // if polygonal cells
  }   // for all loops

  delete helper;

  // Assign output as appropriate
  output->SetVerts(outVerts);
  outVerts->Delete();

  output->SetLines(outLines);
  outLines->Delete();

  output->SetPolys(outPolys);
  outPolys->Delete();

  // Clean up
  output->SetPoints(outPts);
  outPts->Delete();
  loop->Delete();
  poly->Delete();

  return 1;
}

//----------------------------------------------------------------------------
int vtkCookieCutter::RequestUpdateExtent(vtkInformation* vtkNotUsed(request),
  vtkInformationVector** inputVector, vtkInformationVector* outputVector)
{
  // get the info objects
  vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation* loopInfo = inputVector[1]->GetInformationObject(0);
  vtkInformation* outInfo = outputVector->GetInformationObject(0);

  if (loopInfo)
  {
    loopInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER(), 0);
    loopInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES(), 1);
    loopInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS(), 0);
  }
  inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER(),
    outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER()));
  inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES(),
    outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES()));
  inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS(),
    outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS()));
  inInfo->Set(vtkStreamingDemandDrivenPipeline::EXACT_EXTENT(), 1);

  return 1;
}

//----------------------------------------------------------------------------
int vtkCookieCutter::FillInputPortInformation(int port, vtkInformation* info)
{
  if (port == 0)
  {
    info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkPolyData");
    return 1;
  }
  else if (port == 1)
  {
    info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkPolyData");
    return 1;
  }
  return 0;
}

//----------------------------------------------------------------------------
void vtkCookieCutter::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os, indent);

  if (this->Locator)
  {
    os << indent << "Locator: " << this->Locator << "\n";
  }
  else
  {
    os << indent << "Locator: (none)\n";
  }
}

//-----------------------------------------------------------------------------
void vtkCookieCutter::CreateDefaultLocator()
{
  if (this->Locator == nullptr)
  {
    this->Locator = vtkMergePoints::New();
    this->Locator->Register(this);
    this->Locator->Delete();
  }
}
