coloring.cpp

Go to the documentation of this file.

#include <iostream>
#include <iomanip>
#include <limits>
#include <list>
#include <set>
#include <stack>
#include <assert.h>
#include "geometry.hpp"
#include "properties.hpp"
#include "coloring.hpp"
#include "query_point.hpp"
#include <CGAL/IO/Qt_widget.h>
#include <CGAL/IO/Qt_widget_Polygon_2.h>

using namespace Arak;
using namespace Arak::Geometry;

#ifdef SANITY_CHECK
// Always test for validity.
#define TEST_VALIDITY test();
#else
// Test randomly for validity.
#define TEST_VALIDITY if (Arak::Util::default_random.uniform(10000) == 0) test(); 
#endif

double Coloring::Vertex::logSine() {
  Geometry::Point a, b = this->point(), c;
  if (type() == INTERIOR) {
    a = getPrevVertex()->point();
    c = getNextVertex()->point();
  } else if (type() == BOUNDARY) {
    if (hasPrevIntEdge())
      a = getPrevVertex()->point();
    else
      a = getNextVertex()->point();
    c = getNextBdEdge()->getNextVertex()->point();
    if (CGAL::angle(a, b, c) != CGAL::ACUTE) 
      c = getPrevBdEdge()->getPrevVertex()->point();
  } else assert(false);
  return logSineB(a, b, c);
}

bool Coloring::InteriorEdge::crosses(const Geometry::Point& a, 
             const Geometry::Point& b) const {
  return Arak::Geometry::crosses(*this, a, b);
}

Geometry::Point Coloring::Edge::randomPoint(Arak::Util::Random& random) const {
  return random.uniform(this->segment());
}

Geometry::Point 
Coloring::BoundaryEdge::randomPoint(Arak::Util::Random& random) const {
  Geometry::Point p = Edge::randomPoint(random);
#ifndef EXACT_GEOM_CXN
  // If constructions are not exact, then we must explicitly enforce
  // that points sampled from boundary edges remain on the boundary.
  // The technique below works for axis-aligned rectangle windows.
  // For more general windows, we might represent a boundary vertex by
  // its face's corners and a real number between 0 and 1.
  const Coloring::VertexHandle c1 = getPrevCornerVertex();
  const Coloring::VertexHandle c2 = getNextCornerVertex();
  if (c1->x() == c2->x())
    p = Geometry::Point(c1->x(), p.y());
  else if (c1->y() == c2->y())
    p = Geometry::Point(p.x(), c1->y());
  else
    assert(false);
#endif
  return p;
}

void Coloring::test() const {
  // Check the integrity of the interior vertices.
  for (int i = 0; i < numVertices(Vertex::INTERIOR); i++) {
    VertexHandle v = getVertex(Vertex::INTERIOR, i);
    assert(v->index == i);
    assert(v->type() == Vertex::INTERIOR);
    assert(v == v->getNextIntEdge()->getPrevVertex());
    assert(v == v->getPrevIntEdge()->getNextVertex());
#ifdef EXACT_GEOM_CXN
    assert(interior(*v));
#endif
  }
  // Check the integrity of the boundary vertices.
  for (int i = 0; i < numVertices(Vertex::BOUNDARY); i++) {
    VertexHandle v = getVertex(Vertex::BOUNDARY, i);
    assert(v->index == i);
    assert(v->type() == Vertex::BOUNDARY);
    assert(v->hasPrevIntEdge() || v->hasNextIntEdge());
    assert(!v->hasNextIntEdge() || 
     (v == v->getNextIntEdge()->getPrevVertex()));
    assert(!v->hasPrevIntEdge() || 
     (v == v->getPrevIntEdge()->getNextVertex()));
    assert(v == v->getNextBdEdge()->getPrevVertex());
    assert(v == v->getPrevBdEdge()->getNextVertex());
#ifdef EXACT_GEOM_CXN
    assert(boundary(*v));
#endif
  }
  // Check the integrity of the corner vertices.
  for (int i = 0; i < numVertices(Vertex::CORNER); i++) {
    VertexHandle v = getVertex(Vertex::CORNER, i);
    assert(v->index == i);
    assert(v->type() == Vertex::CORNER);
    assert(v == v->getNextBdEdge()->getPrevVertex());
    assert(v == v->getPrevBdEdge()->getNextVertex());
    assert(v->Geometry::Point::operator==(bd[i]));
  }
  // Check the integrity of the interior edges.
  for (int i = 0; i < numInteriorEdges(); i++) {
    IntEdgeHandle e = getInteriorEdge(i);
    assert(e->index == i);
    assert(!(e->getNextVertex() == e->getPrevVertex()));
    assert(e == e->getNextVertex()->getPrevIntEdge());
    assert(e == e->getPrevVertex()->getNextIntEdge());
    assert((e->getNextVertex()->type() == Vertex::INTERIOR) ||
     (e->getNextVertex()->type() == Vertex::BOUNDARY));
    assert((e->getPrevVertex()->type() == Vertex::INTERIOR) ||
     (e->getPrevVertex()->type() == Vertex::BOUNDARY));
    assert(Segment(*(e->getPrevVertex()), *(e->getNextVertex())) == 
     ((Segment)*e));
    assert(e->getNextVertex()->Geometry::Point::operator!=(*(e->getPrevVertex())));
    assert(CGAL::to_double(CGAL::squared_distance(*(e->getPrevVertex()),
              *(e->getNextVertex()))) >=
     minEdgeLengthSq);
#ifdef TRIANGULATION_EDGE_INDEX
    // Check that the edge is indexed in the triangulation.
    TFaceHandle fh;
    int index;
    assert(tri.is_edge(e->prevVertex->tvh, e->nextVertex->tvh, fh, index));
    assert(fh->is_constrained(index));
    assert(tri.is_valid(false, 0));
#endif
  }
  // Check exhaustively for crossing edges.
  for (int i = 0; i < numInteriorEdges(); i++) {
    IntEdgeHandle e = getInteriorEdge(i);
    for (int j = 0; j < numInteriorEdges(); j++) {
      if (i == j) continue;
      IntEdgeHandle f = getInteriorEdge(j);
      assert(!e->crosses(f->source(), f->target()));
    }
  }
  // Check the integrity of the boundary edges.
  for (int i = 0; i < numBoundaryEdges(); i++) {
    BdEdgeHandle e = getBoundaryEdge(i);
    assert(e->index == i);
    assert(!(e->getNextVertex() == e->getPrevVertex()));
    assert(e == e->getNextVertex()->getPrevBdEdge());
    assert(e == e->getPrevVertex()->getNextBdEdge());
    assert((e->getNextVertex()->type() == Vertex::BOUNDARY) ||
     (e->getNextVertex()->type() == Vertex::CORNER));
    assert((e->getPrevVertex()->type() == Vertex::BOUNDARY) ||
     (e->getPrevVertex()->type() == Vertex::CORNER));
    assert(Segment(*(e->getPrevVertex()), *(e->getNextVertex())) == 
     ((Segment)*e));
    assert(e->getNextVertex()->point() != e->getPrevVertex()->point());
    // Check the cached corner pointers.
    if (e->getPrevVertex()->type() == Vertex::CORNER) 
      assert(e->getPrevCornerVertex() == e->getPrevVertex());
    else 
      assert(e->getPrevCornerVertex() == 
       e->getPrevVertex()->getPrevCornerVertex());
    if (e->getNextVertex()->type() == Vertex::CORNER) 
      assert(e->getNextCornerVertex() == e->getNextVertex());
    else 
      assert(e->getNextCornerVertex() == 
       e->getNextVertex()->getNextCornerVertex());
  }
  // Check the color of the query points.
  for (int i = 0; i < queryPoints->size(); i++)
    assert(color(queryPoints->point(i)) == 
     queryPoints->point(i).color());
}

Coloring::IntEdgeHandle Coloring::newInteriorEdge(VertexHandle prev, 
              VertexHandle next) {
  // Check that the new edge is not degerate.
  assert(prev->point() != next->point());
  // Allocate the new edge.
  IntEdgeHandle e = intEdgeStore.newItem();
  e->index = interiorEdges.size();
  interiorEdges.push_back(e);
  e->nextVertex = e->prevVertex = VertexHandle();
  // Assign the edge's new coordinates.
  e->Segment::operator=(Segment(prev->point(), next->point()));
  // Update the vertex and edge pointers.
  e->prevVertex = prev;
  e->nextVertex = next;
  // Ensure the vertices are not yet attached.
  assert(!(e->nextVertex->prevIntEdge.valid()));
  assert(!(e->prevVertex->nextIntEdge.valid()));
  e->nextVertex->prevIntEdge = e->prevVertex->nextIntEdge = e;
#ifdef GRID_EDGE_INDEX
  // Index the edge.
  typedef InteriorEdgeIndex::LineCellIterator Iterator; 
  Iterator it(*intEdgeIndex, *prev, *next, Iterator::SEGMENT()), end;
  while (it != end) {
    InteriorEdgeIndex::Cell& cell = *it;
    InteriorEdgeIndex::Cell::Entry entry = cell.add(e); 
    if (e->unused.empty())
      e->entries.push_front(entry);
    else {
      InteriorEdge::CellEntryList::iterator jt = e->unused.begin();
      *jt = entry;
      e->entries.splice(e->entries.begin(), e->unused, jt);
    }
    ++it;
  }
#endif
#ifdef TRIANGULATION_EDGE_INDEX
  // Add the constraint to the triangulation.
  tri.insert_constraint(prev->tvh, next->tvh);
#endif
  // Return the new edge.
  return *e;
}

Coloring::BdEdgeHandle Coloring::newBoundaryEdge(VertexHandle prev, 
             VertexHandle next,
             VertexHandle prevCorner,
             VertexHandle nextCorner) {
  assert(prev->type() != Vertex::INTERIOR);
  assert(next->type() != Vertex::INTERIOR);
  // Check that the edge is not degenerate.
  assert(prev->point() != next->point());
  // Allocate the new edge.
  BdEdgeHandle e = bdEdgeStore.newItem();
  e->index = boundaryEdges.size();
  boundaryEdges.push_back(e);
  e->nextVertex = e->prevVertex = VertexHandle();
  // Assign the edge's new coordinates.
  e->Segment::operator=(Segment(*prev, *next));
  // Update the vertex and edge pointers.
  e->prevVertex = prev;
  e->nextVertex = next;
  // Ensure the vertices are not already attached.
  assert(!(e->nextVertex->prevBdEdge.valid()));
  assert(!(e->prevVertex->nextBdEdge.valid()));
  e->nextVertex->prevBdEdge = e->prevVertex->nextBdEdge = e;
  e->prevCorner = prevCorner;
  e->nextCorner = nextCorner;
  // Return the new edge.
  return *e;
}

void Coloring::initialize(Geometry::Rectangle &boundary, 
        int gridRows, int gridCols,
        double minEdgeLength) {
  clear();
  this->minEdgeLengthSq = minEdgeLength * minEdgeLength;
  bd = boundary;
  intEdgeIndex = new InteriorEdgeIndex(boundary, gridRows, gridCols);
  queryPoints = NULL;
  // Create the original boundary edges.
  VertexHandle firstVertex = newVertex(Vertex::CORNER, bd[0]);
  VertexHandle prevVertex = firstVertex;
  for (int i = 1; i < 4; i++) {
    VertexHandle nextVertex = newVertex(Vertex::CORNER, bd[i]);
    newBoundaryEdge(prevVertex, nextVertex, prevVertex, nextVertex);
    prevVertex = nextVertex;
  }
  newBoundaryEdge(*prevVertex, firstVertex, *prevVertex, firstVertex);
  assert(numVertices(Vertex::CORNER) == 4);
  assert(numBoundaryEdges() == 4);
  // Allocate a single query point in the center of the window.  (We
  // must have at least one query point to determine which of the two
  // colorings is consistent with the current set of edges.)
  queryPoints = new QueryPointIndex(boundary, 1, 1);
  // Initialize all query points to be white.
  for (int i = 0; i < queryPoints->size(); i++) {
    assert(interior(queryPoints->point(i)));
    queryPoints->point(i).setColor(WHITE);
  }
  TEST_VALIDITY;
#ifdef SANITY_CHECK
  std::cerr << "WARNING: compiled with SANITY_CHECK on; expect poor performance" << std::endl;
#endif
}

Coloring::Coloring() 
  : intEdgeIndex(NULL), queryPoints(NULL) {
  Geometry::Rectangle r(Geometry::Kernel::RT(0), Geometry::Kernel::RT(0), 
      Geometry::Kernel::RT(1), Geometry::Kernel::RT(1));
  initialize(r, 10, 10);
}

