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

  Program:   Visualization Toolkit
  Module:    $RCSfile: vtkDelaunay2D.cxx,v $

  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 "vtkDelaunay2D.h"

#include "vtkAbstractTransform.h"
#include "vtkCellArray.h"
#include "vtkDoubleArray.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkMath.h"
#include "vtkObjectFactory.h"
#include "vtkPlane.h"
#include "vtkPointData.h"
#include "vtkPolyData.h"
#include "vtkPolygon.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkTriangle.h"
#include "vtkTransform.h"

vtkCxxRevisionMacro(vtkDelaunay2D, "$Revision: 1.68.10.1 $");
vtkStandardNewMacro(vtkDelaunay2D);
vtkCxxSetObjectMacro(vtkDelaunay2D,Transform,vtkAbstractTransform);

// Construct object with Alpha = 0.0; Tolerance = 0.00001; Offset = 1.25;
// BoundingTriangulation turned off.
vtkDelaunay2D::vtkDelaunay2D()
{
  this->Alpha = 0.0;
  this->Tolerance = 0.00001;
  this->BoundingTriangulation = 0;
  this->Offset = 1.0;
  this->Transform = NULL;
  this->ProjectionPlaneMode = VTK_DELAUNAY_XY_PLANE;

  // optional 2nd input
  this->SetNumberOfInputPorts(2);
}

vtkDelaunay2D::~vtkDelaunay2D()
{
  if (this->Transform)
    {
    this->Transform->UnRegister(this);
    }
}

//----------------------------------------------------------------------------
// Specify the input data or filter. Old style.
void vtkDelaunay2D::SetSource(vtkPolyData *input)
{
  this->Superclass::SetInput(1, input);
}

//----------------------------------------------------------------------------
// Specify the input data or filter. New style.
void vtkDelaunay2D::SetSourceConnection(vtkAlgorithmOutput *algOutput)
{
  this->Superclass::SetInputConnection(1, algOutput);
}

vtkPolyData *vtkDelaunay2D::GetSource()
{
  if (this->GetNumberOfInputConnections(1) < 1)
    {
    return NULL;
    }
  return vtkPolyData::SafeDownCast(
    this->GetExecutive()->GetInputData(1, 0));
}

// Determine whether point x is inside of circumcircle of triangle
// defined by points (x1, x2, x3). Returns non-zero if inside circle.
// (Note that z-component is ignored.)
int vtkDelaunay2D::InCircle (double x[3], double x1[3], double x2[3], 
                             double x3[3])
{
  double radius2, center[2], dist2;

  radius2 = vtkTriangle::Circumcircle(x1,x2,x3,center);

  // check if inside/outside circumcircle
  dist2 = (x[0]-center[0]) * (x[0]-center[0]) + 
          (x[1]-center[1]) * (x[1]-center[1]);

  if ( dist2 < (0.999999999999*radius2) )
    {
    return 1;
    }
  else
    {
    return 0;
    }
}

#define VTK_DEL2D_TOLERANCE 1.0e-014

// Recursive method to locate triangle containing point. Starts with arbitrary
// triangle (tri) and "walks" towards it. Influenced by some of Guibas and 
// Stolfi's work. Returns id of enclosing triangle, or -1 if no triangle
// found. Also, the array nei[3] is used to communicate info about points
// that lie on triangle edges: nei[0] is neighboring triangle id, and nei[1]
// and nei[2] are the vertices defining the edge.
vtkIdType vtkDelaunay2D::FindTriangle(double x[3], vtkIdType ptIds[3],
                                      vtkIdType tri, double tol,
                                      vtkIdType nei[3], vtkIdList *neighbors)
{
  int i, j, ir, ic, inside, i2, i3;
  vtkIdType *pts, npts, newNei;
  double p[3][3], n[2], vp[2], vx[2], dp, minProj;

  // get local triangle info
  this->Mesh->GetCellPoints(tri,npts,pts);
  for (i=0; i<3; i++) 
    {
    ptIds[i] = pts[i];
    this->GetPoint(ptIds[i], p[i]);
    }

  // Randomization (of find edge neighbora) avoids walking in 
  // circles in certain weird cases
  srand(tri);
  ir = rand() % 3;
  // evaluate in/out of each edge
  for (inside=1, minProj=0.0, ic=0; ic<3; ic++)
    {
    i  = (ir+ic) % 3;
    i2 = (i+1) % 3;
    i3 = (i+2) % 3;

    // create a 2D edge normal to define a "half-space"; evaluate points (i.e.,
    // candiate point and other triangle vertex not on this edge).
    n[0] = -(p[i2][1] - p[i][1]);
    n[1] = p[i2][0] - p[i][0];
    vtkMath::Normalize2D(n);

    // compute local vectors
    for (j=0; j<2; j++)
      {
      vp[j] = p[i3][j] - p[i][j];
      vx[j] = x[j] - p[i][j];
      }

    //check for duplicate point
    vtkMath::Normalize2D(vp);
    if ( vtkMath::Normalize2D(vx) <= tol ) 
      {
      this->NumberOfDuplicatePoints++;
      return -1;
      }

    // see if two points are in opposite half spaces
    dp = vtkMath::Dot2D(n,vx) * (vtkMath::Dot2D(n,vp) < 0 ? -1.0 : 1.0);
    if ( dp < VTK_DEL2D_TOLERANCE )
      {
      if ( dp < minProj ) //track edge most orthogonal to point direction
        {
        inside = 0;
        nei[1] = ptIds[i];
        nei[2] = ptIds[i2];
        minProj = dp;
        }
      }//outside this edge
    }//for each edge

  if ( inside ) // all edges have tested positive
    {
    nei[0] = (-1);
    return tri;
    }

  else if ( !inside && (fabs(minProj) < VTK_DEL2D_TOLERANCE) ) // on edge
    {
    this->Mesh->GetCellEdgeNeighbors(tri,nei[1],nei[2],neighbors);
    nei[0] = neighbors->GetId(0);
    return tri;
    }

  else //walk towards point
    {
    this->Mesh->GetCellEdgeNeighbors(tri,nei[1],nei[2],neighbors);
    if ( (newNei=neighbors->GetId(0)) == nei[0] )
      {
      this->NumberOfDegeneracies++;
      return -1;
      }
    else
      {
      nei[0] = tri;
      return this->FindTriangle(x,ptIds,newNei,tol,nei,neighbors);
      }
    }
}

