struct triangulateio triMesh;
struct triangulateio in;

Tile::Tile()
{
  InitMesh(triMesh);
  InitMesh(in);
}

VOID Tile::CalculateIndexedPoints()
{
  UINT_32 i, uiPoints, uiOuterPoints, uiInnerPoints;
  CListIterInt iter;
  TilePoint *tpt;

  if( m_ppIndexedTilePoints != NULL )
    {
      delete m_ppIndexedTilePoints;
    }

  uiOuterPoints = m_listSphereOutlinePoints.GetLength();
  uiInnerPoints = m_listSphereInnerPoints.GetLength();
  uiPoints = uiOuterPoints + uiInnerPoints;

  m_uiPoints = uiPoints;
  m_ppIndexedTilePoints = new TilePoint*[uiPoints];

  iter.Init(&m_listSphereOutlinePoints);
  for( tpt = (TilePoint *)(iter.PeekFirst()), i=0; 
       tpt != NULL;
       tpt = (TilePoint *)(iter.PeekNext()), i++ )
    {
      ASSERT( i < uiPoints );
      m_ppIndexedTilePoints[i] = tpt;
    }
  ASSERT( i == uiOuterPoints );

  iter.Init(&m_listSphereInnerPoints);
  for( tpt = (TilePoint *)(iter.PeekFirst()); 
       tpt != NULL;
       tpt = (TilePoint *)(iter.PeekNext()), i++ )
    {
      ASSERT( i < uiPoints );
      m_ppIndexedTilePoints[i] = tpt;
    }
  ASSERT( i == uiPoints );
}

TilePoint* Tile::GetIndexedPoint(UINT_32 index)
{
  //ASSERT( index < m_uiPoints );

  if( index  >= m_uiPoints )
    {
      printf( "Triangulation failed because of bad outer contour.\n");
      printf( "Please uncross all edges.\n");
      return NULL;
    }

  return (m_ppIndexedTilePoints[index]);
}

UINT_32 Tile::GetPointIndex(TilePoint *tpt)
{
  UINT_32 i;

  for( i=0; i < m_uiPoints; i++ )
    {
      if( tpt == m_ppIndexedTilePoints[i] )
	{
	  return i;
	}
    }

  ASSERT(FALSE);
  return 0;
}

VOID Tile::CalculateHoleEdges()
{
  UINT_32 i, index;
  TilePoint *tpt1;

  m_uiHoleEdges = 0;

  if( triMesh.edgemarkerlist != NULL )
    {
      for( i = 0; i < (UINT_32)triMesh.numberofedges; i++ )
	{
	  if( 1 == triMesh.edgemarkerlist[i] )             // Boundary
	    {
	      index = triMesh.edgelist[2*i];
	      tpt1 = GetIndexedPoint(index);
		  
	      if( !SphereOutlinePoint(tpt1) ) // Hole
		{
		  m_uiHoleEdges++;
		}	  
	    }
	}
    }
  //printf( "HoleEdges: %lu\n", m_uiHoleEdges );
}