Coloring::Coloring(Geometry::Rectangle &boundary, 
       int gridRows, int gridCols, double minEdgeLength) 
  : intEdgeIndex(NULL), queryPoints(NULL) {
  initialize(boundary, gridRows, gridCols, minEdgeLength);
}

Coloring::Coloring(const Arak::Util::PropertyMap& props) 
  : intEdgeIndex(NULL), queryPoints(NULL) {
  using namespace Arak::Util;
  double x, y, width, height;
  int rows, cols;
  assert(parse(getp(props, "arak.coloring.xmin"), x));
  assert(parse(getp(props, "arak.coloring.ymin"), y));
  assert(parse(getp(props, "arak.coloring.width"), width));
  assert(parse(getp(props, "arak.coloring.height"), height));
  assert(parse(getp(props, "arak.coloring.rows"), rows));
  assert(parse(getp(props, "arak.coloring.cols"), cols));
  Geometry::Rectangle boundary(x, y, x + width, y + height);
  double mel = 0.0;
  if (hasp(props, "arak.coloring.min_edge_length"))
    assert(parse(getp(props, "arak.coloring.min_edge_length"), mel));
  initialize(boundary, rows, cols, mel);
}

void Coloring::clear() {
  // Delete the grid.
  if (intEdgeIndex != NULL) {
    delete intEdgeIndex;
    intEdgeIndex = NULL;
  }
  // Delete the query point index.
  if (queryPoints != NULL) {
    delete queryPoints;
    queryPoints = NULL;
  }
  // Delete all vertex objects in use.
  for (unsigned int i = 0; i < interiorVertices.size(); i++)
    delete &*(interiorVertices[i]);
  interiorVertices.clear();
  for (unsigned int i = 0; i < boundaryVertices.size(); i++)
    delete &*(boundaryVertices[i]);
  boundaryVertices.clear();
  for (unsigned int i = 0; i < cornerVertices.size(); i++)
    delete &*(cornerVertices[i]);
  cornerVertices.clear();
  // Delete all edge objects in use.
  for (unsigned int i = 0; i < interiorEdges.size(); i++)
    delete &*(interiorEdges[i]);
  interiorEdges.clear();
  for (unsigned int i = 0; i < boundaryEdges.size(); i++)
    delete &*(boundaryEdges[i]);
  boundaryEdges.clear();
}

Coloring::~Coloring() {
  clear();
}

void Coloring::addQueryPoints(const QueryPointIndex& index) const {
  // First set the colors of the new points.
  for (int i = 0; i < index.size(); i++) {
    const QueryPoint q = index.point(i);
    q.setColor(color(q));
  }
  // Then make a new query point index containing all points.
  QueryPointIndex *newQueryPoints = 
    new QueryPointIndex(*queryPoints, index);
  // Replace the current index with the new one.
  delete queryPoints;
  queryPoints = newQueryPoints;
  TEST_VALIDITY;
}

Coloring::VertexHandle Coloring::newVertex(Vertex::Type type, 
             const Geometry::Point& p) {
  VertexHandle v = vertexStore.newItem(); 
  // Assign the type, position, and insert the vertex in the
  // triangulation.
  v->t = type;
  v->Geometry::Point::operator=(p);
#ifdef TRIANGULATION_EDGE_INDEX
  v->tvh = tri.insert(p);
  v->tvh->vertex = v;
#endif
  switch (type) {
  case Vertex::INTERIOR:
    v->index = interiorVertices.size();
    interiorVertices.push_back(v);
    break;
  case Vertex::BOUNDARY:
    v->index = boundaryVertices.size();
    boundaryVertices.push_back(v);
    break;
  case Vertex::CORNER:
    v->index = cornerVertices.size();
    cornerVertices.push_back(v);
    break;
  default:
    assert(false);
  }
  // Initialize the other fields of the Vertex object.
  v->nextIntEdge = v->prevIntEdge = IntEdgeHandle();
  v->nextBdEdge = v->prevBdEdge = BdEdgeHandle();
  return *v;
}

void Coloring::freeVertex(Coloring::VertexHandle& v) {
  // Make sure the vertex is not attached to any edges.
  assert(!(v->nextIntEdge.valid()));
  assert(!(v->prevIntEdge.valid()));
  assert(!(v->nextBdEdge.valid()));
  assert(!(v->prevBdEdge.valid()));
#ifdef TRIANGULATION_EDGE_INDEX
  // Remove the vertex from the triangulation.
  tri.remove(v->tvh);
  v->tvh = TVertexHandle();
#endif
  // Remove the vertex from its random access vector.
  int n;
  switch (v->type()) {
  case Vertex::INTERIOR:
    n = interiorVertices.size();
    if (n > 1) {
      VertexHandle w = interiorVertices[n - 1];
      interiorVertices[n - 1] = interiorVertices[v->index];
      interiorVertices[v->index] = w;
      w->index = v->index;
    }
    interiorVertices.pop_back();
    break;
  case Vertex::BOUNDARY:
    n = boundaryVertices.size();
    if (n > 1) {
      VertexHandle w = boundaryVertices[n - 1];
      boundaryVertices[n - 1] = boundaryVertices[v->index];
      boundaryVertices[v->index] = w;
      w->index = v->index;
    }
    boundaryVertices.pop_back();
    break;
  case Vertex::CORNER:
    n = cornerVertices.size();
    if (n > 1) {
      VertexHandle w = cornerVertices[n - 1];
      cornerVertices[n - 1] = cornerVertices[v->index];
      cornerVertices[v->index] = w;
      w->index = v->index;
    }
    cornerVertices.pop_back();
    break;
  default:
    assert(false);
  }
  v->t = Vertex::INVALID;
  v->index = -1;
  vertexStore.deleteItem(&*v);
  // Invalidate the handle.
  v = VertexHandle();
}

void Coloring::freeEdge(Coloring::IntEdgeHandle& e) {
#ifdef TRIANGULATION_EDGE_INDEX
  // Remove the constraint from the triangulation.
  TFaceHandle fh;
  int i;
  assert(tri.is_edge(e->prevVertex->tvh, e->nextVertex->tvh, fh, i));
  tri.remove_constrained_edge(fh, i);
#endif
#ifdef GRID_EDGE_INDEX
  // If the edge is an interior edge, update the index.
  for (InteriorEdge::CellEntryList::iterator it = e->entries.begin();
       it != e->entries.end(); it++)
    it->remove();
  e->unused.splice(e->unused.begin(), e->entries, 
       e->entries.begin(), 
       e->entries.end());
#endif
  // If the edge has neighboring vertices, drop all links.
  if (e->nextVertex.valid()) {
    e->nextVertex->prevIntEdge = IntEdgeHandle();
    e->nextVertex = VertexHandle();
  }
  if (e->prevVertex.valid()) {
    e->prevVertex->nextIntEdge = IntEdgeHandle();
    e->prevVertex = VertexHandle();
  }
  // Remove the edge from its random access vector.
  int n = interiorEdges.size();
  if (n > 1) {
    IntEdgeHandle f = interiorEdges[n - 1];
    interiorEdges[n - 1] = interiorEdges[e->index];
    interiorEdges[e->index] = f;
    f->index = e->index;
  }
  interiorEdges.pop_back();
  e->index = -1;
  intEdgeStore.deleteItem(&*e);
  // Invalidate the handle.
  e = IntEdgeHandle();
}

void Coloring::freeEdge(Coloring::BdEdgeHandle& e) {
  // If the edge has neighboring vertices, drop all links.
  if (e->nextVertex.valid()) {
    e->nextVertex->prevBdEdge = BdEdgeHandle();
    e->nextVertex = VertexHandle();
  }
  if (e->prevVertex.valid()) {
    e->prevVertex->nextBdEdge = BdEdgeHandle();
    e->prevVertex = VertexHandle();
  }
  // Remove the edge from its random access vector.
  int n = boundaryEdges.size();
  if (n > 1) {
    BdEdgeHandle f = boundaryEdges[n - 1];
    boundaryEdges[n - 1] = boundaryEdges[e->index];
    boundaryEdges[e->index] = f;
    f->index = e->index;
  }
  boundaryEdges.pop_back();
  e->index = -1;
  bdEdgeStore.deleteItem(&*e);
  // Invalidate the handle.
  e = BdEdgeHandle();
}

bool Coloring::deleteVertex(VertexHandle v, TEST_IF_VALID_TAG) const {
  TEST_VALIDITY;
  assert(v->type() == Vertex::INTERIOR);
  // Check that v's neighboring vertices are not adjacent.
  if (v->getNextVertex()->hasNextIntEdge() &&
      v->getNextVertex()->getNextVertex() == v->getPrevVertex())
    return false;
  IntEdgeHandle xv = v->getPrevIntEdge();
  IntEdgeHandle vy = v->getNextIntEdge();
  VertexHandle x = xv->getPrevVertex();
  VertexHandle y = vy->getNextVertex();
  Geometry::Point a = *v, b = *x, c = *y;
  // Check that the minimum edge length constraint (if any) is not violated.
  if ((minEdgeLengthSq > 0.0) &&
      (CGAL::to_double(CGAL::squared_distance(b, c)) < minEdgeLengthSq))
    return false;
  // Check to make sure an edge from x to y would not cross any
  // existing edges.  (It cannot possibly cross the edges xv or vy at
  // a point, so the test should not return truth, even if xv and vy
  // are colinear.)
  if (!compatible(b, c))
    return false;
  return true;
}

bool Coloring::deleteVertex(VertexHandle& v, DO_WITHOUT_TEST_TAG) {
  IntEdgeHandle xv = v->getPrevIntEdge();
  IntEdgeHandle vy = v->getNextIntEdge();
  VertexHandle x = xv->getPrevVertex();
  VertexHandle y = vy->getNextVertex();
  Geometry::Point a = *v, b = *x, c = *y;
  // Notify listeners of the elements that will be removed.
  notifyEdgeWillBeRemoved(xv);
  notifyEdgeWillBeRemoved(vy);
  notifyVertexWillBeRemoved(x);
  notifyVertexWillBeRemoved(v);
  notifyVertexWillBeRemoved(y);
  // Perform the update.
  freeEdge(xv);
  freeEdge(vy);
  freeVertex(v);
  IntEdgeHandle xy = newInteriorEdge(x, y);
  // Notify listeners of the new elements.
  notifyEdgeHasBeenAdded(xy);
  notifyVertexHasBeenAdded(x);
  notifyVertexHasBeenAdded(y);
  // Recolor the query points.
  recolored(a, b, c);
  TEST_VALIDITY;
  return true;
}

bool Coloring::deleteVertex(VertexHandle& v, DO_IF_VALID_TAG) {
  if (!deleteVertex(v, TEST_IF_VALID_TAG())) return false;
  return deleteVertex(v, DO_WITHOUT_TEST_TAG());
}

bool Coloring::splitEdge(IntEdgeHandle e, 
       const Geometry::Point& p, 
       TEST_IF_VALID_TAG) const {
  VertexHandle u = e->getPrevVertex();
  VertexHandle v = e->getNextVertex();
  assert(p != *u);
  assert(p != *v);
  // Check that the minimum edge length constraint (if any) is not violated.
  if ((minEdgeLengthSq > 0.0) &&
      ((CGAL::to_double(CGAL::squared_distance(*u, p)) < minEdgeLengthSq) ||
       (CGAL::to_double(CGAL::squared_distance(p, *v)) < minEdgeLengthSq)))
    return false;
  // Note that the two new edges cannot possibly cross the current
  // edge e, so this test correctly determines if the two new edges
  // would cross any edges except e.
  if (!compatible(*u, p) || !compatible(p, *v))
    return false;
  return true;
}

bool Coloring::splitEdge(IntEdgeHandle& e, 
       const Geometry::Point& p, 
       DO_WITHOUT_TEST_TAG) {
  VertexHandle u = e->getPrevVertex();
  VertexHandle v = e->getNextVertex();
  // Notify listeners of the elements that will be removed.
  notifyEdgeWillBeRemoved(e);
  notifyVertexWillBeRemoved(u);
  notifyVertexWillBeRemoved(v);
  // Perform the update
  freeEdge(e);
  VertexHandle x = newVertex(Vertex::INTERIOR, p);
  // Note: the ordering of the vertices is important here.
  IntEdgeHandle ux = newInteriorEdge(u, x);
  IntEdgeHandle xv = newInteriorEdge(x, v);
  // Notify listeners of the new elements.
  notifyVertexHasBeenAdded(x);
  notifyVertexHasBeenAdded(u);
  notifyVertexHasBeenAdded(v);
  notifyEdgeHasBeenAdded(ux);
  notifyEdgeHasBeenAdded(xv);
  // Recolor the query points.
  recolored(*u, *v, p);
  TEST_VALIDITY;
  return true;
}