#undef VTK_DEL2D_TOLERANCE

// Recursive method checks whether edge is Delaunay, and if not, swaps edge.
// Continues until all edges are Delaunay. Points p1 and p2 form the edge in
// question; x is the coordinates of the inserted point; tri is the current
// triangle id.
void vtkDelaunay2D::CheckEdge(vtkIdType ptId, double x[3], vtkIdType p1,
                              vtkIdType p2, vtkIdType tri)
{
  int i;
  vtkIdType *pts, npts, numNei, nei, p3;
  double x1[3], x2[3], x3[3];
  vtkIdList *neighbors;
  vtkIdType swapTri[3];

  this->GetPoint(p1,x1);
  this->GetPoint(p2,x2);

  neighbors = vtkIdList::New();
  neighbors->Allocate(2);

  this->Mesh->GetCellEdgeNeighbors(tri,p1,p2,neighbors);
  numNei = neighbors->GetNumberOfIds();

  if ( numNei > 0 ) //i.e., not a boundary edge
    {
    // get neighbor info including opposite point
    nei = neighbors->GetId(0);
    this->Mesh->GetCellPoints(nei, npts, pts);
    for (i=0; i<2; i++)
      {
      if ( pts[i] != p1 && pts[i] != p2 )
        {
        break;
        }
      }
    p3 = pts[i];
    this->GetPoint(p3,x3);

    // see whether point is in circumcircle
    if ( this->InCircle (x3, x, x1, x2) )
      {// swap diagonal
      this->Mesh->RemoveReferenceToCell(p1,tri);
      this->Mesh->RemoveReferenceToCell(p2,nei);
      this->Mesh->ResizeCellList(ptId,1);
      this->Mesh->AddReferenceToCell(ptId,nei);
      this->Mesh->ResizeCellList(p3,1);
      this->Mesh->AddReferenceToCell(p3,tri);

      swapTri[0] = ptId; swapTri[1] = p3; swapTri[2] = p2;
      this->Mesh->ReplaceCell(tri,3,swapTri);

      swapTri[0] = ptId; swapTri[1] = p1; swapTri[2] = p3;
      this->Mesh->ReplaceCell(nei,3,swapTri);

      // two new edges become suspect
      this->CheckEdge(ptId, x, p3, p2, tri);
      this->CheckEdge(ptId, x, p1, p3, nei);

      }//in circle
    }//interior edge

  neighbors->Delete();
}

