// SPDX-FileCopyrightText: Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
// SPDX-License-Identifier: BSD-3-Clause
#include "vtkCookieCutter.h"

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

#include <memory>
#include <set>
#include <vector>

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

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

// Precision in the parametric coordinate system. Note that nearly coincident
// 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
{

// Helper function to create cutting/trim loops from the second filter input,
// which may consist of polylines and/or polygons. Note that the polylines
// are assumed to close in on themselves, meaning that the first and last
// points are repeated, or if not, implicitly assumed to close and also to
// not self-intersect.  Also, polylines and polygons must consist of three
// points or more (it is assumed that they define a non-zero area; if the
// area is zero the cutting process will likely fail). The valid loops are
// placed into a cell array and returned. If the returned value is a nullptr,
// then no trim lines/polygons are available.
vtkSmartPointer<vtkCellArray> BuildTrimLoops(vtkPolyData* loops)
{
  // See what we have
  vtkCellArray* loopLines = loops->GetLines();
  vtkCellArray* loopPolys = loops->GetPolys();
  vtkIdType numLoopLines = loopLines->GetNumberOfCells();
  vtkIdType numLoopPolys = loopPolys->GetNumberOfCells();

  // If no lines, just return the polygon array
  if (numLoopLines < 1 && numLoopPolys < 1)
  {
    return nullptr;
  }

  else if (numLoopLines > 0 && numLoopPolys < 1)
  {
    return loopLines;
  } // Just has trim lines

  else if (numLoopPolys > 0 && numLoopLines < 1)
  {
    return loopPolys;
  } // Just has trim polygons

  else // Have to create a combined collection of cut polygons and lines
  {
    vtkNew<vtkCellArray> cutLoops;
    vtkIdType npts;
    const vtkIdType* pts;

    for (loopLines->InitTraversal(); loopLines->GetNextCell(npts, pts);)
    {
      // The polyline closes in on itself, to make it a polygon, we eliminate
      // the last point. This is far from foolproof, especially if the polyline
      // self-intersects etc.
      npts = (pts[0] == pts[npts - 1] ? npts - 1 : npts);
      if (npts > 2)
      {
        cutLoops->InsertNextCell(npts, pts);
      }
    }

    // Process polygons
    for (loopPolys->InitTraversal(); loopPolys->GetNextCell(npts, pts);)
    {
      if (npts > 2)
      {
        cutLoops->InsertNextCell(npts, pts);
      }
    }

    return cutLoops;
  } // Combination of trim lines and polygons
} // BuildTrimLoops

// 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 of intersection points via parametric coordinate T---
bool PointSorter(SortPoint const& lhs, SortPoint const& rhs)
{
  return lhs.T < rhs.T;
}

// Vectors are used to hold points, eventually sorted around the perimeter of
// polygons, or along lines.
using SortedPointsType = std::vector<SortPoint>;

// This is an adjunct to the SortPoint struct. It is used when passing point
// data through to the output, and used by vtkAttributeManager. It collects
// information about the point data so that when points are inserted into the
// final polydata output, the appropriate interpolation (or CopyData())
// process can be used.  These point attributes (which are referenced via
// SortPoint::Id in the vtkAttributeManager) may originate from a mesh vertex
// (MeshAttribute); a loop vertex (LoopAttribute); interpolation along a mesh
// edge (MeshEdgeAttribute); or interpolation along a loop edge
// (LoopEdgeAttribute).
struct PointAttribute
{
  enum PointAttributeEnum
  {
    MeshVertex = 0,
    LoopVertex = 1,
    Intersection = 2,
    MeshEdge = 3,
    LoopEdge = 4
  };

  // Note that V0 may be the originating point is from which to copy
  // data, or (V0,V1) may be the edge from which to interpolate data.
  int AttributeType; // one of PointAttributeEnum
  vtkIdType V0;      // Vertex id(s)
  vtkIdType V1;
  double T; // the parametric coordinate along the interpolating edge

  PointAttribute(int attType, vtkIdType v0, vtkIdType v1, double t)
    : AttributeType(attType)
    , V0(v0)
    , V1(v1)
    , T(t)
  {
  }
};

// Vectors are used to hold point attributes, eventually processed during the
// crop operation to produce output point data.
using PointAttributesType = std::vector<PointAttribute>;

// 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
using MergeRangeType = std::vector<MergeRange>;

// 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.emplace_back(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.emplace_back(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.emplace_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.emplace_back(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.emplace_back(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.emplace_back(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.emplace_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.emplace_back(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 manage point and cell attributes, and simplify the
// overall cookie cutting code.
struct vtkAttributeManager
{
  // Support cell data processing
  vtkCellData* InCellData;
  vtkCellData* OutCellData;

  // Support point data processing
  int PointInterpolation;              // indicate how to interpolate point data
  PointAttributesType PointAttributes; // point attributes for this cell
  SortedPointsType* LoopPts;
  SortedPointsType* PolyPts;
  vtkSmartPointer<vtkPointData> MeshPtData;
  vtkSmartPointer<vtkPointData> LoopPtData;
  vtkPointData* OutPtData;

  // Helper function to determine if two vtkDataSetAttributes are equivalent and may
  // be used together for copying and interpolating point data.
  bool HaveEquivalentAttributes(vtkDataSetAttributes* dsa0, vtkDataSetAttributes* dsa1)
  {
    // Check that the number of arrays is the same
    int numArrays0 = dsa0->GetNumberOfArrays();
    int numArrays1 = dsa1->GetNumberOfArrays();
    if (numArrays0 != numArrays1)
    {
      return false;
    }

    // Check that identical arrays exist in both (name, number of components, order)
    for (auto arrayNum = 0; arrayNum < numArrays0; ++arrayNum)
    {
      const char *name0, *name1;
      vtkDataArray* da0 = dsa0->GetArray(arrayNum);
      vtkDataArray* da1 = dsa1->GetArray(arrayNum);
      if (!da0 || !da1 || da0->GetDataType() != da1->GetDataType() ||
        da0->GetNumberOfComponents() != da1->GetNumberOfComponents() ||
        !(name0 = dsa0->GetArrayName(arrayNum)) || !(name1 = dsa1->GetArrayName(arrayNum)) ||
        strcmp(name0, name1) != 0)
      {
        return false;
      }
    }

    return true;
  }

  // Build point data attributes that are the set intersection of the mesh
  // point data attributes, and the loop point data attributes.  The
  // attributes must be equivalent because of the two-way interpolation/copy
  // data process that occurs.
  void IntersectAttributes(vtkPointData* meshPD, vtkPointData* loopPD)
  {
    this->MeshPtData.TakeReference(vtkPointData::New());
    this->LoopPtData.TakeReference(vtkPointData::New());

    // Loop over the the first attributes determining what data arrays are
    // common.
    int numMeshPDArrays = meshPD->GetNumberOfArrays();
    for (auto arrayNum = 0; arrayNum < numMeshPDArrays; ++arrayNum)
    {
      vtkDataArray *da0 = meshPD->GetArray(arrayNum), *da1;
      const char* name = meshPD->GetArrayName(arrayNum);
      if ((da1 = loopPD->GetArray(name)) != nullptr && da1->GetDataType() == da0->GetDataType() &&
        da1->GetNumberOfComponents() == da0->GetNumberOfComponents())
      {
        this->MeshPtData->AddArray(da0);
        this->LoopPtData->AddArray(da1);
      }
    }
  }

  // Constructor: Initialize the process of attribute processing
  vtkAttributeManager(vtkPolyData* input, vtkPolyData* loops, vtkPolyData* output, bool passCD,
    bool passPD, int ptInterpolation)
  {
    // Process cell data if requested
    this->InCellData = nullptr;
    this->OutCellData = nullptr;
    if (passCD)
    {
      this->InCellData = input->GetCellData();
      this->OutCellData = output->GetCellData();
      this->OutCellData->CopyAllocate(this->InCellData);
    }

    // Process point data if requested. Note the complication due to how the
    // new point data is interpolated from line-line intersections.
    this->PointInterpolation = ptInterpolation;
    this->MeshPtData = nullptr;
    this->LoopPtData = nullptr;
    this->OutPtData = nullptr;
    if (passPD)
    {
      vtkPointData* inPD = input->GetPointData();
      vtkPointData* loopPD = loops->GetPointData();
      // If trim edge interpoltion is specified, make sure that the trim loop point
      // attributes match the mesh point attributes. If not, point data is not passed
      // to the output.
      if (this->HaveEquivalentAttributes(inPD, loopPD))
      {
        this->MeshPtData = inPD;
        this->LoopPtData = loopPD;
      }
      else
      {
        vtkLog(WARNING,
          "Mesh and trim loop point data attributes are different, only common point data arrays "
          "will be processed");
        this->IntersectAttributes(inPD, loopPD);
      } // if point attributes identical
      this->OutPtData = output->GetPointData();
      this->OutPtData->CopyAllocate(this->MeshPtData);
    } // if passing point data
  }   // constructor

  // Copy cell data from the input to the output.
  void CopyCellData(vtkIdType inCellId, vtkIdType outCellId)
  {
    if (this->InCellData)
    {
      this->OutCellData->CopyData(this->InCellData, inCellId, outCellId);
    }
  }

  // Processing point data is tricky because two sets of point attributes may
  // be in play. Point which are passed through, can just use the
  // conventional vtkDataSetAttributes::CopyData() approach (copy from input
  // mesh to output mesh). However, intersection points may be interpolated
  // from the mesh edges, or from the trim loop edges. This means using the
  // correct point data to interpolate from (via InterpolateEdge()).

  // In some simple cases, a straightforward copying of point data from the
  // input mesh to the output mesh can be used.
  void CopyPointData(vtkIdType inPtId, vtkIdType outPtId)
  {
    if (this->MeshPtData)
    {
      this->OutPtData->CopyData(this->MeshPtData, inPtId, outPtId);
    }
  }

  // Because the cookie cutter algorithm processes cells on a cell-by-cell
  // basis, before processing each cell the point data management process
  // must be reset. Access to the attribute data is through the
  // sorted points list.
  void InitializePointAttributeCollection(SortedPointsType* polyPts, SortedPointsType* loopPts)
  {
    this->PolyPts = polyPts;
    this->LoopPts = loopPts;
    this->PointAttributes.clear();
  }

  // Add point attribute data to the vector of attributes. Note that there is
  // an implicit ordering of points consisted with SortPoints: the
  // SortPoint::Id refers to the position of the point attribute in the
  // PointAttributes vector.
  void AddPointAttribute(int attType, vtkIdType m0 = (-1), vtkIdType m1 = (-1), vtkIdType l0 = (-1),
    vtkIdType l1 = (-1), double tm = 0.0, double tl = 0.0)
  {
    // See if we are collecting point attribute data
    if (this->OutPtData == nullptr)
    {
      return;
    }

    // See if it's a vertex (CopyData)
    if (attType < PointAttribute::Intersection)
    {
      this->PointAttributes.emplace_back(attType, m0, m1, tm);
      return;
    }

    // See if it's an edge (InterpolateData)
    if (this->PointInterpolation == vtkCookieCutter::USE_MESH_EDGES)
    {
      this->PointAttributes.emplace_back(PointAttribute::MeshEdge, m0, m1, tm);
    }
    else // if ( this->PointInterpolation == vtkCookieCutter::USE_LOOP_EDGES )
    {
      this->PointAttributes.emplace_back(PointAttribute::LoopEdge, l0, l1, tl);
    }
  }

  // Interpolate or pass point data from either the input mesh or trim loop
  // edges, to the output mesh.
  void InterpolatePointData(vtkIdType ptId, vtkIdType outPtId)
  {
    // Make sure we actually want to process point data
    if (!this->OutPtData)
    {
      return;
    }

    // Okay, good to process, grab the point attributes.
    PointAttribute& ptAttr = this->PointAttributes[ptId];

    // Maybe just copying data from mesh or trim loop to output
    if (ptAttr.AttributeType == PointAttribute::MeshVertex)
    {
      this->OutPtData->CopyData(this->MeshPtData, ptAttr.V0, outPtId);
    }
    else if (ptAttr.AttributeType == PointAttribute::LoopVertex)
    {
      this->OutPtData->CopyData(this->LoopPtData, ptAttr.V0, outPtId);
    }

    // Otherwise interpolating along an edge
    else
    {
      if (ptAttr.AttributeType == PointAttribute::MeshEdge)
      {
        this->OutPtData->InterpolateEdge(this->MeshPtData, outPtId, ptAttr.V0, ptAttr.V1, ptAttr.T);
      }
      else // if ( ptAttr.AttributeType == PointAttribute::LoopEdge )
      {
        this->OutPtData->InterpolateEdge(this->LoopPtData, outPtId, ptAttr.V0, ptAttr.V1, ptAttr.T);
      }
    }
  }

}; // vtkAttributeManager

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

  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(
    vtkAttributeManager* attrMgr, vtkIdType ptId, 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;
  vtkPoints* OutPoints;
  vtkCellArray* OutLines;
  vtkCellArray* OutPolys;
  vtkAttributeManager* AttributeManager;
}; // vtkCookieCutterHelper

// vtkCookieCutterHelper methods follow
vtkCookieCutterHelper::vtkCookieCutterHelper(vtkIncrementalPointLocator* locator, vtkPoints* inPts,
  vtkPoints* outPts, vtkCellArray* outLines, vtkCellArray* outPolys, vtkAttributeManager* attMgr)
  : Locator(locator)
  , InPoints(inPts)
  , OutPoints(outPts)
  , OutLines(outLines)
  , OutPolys(outPolys)
  , AttributeManager(attMgr)
{
}

vtkCookieCutterHelper::~vtkCookieCutterHelper() = default;

void vtkCookieCutterHelper::InsertPoint(
  vtkAttributeManager* attrMgr, vtkIdType attrId, double x[3], vtkCellArray* ca)
{
  vtkIdType ptId;
  if (this->Locator->InsertUniquePoint(x, ptId))
  {
    attrMgr->InterpolatePointData(attrId, 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;

  // Begin gathering attribute information for this line.
  vtkAttributeManager* attrMgr = this->AttributeManager;
  attrMgr->InitializePointAttributeCollection(&sortedPoints, nullptr);

  // Insert all polyline points
  for (i = 0; i < npts; ++i)
  {
    t = static_cast<double>(i);
    this->InPoints->GetPoint(pts[i], x);
    sortedPoints.emplace_back(t, SortPoint::VERTEX, -1, -1, x);
    attrMgr->AddPointAttribute(PointAttribute::MeshVertex, pts[i]);
  }

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

    // Traverse polygon loop intersecting each polygon segment
    for (j = 0; j < numLoopPts; ++j)
    {
      vtkIdType l0 = j;
      vtkIdType l1 = (j + 1) % numLoopPts;
      loop->Points->GetPoint(l0, y0);
      loop->Points->GetPoint(l1, y1);
      if ((result = vtkLine::Intersection(x0, x1, y0, y1, u, v)) == vtkLine::Intersect)
      {
        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]);
        attrMgr->AddPointAttribute(PointAttribute::Intersection, m0, m1, l0, l1, u, v);
        u += static_cast<double>(i);
        v += static_cast<double>(j);
        sortedPoints.emplace_back(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.emplace_back(static_cast<double>(i), SortPoint::ON, numInts, onId, c);
            attrMgr->AddPointAttribute(PointAttribute::MeshVertex, m0);
            numInts++;
          }
          vtkLine::DistanceToLine(x1, y0, y1, u, c);
          if (-VTK_DEGENERATE_TOL <= u && u <= (1.0 + VTK_DEGENERATE_TOL))
          {
            sortedPoints.emplace_back(
              static_cast<double>(i) + 1.0, SortPoint::ON, numInts, onId, c);
            attrMgr->AddPointAttribute(PointAttribute::MeshVertex, m1);
            numInts++;
          }
          vtkLine::DistanceToLine(y0, x0, x1, u, c);
          if (-VTK_DEGENERATE_TOL <= u && u <= (1.0 + VTK_DEGENERATE_TOL))
          {
            sortedPoints.emplace_back(static_cast<double>(i) + u, SortPoint::ON, numInts, onId, c);
            attrMgr->AddPointAttribute(PointAttribute::LoopVertex, l0);
            numInts++;
          }
          vtkLine::DistanceToLine(y1, x0, x1, u, c);
          if (-VTK_DEGENERATE_TOL <= u && u <= (1.0 + VTK_DEGENERATE_TOL))
          {
            sortedPoints.emplace_back(static_cast<double>(i) + u, SortPoint::ON, numInts, onId, c);
            attrMgr->AddPointAttribute(PointAttribute::LoopVertex, l1);
            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->AttributeManager->CopyCellData(cellId, newCellId);
    for (i = startIdx; i <= endIdx; ++i)
    {
      InsertPoint(attrMgr, i, 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
  double polyBds[6];
  vtkBoundingBox bbox(poly->GetBounds());
  bbox.Inflate(0, 0, 0.01);
  bbox.GetBounds(polyBds);
  if (loopBds[0] > polyBds[1] || loopBds[1] < polyBds[0] || loopBds[2] > polyBds[3] ||
    loopBds[3] < polyBds[2] || loopBds[4] > polyBds[5] || loopBds[5] < polyBds[4])
  {
    return;
  }

  // Run around the polygon and loop 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.
  double* p = vtkDoubleArray::FastDownCast(poly->Points->GetData())->GetPointer(0);
  vtkIdType i, j, newCellId, numPts = 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();

  // Begin gathering attribute information for this polygon
  vtkAttributeManager* attrMgr = this->AttributeManager;
  attrMgr->InitializePointAttributeCollection(&loopPoints, &polyPoints);

  // The polygon points
  for (i = 0; i < npts; ++i)
  {
    t = static_cast<double>(i);
    this->InPoints->GetPoint(pts[i], x);
    polyPoints.emplace_back(t, SortPoint::VERTEX, numPts, -1, x);
    attrMgr->AddPointAttribute(PointAttribute::MeshVertex, pts[i]);
    numPts++;
  }

  // The loop points
  for (i = 0; i < numLoopPts; ++i)
  {
    t = static_cast<double>(i);
    loop->Points->GetPoint(i, x);
    loopPoints.emplace_back(t, SortPoint::VERTEX, numPts, -1, x);
    attrMgr->AddPointAttribute(PointAttribute::LoopVertex, i);
    numPts++;
  }

  // Now insert intersection points. Note that the attribute data is collected
  // with an implicit attribute id in lock step with numPts.
  for (i = 0; i < npts; ++i)
  {
    // Define an edge from the polygon (m0,m1)
    vtkIdType m0 = pts[i];
    vtkIdType m1 = pts[(i + 1) % npts];
    this->InPoints->GetPoint(m0, x0);
    this->InPoints->GetPoint(m1, x1);

    // Traverse polygon loop intersecting each polygon segment
    for (j = 0; j < numLoopPts; ++j)
    {
      // Define an edge from the loop (j,j+1)
      vtkIdType l0 = j;
      vtkIdType l1 = (j + 1) % numLoopPts;
      loop->Points->GetPoint(l0, y0);
      loop->Points->GetPoint(l1, y1);
      if ((result = vtkLine::Intersection(x0, x1, y0, y1, u, v)) == vtkLine::Intersect)
      {
        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]);
        attrMgr->AddPointAttribute(PointAttribute::Intersection, m0, m1, l0, l1, u, v);
        u += static_cast<double>(i);
        v += static_cast<double>(j);
        polyPoints.emplace_back(u, SortPoint::INTERSECTION, numPts, -1, x);
        loopPoints.emplace_back(v, SortPoint::INTERSECTION, numPts, -1, x);
        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.emplace_back(static_cast<double>(i), SortPoint::ON, numPts, onId, c);
            loopPoints.emplace_back(static_cast<double>(j) + u, SortPoint::ON, numPts, onId, c);
            attrMgr->AddPointAttribute(PointAttribute::MeshVertex, m0);
            numPts++;
          }
          vtkLine::DistanceToLine(x1, y0, y1, u, c);
          if (-VTK_DEGENERATE_TOL <= u && u <= (1.0 + VTK_DEGENERATE_TOL))
          {
            polyPoints.emplace_back(static_cast<double>(i) + 1.0, SortPoint::ON, numPts, onId, c);
            loopPoints.emplace_back(static_cast<double>(j) + u, SortPoint::ON, numPts, onId, c);
            attrMgr->AddPointAttribute(PointAttribute::MeshVertex, m1);
            numPts++;
          }
          vtkLine::DistanceToLine(y0, x0, x1, u, c);
          if (-VTK_DEGENERATE_TOL <= u && u <= (1.0 + VTK_DEGENERATE_TOL))
          {
            polyPoints.emplace_back(static_cast<double>(i) + u, SortPoint::ON, numPts, onId, c);
            loopPoints.emplace_back(static_cast<double>(j), SortPoint::ON, numPts, onId, c);
            attrMgr->AddPointAttribute(PointAttribute::LoopVertex, l0);
            numPts++;
          }
          vtkLine::DistanceToLine(y1, x0, x1, u, c);
          if (-VTK_DEGENERATE_TOL <= u && u <= (1.0 + VTK_DEGENERATE_TOL))
          {
            polyPoints.emplace_back(static_cast<double>(i) + u, SortPoint::ON, numPts, onId, c);
            loopPoints.emplace_back(static_cast<double>(j) + 1.0, SortPoint::ON, numPts, onId, c);
            attrMgr->AddPointAttribute(PointAttribute::LoopVertex, l1);
            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);
  pData->EditableOn();

  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 assumed 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. The variable
  // npts is the number of points, mesh+loop+intersection points, inserted
  // into the loops.
  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->AttributeManager->CopyCellData(cellId, newCellId);

      do
      {
        visited[nextId] = 1;
        numInsertedPts++;
        InsertPoint(attrMgr, nextId, 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->PassCellData = true;
  this->PassPointData = true;
  this->PointInterpolation = USE_MESH_EDGES;

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

  // Process loops from second input. Make sure valid trim loops
  // are available.
  vtkSmartPointer<vtkCellArray> loopPolys = BuildTrimLoops(loops);
  if (loopPolys == nullptr)
  {
    vtkErrorMacro("No trim lines/polygons defined");
    return 1;
  }

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

  // This class assists with the copying and interpolation of point and cell
  // data. While cell data is easy to process, point data can be complex
  // because of the merging of coincident points, and the interpolation of
  // data across mesh edges and/or tril loop edges.
  vtkAttributeManager attrMgr(
    input, loops, output, this->PassCellData, this->PassPointData, this->PointInterpolation);

  vtkCellArray* inVerts = input->GetVerts();
  vtkNew<vtkCellArray> outVerts;
  outVerts->AllocateCopy(inVerts);

  vtkCellArray* inLines = input->GetLines();
  vtkNew<vtkCellArray> outLines;
  outLines->AllocateCopy(inLines);

  vtkCellArray* inPolys = input->GetPolys();
  vtkCellArray* inStrips = input->GetStrips();
  vtkNew<vtkCellArray> outPolys;
  outPolys->AllocateCopy(inPolys);

  // Locator used to merge potentially duplicate points. Make sure that
  // the bounding box provides a little breathing room for proper
  // bounds-bounds intersection checks.
  vtkBoundingBox bbox(loops->GetBounds());
  double tol = 0.001 * bbox.GetDiagonalLength();
  bbox.Inflate(tol);
  if (this->Locator == nullptr)
  {
    this->CreateDefaultLocator();
  }
  double bds[6];
  bbox.GetBounds(bds);
  this->Locator->InitPointInsertion(outPts, bds);

  // Initialize and create polygon representing the loop
  double n[3];
  vtkNew<vtkPolygon> loop;
  vtkNew<vtkPolygon> poly;

  // Setup the lines and polys cropping helper
  std::unique_ptr<vtkCookieCutterHelper> helper(
    new vtkCookieCutterHelper(this->Locator, inPts, outPts, outLines, outPolys, &attrMgr));

  // For each trim loop, crop the input mesh cells. Cells are cropped
  // repeatedly one-by-one in a loop.
  vtkIdType npts;
  const vtkIdType* pts;
  for (loopPolys->InitTraversal(); loopPolys->GetNextCell(npts, pts);)
  {
    if (this->CheckAbort())
    {
      break;
    }
    vtkIdType cellId = 0;
    vtkIdType numLoopPts = npts;
    loop->Initialize(npts, pts, loopPts);
    bbox.SetBounds(loop->GetBounds());
    bbox.Inflate(tol); // Handle exactly aligned x-y plane numerical fuzz
    bbox.GetBounds(bds);
    vtkPolygon::ComputeNormal(loop->Points, n);
    double* l = vtkDoubleArray::FastDownCast(loop->Points->GetData())->GetPointer(0);

    // Start by processing the verts. A simple in/out check.
    if (inVerts->GetNumberOfCells() > 0)
    {
      for (inVerts->InitTraversal(); inVerts->GetNextCell(npts, pts); ++cellId)
      {
        for (vtkIdType i = 0; i < npts; ++i)
        {
          double x[3];
          inPts->GetPoint(pts[i], x);
          if (vtkPolygon::PointInPolygon(x, numLoopPts, l, bds, n) == 1)
          {
            vtkIdType ptId;
            this->Locator->InsertUniquePoint(x, ptId);
            attrMgr.CopyPointData(pts[i], ptId);
            vtkIdType newCellId = outVerts->InsertNextCell(1, &ptId);
            attrMgr.CopyCellData(cellId, newCellId);
          } // if point inside trim loop
        }   // for all points composing vert cell
      }     // for all vert cells
    }       // 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

  // Assign output as appropriate
  output->SetPoints(outPts);
  output->SetVerts(outVerts);
  output->SetLines(outLines);
  output->SetPolys(outPolys);

  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::CreateDefaultLocator()
{
  if (this->Locator == nullptr)
  {
    this->Locator = vtkMergePoints::New();
    this->Locator->Register(this);
    this->Locator->Delete();
  }
}

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

  os << indent << "Pass Cell Data: " << (this->PassCellData ? "On\n" : "Off\n");
  os << indent << "Pass Point Data: " << (this->PassPointData ? "On\n" : "Off\n");
  os << indent << "Point Interpolation: " << this->PointInterpolation << "\n";

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