bool Coloring::splitEdge(IntEdgeHandle& e, 
       const Geometry::Point& p, 
       DO_IF_VALID_TAG) {
  if (!splitEdge(e, p, TEST_IF_VALID_TAG())) return false;
  return splitEdge(e, p, DO_WITHOUT_TEST_TAG());
}

bool Coloring::newInteriorTriangle(const Geometry::Point& xp, 
           const Geometry::Point& yp, 
           const Geometry::Point& zp, 
           TEST_IF_VALID_TAG) const {
#ifdef EXACT_GEOM_CXN
  // Check that all points are not outside the region.
  assert(interior(xp) && interior(yp) && interior(zp));
#endif
  // Check that the minimum edge length constraint (if any) is not violated.
  if ((minEdgeLengthSq > 0.0) &&
      ((CGAL::to_double(CGAL::squared_distance(xp, yp)) < minEdgeLengthSq) ||
       (CGAL::to_double(CGAL::squared_distance(yp, zp)) < minEdgeLengthSq) ||
       (CGAL::to_double(CGAL::squared_distance(zp, xp)) < minEdgeLengthSq)))
    return false;
  // Check that all edges are compatible with the current coloring.
  if (!compatible(xp, yp) ||
      !compatible(yp, zp) || 
      !compatible(zp, xp))
    return false;
  return true;
}

Coloring::VertexHandle 
Coloring::newInteriorTriangle(const Geometry::Point& xp, 
            const Geometry::Point& yp, 
            const Geometry::Point& zp, 
            DO_WITHOUT_TEST_TAG) {
  // Create the vertices.
  VertexHandle x = newVertex(Vertex::INTERIOR, xp);
  VertexHandle y = newVertex(Vertex::INTERIOR, yp);
  VertexHandle z = newVertex(Vertex::INTERIOR, zp);
  // Create the edges.
  IntEdgeHandle xy = newInteriorEdge(x, y);
  IntEdgeHandle yz = newInteriorEdge(y, z);
  IntEdgeHandle zx = newInteriorEdge(z, x);
  // Notify the listeners of the new vertices and edges.
  notifyVertexHasBeenAdded(x);
  notifyVertexHasBeenAdded(y);
  notifyVertexHasBeenAdded(z);
  notifyEdgeHasBeenAdded(xy);
  notifyEdgeHasBeenAdded(yz);
  notifyEdgeHasBeenAdded(zx);
  // Recolor the query points.
  recolored(xp, yp, zp);
  TEST_VALIDITY;
  return x;
}

Coloring::VertexHandle 
Coloring::newInteriorTriangle(const Geometry::Point& xp, 
            const Geometry::Point& yp, 
            const Geometry::Point& zp, 
            DO_IF_VALID_TAG) {
  if (!newInteriorTriangle(xp, yp, zp, TEST_IF_VALID_TAG())) 
    return VertexHandle();
  else return newInteriorTriangle(xp, yp, zp, DO_WITHOUT_TEST_TAG());
}

bool Coloring::newBoundaryTriangle(Coloring::BdEdgeHandle e,
           const Geometry::Point& up, 
           const Geometry::Point& vp, 
           const Geometry::Point& wp, 
           TEST_IF_VALID_TAG) const {
#ifdef EXACT_GEOM_CXN
  // Check that the interior point is inside the region.
  assert(interior(wp));
  // Check that up and vp lie on e.
  assert(e->has_on(up) && e->has_on(vp));
#endif
  // Check that the boundary points are not vertices (i.e., if they
  // are on e, then they are in the interior of e).
  assert(e->source() != up);
  assert(e->source() != vp);
  assert(e->target() != up);
  assert(e->target() != vp);
  // Check that the minimum edge length constraint (if any) is not violated.
  if ((minEdgeLengthSq > 0.0) &&
      ((CGAL::to_double(CGAL::squared_distance(vp, wp)) < minEdgeLengthSq) ||
       (CGAL::to_double(CGAL::squared_distance(wp, up)) < minEdgeLengthSq)))
    return false;
  // Check that both interior edges are compatible with the current coloring.
  if (!compatible(vp, wp) || !compatible(wp, up))
    return false;
  return true;
}

Coloring::VertexHandle 
Coloring::newBoundaryTriangle(Coloring::BdEdgeHandle e,
            const Geometry::Point& up, 
            const Geometry::Point& vp, 
            const Geometry::Point& wp, 
            DO_WITHOUT_TEST_TAG) {
  const int oldNumBdEdges = numBoundaryEdges();
  // Create the vertices.
  VertexHandle u = newVertex(Vertex::BOUNDARY, up);
  VertexHandle v = newVertex(Vertex::BOUNDARY, vp);
  VertexHandle w = newVertex(Vertex::INTERIOR, wp);
  // Create the interior edges.
  IntEdgeHandle uw = newInteriorEdge(u, w);
  IntEdgeHandle wv = newInteriorEdge(w, v);
  // Split the boundary edge e into three parts.
  VertexHandle p = e->getPrevVertex();
  VertexHandle q = e->getNextVertex();
  VertexHandle prevCorner = e->getPrevCornerVertex();
  VertexHandle nextCorner = e->getNextCornerVertex();
  freeEdge(e);
  if (squared_distance(*p, *u) < squared_distance(*p, *v)) {
    newBoundaryEdge(p, u, prevCorner, nextCorner);
    newBoundaryEdge(u, v, prevCorner, nextCorner);
    newBoundaryEdge(v, q, prevCorner, nextCorner);
  } else {
    newBoundaryEdge(p, v, prevCorner, nextCorner);
    newBoundaryEdge(v, u, prevCorner, nextCorner);
    newBoundaryEdge(u, q, prevCorner, nextCorner);
  }
  // Notify the listeners of the new vertices and edges.
  notifyVertexHasBeenAdded(u);
  notifyVertexHasBeenAdded(v);
  notifyVertexHasBeenAdded(w);
  notifyEdgeHasBeenAdded(uw);
  notifyEdgeHasBeenAdded(wv);
  // Recolor the query points.
  recolored(up, vp, wp);
  assert(numBoundaryEdges() - 2 == oldNumBdEdges);
  TEST_VALIDITY;
  return u;
}

Coloring::VertexHandle 
Coloring::newBoundaryTriangle(Coloring::BdEdgeHandle e,
            const Geometry::Point& up, 
            const Geometry::Point& vp, 
            const Geometry::Point& wp, 
            DO_IF_VALID_TAG) {
  if (!newBoundaryTriangle(e, up, vp, wp, TEST_IF_VALID_TAG()))
    return VertexHandle();
  else return newBoundaryTriangle(e, up, vp, wp, DO_WITHOUT_TEST_TAG());
}

bool Coloring::newCornerTriangle(VertexHandle corner, 
         const Geometry::Point& u, 
         const Geometry::Point& v, 
         TEST_IF_VALID_TAG) const {
  TEST_VALIDITY;
  assert(u != *corner);
  assert(v != *corner);
  // Check that the minimum edge length constraint (if any) is not violated.
  if ((minEdgeLengthSq > 0.0) &&
      (CGAL::to_double(CGAL::squared_distance(u, v)) < minEdgeLengthSq))
    return false;
  return compatible(u, v);
}

Coloring::VertexHandle Coloring::newCornerTriangle(VertexHandle corner, 
               const Geometry::Point& u, 
               const Geometry::Point& v, 
               DO_WITHOUT_TEST_TAG) {
  TEST_VALIDITY;
  const int oldNumBdEdges = numBoundaryEdges();
  // Insert the new boundary vertex before the corner.
  BdEdgeHandle prevBdEdge = corner->getPrevBdEdge();
  VertexHandle prevCorner = prevBdEdge->getPrevCornerVertex();
#ifdef EXACT_GEOM_CXN
  assert(prevBdEdge->has_on(u));
#endif
  VertexHandle prevBdVertex = prevBdEdge->getPrevVertex();
  VertexHandle newPrevBdVertex = newVertex(Vertex::BOUNDARY, u);
  freeEdge(prevBdEdge);
  prevBdEdge = newBoundaryEdge(newPrevBdVertex, corner, prevCorner, corner);
  newBoundaryEdge(prevBdVertex, newPrevBdVertex, prevCorner, corner);
  // Insert the new boundary vertex after the corner.
  BdEdgeHandle nextBdEdge = corner->getNextBdEdge();
  VertexHandle nextCorner = nextBdEdge->getNextCornerVertex();
#ifdef EXACT_GEOM_CXN
  assert(nextBdEdge->has_on(v));
#endif
  VertexHandle nextBdVertex = nextBdEdge->getNextVertex();
  VertexHandle newNextBdVertex = newVertex(Vertex::BOUNDARY, v);
  freeEdge(nextBdEdge);
  nextBdEdge = newBoundaryEdge(corner, newNextBdVertex, corner, nextCorner);
  newBoundaryEdge(newNextBdVertex, nextBdVertex, corner, nextCorner);
  // Join them.
  IntEdgeHandle interiorEdge = 
    newInteriorEdge(newPrevBdVertex, newNextBdVertex);

  // Notify the listeners of the new vertices and edges.
  notifyVertexHasBeenAdded(newPrevBdVertex);
  notifyVertexHasBeenAdded(newNextBdVertex);
  notifyEdgeHasBeenAdded(interiorEdge);
  // Recolor the query points.
  recolored(u, v, *corner);
  TEST_VALIDITY;
  assert(numBoundaryEdges() - 2 == oldNumBdEdges);
  return newPrevBdVertex;
}

Coloring::VertexHandle Coloring::newCornerTriangle(VertexHandle corner, 
               const Geometry::Point& u, 
               const Geometry::Point& v, 
               DO_IF_VALID_TAG) {
  if (!newCornerTriangle(corner, u, v, TEST_IF_VALID_TAG())) return false;
  else return newCornerTriangle(corner, u, v, DO_WITHOUT_TEST_TAG());
}

bool Coloring::deleteIntTriangle(VertexHandle vh, TEST_IF_VALID_TAG) const {
  // Ensure this vertex is part of an interior triangle.
  if (vh->type() != Vertex::INTERIOR) return false;
  if (!vh->hasNextIntEdge()) return false;
  VertexHandle vi = vh->getNextVertex();
  if (vi->type() != Vertex::INTERIOR) return false;
  if (!vi->hasNextIntEdge()) return false;
  VertexHandle vj = vi->getNextVertex();
  if (vj->type() != Vertex::INTERIOR) return false;
  if (!vj->hasNextIntEdge()) return false;
  if (vj->getNextVertex() != vh) return false;
  return true;
}

bool Coloring::deleteIntTriangle(VertexHandle& vh, DO_WITHOUT_TEST_TAG) {
  // Ensure this vertex is part of a triangle, and extract the vertices.
  VertexHandle va = vh;
  IntEdgeHandle ab = va->getNextIntEdge();
  VertexHandle vb = ab->getNextVertex();
  IntEdgeHandle bc = vb->getNextIntEdge();
  VertexHandle vc = bc->getNextVertex();
  IntEdgeHandle ca = vc->getNextIntEdge();
  Geometry::Point pa = va->point();
  Geometry::Point pb = vb->point();
  Geometry::Point pc = vc->point();
  // Notify the listeners of the removed vertices and edges.
  notifyEdgeWillBeRemoved(ab);
  notifyEdgeWillBeRemoved(bc);
  notifyEdgeWillBeRemoved(ca);
  notifyVertexWillBeRemoved(va);
  notifyVertexWillBeRemoved(vb);
  notifyVertexWillBeRemoved(vc);
  // Now free the vertices and edges.
  freeEdge(ab);
  freeEdge(bc);
  freeEdge(ca);
  freeVertex(va);
  freeVertex(vb);
  freeVertex(vc);
  vh = VertexHandle();
  // Recolor the query points.
  recolored(pa, pb, pc);
  TEST_VALIDITY;
  return true;
}

bool Coloring::deleteIntTriangle(VertexHandle& vh, DO_IF_VALID_TAG) {
  if (!deleteIntTriangle(vh, TEST_IF_VALID_TAG())) return false;
  else return deleteIntTriangle(vh, DO_WITHOUT_TEST_TAG());
}