// 2D Delaunay triangulation. Steps are as follows:
//   1. For each point
//   2. Find triangle point is in
//   3. Create 3 triangles from each edge of triangle that point is in
//   4. Recursively evaluate Delaunay criterion for each edge neighbor
//   5. If criterion not satisfied; swap diagonal
// 
int vtkDelaunay2D::RequestData(
  vtkInformation *vtkNotUsed(request),
  vtkInformationVector **inputVector,
  vtkInformationVector *outputVector)
{
  // get the info objects
  vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation *sourceInfo = inputVector[1]->GetInformationObject(0);
  vtkInformation *outInfo = outputVector->GetInformationObject(0);

  // get the input and ouptut
  vtkPointSet *input = vtkPointSet::SafeDownCast(
    inInfo->Get(vtkDataObject::DATA_OBJECT()));
  vtkPolyData *source = 0;
  if (sourceInfo)
    {
    source =
      vtkPolyData::SafeDownCast(sourceInfo->Get(vtkDataObject::DATA_OBJECT()));
    }
  vtkPolyData *output = vtkPolyData::SafeDownCast(
    outInfo->Get(vtkDataObject::DATA_OBJECT()));

  vtkIdType numPoints, i;
  vtkIdType numTriangles = 0;
  vtkIdType ptId, tri[4], nei[3];
  vtkIdType p1 = 0;
  vtkIdType p2 = 0;
  vtkPoints *inPoints;
  vtkPoints *points;
  vtkPoints *tPoints = NULL;
  vtkCellArray *triangles;
  int ncells;
  vtkIdType nodes[4][3], *neiPts;
  vtkIdType *triPts = 0;
  vtkIdType numNeiPts;
  vtkIdType npts = 0;
  vtkIdType pts[3];
  vtkIdList *neighbors, *cells;
  double center[3], radius, tol, x[3];
  int *triUse = NULL;
  double *bounds;

  vtkDebugMacro(<<"Generating 2D Delaunay triangulation");

  if (this->Transform && this->BoundingTriangulation)
    {
    vtkWarningMacro(<<"Bounding triangulation cannot be used when an input transform is specified.  Output will not contain bounding triangulation.");
    }
  
  if (this->ProjectionPlaneMode == VTK_BEST_FITTING_PLANE && this->BoundingTriangulation)
    {
    vtkWarningMacro(<<"Bounding triangulation cannot be used when the best fitting plane option is on.  Output will not contain bounding triangulation.");
    }
  
  // Initialize; check input
  //
  if ( (inPoints=input->GetPoints()) == NULL )
    {
    vtkErrorMacro("<<Cannot triangulate; no input points");
    return 1;
    }

  if ( (numPoints=inPoints->GetNumberOfPoints()) <= 2 )
    {
    vtkErrorMacro("<<Cannot triangulate; need at least 3 input points");
    return 1;
    }
  
  neighbors = vtkIdList::New(); neighbors->Allocate(2);
  cells = vtkIdList::New(); cells->Allocate(64);
  
  this->NumberOfDuplicatePoints = 0;
  this->NumberOfDegeneracies = 0;

  this->Mesh = vtkPolyData::New();
  

  // If the user specified a transform, apply it to the input data.
  //
  // Only the input points are transformed.  We do not bother
  // transforming the source points (if specified).  The reason is
  // that only the topology of the Source is used during the constrain
  // operation.  The point ids in the Source topology are assumed to
  // reference points in the input. So, when an input transform is
  // used, only the input points are transformed.  We do not bother
  // with transforming the Source points since they are never
  // referenced.
  if (this->Transform)
    {
    tPoints = vtkPoints::New();
    this->Transform->TransformPoints(inPoints, tPoints);
    }
  else
    {
    // If the user asked this filter to compute the best fitting plane,
    // proceed to compute the plane and generate a transform that will
    // map the input points into that plane.
    if(this->ProjectionPlaneMode == VTK_BEST_FITTING_PLANE)
      {
      this->SetTransform( this->ComputeBestFittingPlane(input) );
      tPoints = vtkPoints::New();
      this->Transform->TransformPoints(inPoints, tPoints);
      }
    }
  
  // Create initial bounding triangulation. Have to create bounding points.
  // Initialize mesh structure.
  //
  points = vtkPoints::New(); 
  // This will copy doubles to doubles if the input is double. 
  points->SetDataTypeToDouble();
  points->SetNumberOfPoints(numPoints);
  if (!this->Transform)
    {
    points->DeepCopy(inPoints);
    }
  else
    {
    points->DeepCopy(tPoints);
    tPoints->Delete();
    tPoints = NULL;
    }

  bounds = points->GetBounds();
  center[0] = (bounds[0]+bounds[1])/2.0;
  center[1] = (bounds[2]+bounds[3])/2.0;
  center[2] = (bounds[4]+bounds[5])/2.0;
  tol = input->GetLength();
  radius = this->Offset * tol;
  tol *= this->Tolerance;

  for (ptId=0; ptId<8; ptId++)
    {
    x[0] = center[0]
      + radius*cos((double)(45.0*ptId)*vtkMath::DegreesToRadians());
    x[1] = center[1]
      + radius*sin((double)(45.0*ptId)*vtkMath::DegreesToRadians());
    x[2] = center[2];
    points->InsertPoint(numPoints+ptId,x);
    }
  // We do this for speed accessing points
  this->Points = ((vtkDoubleArray *)points->GetData())->GetPointer(0);

  triangles = vtkCellArray::New();
  triangles->Allocate(triangles->EstimateSize(2*numPoints,3));

  //create bounding triangles (there are six)
  pts[0] = numPoints; pts[1] = numPoints + 1; pts[2] = numPoints + 2;
  triangles->InsertNextCell(3,pts);
  pts[0] = numPoints + 2; pts[1] = numPoints + 3; pts[2] = numPoints + 4;
  triangles->InsertNextCell(3,pts);
  pts[0] = numPoints + 4; pts[1] = numPoints + 5; pts[2] = numPoints + 6;
  triangles->InsertNextCell(3,pts);
  pts[0] = numPoints + 6; pts[1] = numPoints + 7; pts[2] = numPoints + 0;
  triangles->InsertNextCell(3,pts);
  pts[0] = numPoints + 0; pts[1] = numPoints + 2; pts[2] = numPoints + 6;
  triangles->InsertNextCell(3,pts);
  pts[0] = numPoints + 2; pts[1] = numPoints + 4; pts[2] = numPoints + 6;
  triangles->InsertNextCell(3,pts);
  tri[0] = 0;

  this->Mesh->SetPoints(points);
  this->Mesh->SetPolys(triangles);
  this->Mesh->BuildLinks(); //build cell structure

  // For each point; find triangle containing point. Then evaluate three 
  // neighboring triangles for Delaunay criterion. Triangles that do not 
  // satisfy criterion have their edges swapped. This continues recursively 
  // until all triangles have been shown to be Delaunay.
  //
  for (ptId=0; ptId < numPoints; ptId++)
    {
    this->GetPoint(ptId,x); 
    nei[0] = (-1); //where we are coming from...nowhere initially

    if ( (tri[0] = this->FindTriangle(x,pts,tri[0],tol,nei,neighbors)) >= 0 )
      {
      if ( nei[0] < 0 ) //in triangle
        {
        //delete this triangle; create three new triangles
        //first triangle is replaced with one of the new ones
        nodes[0][0] = ptId; nodes[0][1] = pts[0]; nodes[0][2] = pts[1];
        this->Mesh->RemoveReferenceToCell(pts[2], tri[0]);
        this->Mesh->ReplaceCell(tri[0], 3, nodes[0]);
        this->Mesh->ResizeCellList(ptId,1);
        this->Mesh->AddReferenceToCell(ptId,tri[0]);

        //create two new triangles
        nodes[1][0] = ptId; nodes[1][1] = pts[1]; nodes[1][2] = pts[2];
        tri[1] = this->Mesh->InsertNextLinkedCell(VTK_TRIANGLE, 3, nodes[1]);

        nodes[2][0] = ptId; nodes[2][1] = pts[2]; nodes[2][2] = pts[0];
        tri[2] = this->Mesh->InsertNextLinkedCell(VTK_TRIANGLE, 3, nodes[2]);

        // Check edge neighbors for Delaunay criterion. If not satisfied, flip
        // edge diagonal. (This is done recursively.)
        this->CheckEdge(ptId, x, pts[0], pts[1], tri[0]);
        this->CheckEdge(ptId, x, pts[1], pts[2], tri[1]);
        this->CheckEdge(ptId, x, pts[2], pts[0], tri[2]);
        }

      else // on triangle edge
        {
        //update cell list
        this->Mesh->GetCellPoints(nei[0],numNeiPts,neiPts);
        for (i=0; i<3; i++)
          {
          if ( neiPts[i] != nei[1] && neiPts[i] != nei[2] ) 
            {
            p1 = neiPts[i];
            }
          if ( pts[i] != nei[1] && pts[i] != nei[2] ) 
            {
            p2 = pts[i];
            }
          }
        this->Mesh->ResizeCellList(p1,1);
        this->Mesh->ResizeCellList(p2,1);

        //replace two triangles
        this->Mesh->RemoveReferenceToCell(nei[2],tri[0]);
        this->Mesh->RemoveReferenceToCell(nei[2],nei[0]);
        nodes[0][0] = ptId; nodes[0][1] = p2; nodes[0][2] = nei[1];
        this->Mesh->ReplaceCell(tri[0], 3, nodes[0]);
        nodes[1][0] = ptId; nodes[1][1] = p1; nodes[1][2] = nei[1];
        this->Mesh->ReplaceCell(nei[0], 3, nodes[1]);
        this->Mesh->ResizeCellList(ptId, 2);
        this->Mesh->AddReferenceToCell(ptId,tri[0]);
        this->Mesh->AddReferenceToCell(ptId,nei[0]);

        tri[1] = nei[0];

        //create two new triangles
        nodes[2][0] = ptId; nodes[2][1] = p2; nodes[2][2] = nei[2];
        tri[2] = this->Mesh->InsertNextLinkedCell(VTK_TRIANGLE, 3, nodes[2]);

        nodes[3][0] = ptId; nodes[3][1] = p1; nodes[3][2] = nei[2];
        tri[3] = this->Mesh->InsertNextLinkedCell(VTK_TRIANGLE, 3, nodes[3]);

        // Check edge neighbors for Delaunay criterion.
        for ( i=0; i<4; i++ )
          {
          this->CheckEdge (ptId, x, nodes[i][1], nodes[i][2], tri[i]);
          }
        }
      }//if triangle found

    else
      {
      tri[0] = 0; //no triangle found
      }

    if ( ! (ptId % 1000) ) 
      {
      vtkDebugMacro(<<"point #" << ptId);
      this->UpdateProgress ((double)ptId/numPoints);
      if (this->GetAbortExecute()) 
        {
        break;
        }
      }
    
    }//for all points

  vtkDebugMacro(<<"Triangulated " << numPoints <<" points, " 
                << this->NumberOfDuplicatePoints
                << " of which were duplicates");

  if ( this->NumberOfDegeneracies > 0 )
    {
    vtkDebugMacro(<< this->NumberOfDegeneracies 
    << " degenerate triangles encountered, mesh quality suspect");
    }

  // Finish up by recovering the boundary, or deleting all triangles connected 
  // to the bounding triangulation points or not satisfying alpha criterion,
  if ( !this->BoundingTriangulation || this->Alpha > 0.0 || source )
    {
    numTriangles = this->Mesh->GetNumberOfCells();
    if ( source ) 
      {
      triUse = this->RecoverBoundary(source);
      }
    else
      {
      triUse = new int[numTriangles];
      for (i=0; i<numTriangles; i++) 
        {
        triUse[i] = 1;
        }
      }
    }

  // Delete triangles connected to boundary points (if not desired)
  if ( ! this->BoundingTriangulation )
    {
    for (ptId=numPoints; ptId < (numPoints+8); ptId++)
      {
      this->Mesh->GetPointCells(ptId, cells);
      ncells = cells->GetNumberOfIds();
      for (i=0; i < ncells; i++)
        {
        triUse[cells->GetId(i)] = 0; //mark as deleted
        }
      }
    }

  // If non-zero alpha value, then figure out which parts of mesh are
  // contained within alpha radius.
  //
  if ( this->Alpha > 0.0 )
    {
    double alpha2 = this->Alpha * this->Alpha;
    double x1[3], x2[3], x3[3];
    double xx1[3], xx2[3], xx3[3];
    vtkIdType cellId, numNei, ap1, ap2, neighbor;

    vtkCellArray *alphaVerts = vtkCellArray::New();
    alphaVerts->Allocate(numPoints);
    vtkCellArray *alphaLines = vtkCellArray::New();
    alphaLines->Allocate(numPoints);

    char *pointUse = new char[numPoints+8];
    for (ptId=0; ptId < (numPoints+8); ptId++)
      {
      pointUse[ptId] = 0;
      }

    //traverse all triangles; evaluating Delaunay criterion
    for (i=0; i < numTriangles; i++)
      {
      if ( triUse[i] == 1 )
        {
        this->Mesh->GetCellPoints(i, npts, triPts);

        // if any point is one of the bounding points that was added
        // at the beginning of the algorithm, then grab the points
        // from the variable "points" (this list has the boundary
        // points and the original points have been transformed by the
        // input transform).  if none of the points are bounding points,
        // then grab the points from the variable "inPoints" so the alpha
        // criterion is applied in the nontransformed space.
        if (triPts[0]<numPoints && triPts[1]<numPoints && triPts[2]<numPoints)
          {
          inPoints->GetPoint(triPts[0],x1);
          inPoints->GetPoint(triPts[1],x2);
          inPoints->GetPoint(triPts[2],x3);
          }
        else
          {
          points->GetPoint(triPts[0],x1);
          points->GetPoint(triPts[1],x2);
          points->GetPoint(triPts[2],x3);
          }

        // evaluate the alpha criterion in 3D
        vtkTriangle::ProjectTo2D(x1, x2, x3, xx1, xx2, xx3);
        if ( vtkTriangle::Circumcircle(xx1,xx2,xx3,center) > alpha2 )
          {
          triUse[i] = 0;
          }
        else
          {
          for (int j=0; j<3; j++)
            {
            pointUse[triPts[j]] = 1; 
            }
          }
        }//if non-deleted triangle
      }//for all triangles

    //traverse all edges see whether we need to create some
    for (cellId=0, triangles->InitTraversal(); 
    triangles->GetNextCell(npts,triPts); cellId++)
      {
      if ( ! triUse[cellId] )
        {
        for (i=0; i < npts; i++) 
          {
          ap1 = triPts[i];
          ap2 = triPts[(i+1)%npts];

          if (this->BoundingTriangulation || (ap1<numPoints && ap2<numPoints))
            {
            this->Mesh->GetCellEdgeNeighbors(cellId,ap1,ap2,neighbors);
            numNei = neighbors->GetNumberOfIds();

            if ( numNei < 1 || ((neighbor=neighbors->GetId(0)) > cellId 
                                && !triUse[neighbor]) )
              {//see whether edge is shorter than Alpha

              // same argument as above, if one is a boundary point, get
              // it using this->GetPoint() which are transformed points. if
              // neither of the points are boundary points, get the from
              // inPoints (untransformed points) so alpha comparison done
              // untransformed space
              if (ap1 < numPoints && ap2 < numPoints)
                {
                inPoints->GetPoint(ap1,x1);
                inPoints->GetPoint(ap2,x2);
                }
              else
                {
                this->GetPoint(ap1,x1);
                this->GetPoint(ap2,x2);
                }
              if ( (vtkMath::Distance2BetweenPoints(x1,x2)*0.25) <= alpha2 )
                {
                pointUse[ap1] = 1; pointUse[ap2] = 1;
                pts[0] = ap1;
                pts[1] = ap2;
                alphaLines->InsertNextCell(2,pts);
                }//if passed test
              }//test edge
            }//if valid edge
          }//for all edges of this triangle
        }//if triangle not output
      }//for all triangles

    //traverse all points, create vertices if none used
    for (ptId=0; ptId<(numPoints+8); ptId++)
      {
      if ( !pointUse[ptId]
           &&  (ptId < numPoints || this->BoundingTriangulation) )
        {
        pts[0] = ptId;
        alphaVerts->InsertNextCell(1,pts);
        }
      }

    // update output
    delete [] pointUse;
    output->SetVerts(alphaVerts);
    alphaVerts->Delete();
    output->SetLines(alphaLines);
    alphaLines->Delete();
    }

  // Update output; free up supporting data structures.
  //
  if ( this->BoundingTriangulation && !this->Transform)
    {
    output->SetPoints(points);
    }  
  else
    {
    output->SetPoints(inPoints);
    output->GetPointData()->PassData(input->GetPointData());
    }

  if ( this->Alpha <= 0.0 && this->BoundingTriangulation && !source )
    {
    output->SetPolys(triangles);
    }
  else
    {
    vtkCellArray *alphaTriangles = vtkCellArray::New();
    alphaTriangles->Allocate(numTriangles);
    vtkIdType *alphaTriPts;

    for (i=0; i<numTriangles; i++)
      {
      if ( triUse[i] )
        {
        this->Mesh->GetCellPoints(i,npts,alphaTriPts);
        alphaTriangles->InsertNextCell(3,alphaTriPts);
        }
      }
    output->SetPolys(alphaTriangles);
    alphaTriangles->Delete();
    delete [] triUse;
    }

  points->Delete();
  triangles->Delete();
  this->Mesh->Delete();
  neighbors->Delete();
  cells->Delete();

  // If the best fitting option was ON, then the current transform
  // is the one that was computed internally. We must now destroy it.
  if( this->ProjectionPlaneMode == VTK_BEST_FITTING_PLANE )
    {
    if (this->Transform)
      {
      this->Transform->UnRegister(this);
      this->Transform = NULL;
      }
    }

  output->Squeeze();

  return 1;
}