VOID Tile::Foo()
{
  // This produces the correct results

  //struct triangulateio in, vorout, out;
  struct triangulateio vorout;

  // Flatten points to 2D and calculate indices
  CalculateStereographicPoints();
  CalculateIndexedPoints();

  if( triMesh.pointlist != 0 )
    {
      safefree(triMesh.pointlist);
      safefree(triMesh.pointattributelist);
      safefree(triMesh.pointmarkerlist);
      safefree(triMesh.trianglelist);
      safefree(triMesh.triangleattributelist);
      safefree(triMesh.trianglearealist);
      safefree(triMesh.neighborlist);
      safefree(triMesh.segmentlist);
      safefree(triMesh.segmentmarkerlist);
      safefree(triMesh.edgelist);
      safefree(triMesh.edgemarkerlist);
    }

  in.pointlist = NULL;
  in.pointattributelist = NULL;
  in.pointmarkerlist = NULL;
  in.trianglelist = NULL;
  in.triangleattributelist = NULL;
  in.trianglearealist = NULL;
  in.neighborlist = NULL;
  in.segmentlist = NULL;
  in.segmentmarkerlist = NULL;
  in.holelist = NULL;
  in.regionlist = NULL;
  in.edgelist = NULL;
  in.edgemarkerlist = NULL;
  in.normlist = NULL;

  /*
  out.pointlist = NULL;
  out.pointattributelist = NULL;
  out.pointmarkerlist = NULL;
  out.trianglelist = NULL;
  out.triangleattributelist = NULL;
  out.trianglearealist = NULL;
  out.neighborlist = NULL;
  out.segmentlist = NULL;
  out.segmentmarkerlist = NULL;
  out.holelist = NULL;
  out.regionlist = NULL;
  out.edgelist = NULL;
  out.edgemarkerlist = NULL;
  out.normlist = NULL;
  */

  ////////////////////////////////////////////////////////////////////////////
  in.numberofpoints = 9;
  in.numberofpointattributes = 1;
  in.pointlist = (REAL *) malloc(in.numberofpoints * 2 * sizeof(REAL));
  in.pointlist[0] = 0;
  in.pointlist[1] = -2.98858;
  in.pointlist[2] = 4.24264;
  in.pointlist[3] = -0.816496;
  in.pointlist[4] = 0;
  in.pointlist[5] = 0.800788; 
  in.pointlist[6] = -4.24264;
  in.pointlist[7] = -0.816497;
  in.pointlist[8] = -0.202093;
  in.pointlist[9] = -0.578931;
  in.pointlist[10] = -1.33512;
  in.pointlist[11] = -1.29447;
  in.pointlist[12] = -0.0570442;
  in.pointlist[13] = -2.16679;
  in.pointlist[14] = 1.68719;
  in.pointlist[15] = -1.19395;
  in.pointlist[16] = -0.202093;
  in.pointlist[17] = -0.578931;

  // Must have attributelist and pointmarkerlist or segfaults
  in.pointattributelist = (REAL *) malloc(in.numberofpoints * in.numberofpointattributes *
                                          sizeof(REAL));
  in.pointattributelist[0] = 0;
  in.pointattributelist[1] = 0;
  in.pointattributelist[2] = 0;
  in.pointattributelist[3] = 0;
  in.pointattributelist[4] = 0;
  in.pointattributelist[5] = 0;
  in.pointattributelist[6] = 0;
  in.pointattributelist[7] = 0;
  in.pointattributelist[8] = 0;

  in.pointmarkerlist = (int *) malloc(in.numberofpoints * sizeof(int));
  in.pointmarkerlist[0] = 0;
  in.pointmarkerlist[1] = 0;
  in.pointmarkerlist[2] = 0;
  in.pointmarkerlist[3] = 0;
  in.pointmarkerlist[4] = 0;
  in.pointmarkerlist[5] = 0;
  in.pointmarkerlist[6] = 0;
  in.pointmarkerlist[7] = 0;
  in.pointmarkerlist[8] = 0;

  in.numberofsegments = 8;
  in.segmentlist = (int *) malloc(in.numberofsegments * 2 * sizeof(int));
  in.segmentlist[0] = 0;
  in.segmentlist[1] = 1;
  in.segmentlist[2] = 1;
  in.segmentlist[3] = 2;
  in.segmentlist[4] = 2;
  in.segmentlist[5] = 3;
  in.segmentlist[6] = 3;
  in.segmentlist[7] = 0;
  in.segmentlist[8] = 4;
  in.segmentlist[9] = 5;
  in.segmentlist[10] = 5;
  in.segmentlist[11] = 6;
  in.segmentlist[12] = 6;
  in.segmentlist[13] = 7;
  in.segmentlist[14] = 7;
  in.segmentlist[15] = 8;

  in.numberofholes = 1;
  in.holelist = (REAL *) malloc(in.numberofholes * 2 * sizeof(REAL));
  in.holelist[0] = 0.010;
  in.holelist[1] = -1.5;

  in.numberofregions = 0;

  // Don't need region
  /*
  in.numberofregions = 1;
  in.regionlist = (REAL *) malloc(in.numberofregions * 4 * sizeof(REAL));
  in.regionlist[0] = 0.0;
  in.regionlist[1] = 0.0;
  in.regionlist[2] = 0.0;
  in.regionlist[3] = 0.0;
*/

  printf("Input point set:\n\n");
  report(&in, 1, 0, 0, 1, 1, 0, 1);

  /* Triangulate the points.  Switches are chosen to read and write a  */
  /*   PSLG (p), preserve the convex hull (c), number everything from  */
  /*   zero (z), assign a regional attribute to each element (A), and  */
  /*   produce an edge list (e), a Voronoi diagram (v), and a triangle */
  /*   neighbor list (n).                                              */

  /* triangulate("pczAevn", &in, &out, &vorout); */
  triangulate("pze", &in, &triMesh, &vorout);

  printf("Initial triangulation:\n\n");
  report(&triMesh, 1, 1, 0, 1, 1, 0, 1);

  safefree(in.pointlist);
  safefree(in.pointattributelist);
  safefree(in.pointmarkerlist);
  safefree(in.regionlist);
}