bool Coloring::deleteBdTriangle(VertexHandle vh, TEST_IF_VALID_TAG) const {
  // Try to back up the vertex handle to the root of a chain.
  if (vh->hasPrevIntEdge()) vh = vh->getPrevVertex();
  if (vh->hasPrevIntEdge()) vh = vh->getPrevVertex();
  // If we are not at a root of the chain, this is not a boundary
  // triangle.
  if (vh->type() != Vertex::BOUNDARY) return false;
  // Check to see if this is a boundary triangle.
  VertexHandle vi = vh->getNextVertex();
  if (vi->type() != Vertex::INTERIOR) return false;
  if (!vi->hasNextIntEdge()) return false;
  VertexHandle vj = vi->getNextVertex();
  if (vj->type() != Vertex::BOUNDARY) return false;
  // Check to make sure the boundary vertices are on the same face.
  if (vh->getNextCornerVertex() != vj->getNextCornerVertex()) 
    return false;
  return true;
}

void Coloring::disconnectFromBd(VertexHandle v) {
  assert(v->type() == Vertex::BOUNDARY);
  BdEdgeHandle prevBoundaryEdge = v->getPrevBdEdge();
  BdEdgeHandle nextBoundaryEdge = v->getNextBdEdge();
  VertexHandle prevVertex = prevBoundaryEdge->getPrevVertex();
  VertexHandle nextVertex = nextBoundaryEdge->getNextVertex();
  VertexHandle prevCorner = prevBoundaryEdge->getPrevCornerVertex();
  VertexHandle nextCorner = nextBoundaryEdge->getNextCornerVertex();
  freeEdge(prevBoundaryEdge);
  freeEdge(nextBoundaryEdge);
  newBoundaryEdge(prevVertex, nextVertex, prevCorner, nextCorner);
}

bool Coloring::deleteBdTriangle(VertexHandle& vh, DO_WITHOUT_TEST_TAG) {
  // Try to back up the vertex handle to the root of a chain.
  if (vh->hasPrevIntEdge()) vh = vh->getPrevVertex();
  if (vh->hasPrevIntEdge()) vh = vh->getPrevVertex();
  VertexHandle va = vh;
  IntEdgeHandle ab = va->getNextIntEdge();
  VertexHandle vb = ab->getNextVertex();
  IntEdgeHandle bc = vb->getNextIntEdge();
  VertexHandle vc = bc->getNextVertex();
  Geometry::Point pa = va->point();
  Geometry::Point pb = vb->point();
  Geometry::Point pc = vc->point();
  // Notify the listeners of the removed vertices and edges.
  notifyEdgeWillBeRemoved(ab);
  notifyEdgeWillBeRemoved(bc);
  notifyVertexWillBeRemoved(va);
  notifyVertexWillBeRemoved(vb);
  notifyVertexWillBeRemoved(vc);
  // Now free the vertices and edges.
  freeEdge(ab);
  freeEdge(bc);
  disconnectFromBd(va);
  disconnectFromBd(vc);
  freeVertex(va);
  freeVertex(vb);
  freeVertex(vc);
  vh = VertexHandle();
  // Recolor the query points.
  recolored(pa, pb, pc);
  TEST_VALIDITY;
  return true;
}

bool Coloring::deleteBdTriangle(VertexHandle& vh, DO_IF_VALID_TAG) {
  if (!deleteBdTriangle(vh, TEST_IF_VALID_TAG())) return false;
  else return deleteBdTriangle(vh, DO_WITHOUT_TEST_TAG());
}

bool Coloring::deleteCornerTriangle(VertexHandle vh, TEST_IF_VALID_TAG) const {
  // Try to back up the vertex handle to the root of a chain.
  if (vh->hasPrevIntEdge()) vh = vh->getPrevVertex();
  // If we are not at a root of the chain, this is not a corner
  // triangle.
  if (vh->type() != Vertex::BOUNDARY) return false;
  // Check to see if this is a boundary triangle.
  VertexHandle vi = vh->getNextVertex();
  if (vi->type() != Vertex::BOUNDARY) return false;
  return true;
}

bool Coloring::deleteCornerTriangle(VertexHandle& vh, DO_WITHOUT_TEST_TAG) {
  // Back up the vertex handle to the root of a chain.
  if (vh->hasPrevIntEdge()) vh = vh->getPrevVertex();
  VertexHandle va = vh;
  IntEdgeHandle ab = va->getNextIntEdge();
  VertexHandle vb = ab->getNextVertex();
  VertexHandle vaNextCorner = va->getNextCornerVertex();
  VertexHandle vbNextCorner = vb->getNextCornerVertex();
  Geometry::Point pa = va->point();
  Geometry::Point pb = vb->point();
  Geometry::Point pc;
  if (vaNextCorner->getNextCornerVertex() == vbNextCorner)
    pc = vaNextCorner->point();
  else if (vbNextCorner->getNextCornerVertex() == vaNextCorner)
    pc = vbNextCorner->point();
  // Notify the listeners of the removed vertices and edges.
  notifyEdgeWillBeRemoved(ab);
  notifyVertexWillBeRemoved(va);
  notifyVertexWillBeRemoved(vb);
  // Now free the vertices and edges.
  freeEdge(ab);
  disconnectFromBd(va);
  disconnectFromBd(vb);
  freeVertex(va);
  freeVertex(vb);
  vh = VertexHandle();
  // Recolor the query points.
  recolored(pa, pb, pc);
  TEST_VALIDITY;
  return true;
}

bool Coloring::deleteCornerTriangle(VertexHandle& vh, DO_IF_VALID_TAG) {
  if (!deleteCornerTriangle(vh, TEST_IF_VALID_TAG())) return false;
  else return deleteCornerTriangle(vh, DO_WITHOUT_TEST_TAG());
}

bool Coloring::moveIntVertex(Coloring::VertexHandle v, 
           const Geometry::Point& p, 
           TEST_IF_VALID_TAG) const {
  assert(v->type() == Vertex::INTERIOR);
#ifdef EXACT_GEOM_CXN
  assert(interior(p));
#endif
  if (v->point() == p) return true;
  // Check that the minimum edge length constraint (if any) is not violated.
  if ((minEdgeLengthSq > 0.0) &&
      ((CGAL::to_double(CGAL::squared_distance(p, v->getNextVertex()->point())) < minEdgeLengthSq) ||
       (CGAL::to_double(CGAL::squared_distance(p, v->getPrevVertex()->point())) < minEdgeLengthSq)))
    return false;
  if (!compatible(p, v->getNextVertex()->point()) ||
      !compatible(p, v->getPrevVertex()->point()))
    return false;
  return true;
}

bool Coloring::moveIntVertex(Coloring::VertexHandle& v, 
           const Geometry::Point& p, 
           DO_WITHOUT_TEST_TAG) {
  IntEdgeHandle uv = v->getPrevIntEdge();
  VertexHandle u = uv->getPrevVertex();
  IntEdgeHandle vw = v->getNextIntEdge();
  VertexHandle w = vw->getNextVertex();
  // Notify the listeners of the removed vertices and edges.
  notifyEdgeWillBeRemoved(uv);
  notifyEdgeWillBeRemoved(vw);
  notifyVertexWillBeRemoved(u);
  notifyVertexWillBeRemoved(v);
  notifyVertexWillBeRemoved(w);
  // Perform the update.
  Geometry::Point oldLoc = *v;
  freeEdge(uv);
  freeEdge(vw);
  freeVertex(v);
  v = newVertex(Vertex::INTERIOR, p);
  uv = newInteriorEdge(u, v);
  vw = newInteriorEdge(v, w);
  // Notify the listeners of the new vertices and edges.
  notifyVertexHasBeenAdded(u);
  notifyVertexHasBeenAdded(v);
  notifyVertexHasBeenAdded(w);
  notifyEdgeHasBeenAdded(uv);
  notifyEdgeHasBeenAdded(vw);
  // Recolor the query points.
  recolored(*u, p, *w, oldLoc);
  TEST_VALIDITY;
  return true;
}

bool Coloring::moveIntVertex(Coloring::VertexHandle& v, 
           const Geometry::Point& p, 
           DO_IF_VALID_TAG) {
  if (!moveIntVertex(v, p, TEST_IF_VALID_TAG())) return false;
  return moveIntVertex(v, p, DO_WITHOUT_TEST_TAG());
}

bool Coloring::moveVertexAlongBd(Coloring::VertexHandle v, 
         const Geometry::Point& p, 
         TEST_IF_VALID_TAG) const {
  TEST_VALIDITY;
  // Check that the vertex is a boundary vertex and that the new
  // location is on one or the other incident boundary edge.
  assert(v->type() == Vertex::BOUNDARY);
#ifdef EXACT_GEOM_CXN
  assert(v->getPrevBdEdge()->has_on(p) || v->getNextBdEdge()->has_on(p));
#endif
  if (v->point() == p) return true;
  // Check that the interior edge incident to v will not cross any
  // other edges in the coloring and that the minimum edge length
  // constraint (if any) is not violated.
  bool isRoot = v->hasNextIntEdge();
  IntEdgeHandle e = isRoot ? v->getNextIntEdge() : v->getPrevIntEdge();
  VertexHandle u = isRoot ? e->getNextVertex() : e->getPrevVertex();
  if ((minEdgeLengthSq > 0.0) &&
      (CGAL::to_double(CGAL::squared_distance(*u, p)) < minEdgeLengthSq))
    return false;
  if (!compatible(*u, p)) return false;
  return true;
}
  
bool Coloring::moveVertexAlongBd(Coloring::VertexHandle& v, 
         const Geometry::Point& p, 
         DO_WITHOUT_TEST_TAG) {
  const int oldNumBdEdges = numBoundaryEdges();
  // Check that the interior edge incident to v will not cross any
  // other edges in the coloring.
  bool isRoot = v->hasNextIntEdge();
  IntEdgeHandle e = isRoot ? v->getNextIntEdge() : v->getPrevIntEdge();
  VertexHandle u = isRoot ? e->getNextVertex() : e->getPrevVertex();
  // Store some information.
  const Geometry::Point oldLoc = v->point();
  BdEdgeHandle prevBdEdge = v->getPrevBdEdge();
  BdEdgeHandle nextBdEdge = v->getNextBdEdge();
  VertexHandle prevCorner = nextBdEdge->getPrevCornerVertex();
  VertexHandle nextCorner = nextBdEdge->getNextCornerVertex();
  VertexHandle nextBdVertex = nextBdEdge->getNextVertex();
  VertexHandle prevBdVertex = prevBdEdge->getPrevVertex();
  // Notify the listeners of the removed vertices and edges.
  notifyEdgeWillBeRemoved(e);
  notifyVertexWillBeRemoved(u);
  notifyVertexWillBeRemoved(v);
  // Perform the update.  Free the old edges and vertex.
  freeEdge(e);
  freeEdge(prevBdEdge);
  freeEdge(nextBdEdge);
  freeVertex(v);
  // Create the new edges and vertex.
  v = newVertex(Vertex::BOUNDARY, p);
  e = isRoot ? newInteriorEdge(v, u) : newInteriorEdge(u, v);
  prevBdEdge = newBoundaryEdge(prevBdVertex, v, prevCorner, nextCorner);
  nextBdEdge = newBoundaryEdge(v, nextBdVertex, prevCorner, nextCorner);
  // Notify the listeners of the new vertices and edges.
  notifyVertexHasBeenAdded(u);
  notifyVertexHasBeenAdded(v);
  notifyEdgeHasBeenAdded(e);
  // Recolor the query points.
  recolored(*u, oldLoc, p);
  assert(oldNumBdEdges == numBoundaryEdges());
  TEST_VALIDITY;  
  return true;
}

bool Coloring::moveVertexAlongBd(Coloring::VertexHandle& v, 
         const Geometry::Point& p, 
         DO_IF_VALID_TAG) {
  if (!moveVertexAlongBd(v, p, TEST_IF_VALID_TAG())) return false;
  return moveVertexAlongBd(v, p, DO_WITHOUT_TEST_TAG());
}