// Methods used to recover edges. Uses lines and polygons to determine boundary
// and inside/outside.
//
// Only the topology of the Source is used during the constrain operation.
// The point ids in the Source topology are assumed to reference points
// in the input. So, when an input transform is used, only the input points
// are transformed.  We do not bother with transforming the Source points
// since they are never referenced.
int *vtkDelaunay2D::RecoverBoundary(vtkPolyData *source)
{
  vtkCellArray *lines=source->GetLines();
  vtkCellArray *polys=source->GetPolys();
  vtkIdType *pts = 0;
  vtkIdType npts = 0;
  vtkIdType i, p1, p2;
  int *triUse;
  
  // Recover the edges of the mesh
  for ( lines->InitTraversal(); lines->GetNextCell(npts,pts); )
    {
    for (i=0; i<(npts-1); i++)
      {
      p1 = pts[i];
      p2 = pts[i+1];
      if ( ! this->Mesh->IsEdge(p1,p2) )
        {
        this->RecoverEdge(p1, p2);
        }
      }
    }

  // Recover the enclosed regions (polygons) of the mesh
  for ( polys->InitTraversal(); polys->GetNextCell(npts,pts); )
    {
    for (i=0; i<npts; i++)
      {
      p1 = pts[i];
      p2 = pts[(i+1)%npts];
      if ( ! this->Mesh->IsEdge(p1,p2) )
        {
        this->RecoverEdge(p1, p2);
        }
      }
    }
  
  // Generate inside/outside marks on mesh
  int numTriangles = this->Mesh->GetNumberOfCells();
  triUse = new int[numTriangles];
  for (i=0; i<numTriangles; i++) 
    {
    triUse[i] = 1;
    }
  
  // Use any polygons to mark inside and outside. (Note that if an edge was not
  // recovered, we're going to have a problem.) The first polygon is assumed to
  // define the outside of the polygon; additional polygons carve out inside
  // holes.
  this->FillPolygons(polys, triUse);

  return triUse;
}

