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

  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 "vtkCellArray.h"
#include "vtkExecutive.h"
#include "vtkGarbageCollector.h"
#include "vtkLine.h"
#include "vtkMath.h"
#include "vtkObjectFactory.h"
#include "vtkCellData.h"
#include "vtkPointData.h"
#include "vtkPoints.h"
#include "vtkPolyData.h"
#include "vtkPolygon.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkNew.h"

#include <vector>
#include <set>

vtkStandardNewMacro(vtkCookieCutter);

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

    return;
  }

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

  // Convenience method------------------------------------------------------
  inline void InsertPoint(double x[3], vtkPoints *outPts, vtkCellArray *ca)
  {
    vtkIdType newPtId = outPts->InsertNextPoint(x);
    ca->InsertCellPoint(newPtId);
  }

  // Process a polyline-------------------------------------------------------
  void CropLine(vtkIdType cellId, vtkIdType cellOffset, vtkIdType npts,
                vtkIdType *pts, vtkPoints *inPts,
                vtkPolygon *loop, double *l, double loopBds[6], double n[3],
                vtkCellData *inCellData, vtkPoints *outPts,
                vtkCellArray *outLines, vtkCellData *outCellData)
  {
    // 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);
      inPts->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)
    {
      inPts->GetPoint(pts[i],x0);
      inPts->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 = outLines->InsertNextCell(numInsertedPts) + cellOffset;
      outCellData->CopyData(inCellData,cellId,newCellId);
      for (i=startIdx; i <= endIdx; ++i)
      {
        InsertPoint(sortedPoints[i].X, outPts, outLines);
      }
      startIdx = endIdx;
    }//over all sorted points
  }//CropLine

  // 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)
  {
    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;
    unsigned short 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;
  }

  // Process a polygon-------------------------------------------------------
  void CropPoly(vtkIdType cellId, vtkIdType cellOffset,
                vtkPolygon *poly, vtkIdType npts,
                vtkIdType *pts, vtkPoints *inPts, vtkPolygon *loop,
                double *l, double loopBds[6], double n[3],
                vtkCellData *inCellData, vtkPoints *outPts,
                vtkCellArray *outPolys, vtkCellData *outCellData)
  {
    // Make sure that this is a valid polygon
    if ( npts < 3 )
    {
      return;
    }

    // Make sure that the polyons 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);
      inPts->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)
    {
      inPts->GetPoint(pts[i],x0);
      inPts->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
    // everythig 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.GetPointer());
    pData->SetLines(pDataLines.GetPointer());

    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 ( ! ResolveTopology(pData.GetPointer()) )
    {
      return;
    }

    // Build loops (i.e., output polygons)
    vtkIdType numInsertedPts, *cells, thisCell, nextId, startId;
    unsigned short ncells;
    vtkIdType thisNPts, *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 = outPolys->InsertNextCell(numInsertedPts) + cellOffset;
        outCellData->CopyData(inCellData,cellId,newCellId);

        do
        {
          visited[nextId] = 1;
          numInsertedPts++;
          InsertPoint(pDataPts->GetPoint(nextId), outPts, 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 );
        outPolys->UpdateCellCount(numInsertedPts);
      }
    }
  }//That was easy! CropPoly

} //anonymous namespace

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

//----------------------------------------------------------------------------
vtkCookieCutter::~vtkCookieCutter()
{
}

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

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

  vtkIdType numPts, numLoopPts;
  if ( (numPts = 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->Allocate(inVerts->GetNumberOfCells(),1);

  vtkCellArray *inLines = input->GetLines();
  vtkCellArray *outLines = vtkCellArray::New();
  outLines->Allocate(inLines->GetNumberOfCells(),1);

  vtkCellArray *inPolys = input->GetPolys();
  vtkCellArray *inStrips = input->GetStrips();
  vtkCellArray *outPolys = vtkCellArray::New();
  outPolys->Allocate(inPolys->GetNumberOfCells(),1);

  // 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 i, npts, *pts, cellId=0, newPtId, newCellId, cellOffset=0;
  vtkCellArray *loopPolys = loops->GetPolys();
  for (loopPolys->InitTraversal(); loopPolys->GetNextCell(npts,pts); )
  {
    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 ( i=0; i<npts; ++i)
        {
          inPts->GetPoint(pts[i], x);
          if ( vtkPolygon::PointInPolygon(x, numLoopPts, l, bds, n) == 1 )
          {
            newPtId = outPts->InsertNextPoint(x);
            newCellId = outVerts->InsertNextCell(1,&newPtId);
            outCellData->CopyData(inCellData,cellId,newCellId);
          }
        }
      }//for all verts
    }//if vert cells

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

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

    // Now process triangle strips
    if ( inStrips->GetNumberOfCells() > 0 )
    {
      cellOffset = outVerts->GetNumberOfCells() + outLines->GetNumberOfCells();
      for (inStrips->InitTraversal(); inStrips->GetNextCell(npts,pts); ++cellId)
      {
        vtkIdType numTriPts=3, triPts[3];
        for ( 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);
          CropPoly(cellId,cellOffset, poly,numTriPts,triPts,inPts, loop,l,bds,n,inCellData,
                   outPts,outPolys,outCellData);
        }
      }
    }//if polygonal cells
  }//for all loops

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

}