bool Coloring::moveBdVertexPastCorner(Coloring::BdEdgeHandle bdEdge, 
              const Geometry::Point& p, 
              TEST_IF_VALID_TAG) const {
  TEST_VALIDITY;
  bool prevIsCorner = 
    (bdEdge->getPrevVertex()->type() == Coloring::Vertex::CORNER);
  Coloring::VertexHandle corner = 
    (prevIsCorner ? bdEdge->getPrevVertex() : bdEdge->getNextVertex());
  Coloring::VertexHandle bdVertex = 
    (prevIsCorner ? bdEdge->getNextVertex() : bdEdge->getPrevVertex());
  assert(bdVertex->type() == Coloring::Vertex::BOUNDARY);
#ifdef EXACT_GEOM_CXN
  Coloring::BdEdgeHandle otherEdge =
    (prevIsCorner ? corner->getPrevBdEdge() : corner->getNextBdEdge());
  assert(otherEdge->has_on(p));
#endif
  bool isRoot = bdVertex->hasNextIntEdge();
  Coloring::IntEdgeHandle intEdge = (isRoot ? 
             bdVertex->getNextIntEdge() : 
             bdVertex->getPrevIntEdge());
  Coloring::VertexHandle adjVertex = (isRoot ? 
         intEdge->getNextVertex() : 
         intEdge->getPrevVertex());
  // If the adjacent edge is a boundary vertex on the target face,
  // abort the update.
  if ((adjVertex->type() == Vertex::BOUNDARY) &&
      ((prevIsCorner && (adjVertex->getNextCornerVertex() == corner)) ||
       (!prevIsCorner && (adjVertex->getPrevCornerVertex() == corner))))
    return false;
  // Check that the minimum edge length constraint (if any) is not violated.
  if ((minEdgeLengthSq > 0.0) &&
      (CGAL::to_double(CGAL::squared_distance(p, *adjVertex)) < minEdgeLengthSq))
    return false;
  // Check that the new edge would not cross any existing edges in the
  // coloring.
  if (!compatible(p, *adjVertex)) return false;
  return true;
}

bool Coloring::moveBdVertexPastCorner(Coloring::BdEdgeHandle& bdEdge, 
              const Geometry::Point& p, 
              DO_WITHOUT_TEST_TAG) {
  const int oldNumBdEdges = numBoundaryEdges();
  bool prevIsCorner = 
    (bdEdge->getPrevVertex()->type() == Coloring::Vertex::CORNER);
  Coloring::VertexHandle corner = 
    (prevIsCorner ? bdEdge->getPrevVertex() : bdEdge->getNextVertex());
  Coloring::VertexHandle bdVertex = 
    (prevIsCorner ? bdEdge->getNextVertex() : bdEdge->getPrevVertex());
  assert(bdVertex->type() == Coloring::Vertex::BOUNDARY);
  bool isRoot = bdVertex->hasNextIntEdge();
  Coloring::IntEdgeHandle intEdge = (isRoot ? 
             bdVertex->getNextIntEdge() : 
             bdVertex->getPrevIntEdge());
  Coloring::VertexHandle adjVertex = (isRoot ? 
         intEdge->getNextVertex() : 
         intEdge->getPrevVertex());
  // Notify the listeners of the removed vertices and edges.
  notifyEdgeWillBeRemoved(intEdge);
  notifyVertexWillBeRemoved(bdVertex);
  notifyVertexWillBeRemoved(adjVertex);
  // Perform the update.  
  Geometry::Point oldLoc = *bdVertex;
  Coloring::BdEdgeHandle nextBdEdge = bdEdge->getNextVertex()->getNextBdEdge();
  Coloring::VertexHandle nextBdVertex = nextBdEdge->getNextVertex();
  Coloring::BdEdgeHandle prevBdEdge = bdEdge->getPrevVertex()->getPrevBdEdge();
  Coloring::VertexHandle prevBdVertex = prevBdEdge->getPrevVertex();
  Coloring::VertexHandle prevCorner = corner->getPrevCornerVertex();
  Coloring::VertexHandle nextCorner = corner->getNextCornerVertex();
  // Remove the existing edges and vertex.
  freeEdge(prevBdEdge);
  freeEdge(bdEdge);
  freeEdge(nextBdEdge);
  freeEdge(intEdge);
  freeVertex(bdVertex);
  // Insert the new edges and vertex.
  bdVertex = newVertex(Vertex::BOUNDARY, p);
  intEdge = isRoot ? newInteriorEdge(bdVertex, adjVertex) :
    newInteriorEdge(adjVertex, bdVertex);
  if (prevIsCorner) {
    newBoundaryEdge(prevBdVertex, bdVertex, prevCorner, corner);
    bdEdge = newBoundaryEdge(bdVertex, corner, prevCorner, corner);
    newBoundaryEdge(corner, nextBdVertex, corner, nextCorner);
  } else {
    newBoundaryEdge(prevBdVertex, corner, prevCorner, corner);
    bdEdge = newBoundaryEdge(corner, bdVertex, corner, nextCorner);
    newBoundaryEdge(bdVertex, nextBdVertex, corner, nextCorner);
  }
  // Notify the listeners of the new vertices and edges.
  notifyVertexHasBeenAdded(bdVertex);
  notifyVertexHasBeenAdded(adjVertex);
  notifyEdgeHasBeenAdded(intEdge);
  // Recolor the query points.
  recolored(*adjVertex, p, *corner, oldLoc);
  TEST_VALIDITY;
  assert(oldNumBdEdges == numBoundaryEdges());
  return true;
}

bool Coloring::moveBdVertexPastCorner(Coloring::BdEdgeHandle& bdEdge, 
              const Geometry::Point& p, 
              DO_IF_VALID_TAG) {
  if (!moveBdVertexPastCorner(bdEdge, p, TEST_IF_VALID_TAG())) return false;
  return moveBdVertexPastCorner(bdEdge, p, DO_WITHOUT_TEST_TAG());
}

bool Coloring::recolorQuadrilateral(IntEdgeHandle ab, 
            IntEdgeHandle cd, 
            bool ac, 
            TEST_IF_VALID_TAG) const {
  VertexHandle a = ab->getPrevVertex();
  VertexHandle b = ab->getNextVertex();
  VertexHandle c = cd->getPrevVertex();
  VertexHandle d = cd->getNextVertex();
  /*
   * Test to make sure the update would yield a valid coloring.  The
   * edges we will create are {w,x} and {y,z}.
   */
  VertexHandle w = a;
  VertexHandle x = (ac ? c : d);
  VertexHandle y = b;
  VertexHandle z = (ac ? d : c);
  // Self-edges are not allowed. (This can happen when the edges are
  // adjacent.)
  if ((w == x) || (y == z)) return false;
  // The two new edges may not cross.
  if (CGAL::do_intersect(Segment(*w, *x), Segment(*y, *z))) return false;
  // The new edges must not duplicate existing edges.  This can occur
  // when the original edges are separated by one edge.
  if ((d->hasNextIntEdge() && (d->getNextVertex() == a)) 
      ||
      (b->hasNextIntEdge() && (b->getNextVertex() == c)))
    return false;
  // The new edges must not join two boundary vertices on the same
  // face.
  if ((w->type() == Vertex::BOUNDARY) && 
      (x->type() == Vertex::BOUNDARY) &&
      (x->getNextCornerVertex() == w->getNextCornerVertex()))
    return false;
  if ((y->type() == Vertex::BOUNDARY) && 
      (z->type() == Vertex::BOUNDARY) &&
      (y->getNextCornerVertex() == z->getNextCornerVertex()))
    return false;
  // Check that the minimum edge length constraint (if any) is not violated.
  if ((minEdgeLengthSq > 0.0) &&
      ((CGAL::to_double(CGAL::squared_distance(*w, *x)) < minEdgeLengthSq) ||
       (CGAL::to_double(CGAL::squared_distance(*y, *z)) < minEdgeLengthSq)))
    return false;
  // The two new edges must be compatible with the current coloring.
  // (They cannot possibly cross the two removed edges because of the
  // shared vertices.)
  if (!compatible(*w, *x) || !compatible(*y, *z)) return false;
  return true;
}

bool Coloring::recolorQuadrilateral(IntEdgeHandle& ab, 
            IntEdgeHandle& cd, 
            bool ac, 
            DO_WITHOUT_TEST_TAG) {
  VertexHandle a = ab->getPrevVertex();
  VertexHandle b = ab->getNextVertex();
  VertexHandle c = cd->getPrevVertex();
  VertexHandle d = cd->getNextVertex();
  // Notify the listeners of the removed vertices and edges.
  notifyEdgeWillBeRemoved(ab);
  notifyEdgeWillBeRemoved(cd);
  notifyVertexWillBeRemoved(a);
  notifyVertexWillBeRemoved(b);
  notifyVertexWillBeRemoved(c);
  notifyVertexWillBeRemoved(d);
  if (ac) {
    // The two sources (a and c) are to be joined, which means we must
    // reverse some edges.  To determine which must be reversed, we
    // must first figure out if ab and cd are part of the same path.
    // The first step is to remove the current edges.
    freeEdge(ab);
    freeEdge(cd);
    // Determine if b leads to c.
    VertexHandle v = b;
    while ((v != c) && v->hasNextIntEdge()) v = v->getNextVertex();
    const bool b2c = (v == c);
    // Determine if d leads to a.
    v = d;
    while ((v != a) && v->hasNextIntEdge()) v = v->getNextVertex();
    const bool d2a = (v == a);
    /* There are four cases to consider: 
     * 
     * 1. If b leads to c and d leads to a, then both edges are part
     * of a cycle.  In this case, we reverse all edges on the path
     * from b to c, connect a to c, and connect b to d.
     *
     * 2. If b leads to c but d does not lead to a, then both edges
     * are part of a chain in which ab precedes cd.  In this case, we
     * reverse all edges on the path from b to c, connect a to c, and
     * connect b to d.
     * 
     * 3. If d leads to a but b does not lead to c, then both edges
     * are part of a chain in which cd precedes ab.  In this case, we
     * reverse all edges on the path from d to a, connect c to a, and
     * connect d to b.
     *
     * 4. If b does not lead to c and d does not lead to a then ab and
     * cd are not part of the same path.  In this case we reverse all
     * edges before c and after d and then connect a to c and d to b.
     */
    // In cases 1, 2, and 4 we reverse all edges from b to c.
    if (!(d2a && !b2c)) {
      // Walk backwards from c, reversing all edges until we run out of
      // edges.  (Because we've already removed the edges entering b,
      // this stops correctly)
      v = c;
      VertexHandle u = 
  v->hasPrevIntEdge() ? v->getPrevVertex() : VertexHandle();
      if (u.valid()) {
  // Free the edge to v from u.
  IntEdgeHandle uv = v->getPrevIntEdge();
  freeEdge(uv);
      }
      // We have now set up the loop invariant.  If u is valid, then it
      // was the previous vertex before v, but we've removed that edge.
      while (u.valid()) {
  // Get the vertex before u.
  VertexHandle t = 
    u->hasPrevIntEdge() ? u->getPrevVertex() : VertexHandle();
  if (t.valid()) {
    // Free the edge to u from t.
    IntEdgeHandle tu = u->getPrevIntEdge();
    freeEdge(tu);
  }
  // Add an edge to u from v.
  newInteriorEdge(v, u);
  // Regress the handles.
  v = u;
  u = t;
      }
    }
    // In cases 3 and 4 we reverse the edges from d to a.
    if (!b2c) {
      // Walk forwards from d, reversing all edges until we run out of
      // edges.  Because we have already removed the edge from a to c
      // this will stop correctly.
      v = d;
      VertexHandle u = 
  v->hasNextIntEdge() ? v->getNextVertex() : VertexHandle();
      if (u.valid()) {
  // Free the edge from v to u.
  IntEdgeHandle vu = v->getNextIntEdge();
  freeEdge(vu);
      }
      // We have now set up the loop invariant.  If u is valid, then it
      // was the next vertex after v, but we've removed that edge.
      while (u.valid()) {
  // Get the vertex after u.
  VertexHandle t = 
    u->hasNextIntEdge() ? u->getNextVertex() : VertexHandle();
  if (t.valid()) {
    // Free the edge from u to t.
    IntEdgeHandle ut = u->getNextIntEdge();
    freeEdge(ut);
  }
  // Add an edge from u back to v.
  newInteriorEdge(u, v);
  // Advance the handles.
  v = u;
  u = t;
      }
    }
    // Now we're ready to install the new edges.
    if (b2c && d2a) { // case 1
      ab = newInteriorEdge(a, c);
      cd = newInteriorEdge(b, d);
    } else if (b2c && !d2a) { // case 2
      ab = newInteriorEdge(a, c);
      cd = newInteriorEdge(b, d);
    } else if (!b2c && d2a) { // case 3
      ab = newInteriorEdge(c, a);
      cd = newInteriorEdge(d, b);
    } else { // case 4
      ab = newInteriorEdge(a, c);
      cd = newInteriorEdge(d, b);
    }
  } else {
    // We are joining sources to targets, so no reversing is required.
    freeEdge(ab);
    freeEdge(cd);
    ab = newInteriorEdge(a, d);
    cd = newInteriorEdge(c, b);
  }
  // Notify the listeners of the new vertices and edges.
  notifyVertexHasBeenAdded(a);
  notifyVertexHasBeenAdded(b);
  notifyVertexHasBeenAdded(c);
  notifyVertexHasBeenAdded(d);
  notifyEdgeHasBeenAdded(ab);
  notifyEdgeHasBeenAdded(cd);
  // Recolor the query points.
  if (ac) 
    recolored(*b, *a, *c, *d);
  else 
    recolored(*b, *a, *d, *c);
  TEST_VALIDITY;
  return true;
}