// Method attempts to recover an edge by retriangulating mesh around the edge.
// What we do is identify a "submesh" of triangles that includes the edge to recover.
// Then we split the submesh in two with the recovered edge, and triangulate each of
// the two halves. If any part of this fails, we leave things alone.
int vtkDelaunay2D::RecoverEdge(vtkIdType p1, vtkIdType p2)
{
  vtkIdType cellId = 0;
  int i, j, k;
  double p1X[3], p2X[3], xyNormal[3], splitNormal[3], p21[3];
  double x1[3], x2[3], sepNormal[3], v21[3];
  int ncells, v1=0, v2=0, signX1=0, signX2, signP1, signP2;
  vtkIdType *pts, *leftTris, *rightTris, npts, numRightTris, numLeftTris;
  int success=0;

  vtkIdList *cells=vtkIdList::New(); cells->Allocate(64);  
  vtkIdList *tris=vtkIdList::New(); tris->Allocate(64);
  vtkPolygon *rightPoly=vtkPolygon::New();
  vtkPolygon *leftPoly=vtkPolygon::New();
  vtkIdList *leftChain=leftPoly->GetPointIds();
  vtkIdList *rightChain=rightPoly->GetPointIds();
  vtkPoints *leftChainX=leftPoly->GetPoints();
  vtkPoints *rightChainX=rightPoly->GetPoints();
  vtkIdList *neis=vtkIdList::New(); neis->Allocate(4);
  vtkIdList *rightPtIds=vtkIdList::New(); rightPtIds->Allocate(64);
  vtkIdList *leftPtIds=vtkIdList::New(); leftPtIds->Allocate(64);
  vtkPoints *rightTriPts=vtkPoints::New(); rightTriPts->Allocate(64);
  vtkPoints *leftTriPts=vtkPoints::New(); leftTriPts->Allocate(64);
  
  // Compute a split plane along (p1,p2) and parallel to the z-axis.
  //
  this->GetPoint(p1,p1X); p1X[2] = 0.0; //split plane point
  this->GetPoint(p2,p2X); p2X[2] = 0.0; //split plane point
  xyNormal[0] = xyNormal[1] = 0.0; xyNormal[2] = 1.0;
  for (i=0; i<3; i++ ) 
    {
    p21[i] = p2X[i] - p1X[i]; //working in x-y plane
    }

  vtkMath::Cross (p21,xyNormal,splitNormal);
  if ( vtkMath::Normalize(splitNormal) == 0.0 ) 
    {//Usually means coincident points
    goto FAILURE;
    }

  // Identify a triangle connected to the point p1 containing a portion of the edge.
  //
  this->Mesh->GetPointCells(p1, cells);
  ncells = cells->GetNumberOfIds();
  for (i=0; i < ncells; i++)
    {
    cellId = cells->GetId(i);
    this->Mesh->GetCellPoints(cellId, npts, pts);
    for (j=0; j<3; j++)
      {
      if ( pts[j] == p1 ) 
        {
        break;
        }
      }
    v1 = pts[(j+1)%3];
    v2 = pts[(j+2)%3];
    this->GetPoint(v1,x1); x1[2] = 0.0;
    this->GetPoint(v2,x2); x2[2] = 0.0;
    signX1 = (vtkPlane::Evaluate(splitNormal, p1X, x1) > 0.0 ? 1 : -1);
    signX2 = (vtkPlane::Evaluate(splitNormal, p1X, x2) > 0.0 ? 1 : -1);
    if ( signX1 != signX2 ) //points of triangle on either side of edge
      {
      // determine if edge separates p1 from p2 - then we've found triangle
      v21[0] = x2[0] - x1[0]; //working in x-y plane
      v21[1] = x2[1] - x1[1];
      v21[2] = 0.0;

      vtkMath::Cross (v21,xyNormal,sepNormal);
      if ( vtkMath::Normalize(sepNormal) == 0.0 )
        { //bad mesh
        goto FAILURE;
        }

      signP1 = (vtkPlane::Evaluate(sepNormal, x1, p1X) > 0.0 ? 1 : -1);
      signP2 = (vtkPlane::Evaluate(sepNormal, x1, p2X) > 0.0 ? 1 : -1);
      if ( signP1 != signP2 ) //is a separation line
        {
        break;
        }
      }
    } //for all cells

  if ( i >= ncells ) 
    {//something is really screwed up
    goto FAILURE;
    }

  // We found initial triangle; begin to track triangles containing edge. Also,
  // the triangle defines the beginning of two "chains" which form a boundary of
  // enclosing triangles around the edge. Create the two chains (from p1 to p2).
  // (The chains are actually defining two polygons on either side of the edge.)
  // 
  tris->InsertId(0, cellId);
  rightChain->InsertId(0, p1); rightChainX->InsertPoint(0, p1X);
  leftChain->InsertId(0, p1); leftChainX->InsertPoint(0, p1X);
  if ( signX1 > 0 )
    {
    rightChain->InsertId(1, v1); rightChainX->InsertPoint(1, x1);
    leftChain->InsertId(1, v2); leftChainX->InsertPoint(1, x2);
    }
  else
    {
    leftChain->InsertId(1, v1); leftChainX->InsertPoint(1, x1);
    rightChain->InsertId(1, v2); rightChainX->InsertPoint(1, x2);
    }
  
  // Walk along triangles (edge neighbors) towards point p2.
  while ( v1 != p2 )
    {
    this->Mesh->GetCellEdgeNeighbors(cellId, v1, v2, neis);
    if ( neis->GetNumberOfIds() != 1 )
      {//Mesh is folded or degenerate
      goto FAILURE;
      }
    cellId = neis->GetId(0);
    tris->InsertNextId(cellId);
    this->Mesh->GetCellPoints(cellId, npts, pts);
    for (j=0; j<3; j++)
      {
      if ( pts[j] != v1 && pts[j] != v2 )
        {//found point opposite current edge (v1,v2)
        if ( pts[j] == p2 )
          {
          v1 = p2; //this will cause the walk to stop
          rightChain->InsertNextId(p2); rightChainX->InsertNextPoint(p2X);
          leftChain->InsertNextId(p2); leftChainX->InsertNextPoint(p2X);
          }
        else
          {//keep walking
          this->GetPoint(pts[j], x1); x1[2] = 0.0;
          if ( vtkPlane::Evaluate(splitNormal, p1X, x1) > 0.0 )
            {
            v1 = pts[j];
            rightChain->InsertNextId(v1); rightChainX->InsertNextPoint(x1);
            }
          else
            {
            v2 = pts[j];
            leftChain->InsertNextId(v2); leftChainX->InsertNextPoint(x1);
            }
          }
        break;
        }//else found opposite point
      }//for all points in triangle
    }//while walking
  
  // Now that the to chains are formed, each chain forms a polygon (along with
  // the edge (p1,p2)) that requires triangulation. If we can successfully
  // triangulate the two polygons, we will delete the triangles contained within
  // the chains and replace them with the new triangulation.
  //
  success = 1;
  success &= (rightPoly->Triangulate(0, rightPtIds, rightTriPts));
  numRightTris = rightPtIds->GetNumberOfIds() / 3;

  success &= (leftPoly->Triangulate(0, leftPtIds, leftTriPts));
  numLeftTris = leftPtIds->GetNumberOfIds() / 3;
  
  if ( ! success ) 
    {//polygons on either side of edge are poorly shaped
    goto FAILURE;
    }
  
  // Okay, delete the old triangles and replace them with new ones. There should be
  // the same number of new triangles as old ones.
  leftTris = leftPtIds->GetPointer(0);
  for ( j=i=0; i<numLeftTris; i++, j++, leftTris+=3)
    {
    cellId = tris->GetId(j);
    this->Mesh->RemoveCellReference(cellId);
    for (k=0; k<3; k++)
      {//allocate new space for cell lists
      this->Mesh->ResizeCellList(leftTris[k],1);
      }
    this->Mesh->ReplaceLinkedCell(cellId, 3, leftTris);
    }

  rightTris = rightPtIds->GetPointer(0);
  for ( i=0; i<numRightTris; i++, j++, rightTris+=3)
    {
    cellId = tris->GetId(j);
    this->Mesh->RemoveCellReference(cellId);
    for (k=0; k<3; k++)
      {//allocate new space for cell lists
      this->Mesh->ResizeCellList(rightTris[k],1);
      }
    this->Mesh->ReplaceLinkedCell(cellId, 3, rightTris);
    }

  FAILURE:
    tris->Delete(); cells->Delete();
    leftPoly->Delete(); rightPoly->Delete(); neis->Delete();
    rightPtIds->Delete(); leftPtIds->Delete();
    rightTriPts->Delete(); leftTriPts->Delete();
    return success;
}