VOID Tile::ClearMesh(struct triangulateio &mesh)
{
  safefree(mesh.pointlist);
  safefree(mesh.pointattributelist);
  safefree(mesh.pointmarkerlist);
  safefree(mesh.trianglelist);
  safefree(mesh.triangleattributelist);
  safefree(mesh.trianglearealist);
  safefree(mesh.neighborlist);
  safefree(mesh.segmentlist);
  safefree(mesh.segmentmarkerlist);
  safefree(mesh.holelist);
  safefree(mesh.regionlist);
  safefree(mesh.edgelist);
  safefree(mesh.edgemarkerlist);
  safefree(mesh.normlist);

  InitMesh(mesh);
}

VOID Tile::InitMesh(struct triangulateio &mesh)
{
  mesh.pointlist = NULL;
  mesh.pointattributelist = NULL;
  mesh.pointmarkerlist = NULL;
  mesh.trianglelist = NULL;
  mesh.triangleattributelist = NULL;
  mesh.trianglearealist = NULL;
  mesh.neighborlist = NULL;
  mesh.segmentlist = NULL;
  mesh.segmentmarkerlist = NULL;
  mesh.holelist = NULL;
  mesh.regionlist = NULL;
  mesh.edgelist = NULL;
  mesh.edgemarkerlist = NULL;
  mesh.normlist = NULL;

  mesh.numberofpoints = 0;
  mesh.numberofpointattributes = 0;
  mesh.numberoftriangles = 0;
  mesh.numberofcorners = 0;
  mesh.numberoftriangleattributes = 0;
  mesh.numberofsegments = 0;
  mesh.numberofedges = 0;
}