bool Coloring::recolorQuadrilateral(IntEdgeHandle& ab, 
            IntEdgeHandle& cd, 
            bool ac, 
            DO_IF_VALID_TAG) {
  if (!recolorQuadrilateral(ab, cd, ac, TEST_IF_VALID_TAG())) return false;
  return recolorQuadrilateral(ab, cd, ac, DO_WITHOUT_TEST_TAG());
}

#ifdef TRIANGULATION_EDGE_INDEX
Coloring::Tracer::Tracer(const Geometry::Point& source,
       const Geometry::Point& target,
       const Coloring::TriangulationIndex& tri,
       bool isRay,
       const Coloring::TFaceHandle hint) 
  : source(source), target(target), tri(&tri), circ(),
    edge(Coloring::TFaceHandle(), 0), 
    forwards(true), valid(true), isRay(isRay) { 
  // Initialize the face handle to a face containing the source.
  TFaceHandle fh = tri.locate(this->source, hint);
  // Initialize the face circulator at this face handle.
  Coloring::TFaceCirculator start = circ = tri.line_walk(source, target, fh);
  assert(tri.is_infinite(start) == false);
  // Check to see if the object is a segment and this face contains
  // the target.  If it does, then this tracer is initialized to the
  // end iterator.
  if (!isRay && 
      (tri.triangle(start).has_on_unbounded_side(target) == false)) {
    valid = false;
    return;
  }
  // Test to see if we move towards the target.
  ++circ;
  int nbr;
  assert(start->has_neighbor(circ, nbr));
  Geometry::Segment e = tri.segment(start, nbr);
  if (tri.is_infinite(circ) == false) {
    Geometry::Point p = start->vertex(tri.cw(nbr))->point();
    Geometry::Point q = start->vertex(tri.ccw(nbr))->point();
    Geometry::Point s = source - (target - source);
    Geometry::Point t = target + (target - source);
    CGAL::Orientation so = CGAL::orientation(p, q, s);
    CGAL::Orientation to = CGAL::orientation(p, q, t);
    assert(so != CGAL::COLLINEAR);
    assert(to != CGAL::COLLINEAR);
    forwards = (so != to);
  } else 
    forwards = false;
  // Now that we know which way it will walk, reset the circulator.
  circ = tri.line_walk(source, target, fh);
  // Increment the iterator.
  this->operator++();
}

Coloring::Tracer& Coloring::Tracer::operator++() {
  assert(valid);
  // Advance the underlying circulator until a constrained edge is
  // crossed.
  while (true) {
    assert(tri->is_infinite(circ) == false);
    Coloring::TFaceCirculator prev = circ;
    if (forwards) ++circ; else --circ;
    // If the shared edge between the current and next faces is a
    // constrained edge that crosses the segment, stop the search.
    int nbr;
    if (prev->has_neighbor(circ, nbr)) {
      // Check to see if the shared edge is a constrained edge that
      // crosses the segment.
      if (prev->is_constrained(nbr) && 
    crosses(tri->segment(prev, nbr), source, target)) {
  edge = TEdge(prev, nbr);
  return *this;
      }
      // Check if we are done with the traversal.
      if (!isRay) {
  Geometry::Point p = prev->vertex(tri->cw(nbr))->point();
  Geometry::Point q = prev->vertex(tri->ccw(nbr))->point();
  CGAL::Orientation so = CGAL::orientation(p, q, source);
  CGAL::Orientation to = CGAL::orientation(p, q, target);
  if ((to == CGAL::COLLINEAR) || (so == to)) {
    valid = false;
    return *this;
  }
      }
    } else {
      // The current and next faces do not share an edge.  In this
      // case, they share a vertex that is on the search line.  
    }
    // Check if we have passed out of the convex hull.
    if (tri->is_infinite(circ)) {
      valid = false;
      return *this;
    }
  }
}
#endif // #ifdef TRIANGULATION_EDGE_INDEX

#ifdef GRID_EDGE_INDEX
bool Coloring::compatible(const Geometry::Point &xp, 
        const Geometry::Point &yp) const {
  typedef InteriorEdgeIndex::ConstLineCellIterator Iterator; 
  Iterator it(*intEdgeIndex, xp, yp, Iterator::SEGMENT()), end;
  typedef InteriorEdgeIndex::Cell Cell;
  while (it != end) {
    const InteriorEdgeIndex::Cell& cell = *it;
    const InteriorEdgeIndex::Cell::ItemList& edges = cell.getItemList();
    typedef InteriorEdgeIndex::Cell::ItemList::const_iterator EdgeIterator;
    for (EdgeIterator jt = edges.begin(); jt != edges.end(); jt++) {
      const IntEdgeHandle e = *jt;
      if (e->crosses(xp, yp)) {
  return false;
      }
    }
    ++it;
  }
  return true;
}
#endif // #ifdef GRID_EDGE_INDEX

#ifdef TRIANGULATION_EDGE_INDEX
bool Coloring::compatible(const Geometry::Point &xp, 
        const Geometry::Point &yp,
        TFaceHandle fh) const {
  Tracer trace(xp, yp, tri, outside(yp), fh);
  bool found = (trace != Tracer());
  return !found;
}
#endif // #ifdef TRIANGULATION_EDGE_INDEX

bool Coloring::trace(const Geometry::Point &xp, 
         const Geometry::Point &yp,
         IntEdgeHandle& edge,
         Geometry::Point& ipt,
         Geometry::Kernel::FT& sd) const {
#ifdef TRIANGULATION_EDGE_INDEX
  Tracer trace(xp, yp, tri);
  if (trace != Tracer()) {
    TFaceHandle fh = trace->first;
    int i = trace->second;
    TVertexHandle tva = fh->vertex(tri.cw(i));
    TVertexHandle tvb = fh->vertex(tri.ccw(i));
    VertexHandle va = tva->vertex;
    VertexHandle vb = tvb->vertex;
    if (va->hasPrevIntEdge() && (va->getPrevVertex() == vb))
      edge = va->getPrevIntEdge();
    else if (va->hasNextIntEdge() && (va->getNextVertex() == vb))
      edge = va->getNextIntEdge();
    else 
      assert(false);
    CGAL::Object result = CGAL::intersection(edge->segment(), 
               Geometry::Segment(xp, yp));
    assert(CGAL::assign(ipt, result));
    sd = squared_distance(ipt, xp);
    return true;
  } else
    return false;
#endif
#ifdef GRID_EDGE_INDEX
  typedef InteriorEdgeIndex::ConstLineCellIterator Iterator; 
  Iterator it(*intEdgeIndex, xp, yp, Iterator::SEGMENT()), end;
  typedef const InteriorEdgeIndex::Cell Cell;
  const Segment segment(xp, yp);
  edge = IntEdgeHandle();
  while (it != end) {
    Cell& cell = *it;
    const Cell::ItemList& edges = cell.getItemList();
    typedef Cell::ItemList::const_iterator EdgeIterator;
    for (EdgeIterator jt = edges.begin(); jt != edges.end(); jt++) {
      const IntEdgeHandle e = *jt;
      Geometry::Point intersection;
      CGAL::Object result = CGAL::intersection(*e, segment);
      if (CGAL::assign(intersection, result)) {
  Geometry::Kernel::FT d = squared_distance(intersection, xp);
  if ((!edge.valid()) || (d < sd)) {
    edge = e;
    sd = d;
    ipt = intersection;
  }
      }
    }
    if (edge.valid()) break;
    ++it;
  }
  return edge.valid();
#endif
}

Geometry::Point Coloring::randomPoint(Arak::Util::Random& random) const {
  return Geometry::Point(random.uniform(CGAL::to_double(bd.xmin()), 
          CGAL::to_double(bd.xmax())),
       random.uniform(CGAL::to_double(bd.ymin()), 
          CGAL::to_double(bd.ymax())));
}

Color Coloring::color(const Geometry::Point& point) const {
  // Find the query point closest to the supplied point.
  QueryPoint q = queryPoints->closest(point);
  Color c = q.color();
  assert(c != INVALID_COLOR);
  if (q == point) return c;
#ifndef GRID_EDGE_INDEX
  // Trace the segment from the query point to the test point,
  // inverting the color for every edge crossed.
  Tracer trace(q, point, tri);
  while (trace != Tracer()) {
    ++trace;
    c = opposite(c);
  }
  return c;
#else
  // We must collect a list of crossing edges to make sure they are
  // not counted twice.
  std::set<IntEdgeHandle> counted;
  typedef InteriorEdgeIndex::ConstLineCellIterator Iterator; 
  Iterator it(*intEdgeIndex, q, point, Iterator::SEGMENT()), end;
  typedef InteriorEdgeIndex::Cell Cell;
  while (it != end) {
    const Cell& cell = *it;
    const Cell::ItemList& edges = cell.getItemList();
    typedef Cell::ItemList::const_iterator EdgeIterator;
    for (EdgeIterator jt = edges.begin(); jt != edges.end(); jt++) {
      const IntEdgeHandle e = *jt;
      if ((counted.find(e) == counted.end()) && e->crosses(q, point)) {
  counted.insert(e);
  c = opposite(c);
      }
    }
    ++it;
  }
  return c;
#endif
}

bool Coloring::solid(const Geometry::Polygon& poly, Arak::Color color) const {
#ifdef GRID_EDGE_INDEX
  // First check that a point in the closure of poly is colored
  // correctly.  
  if (this->color(poly[0]) != color) return false;
  // Check that the segments of the boundary do not cross a color
  // discontinuity.
  int size = poly.size();
  for (int i = 0; i < size; i++)
    if (!compatible(poly[i], poly[(i + 1) % size])) return false;
  // Finally, check that no vertex of the coloring falls in the
  // polygon.
  typedef InteriorEdgeIndex::ConstConvexPolygonCellIterator CellIterator;
  typedef InteriorEdgeIndex::Cell Cell;
  typedef Cell::ItemList::const_iterator EdgeIterator;
  CellIterator end;
  for (CellIterator it(*intEdgeIndex, poly); it != end; ++it) {
    const Cell& cell = *it;
    const Cell::ItemList& edges = cell.getItemList();
    for (EdgeIterator jt = edges.begin(); jt != edges.end(); jt++) {
      const IntEdgeHandle e = *jt;
      const VertexHandle u = e->getPrevVertex();
      const VertexHandle v = e->getNextVertex();
      if (!cell.has_on_unbounded_side(u->point()) && // pre-filter
    !poly.has_on_unbounded_side(u->point()))
  return false;
      if (!cell.has_on_unbounded_side(v->point()) && // pre-filter
    !poly.has_on_unbounded_side(v->point()))
  return false;
    }
  }
  // The closure of the polygon is colored correctly.
  return true;
#endif
#ifdef TRIANGULATION_EDGE_INDEX
  // Traverse all triangles of the triangulation that intersect the
  // polygon and ensure that they are all colored appropriately.
  assert(false); // TODO
#endif 
}

void Coloring::addListener(Listener& listener) const {
  // Make sure not to add the same listener twice.
  if (std::find(listeners.begin(), listeners.end(), &listener) == 
      listeners.end())
    listeners.push_front(&listener);
}

void Coloring::recolored(const Geometry::Point& a,
       const Geometry::Point& b,
       const Geometry::Point& c) {
  queryPoints->recolor(a, b, c);
  typedef std::list<Listener*>::iterator ListenerIterator;
  for (ListenerIterator it = listeners.begin(); it != listeners.end(); it++) {
    Listener *listener = *it;
    listener->recolored(a, b, c);
  }
}

void Coloring::recolored(const Geometry::Point& a,
       const Geometry::Point& b,
       const Geometry::Point& c,
       const Geometry::Point& d) {
  queryPoints->recolor(a, b, c, d);
  typedef std::list<Listener*>::iterator ListenerIterator;
  for (ListenerIterator it = listeners.begin(); it != listeners.end(); it++) {
    Listener *listener = *it;
    listener->recolored(a, b, c, d);
  }
}