void vtkDelaunay2D::FillPolygons(vtkCellArray *polys, int *triUse)
{
  vtkIdType p1, p2, j, kk;
  int i, k;
  vtkIdType *pts = 0;
  vtkIdType *triPts;
  vtkIdType npts = 0;
  vtkIdType numPts;
  static double xyNormal[3]={0.0,0.0,1.0};
  double negDir[3], x21[3], x1[3], x2[3], x[3];
  vtkIdList *neis=vtkIdList::New();
  vtkIdType cellId, numNeis;
  vtkIdList *currentFront = vtkIdList::New(), *tmpFront;
  vtkIdList *nextFront = vtkIdList::New();
  vtkIdType numCellsInFront, neiId;
  vtkIdType numTriangles=this->Mesh->GetNumberOfCells();

  // Loop over edges of polygon, marking triangles on "outside" of polygon as outside.
  // Then perform a fill.
  for ( polys->InitTraversal(); polys->GetNextCell(npts,pts); )
    {
    currentFront->Reset();
    for (i=0; i<npts; i++)
      {
      p1 = pts[i];
      p2 = pts[(i+1)%npts];
      if ( ! this->Mesh->IsEdge(p1,p2) )
        {
        vtkWarningMacro(<<"Edge not recovered, polygon fill suspect");
        }
      else //Mark the "outside" triangles
        {
        neis->Reset();
        this->GetPoint(p1,x1);
        this->GetPoint(p2,x2);
        for (j=0; j<3; j++)
          {
          x21[j] = x2[j] - x1[j];
          }
        vtkMath::Cross (x21,xyNormal,negDir);
        this->Mesh->GetCellEdgeNeighbors(-1, p1, p2, neis); //get both triangles
        numNeis = neis->GetNumberOfIds();
        for (j=0; j<numNeis; j++)
          {//find the vertex not on the edge; evaluate it (and the cell) in/out
          cellId = neis->GetId(j);
          this->Mesh->GetCellPoints(cellId, numPts, triPts);
          for (k=0; k<3; k++)
            {
            if ( triPts[k] != p1 && triPts[k] != p2 ) 
              {
              break;
              }
            }
          this->GetPoint(triPts[k],x); x[2] = 0.0;
          if ( vtkPlane::Evaluate(negDir, x1, x) > 0.0 )
            {
            triUse[cellId] = 0;
            currentFront->InsertNextId(cellId);
            }
          else
            {
            triUse[cellId] = -1;
            }
          }
        }//edge was recovered
      }//for all edges in polygon

    // Okay, now perform a fill operation (filling "outside" values).
    //
    while ( (numCellsInFront = currentFront->GetNumberOfIds()) > 0 )
      {
      for (j=0; j < numCellsInFront; j++)
        {
        cellId = currentFront->GetId(j);

        this->Mesh->GetCellPoints(cellId, numPts, triPts);
        for (k=0; k<3; k++)
          {
          p1 = triPts[k];
          p2 = triPts[(k+1)%3];

          this->Mesh->GetCellEdgeNeighbors(cellId, p1, p2, neis);
          numNeis = neis->GetNumberOfIds();
          for (kk=0; kk<numNeis; kk++)
            {
            neiId = neis->GetId(kk);
            if ( triUse[neiId] == 1 ) //0 is what we're filling with
              {
              triUse[neiId] = 0;
              nextFront->InsertNextId(neiId);
              }
            }//mark all neigbors
          }//for all edges of cell
        } //all cells in front

      tmpFront = currentFront;
      currentFront = nextFront;
      nextFront = tmpFront;
      nextFront->Reset();
      } //while still advancing

    }//for all polygons

  //convert all unvisited to inside
  for (i=0; i<numTriangles; i++) 
    {
    if ( triUse[i] == -1 )
      {
      triUse[i] = 1;
      }
    }

  currentFront->Delete();
  nextFront->Delete();
  neis->Delete();
}

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