VOID Tile::Triangulate()
{
  CListIterInt iter;
  TilePoint *tpt;
  UINT_32 i, j, uiPoints, uiOuterPoints, uiInnerPoints, uiSegments, uiHolePoints;
  UINT_32 uiPointAttributes;
  Point pt;

  // Flatten points to 2D and calculate indices
  CalculateStereographicPoints();
  CalculateIndexedPoints();

  ////////////////////////////////////////////////////////////
  // Free all allocated arrays, including those allocated by Triangle.
  // Initialize in and out

  //if( triMesh.numberofpoints != 0 )
  ClearMesh(triMesh);
  ClearMesh(in);

  ////////////////////////////////////////////////////////////
  // Define points for triangulation
  uiOuterPoints = m_listSphereOutlinePoints.GetLength();
  uiInnerPoints = m_listSphereInnerPoints.GetLength();
  uiHolePoints = m_listSphereHolePoints.GetLength();
  uiPoints = uiOuterPoints + uiInnerPoints;

  in.numberofpoints = uiPoints;
  in.pointlist = (REAL *) malloc(in.numberofpoints * 2 * sizeof(REAL));

  ////////////////////////////////////////////////////////////
  // points (markers and attributes)
  iter.Init(&m_listSphereOutlinePoints);
  for( tpt = (TilePoint *)(iter.PeekFirst()), i=0; 
       tpt != NULL;
       tpt = (TilePoint *)(iter.PeekNext()), i++ )
    {
      pt = tpt->GetStereographicPoint();
      pt = ScaleStereographicPoint(pt);

      in.pointlist[(2*i)]   = pt[X];
      in.pointlist[(2*i)+1] = pt[Y];
    }
  ASSERT( i == uiOuterPoints );

  iter.Init(&m_listSphereInnerPoints);
  for( tpt = (TilePoint *)(iter.PeekFirst()); 
       tpt != NULL;
       tpt = (TilePoint *)(iter.PeekNext()), i++ )
    {
      pt = tpt->GetStereographicPoint();
      pt = ScaleStereographicPoint(pt);

      in.pointlist[(2*i)]   = pt[X];
      in.pointlist[(2*i)+1] = pt[Y];
    }
  ASSERT( i == uiPoints );

  ////////////////////////////////////////////////////////////
  // attrributelist and pointmarkerlist
  in.numberofpointattributes = 1; 
  uiPointAttributes = in.numberofpoints * in.numberofpointattributes;
  in.pointattributelist = (REAL *) malloc(in.numberofpoints * in.numberofpointattributes *
                                          sizeof(REAL));
  
  for( i = 0; i < uiPointAttributes; i++ )
    {
      in.pointattributelist[i] = 0;
    }

  in.pointmarkerlist = (int *) malloc(in.numberofpoints * sizeof(int));

  for( i = 0; i < (UINT_32)in.numberofpoints; i++ )
    {
      in.pointmarkerlist[i] = 0;
    }

  ////////////////////////////////////////////////////////////
  // segments
  uiSegments = uiOuterPoints;
  iter.Init(&m_listSphereInnerPoints);
  for( tpt = (TilePoint *)(iter.PeekFirst()); 
       tpt != NULL;
       tpt = (TilePoint *)(iter.PeekNext()) )
    {
      if( tpt->GetConnected() )
	{
	  uiSegments++;
	}
    }
  in.numberofsegments = uiSegments;
  in.segmentlist = (int *) malloc(in.numberofsegments * 2 * sizeof(int));

  // make simple polygon to triangulate (make sure closed or problems)
  iter.Init(&m_listSphereOutlinePoints);
  for( tpt = (TilePoint *)(iter.PeekFirst()), i=0; 
       tpt != NULL;
       tpt = (TilePoint *)(iter.PeekNext()), i++ )
    {
      in.segmentlist[(2*i)] = i;

      if( i+1 > uiOuterPoints-1 )
	in.segmentlist[(2*i)+1] = 0;
      else
	in.segmentlist[(2*i)+1] = i+1;
    }
  ASSERT( uiOuterPoints == i );

  iter.Init(&m_listSphereInnerPoints);
  for( tpt = (TilePoint *)(iter.PeekFirst()), j=0; 
       tpt != NULL;
       tpt = (TilePoint *)(iter.PeekNext()), j++ )
    {
      if( tpt->GetConnected() )
	{
	  in.segmentlist[(2*i)]   = uiOuterPoints+j;
	  in.segmentlist[(2*i)+1] = uiOuterPoints+j+1;

	  i++;
	}
    }
  ASSERT( uiSegments == i );

  in.numberofholes = uiHolePoints;
  if( uiHolePoints > 0 )
    {
      in.holelist = (REAL *) malloc(in.numberofholes * 2 * sizeof(REAL));

      iter.Init(&m_listSphereHolePoints);
      for( tpt = (TilePoint *)(iter.PeekFirst()), i=0; 
	   tpt != NULL;
	   tpt = (TilePoint *)(iter.PeekNext()), i++ )
	{
	  pt = tpt->GetStereographicPoint();
	  pt = ScaleStereographicPoint(pt);

	  in.holelist[(2*i)+0] = pt[X]; 
	  in.holelist[(2*i)+1] = pt[Y]; 
	}
      ASSERT( uiHolePoints == i);
    }

  in.numberofregions = 0;  


  printf("Input point set:\n\n");
  report(&in, 1, 0, 0, 1, 1, 0, 1);

  ////////////////////////////////////////////////////////////
  // Triangulate points. Switch -pze (c)
  triangulate("pze", &in, &triMesh, (struct triangulateio *) NULL);

  printf("Triangulation:\n\n");
  report(&triMesh, 1, 1, 0, 1, 1, 0, 1);

  printf("------------------------------------------------------\n");

  ////////////////////////////////////////////////////////////
  // Free all allocated arrays, including those allocated by Triangle.

  //free(in.pointlist);
  //free(in.segmentlist);
  //ClearMesh(in);

  printf("------------------------------------------------------\n\n");


  ////////////////////////////////////////////////
  CalculateHoleEdges();
}