void Coloring::getPointWithColor(Geometry::Point& p, Color& c) const {
  const QueryPoint& q = queryPoints->point(0);
  p = q;
  c = q.color();
}

void Coloring::visualize(CGAL::Qt_widget& widget,
       bool drawEdges,
       bool drawVertices,
       bool drawRegions,
       bool drawQueries) const {
  using namespace CGAL;
  if (drawRegions) {
    // Compute a triangulation representation of the coloring.
    typedef TriangulationTraits<FaceColorInfo>::Triangulation Triangulation;
    Triangulation tri;
    this->toTriangulation(tri);
    // Draw the faces.
    widget << LineWidth(0);
    Triangulation::Finite_faces_iterator it = tri.finite_faces_begin();
    Triangulation::Finite_faces_iterator end = tri.finite_faces_end();
    while (it != end) {
      Geometry::Triangle triangle = tri.triangle(it);
      bool white = (it->info().color == Arak::WHITE);
      // widget << CGAL::GRAY; // Uncomment to show triangulation edges
      widget << (white ? CGAL::WHITE : CGAL::BLACK);
      widget << CGAL::FillColor(white ? CGAL::WHITE : CGAL::BLACK);
      widget << LineWidth(1) << triangle;
      it++;
    }
    // Draw the window.
    widget << CGAL::BLACK << noFill << LineWidth(1) << bd;
  }
  if (drawQueries) {
    // Draw the query points.
    widget << PointSize(4) << PointStyle(DISC);
    const QueryPointIndex& queryPoints = getQueryPoints();
    for (int i = 0; i < queryPoints.size(); i++) {
      if (queryPoints.point(i).color() == BLACK) 
  widget << CGAL::BLACK << queryPoints.point(i);
      else
  widget << CGAL::BLUE << queryPoints.point(i);
    }
  }
  if (drawEdges) {
    /*
     * Draw the edges.
     */
    // Draw the interior edges.
    widget << BLUE << LineWidth(2);
    for (int i = 0; i < numInteriorEdges(); i++) {
      const Coloring::IntEdgeHandle e = getInteriorEdge(i);
      widget << e->segment();
    }
    // Draw the boundary edges.
    widget << RED;
    for (int i = 0; i < numBoundaryEdges(); i++) {
      const Coloring::BdEdgeHandle e = getBoundaryEdge(i);
      widget << e->segment();
    }
  }
  if (drawVertices) {
    // Draw the interior vertices.
    widget << PointSize(6) << PointStyle(DISC) << BLUE;
    for (int i = 0; i < numVertices(Coloring::Vertex::INTERIOR); i++) {
      const Coloring::VertexHandle v = 
  getVertex(Coloring::Vertex::INTERIOR, i);
      widget << v->point();
    }
    // Draw the boundary vertices.
    widget << RED;
    for (int i = 0; i < numVertices(Coloring::Vertex::BOUNDARY); i++) {
      const Coloring::VertexHandle v = 
  getVertex(Coloring::Vertex::BOUNDARY, i);
      widget << v->point();
    }
  }
}


#define COLORING_HEADER "\
#######################################################################\n\
#                                                                      \n\
# This file represents a binary polygonal coloring of a rectangular    \n\
# subset of the plane.  It is described by its boundary, the           \n\
# vertices and edges of the coloring, and the color of a point in the  \n\
# interior of the coloring.                                            \n\
#                                                                      \n\
# The point coordinates are represented to 40 decimal places to        \n\
# prevent truncation and rounding from causing vertices that are close \n\
# to appear at the same location (and related problems).               \n\
#                                                                      \n\
#######################################################################\n"

void Coloring::write(std::ostream& out) const {
  // Write the header.
  out << COLORING_HEADER << std::endl;
  // Write the boundary of the coloring.
  out << "# The coordinates of the boundary rectangle: xmin ymin xmax ymax" 
      << std::endl << bd << std::endl;
  // Write out the vertices.
  out << std::endl << "# Number of boundary vertices:" << std::endl;
  out << numVertices(Vertex::BOUNDARY) << std::endl;
  out << "# Boundary vertices [index x y]:" << std::endl;
  out << std::scientific << std::setprecision(40);
  for (int i = 0; i < numVertices(Vertex::BOUNDARY); i++) {
    VertexHandle v = getVertex(Vertex::BOUNDARY, i);
    out << v->index << " " 
  << CGAL::to_double(v->point().x()) << " " 
  << CGAL::to_double(v->point().y()) << std::endl;
  }
  out << std::endl << "# Number of interior vertices:" << std::endl;
  out << numVertices(Vertex::INTERIOR) << std::endl;
  out << "# Interior vertices [index x y]:" << std::endl;
  for (int i = 0; i < numVertices(Vertex::INTERIOR); i++) {
    VertexHandle v = getVertex(Vertex::INTERIOR, i);
    out << v->index << " " 
  << CGAL::to_double(v->point().x()) << " " 
  << CGAL::to_double(v->point().y()) << std::endl;
  }
  // Write out the boundary edges.
  out << std::endl << "# Number of boundary edges:" << std::endl;
  out << numBoundaryEdges() << std::endl;
  out << "# Boundary edges; each edge is specified by:" << std::endl;
  out << "# a. index " << std::endl;
  out << "# b. previous corner index (0-3, counter-clockwise from lower-left)"
      << std::endl;
  out << "# c. next corner index (0-3, counter-clockwise from lower-left)" 
      << std::endl;
  out << "# d. previous vertex type (1:boundary, 2:corner)" 
      << std::endl;
  out << "# e. previous vertex index" << std::endl;
  out << "# f. next vertex type (1:boundary, 2:corner)" 
      << std::endl;
  out << "# g. next vertex index" << std::endl;
  for (int i = 0; i < numBoundaryEdges(); i++) {
    BdEdgeHandle e = getBoundaryEdge(i);
    out << e->index << " " 
  << e->getPrevCornerVertex()->index << " "
  << e->getNextCornerVertex()->index << " "
  << e->getPrevVertex()->type() << " "
  << e->getPrevVertex()->index << " " 
  << e->getNextVertex()->type() << " "
  << e->getNextVertex()->index << std::endl;
  }
  out << std::endl << "# Number of interior edges:" << std::endl;
  out << numInteriorEdges() << std::endl;
  out << "# Interior edges; each edge is specified by:" << std::endl;
  out << "# a. index " << std::endl;
  out << "# b. previous vertex type (0:interior, 1:boundary)" 
      << std::endl;
  out << "# c. previous vertex index" << std::endl;
  out << "# d. next vertex type (0:interior, 1:boundary)" 
      << std::endl;
  out << "# e. next vertex index" << std::endl;
  for (int i = 0; i < numInteriorEdges(); i++) {
    IntEdgeHandle e = getInteriorEdge(i);
    out << e->index << " " 
  << e->getPrevVertex()->type() << " "
  << e->getPrevVertex()->index << " " 
  << e->getNextVertex()->type() << " "
  << e->getNextVertex()->index << std::endl;
  }
  // Write the colored point.
  out << std::endl << "# Location of query point:" << std::endl;
  Geometry::Point point;
  Color color;
  getPointWithColor(point, color);
  out << point << std::endl;
  out << std::endl << "# Color of query point (0:white, 1:black):" 
      << std::endl;
  out << color << std::endl;
  
  out << std::endl << "# End of coloring" << std::endl;
}

void Coloring::writeBinary(std::ostream& out) const {
  double d;
  unsigned short int n;
  unsigned char c;
  // Write the boundary of the coloring.
  d = CGAL::to_double(bd.xmin());
  out.write(reinterpret_cast<char*>(&d), sizeof(d));
  d = CGAL::to_double(bd.ymin());
  out.write(reinterpret_cast<char*>(&d), sizeof(d));
  d = CGAL::to_double(bd.xmax());
  out.write(reinterpret_cast<char*>(&d), sizeof(d));
  d = CGAL::to_double(bd.ymax());
  out.write(reinterpret_cast<char*>(&d), sizeof(d));
  // Write out the vertices.
  n = numVertices(Vertex::BOUNDARY);
  out.write(reinterpret_cast<char*>(&n), sizeof(n));
  for (int i = 0; i < numVertices(Vertex::BOUNDARY); i++) {
    VertexHandle v = getVertex(Vertex::BOUNDARY, i);
    /*
    n = v->index;
    out.write(reinterpret_cast<char*>(&n), sizeof(n));
    */
    d = CGAL::to_double(v->point().x());
    out.write(reinterpret_cast<char*>(&d), sizeof(d));
    d = CGAL::to_double(v->point().y());
    out.write(reinterpret_cast<char*>(&d), sizeof(d));
  }
  n = numVertices(Vertex::INTERIOR);
  out.write(reinterpret_cast<char*>(&n), sizeof(n));
  for (int i = 0; i < numVertices(Vertex::INTERIOR); i++) {
    VertexHandle v = getVertex(Vertex::INTERIOR, i);
    /*
    n = v->index;
    out.write(reinterpret_cast<char*>(&n), sizeof(n));
    */
    d = CGAL::to_double(v->point().x());
    out.write(reinterpret_cast<char*>(&d), sizeof(d));
    d = CGAL::to_double(v->point().y());
    out.write(reinterpret_cast<char*>(&d), sizeof(d));
  }
  // Write out the boundary edges.
  n = numBoundaryEdges();
  out.write(reinterpret_cast<char*>(&n), sizeof(n));
  for (int i = 0; i < numBoundaryEdges(); i++) {
    BdEdgeHandle e = getBoundaryEdge(i);
    /*
    n = e->index;
    out.write(reinterpret_cast<char*>(&n), sizeof(n));
    */
    n = e->getPrevCornerVertex()->index;
    out.write(reinterpret_cast<char*>(&n), sizeof(n));
    n = e->getNextCornerVertex()->index;
    out.write(reinterpret_cast<char*>(&n), sizeof(n));
    c = e->getPrevVertex()->type();
    out.write(reinterpret_cast<char*>(&c), sizeof(c));
    n = e->getPrevVertex()->index;
    out.write(reinterpret_cast<char*>(&n), sizeof(n));
    c = e->getNextVertex()->type();
    out.write(reinterpret_cast<char*>(&c), sizeof(c));
    n = e->getNextVertex()->index;
    out.write(reinterpret_cast<char*>(&n), sizeof(n));
  }
  // Write out the interior edges.
  n = numInteriorEdges();
  out.write(reinterpret_cast<char*>(&n), sizeof(n));
  for (int i = 0; i < numInteriorEdges(); i++) {
    IntEdgeHandle e = getInteriorEdge(i);
    /*
    n = e->index;
    out.write(reinterpret_cast<char*>(&n), sizeof(n));
    */
    c = e->getPrevVertex()->type();
    out.write(reinterpret_cast<char*>(&c), sizeof(c));
    n = e->getPrevVertex()->index;
    out.write(reinterpret_cast<char*>(&n), sizeof(n));
    c = e->getNextVertex()->type();
    out.write(reinterpret_cast<char*>(&c), sizeof(c));
    n = e->getNextVertex()->index;
    out.write(reinterpret_cast<char*>(&n), sizeof(n));
  }
  // Write the colored point.
  Geometry::Point point;
  Color color;
  getPointWithColor(point, color);
  d = CGAL::to_double(point.x());
  out.write(reinterpret_cast<char*>(&d), sizeof(d));
  d = CGAL::to_double(point.y());
  out.write(reinterpret_cast<char*>(&d), sizeof(d));
  c = color;
  out.write(reinterpret_cast<char*>(&c), sizeof(c));
}