//----------------------------------------------------------------------------
vtkAbstractTransform * vtkDelaunay2D::ComputeBestFittingPlane(
  vtkPointSet *input)
{
  vtkIdType numPts=input->GetNumberOfPoints();
  double m[9], v[3], x[3];
  vtkIdType ptId;
  int i;
  double *c1, *c2, *c3, det;
  double normal[3];
  double origin[3];

  const double tolerance = 1.0e-03;
    
  //  This code was taken from the vtkTextureMapToPlane class
  //  and slightly modified.
  //  
  for (i=0; i<3; i++)
    {
    normal[i] = 0.0;
    }

  //  Compute least squares approximation.  
  //  Compute 3x3 least squares matrix
  v[0] = v[1] = v[2] = 0.0;
  for (i=0; i<9; i++)
    {
    m[i] = 0.0;
    }

  for (ptId=0; ptId < numPts; ptId++) 
    {
    input->GetPoint(ptId, x);

    v[0] += x[0]*x[2];
    v[1] += x[1]*x[2];
    v[2] += x[2];

    m[0] += x[0]*x[0];
    m[1] += x[0]*x[1];
    m[2] += x[0];

    m[3] += x[0]*x[1];
    m[4] += x[1]*x[1];
    m[5] += x[1];

    m[6] += x[0];
    m[7] += x[1];
    }
  m[8] = numPts;

  origin[0] = m[2] / numPts;
  origin[1] = m[5] / numPts;
  origin[2] = v[2] / numPts;

  //  Solve linear system using Kramers rule
  //
  c1 = m; c2 = m+3; c3 = m+6;
  if ( (det = vtkMath::Determinant3x3 (c1,c2,c3)) > tolerance )
    {
    normal[0] =  vtkMath::Determinant3x3 (v,c2,c3) / det;
    normal[1] =  vtkMath::Determinant3x3 (c1,v,c3) / det;
    normal[2] = -1.0; // because of the formulation
    }

  vtkTransform * transform = vtkTransform::New();

  // Set the new Z axis as the normal to the best fitting
  // plane.
  double zaxis[3];
  zaxis[0] = 0;
  zaxis[1] = 0;
  zaxis[2] = 1;

  double rotationAxis[3];
  
  vtkMath::Normalize(normal);
  vtkMath::Cross(normal,zaxis,rotationAxis);
  vtkMath::Normalize(rotationAxis);

  const double rotationAngle = 
    180.0*acos(vtkMath::Dot(zaxis,normal))/3.1415926; 

  transform->PreMultiply();
  transform->Identity();

  transform->RotateWXYZ(rotationAngle,
    rotationAxis[0], rotationAxis[1], rotationAxis[2]);

  // Set the center of mass as the origin of coordinates
  transform->Translate( -origin[0], -origin[1], -origin[2] );

  return transform;
}

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

  os << indent << "Alpha: " << this->Alpha << "\n";
  os << indent << "ProjectionPlaneMode: " 
     << ((this->ProjectionPlaneMode == VTK_BEST_FITTING_PLANE)? "Best Fitting Plane" : "XY Plane") << "\n";
  os << indent << "Transform: " << (this->Transform ? "specified" : "none") << "\n";
  os << indent << "Tolerance: " << this->Tolerance << "\n";
  os << indent << "Offset: " << this->Offset << "\n";
  os << indent << "Bounding Triangulation: "
     << (this->BoundingTriangulation ? "On\n" : "Off\n");
}