void Coloring::writeXfig(std::ostream& out) const {
  int xfigUnitsPerInch = 1200;
  double nativeUnitsPerInch = CGAL::to_double(bd.xmax() - bd.xmin()) / 7.5;
  double xfigUnitsPerNativeUnit = 
    double(xfigUnitsPerInch) / nativeUnitsPerInch;

  out << "#FIG 3.2" << std::endl;
  out << "Portrait" << std::endl;
  out << "Center" << std::endl;
  out << "Inches" << std::endl;
  out << "Letter" << std::endl;
  out << "100.00" << std::endl;
  out << "Single" << std::endl;
  out << "-2" << std::endl;
  out << xfigUnitsPerInch << " 2" << std::endl;

  // Write the window as a box.
  out << "2 " // object code (polyline)
      << "2 " // sub-type (box)
      << "0 " // line style
      << "1 " // thickness (1/80 inch)
      << "0 " // pen color (black)
      << "7 " // fill color (none)
      << "50 " // depth 
      << "-1 " // pen style (not used)
      << "-1 " // area fill (-1 = no fill)
      << "0.000 " // style (float)
      << "0 " // join style
      << "0 " // cap style
      << "-1 " // radius (of arc boxes)
      << "0 " // forwards arrow (0:off, 1:on)
      << "0 " // backwards arrow (0:off, 1:on)
      << "5 " // number of points
      << std::endl;
  // Each point is written as "x y ", and the first point is repeated
  // to close the box.  (The Y value is negated because FIG files place
  // the origin at the upper left corner.)
  out << "     ";
  for (int i = 0; i < 5; i++) 
    out << int(CGAL::to_double(bd.vertex(i).x()) * xfigUnitsPerNativeUnit)
  << " "
  << -int(CGAL::to_double(bd.vertex(i).y()) * xfigUnitsPerNativeUnit)
  << " ";
  out << std::endl;

  // Compute a triangulation representation of the coloring.
  typedef TriangulationTraits<FaceColorInfo>::Triangulation Triangulation;
  Triangulation tri;
  this->toTriangulation(tri);
  // Draw the faces.
  Triangulation::Finite_faces_iterator it = tri.finite_faces_begin();
  Triangulation::Finite_faces_iterator end = tri.finite_faces_end();
  while (it != end) {
    Geometry::Triangle triangle = tri.triangle(it);
    bool white = (it->info().color == Arak::WHITE);
    if (!white) {
      // Write the triangle.
      out << "2 " // object code (polyline)
    << "1 " // sub-type (polyline)
    << "0 " // line style
    << "1 " // thickness (1/80 inch)
    << "0 " // pen color (black)
    << "-1 " // fill color (black)
    << "50 " // depth 
    << "-1 " // pen style (not used)
    << "20 " // area fill (-1 = no fill)
    << "0.000 " // style (float)
    << "0 " // join style
    << "0 " // cap style
    << "-1 " // radius (of arc boxes)
    << "0 " // forwards arrow (0:off, 1:on)
    << "0 " // backwards arrow (0:off, 1:on)
    << "4 " // number of points
    << std::endl;
      out << "     ";
      for (int i = 0; i < 4; i++) 
    out << int(CGAL::to_double(triangle.vertex(i).x()) * 
         xfigUnitsPerNativeUnit)
  << " "
  << -int(CGAL::to_double(triangle.vertex(i).y()) * 
    xfigUnitsPerNativeUnit)
  << " ";
      out << std::endl;
    }
    it++;
  }
  out << std::endl;
}

bool Coloring::read(std::istream& in) {
  int rows = intEdgeIndex->numRows();
  int cols = intEdgeIndex->numCols();
  Geometry::Rectangle bd;
  in >> skipcomments;
  in >> bd;
  if (!in.good()) return false;
  initialize(bd, rows, cols);
  // Remove the boundary edges.
  while (numBoundaryEdges() > 0) {
    BdEdgeHandle e = getBoundaryEdge(0);
    freeEdge(e);
  }
  // Read the number of boundary vertices.
  int n;
  Geometry::Point p;
  in >> skipcomments;
  in >> n;
  in >> skipcomments;
  int index;
  if (!in.good()) return false;
  // Read the boundary vertices.
  for (int i = 0; i < n; i++) {
    in >> index;
    assert(index == i);
    in >> p;
    if (!in.good()) return false;
    assert(newVertex(Vertex::BOUNDARY, p)->index == i);
  }
  // Read the number of interior vertices.
  in >> skipcomments;
  in >> n;
  in >> skipcomments;
  if (!in.good()) return false;
  // Read the interior vertices.
  for (int i = 0; i < n; i++) {
    in >> index;
    assert(index == i);
    in >> p;
    if (!in.good()) return false;
    assert(newVertex(Vertex::INTERIOR, p)->index == i);
  }
  // Read the number of boundary edges.
  in >> skipcomments;
  in >> n;
  in >> skipcomments;
  if (!in.good()) return false;
  int prevTypeInt, nextTypeInt;
  Vertex::Type prevType, nextType;
  int prevIndex, nextIndex;
  int prevCornerIndex, nextCornerIndex;
  // Read the boundary edges.
  for (int i = 0; i < n; i++) {
    in >> index;
    assert(index == i);
    in >> prevCornerIndex 
       >> nextCornerIndex 
       >> prevTypeInt
       >> prevIndex 
       >> nextTypeInt
       >> nextIndex;
    if (!in.good()) return false;
    switch (prevTypeInt) {
    case 1: prevType = Vertex::BOUNDARY; break;
    case 2: prevType = Vertex::CORNER; break;
    default: return false;
    };
    switch (nextTypeInt) {
    case 1: nextType = Vertex::BOUNDARY; break;
    case 2: nextType = Vertex::CORNER; break;
    default: return false;
    };
    newBoundaryEdge(getVertex(prevType, prevIndex),
        getVertex(nextType, nextIndex),
        getVertex(Vertex::CORNER, prevCornerIndex),
        getVertex(Vertex::CORNER, nextCornerIndex));
  }
  // Read the number of interior edges.
  in >> skipcomments;
  in >> n;
  in >> skipcomments;
  if (!in.good()) return false;
  // Read the interior edges.
  for (int i = 0; i < n; i++) {
    in >> index;
    assert(index == i);
    in >> prevTypeInt >> prevIndex >> nextTypeInt >> nextIndex;
    if (!in.good()) return false;
    switch (prevTypeInt) {
    case 0: prevType = Vertex::INTERIOR; break;
    case 1: prevType = Vertex::BOUNDARY; break;
    default: return false;
    };
    switch (nextTypeInt) {
    case 0: nextType = Vertex::INTERIOR; break;
    case 1: nextType = Vertex::BOUNDARY; break;
    default: return false;
    };
    newInteriorEdge(getVertex(prevType, prevIndex),
        getVertex(nextType, nextIndex));
  }
  // Read the query point.
  int colorInt;
  Color c;
  in >> skipcomments >> p >> skipcomments >> colorInt;
  if (!in.good()) return false;
  switch (colorInt) {
  case 0: c = WHITE; break;
  case 1: c = BLACK; break;
  default: return false;
  }
  // If the current color of this point is wrong, invert the (only)
  // query point's color.
  assert(queryPoints->size() == 1);
  QueryPoint& q = queryPoints->point(0);
  if (color(p) != c) q.setColor(opposite(q.color()));
  // Check that the coloring is valid.
  test();
  return true;
}

bool Coloring::readBinary(std::istream& in) {
  // Read the boundary.
  double xmin, ymin, xmax, ymax;
  in.read(reinterpret_cast<char*>(&xmin), sizeof(xmin));
  in.read(reinterpret_cast<char*>(&ymin), sizeof(ymin));
  in.read(reinterpret_cast<char*>(&xmax), sizeof(xmax));
  in.read(reinterpret_cast<char*>(&ymax), sizeof(ymax));
  if (!in.good()) return false;
  bd = Geometry::Rectangle(Geometry::Kernel::FT(xmin),
         Geometry::Kernel::FT(ymin),
         Geometry::Kernel::FT(xmax),
         Geometry::Kernel::FT(ymax));
  int rows = intEdgeIndex->numRows();
  int cols = intEdgeIndex->numCols();
  initialize(bd, rows, cols);
  // Remove the boundary edges.
  while (numBoundaryEdges() > 0) {
    BdEdgeHandle e = getBoundaryEdge(0);
    freeEdge(e);
  }
  // Read the number of boundary vertices.
  unsigned short int n;
  in.read(reinterpret_cast<char*>(&n), sizeof(n));
  double x, y;
  if (!in.good()) return false;
  // Read the boundary vertices.
  for (int i = 0; i < n; i++) {
    /*
    in.read(reinterpret_cast<char*>(&index), sizeof(index));
    assert(index == i);
    */
    in.read(reinterpret_cast<char*>(&x), sizeof(x));
    in.read(reinterpret_cast<char*>(&y), sizeof(y));
    if (!in.good()) return false;
    Geometry::Point p(x, y);
    assert(newVertex(Vertex::BOUNDARY, p)->index == i);
  }
  // Read the number of interior vertices.
  in.read(reinterpret_cast<char*>(&n), sizeof(n));
  if (!in.good()) return false;
  // Read the interior vertices.
  for (int i = 0; i < n; i++) {
    /*
    in.read(reinterpret_cast<char*>(&index), sizeof(index));
    assert(index == i);
    */
    in.read(reinterpret_cast<char*>(&x), sizeof(x));
    in.read(reinterpret_cast<char*>(&y), sizeof(y));
    if (!in.good()) return false;
    Geometry::Point p(x, y);
    assert(newVertex(Vertex::INTERIOR, p)->index == i);
  }
  // Read the number of boundary edges.
  in.read(reinterpret_cast<char*>(&n), sizeof(n));
  if (!in.good()) return false;
  unsigned char prevTypeInt, nextTypeInt;
  Vertex::Type prevType, nextType;
  unsigned short int prevIndex, nextIndex;
  unsigned short int prevCornerIndex, nextCornerIndex;
  // Read the boundary edges.
  for (int i = 0; i < n; i++) {
    /*
    in.read(reinterpret_cast<char*>(&index), sizeof(index));
    assert(index == i);
    */
    in.read(reinterpret_cast<char*>(&prevCornerIndex), 
      sizeof(prevCornerIndex));
    in.read(reinterpret_cast<char*>(&nextCornerIndex), 
      sizeof(nextCornerIndex));
    in.read(reinterpret_cast<char*>(&prevTypeInt), sizeof(prevTypeInt));
    in.read(reinterpret_cast<char*>(&prevIndex), sizeof(prevIndex));
    in.read(reinterpret_cast<char*>(&nextTypeInt), sizeof(nextTypeInt));
    in.read(reinterpret_cast<char*>(&nextIndex), sizeof(nextIndex));
    if (!in.good()) return false;
    switch (prevTypeInt) {
    case 1: prevType = Vertex::BOUNDARY; break;
    case 2: prevType = Vertex::CORNER; break;
    default: return false;
    };
    switch (nextTypeInt) {
    case 1: nextType = Vertex::BOUNDARY; break;
    case 2: nextType = Vertex::CORNER; break;
    default: return false;
    };
    newBoundaryEdge(getVertex(prevType, prevIndex),
        getVertex(nextType, nextIndex),
        getVertex(Vertex::CORNER, prevCornerIndex),
        getVertex(Vertex::CORNER, nextCornerIndex));
  }
  // Read the number of interior edges.
  in.read(reinterpret_cast<char*>(&n), sizeof(n));
  if (!in.good()) return false;
  // Read the interior edges.
  for (int i = 0; i < n; i++) {
    /*
    in.read(reinterpret_cast<char*>(&index), sizeof(index));
    assert(index == i);
    */
    in.read(reinterpret_cast<char*>(&prevTypeInt), sizeof(prevTypeInt));
    in.read(reinterpret_cast<char*>(&prevIndex), sizeof(prevIndex));
    in.read(reinterpret_cast<char*>(&nextTypeInt), sizeof(nextTypeInt));
    in.read(reinterpret_cast<char*>(&nextIndex), sizeof(nextIndex));
    if (!in.good()) return false;
    switch (prevTypeInt) {
    case 0: prevType = Vertex::INTERIOR; break;
    case 1: prevType = Vertex::BOUNDARY; break;
    default: return false;
    };
    switch (nextTypeInt) {
    case 0: nextType = Vertex::INTERIOR; break;
    case 1: nextType = Vertex::BOUNDARY; break;
    default: return false;
    };
    newInteriorEdge(getVertex(prevType, prevIndex),
        getVertex(nextType, nextIndex));
  }
  // Read the query point.
  unsigned char colorInt;
  in.read(reinterpret_cast<char*>(&x), sizeof(x));
  in.read(reinterpret_cast<char*>(&y), sizeof(y));
  Geometry::Point p(x, y);
  in.read(reinterpret_cast<char*>(&colorInt), sizeof(colorInt));
  if (!in.good()) return false;
  Color c;
  switch (colorInt) {
  case 0: c = WHITE; break;
  case 1: c = BLACK; break;
  default: return false;
  }
  // If the current color of this point is wrong, invert the (only)
  // query point's color.
  assert(queryPoints->size() == 1);
  QueryPoint& q = queryPoints->point(0);
  if (color(p) != c) q.setColor(opposite(q.color()));
  // Check that the coloring is valid.
  test();
  return true;
}

---