// medial_axis.hpp header file
// derived from
// Boost.Polygon library voronoi_diagram.hpp header file
// which is
//          Copyright Andrii Sydorchuk 2010-2012.
// Distributed under the Boost Software License, Version 1.0.
//    (See accompanying file LICENSE_1_0.txt or copy at
//          http://www.boost.org/LICENSE_1_0.txt)

// See http://www.boost.org for updates, documentation, and revision history.

// Derivative work by Michael E. Sheldrake, Copyright 2013, distributed
// under the same terms as the license above.

// This is essentially boost/polygon/voronoi_diagram.hpp adapted to further 
// process the Voronoi diagram graph structure to make it represent the 
// medial axis of a polygon (with holes) represented by segment input.

#ifndef BOOST_POLYGON_MEDIAL_AXIS
#define BOOST_POLYGON_MEDIAL_AXIS

#include <vector>
#include <utility>
#include <cstdio>
#include <string>
#include <math.h>

#include "boost/lexical_cast.hpp"
#include "boost/polygon/detail/voronoi_ctypes.hpp"
#include "boost/polygon/detail/voronoi_structures.hpp"

#include "boost/polygon/voronoi_geometry_type.hpp"

#define to_str(n) boost::lexical_cast<std::string>( n )

namespace boost {
namespace polygon {

// Forward declarations.
template <typename T>
class medial_axis_edge;
 
// Represents Voronoi cell.
// Data members:
//   1) index of the source within the initial input set
//   2) pointer to the incident edge
//   3) mutable color member
// Cell may contain point or segment site inside.
template <typename T>
class medial_axis_cell {
 public:
  typedef T coordinate_type;
  typedef std::size_t color_type;
  typedef medial_axis_edge<coordinate_type> medial_axis_edge_type;
  typedef std::size_t source_index_type;
  typedef SourceCategory source_category_type;

  medial_axis_cell(source_index_type source_index,
               source_category_type source_category) :
      source_index_(source_index),
      incident_edge_(NULL),
      color_(source_category) {}

  // Returns true if the cell contains point site, false else.
  bool contains_point() const {
    source_category_type source_category = this->source_category();
    return belongs(source_category, GEOMETRY_CATEGORY_POINT);
  }

  // Returns true if the cell contains segment site, false else.
  bool contains_segment() const {
    source_category_type source_category = this->source_category();
    return belongs(source_category, GEOMETRY_CATEGORY_SEGMENT);
  }

  source_index_type source_index() const {
    return source_index_;
  }

  source_category_type source_category() const {
    return static_cast<source_category_type>(color_ & SOURCE_CATEGORY_BITMASK);
  }

  // Degenerate cells don't have any incident edges.
  bool is_degenerate() const { return incident_edge_ == NULL; }

  medial_axis_edge_type* incident_edge() { return incident_edge_; }
  const medial_axis_edge_type* incident_edge() const { return incident_edge_; }
  void incident_edge(medial_axis_edge_type* e) { incident_edge_ = e; }

  color_type color() const { return color_ >> BITS_SHIFT; }
  void color(color_type color) const {
    color_ &= BITS_MASK;
    color_ |= color << BITS_SHIFT;
  }

 private:
  // 5 color bits are reserved.
  enum Bits {
    BITS_SHIFT = 0x5,
    BITS_MASK = 0x1F
  };

  source_index_type source_index_;
  medial_axis_edge_type* incident_edge_;
  mutable color_type color_;
};

// Represents Voronoi vertex.
// Data members:
//   1) vertex coordinates
//   2) radius of a maximal inscribed circle to the polygon at the vertex
//   3) pointer to the incident edge
//   4) mutable color member
template <typename T>
class medial_axis_vertex {
 public:
  typedef T coordinate_type;
  typedef std::size_t color_type;
  typedef medial_axis_edge<coordinate_type> medial_axis_edge_type;

  medial_axis_vertex(const coordinate_type& x, const coordinate_type& y,
                     const coordinate_type& r=0) :
      x_(x),
      y_(y),
      r_(r),
      incident_edge_(NULL),
      color_(0) {}

  const coordinate_type& x() const { return x_; }
  const coordinate_type& y() const { return y_; }
  const coordinate_type& r() const { return r_; }

  bool is_degenerate() const { return incident_edge_ == NULL; }

  medial_axis_edge_type* incident_edge() { return incident_edge_; }
  const medial_axis_edge_type* incident_edge() const { return incident_edge_; }
  void incident_edge(medial_axis_edge_type* e) { incident_edge_ = e; }

  color_type color() const { return color_ >> BITS_SHIFT; }
  void color(color_type color) const {
    color_ &= BITS_MASK;
    color_ |= color << BITS_SHIFT;
  }

 private:
  // 5 color bits are reserved.
  enum Bits {
    BITS_SHIFT = 0x5,
    BITS_MASK = 0x1F
  };

  coordinate_type x_;
  coordinate_type y_;
  coordinate_type r_;
  medial_axis_edge_type* incident_edge_;
  mutable color_type color_;
};

// Half-edge data structure. Represents Voronoi edge.
// Data members:
//   1) pointer to the corresponding cell
//   2) pointer to the vertex that is the starting
//      point of the half-edge
//   3) pointer to the twin edge
//   4) pointer to the CCW next edge
//   5) pointer to the CCW prev edge
//   6) boolean indicating whether foot coordinates have been set
//   7) mutable color member
template <typename T>
class medial_axis_edge {
 public:
  typedef T coordinate_type;
  typedef medial_axis_cell<coordinate_type> medial_axis_cell_type;
  typedef medial_axis_vertex<coordinate_type> medial_axis_vertex_type;
  typedef medial_axis_edge<coordinate_type> medial_axis_edge_type;
  typedef std::size_t color_type;

  medial_axis_edge(bool is_linear, bool is_primary) :
      cell_(NULL),
      vertex_(NULL),
      twin_(NULL),
      next_(NULL),
      prev_(NULL),
      footset_(false),
      color_(0) {
    if (is_linear)
      color_ |= BIT_IS_LINEAR;
    if (is_primary)
      color_ |= BIT_IS_PRIMARY;
  }

  medial_axis_cell_type* cell() { return cell_; }
  const medial_axis_cell_type* cell() const { return cell_; }
  void cell(medial_axis_cell_type* c) { cell_ = c; }

  medial_axis_vertex_type* vertex0() { return vertex_; }
  const medial_axis_vertex_type* vertex0() const { return vertex_; }
  void vertex0(medial_axis_vertex_type* v) { vertex_ = v; }

  medial_axis_vertex_type* vertex1() { return twin_->vertex0(); }
  const medial_axis_vertex_type* vertex1() const { return twin_->vertex0(); }

  medial_axis_edge_type* twin() { return twin_; }
  const medial_axis_edge_type* twin() const { return twin_; }
  void twin(medial_axis_edge_type* e) { twin_ = e; }

  medial_axis_edge_type* next() { return next_; }
  const medial_axis_edge_type* next() const { return next_; }
  void next(medial_axis_edge_type* e) { next_ = e; }

  medial_axis_edge_type* prev() { return prev_; }
  const medial_axis_edge_type* prev() const { return prev_; }
  void prev(medial_axis_edge_type* e) { prev_ = e; }

  // Returns a pointer to the rotation next edge
  // over the starting point of the half-edge.
  medial_axis_edge_type* rot_next() { return prev_->twin(); }
  const medial_axis_edge_type* rot_next() const { return prev_->twin(); }

  // Returns a pointer to the rotation prev edge
  // over the starting point of the half-edge.
  medial_axis_edge_type* rot_prev() { return twin_->next(); }
  const medial_axis_edge_type* rot_prev() const { return twin_->next(); }

  // Returns true if the edge is finite (segment, parabolic arc).
  // Returns false if the edge is infinite (ray, line).
  bool is_finite() const { return vertex0() && vertex1(); }

  // Returns true if the edge is infinite (ray, line).
  // Returns false if the edge is finite (segment, parabolic arc).
  bool is_infinite() const { return !vertex0() || !vertex1(); }

  // Returns true if the edge is linear (segment, ray, line).
  // Returns false if the edge is curved (parabolic arc).
  bool is_linear() const {
    return (color_ & BIT_IS_LINEAR) ? true : false;
  }

  // Returns true if the edge is curved (parabolic arc).
  // Returns false if the edge is linear (segment, ray, line).
  bool is_curved() const {
    return (color_ & BIT_IS_LINEAR) ? false : true;
  }

  // Returns false if edge goes through the endpoint of the segment.
  // Returns true else.
  bool is_primary() const {
    return (color_ & BIT_IS_PRIMARY) ? true : false;
  }

  // Returns true if edge goes through the endpoint of the segment.
  // Returns false else.
  bool is_secondary() const {
    return (color_ & BIT_IS_PRIMARY) ? false : true;
  }

  color_type color() const { return color_ >> BITS_SHIFT; }
  void color(color_type color) const {
    color_ &= BITS_MASK;
    color_ |= color << BITS_SHIFT;
  }
  
  // foot: where radius from vertex0 touches source segment at a 90 degree angle
  const detail::point_2d<default_voronoi_builder::int_type>* foot() const { 
    if (!footset_) {return NULL;}
    return &foot_;
  }
  void foot(coordinate_type x, coordinate_type y) {
    footset_ = true;
    foot_.x(x);
    foot_.y(y);
  }

 private:
  // 5 color bits are reserved.
  enum Bits {
    BIT_IS_LINEAR = 0x1,  // linear is opposite to curved
    BIT_IS_PRIMARY = 0x2,  // primary is opposite to secondary

    BITS_SHIFT = 0x5,
    BITS_MASK = 0x1F
  };

  medial_axis_cell_type* cell_;
  medial_axis_vertex_type* vertex_;
  medial_axis_edge_type* twin_;
  medial_axis_edge_type* next_;
  medial_axis_edge_type* prev_;
  mutable color_type color_;
  mutable detail::point_2d<default_voronoi_builder::int_type> foot_;
  bool footset_;
  mutable detail::point_2d<default_voronoi_builder::int_type> p1_;

};

template <typename T>
struct medial_axis_traits {
  typedef T coordinate_type;
  typedef medial_axis_cell<coordinate_type> cell_type;
  typedef medial_axis_vertex<coordinate_type> vertex_type;
  typedef medial_axis_edge<coordinate_type> edge_type;
  typedef class {
   public:
    enum { ULPS = 128 };
    bool operator()(const vertex_type& v1, const vertex_type& v2) const {
      return (ulp_cmp(v1.x(), v2.x(), ULPS) ==
              detail::ulp_comparison<T>::EQUAL) &&
             (ulp_cmp(v1.y(), v2.y(), ULPS) ==
              detail::ulp_comparison<T>::EQUAL);
    }
   private:
    typename detail::ulp_comparison<T> ulp_cmp;
  } vertex_equality_predicate_type;
};

// Voronoi output data structure.
// CCW ordering is used on the faces perimeter and around the vertices.
template <typename T, typename TRAITS = medial_axis_traits<T> >
class medial_axis {
 public:
  typedef typename TRAITS::coordinate_type coordinate_type;
  typedef typename TRAITS::cell_type cell_type;
  typedef typename TRAITS::vertex_type vertex_type;
  typedef typename TRAITS::edge_type edge_type;

  typedef std::vector<cell_type> cell_container_type;
  typedef typename cell_container_type::const_iterator const_cell_iterator;

  typedef std::vector<vertex_type> vertex_container_type;
  typedef typename vertex_container_type::const_iterator const_vertex_iterator;

  typedef std::vector<edge_type> edge_container_type;
  typedef typename edge_container_type::const_iterator const_edge_iterator;

  medial_axis() {}

  void clear() {
    cells_.clear();
    vertices_.clear();
    edges_.clear();
  }

  const cell_container_type& cells() const {
    return cells_;
  }

  const vertex_container_type& vertices() const {
    return vertices_;
  }

  const edge_container_type& edges() const {
    return edges_;
  }

  const std::string& event_log() const {
    return event_log_;
  }

  std::size_t num_cells() const {
    return cells_.size();
  }

  std::size_t num_edges() const {
    return edges_.size();
  }

  std::size_t num_vertices() const {
    return vertices_.size();
  }

  void _reserve(int num_sites) {
    cells_.reserve(num_sites);
    vertices_.reserve(num_sites << 1);
    edges_.reserve((num_sites << 2) + (num_sites << 1));
  }

  template <typename CT>
  void _process_single_site(const detail::site_event<CT>& site) {
    cells_.push_back(cell_type(site.initial_index(), site.source_category()));
  }

  // Insert a new half-edge into the output data structure.
  // Takes as input left and right sites that form a new bisector.
  // Returns a pair of pointers to a new half-edges.
  template <typename CT>
  std::pair<void*, void*> _insert_new_edge(
      const detail::site_event<CT>& site1,
      const detail::site_event<CT>& site2) {
    //printf("site event\n");
    // Get sites' indexes.
    int site_index1 = site1.sorted_index();
    int site_index2 = site2.sorted_index();

    bool is_linear = is_linear_edge(site1, site2);
    bool is_primary = is_primary_edge(site1, site2);

    // Create a new half-edge that belongs to the first site.
    edges_.push_back(edge_type(is_linear, is_primary));
    edge_type& edge1 = edges_.back();
 
    // Create a new half-edge that belongs to the second site.
    edges_.push_back(edge_type(is_linear, is_primary));
    edge_type& edge2 = edges_.back();

    // Add the initial cell during the first edge insertion.
    if (cells_.empty()) {
      cells_.push_back(cell_type(
          site1.initial_index(), site1.source_category()));
    }

    // The second site represents a new site during site event
    // processing. Add a new cell to the cell records.
    cells_.push_back(cell_type(
        site2.initial_index(), site2.source_category()));

    // Set up pointers to cells.
    edge1.cell(&cells_[site_index1]);
    edge2.cell(&cells_[site_index2]);

    // Set up twin pointers.
    edge1.twin(&edge2);
    edge2.twin(&edge1);
    
// if the edge is curved we can display a parabola between the point and
// segment that define it, and we know the start and end feet for the 
// half edge on the point side - just the point. We might be able to 
// figure the feet on the segment too, but not sure.

// if the edge is straight, feet at start and end are probably just
// site vertices - but which? first or second? and what about seeing
// a vertex multiple times for the three events from a segment?
// Maybe only figure foot when the segment event happens. What about the
// inverted segment? maybe only deal with foot when handling the segment the
// second time around.

    event_log_ += "<g id=\"sites"+to_str((UV) &site1)+"_"+to_str((UV) &site2)+"\" ";
    event_log_ += " class=\"ine1"+to_str(is_linear?"linear":"curved")+to_str(is_primary?"primary":"secondary")+"\">\n";

    bool showfeet = false;

    if (false) {
    if (belongs(site1.source_category(), GEOMETRY_CATEGORY_POINT)) {
      event_log_ += "<circle id=\"site" + to_str((UV) &site1) + "\" ";
      event_log_ += "cx=\""+to_str(site1.point0().x())+"\" cy=\""+to_str(site1.point0().y())+"\" r=\"8000\" class=\"es1p\"/>\n";
    } else if (belongs(site1.source_category(), GEOMETRY_CATEGORY_SEGMENT)) {
      event_log_ += "<line id=\"site" + to_str((UV) &site1) + "\" ";
      event_log_ += " x1=\""+to_str(site1.point0().x())+"\" y1=\""+to_str(site1.point0().y())+"\"";
      event_log_ += " x2=\""+to_str(site1.point1().x())+"\" y2=\""+to_str(site1.point1().y())+"\"/>\n";
    } else {
      event_log_ += "<!-- no site 1 -->\n";
    }
    if (belongs(site2.source_category(), GEOMETRY_CATEGORY_POINT)) {
      event_log_ += "<circle id=\"site" + to_str((UV) &site2) + "\" ";
      event_log_ += "cx=\""+to_str(site2.point0().x())+"\" cy=\""+to_str(site2.point0().y())+"\" r=\"8000\" class=\"es2p\"/>\n";
    } else if (belongs(site2.source_category(), GEOMETRY_CATEGORY_SEGMENT)) {
      event_log_ += "<line id=\"site" + to_str((UV) &site2) + "\" ";
      event_log_ += " x1=\""+to_str(site2.point0().x())+"\" y1=\""+to_str(site2.point0().y())+"\"";
      event_log_ += " x2=\""+to_str(site2.point1().x())+"\" y2=\""+to_str(site2.point1().y())+"\"/>\n";
    } else {
      event_log_ += "<!-- no site 2 -->\n";
    }
    if (!is_linear) {
      event_log_ += "<!-- curved -->\n";
    } else {
      event_log_ += "<!-- linear -->\n";
    }
    }
    
    //set foot
    // this needs work - see note for "set foot" below in the 
    // _insert_new_edge() function for circle events
    
    // Though the foot determined here for curved edges is likely correct
    // it's also likely it's getting reset later when a vertex is available.
    
    if (!is_linear
        //belongs(site1.source_category(), GEOMETRY_CATEGORY_SEGMENT)
        //|| belongs(site2.source_category(), GEOMETRY_CATEGORY_SEGMENT)
       ) {
      if (edge1.cell()->contains_point() 
          //&& edge2.cell()->contains_segment()
         ) {
        // belongs(site1.source_category(), GEOMETRY_CATEGORY_POINT)
        // && belongs(site2.source_category(), GEOMETRY_CATEGORY_SEGMENT)) {
        edge1.foot(site1.point0().x(), site1.point0().y());
        //printf("ine1 foot 1\n");
        if (showfeet) {
        event_log_ += "<circle id=\"foot" + to_str((UV) &site1) + "_f\" ";
        event_log_ += "cx=\""+to_str(site1.point0().x())+"\" cy=\""+to_str(site1.point0().y())+"\" r=\"5000\" class=\"ine1f\"/>\n";
        }
      }
      if (edge2.cell()->contains_point() 
          //&& edge1.cell()->contains_segment()
         ) {
        // belongs(site2.source_category(), GEOMETRY_CATEGORY_POINT)
        // && belongs(site1.source_category(), GEOMETRY_CATEGORY_SEGMENT)) {
        edge2.foot(site2.point0().x(), site2.point0().y());
        //printf("ine1 foot 2\n");
        if (showfeet) {
        event_log_ += "<circle id=\"foot" + to_str((UV) &site2) + "_f\" ";
        event_log_ += "cx=\""+to_str(site2.point0().x())+"\" cy=\""+to_str(site2.point0().y())+"\" r=\"5000\" class=\"ine1f\"/>\n";
        }
      }

    }
    
    event_log_ += "</g>\n";

    // Return a pointer to the new half-edge.
    return std::make_pair(&edge1, &edge2);
  }

  // Insert a new half-edge into the output data structure with the
  // start at the point where two previously added half-edges intersect.
  // Takes as input two sites that create a new bisector, circle event
  // that corresponds to the intersection point of the two old half-edges,
  // pointers to those half-edges. Half-edges' direction goes out of the
  // new Voronoi vertex point. Returns a pair of pointers to a new half-edges.
  template <typename CT1, typename CT2>
  std::pair<void*, void*> _insert_new_edge(
      const detail::site_event<CT1>& site1,
      const detail::site_event<CT1>& site3,
      const detail::circle_event<CT2>& circle,
      void* data12, void* data23) {
    edge_type* edge12 = static_cast<edge_type*>(data12);
    edge_type* edge23 = static_cast<edge_type*>(data23);
    //printf("circle event\n");

    // Add a new Voronoi vertex.
    vertices_.push_back(vertex_type(circle.x(), circle.y(), 
                                    circle.lower_x() - circle.x()));
    vertex_type& new_vertex = vertices_.back();

    event_log_ += "<g id=\"sites"+to_str((UV) &site1)+"_"+to_str((UV) &site3)+"\" ";
    event_log_ += " class=\"ine2"+to_str(is_linear_edge(site1, site3)?"linear":"curved")+to_str(is_primary_edge(site1, site3)?"primary":"secondary")+"\">\n";
    
    if (false) {
    event_log_ += "<!-- edge12 is "+to_str(edge12->is_curved()?"curved":"linear")+" "+to_str(edge12->is_primary()?"primary":"secondary")+" -->\n";
    event_log_ += "<!-- edge23 is "+to_str(edge23->is_curved()?"curved":"linear")+" "+to_str(edge23->is_primary()?"primary":"secondary")+" -->\n";
    event_log_ += "<circle id=\"sites" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"_pt\" ";
    event_log_ += "cx=\""+to_str(new_vertex.x())+"\" cy=\""+to_str(new_vertex.y())+"\" r=\"8000\" class=\"cirevtpt\"/>\n";
    event_log_ += "<circle id=\"sites" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"_cir\" ";
    event_log_ += "cx=\""+to_str(new_vertex.x())+"\" cy=\""+to_str(new_vertex.y())+"\" r=\""+to_str(new_vertex.r())+"\" class=\"cirevtcir\"/>\n";
    }
    
    //printf("whatverts? %d %d %d %d\n",(edge12->vertex0()?1:0),(edge12->vertex1()?1:0),(edge23->vertex0()?1:0),(edge23->vertex1()?1:0));

    // Update vertex pointers of the old edges.
    edge12->vertex0(&new_vertex);
    edge23->vertex0(&new_vertex);
    
// It appears that the "old" edges never have vertex0 set before the above.
// But often one or both of them already have vertex1 set.

if (false) {
    if (edge12->vertex0() && edge12->vertex1()) {
      event_log_ += "<line id=\"e12_" + to_str((UV) &site1) + "_" + to_str((UV) &site3) + "\" class=\"edge12\" ";
      event_log_ += " x1=\""+to_str(edge12->vertex0()->x())+"\" y1=\""+to_str(edge12->vertex0()->y())+"\"";
      event_log_ += " x2=\""+to_str(edge12->vertex1()->x())+"\" y2=\""+to_str(edge12->vertex1()->y())+"\"/>\n";
      event_log_ += "<!-- "+to_str(edge12->is_curved()?"curved":"linear")+" -->\n";
      //printf( "WE ARE HERE1\n");
    }
    if (edge23->vertex0() && edge23->vertex1()) {
      event_log_ += "<line id=\"e12_" + to_str((UV) &site1) + "_" + to_str((UV) &site3) + "\" class=\"edge23\" ";
      event_log_ += " x1=\""+to_str(edge23->vertex0()->x())+"\" y1=\""+to_str(edge23->vertex0()->y())+"\"";
      event_log_ += " x2=\""+to_str(edge23->vertex1()->x())+"\" y2=\""+to_str(edge23->vertex1()->y())+"\"/>\n";
      event_log_ += "<!-- "+to_str(edge23->is_curved()?"curved":"linear")+" -->\n";
      //printf( "WE ARE HERE2\n");
    }
    if (belongs(site1.source_category(), GEOMETRY_CATEGORY_SEGMENT)) {
      event_log_ += "<line id=\"siteseg1_" + to_str((UV) &site1) + "_" + to_str((UV) &site3) + "\" class=\"siteseg1\" ";
      event_log_ += " x1=\""+to_str(site1.point0().x())+"\" y1=\""+to_str(site1.point0().y())+"\"";
      event_log_ += " x2=\""+to_str(site1.point1().x())+"\" y2=\""+to_str(site1.point1().y())+"\"/>\n";
    } else if (belongs(site1.source_category(), GEOMETRY_CATEGORY_POINT)) {
      event_log_ += "<circle id=\"sitepoint1_" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"\" ";
      event_log_ += "cx=\""+to_str(site1.point0().x())+"\" cy=\""+to_str(site1.point0().y())+"\" r=\"40000\" class=\"sitept1\"/>\n";
    }
    if (belongs(site3.source_category(), GEOMETRY_CATEGORY_SEGMENT)) {
      event_log_ += "<line id=\"siteseg3_" + to_str((UV) &site1) + "_" + to_str((UV) &site3) + "\" class=\"siteseg3\" ";
      event_log_ += " x1=\""+to_str(site3.point0().x())+"\" y1=\""+to_str(site3.point0().y())+"\"";
      event_log_ += " x2=\""+to_str(site3.point1().x())+"\" y2=\""+to_str(site3.point1().y())+"\"/>\n";
    } else if (belongs(site3.source_category(), GEOMETRY_CATEGORY_POINT)) {
      event_log_ += "<circle id=\"sitepoint3_" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"\" ";
      event_log_ += "cx=\""+to_str(site3.point0().x())+"\" cy=\""+to_str(site3.point0().y())+"\" r=\"40000\" class=\"sitept3\"/>\n";
    }
}
    bool is_linear = is_linear_edge(site1, site3);
    bool is_primary = is_primary_edge(site1, site3);

    // Add a new half-edge.
    edges_.push_back(edge_type(is_linear, is_primary));
    edge_type& new_edge1 = edges_.back();
    new_edge1.cell(&cells_[site1.sorted_index()]);

    // Add a new half-edge.
    edges_.push_back(edge_type(is_linear, is_primary));
    edge_type& new_edge2 = edges_.back();
    new_edge2.cell(&cells_[site3.sorted_index()]);

    // Update twin pointers.
    new_edge1.twin(&new_edge2);
    new_edge2.twin(&new_edge1);

    // Update vertex pointer.
    new_edge2.vertex0(&new_vertex);

    // Update Voronoi prev/next pointers.
    edge12->prev(&new_edge1);
    new_edge1.next(edge12);
    edge12->twin()->next(edge23);
    edge23->prev(edge12->twin());
    edge23->twin()->next(&new_edge2);
    new_edge2.prev(edge23->twin());

    //set foot
    // It's possible that we can do all foot-finding in this event processing
    // (here in the circle event, and in the other site even code above).
    // But we haven't completely understood or diagramed exactly what edges and 
    // vertices are available during these events. With a rough mental sketch
    // of whats going on, we've been able to quickly work up this code to 
    // calculate enough feet, and infer others later, to handle most cases, and
    // demonstrate that this medial axis refinement of the Voronoi diagram
    // should work.
    // 
    // What's needed next is to properly analyze-diagram-understand what's
    // happening during site/circle events, so the code here can be extended
    // a bit to really cover all cases of calculating the foot, or as many
    // cases as possible.
    

    bool showfeet = false;


    // for edges going into corners

    // edge12 into corner

    // note the cast from float to int for these ==s : LHS is an int type, RHS is float
    // nope: that didn't fix it, and that's what we want to avoid anyway
    if (edge12->vertex1()
        //&& (  (site1.point1().x() == (coordinate_type) edge12->vertex1()->x() && site1.point1().y() == (coordinate_type) edge12->vertex1()->y())
        //   || (site1.point0().x() == (coordinate_type) edge12->vertex1()->x() && site1.point0().y() == (coordinate_type) edge12->vertex1()->y())
        //   )
// If this really works, do same/similar for next if for edge23
// ... catches more cases than the target case, putting some feet not on polygon segments
//     but see if something similarly topological can do the job here
        && edge12->is_primary() 
        && edge12->next() && !edge12->next()->is_primary()
        && (  edge12->next()->cell()->contains_point()
           || edge12->next()->twin()->cell()->contains_point()
           )
        && belongs(site1.source_category(), GEOMETRY_CATEGORY_SEGMENT)
       ) {
      
      double x0 = site1.point0().x();
      double y0 = site1.point0().y();
      double x1 = site1.point1().x();
      double y1 = site1.point1().y();
      double x = new_vertex.x();
      double y = new_vertex.y();
      makefoot(x, y, x0, y0, x1, y1);
      edge12->foot(x, y);
      if (showfeet) {
      event_log_ += "<circle id=\"footcornert_" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"\" ";
      event_log_ += "cx=\""+to_str(x)+"\" cy=\""+to_str(y)+"\" r=\"220000\" class=\"evtfootc\"/>\n";
      }

      // and probably:
      //edge12->twin()->foot(edge12->vertex1()->x(), edge12->vertex1()->y());

      //event_log_ += "<circle id=\"footcorner_" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"\" ";
      //event_log_ += "cx=\""+to_str(edge12->vertex1()->x())+"\" cy=\""+to_str(edge12->vertex1()->y())+"\" r=\"50000\" class=\"evtfootc\"/>\n";


      // and probably too:

      // this had bad effect in one place in t/test, similar to case below for edge23
      //reflect(x, y, edge12->vertex0()->x(), edge12->vertex0()->y(),
      //              edge12->vertex1()->x(), edge12->vertex1()->y());
      //edge23->foot(x, y);
      if (showfeet) {
      //event_log_ += "<circle id=\"footcornerm_" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"\" ";
      //event_log_ += "cx=\""+to_str(x)+"\" cy=\""+to_str(y)+"\" r=\"220000\" class=\"evtfootc\"/>\n";
      }
      //printf("pre3_1\n");
      //printf("AH 1\n");
    }

    // edge23 into corner
    if (edge23->vertex1()
        //&& 
        //(  (edge23->vertex1()->x() == site3.point1().x() && edge23->vertex1()->y() == site3.point1().y())
        //|| (edge23->vertex1()->x() == site3.point0().x() && edge23->vertex1()->y() == site3.point0().y())
        //)
        //&& edge23->is_primary() && edge23->next() && !edge23->next()->is_primary()

        && edge23->is_primary() 
        && edge23->next() && !edge23->next()->is_primary()
        && (
           edge23->next()->cell()->contains_point()
           || 
           edge23->next()->twin()->cell()->contains_point()
           )


       && belongs(site3.source_category(), GEOMETRY_CATEGORY_SEGMENT)) {
      
      // this one for edge12 made extraneous feet
      // so commenting out here too - but without demonstrating it's wrong here too
      //edge23->twin()->foot(edge23->vertex1()->x(), edge23->vertex1()->y());
      //event_log_ += "<circle id=\"footcorneraa_" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"\" ";
      //event_log_ += "cx=\""+to_str(edge23->vertex1()->x())+"\" cy=\""+to_str(edge23->vertex1()->y())+"\" r=\"110000\" class=\"evtfootc\"/>\n";
      
      // and :
      double x0 = site3.point0().x();
      double y0 = site3.point0().y();
      double x1 = site3.point1().x();
      double y1 = site3.point1().y();
      double x = new_vertex.x();
      double y = new_vertex.y();
      makefoot(x, y, x0, y0, x1, y1);
      new_edge2.foot(x, y);
      
      if (showfeet) {
      event_log_ += "<circle id=\"footcornerbb_" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"\" ";
      event_log_ += "cx=\""+to_str(x)+"\" cy=\""+to_str(y)+"\" r=\"110000\" class=\"evtfootc\"/>\n";
      }
      // this had bad effect in couple places in t/test
      //reflect(x, y, edge23->vertex0()->x(), edge23->vertex0()->y(),
      //              edge23->vertex1()->x(), edge23->vertex1()->y());
      //edge23->foot(x, y);
      
      if (showfeet) {
      //event_log_ += "<circle id=\"footcornercc_" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"\" ";
      //event_log_ += "cx=\""+to_str(x)+"\" cy=\""+to_str(y)+"\" r=\"110000\" class=\"evtfootc\"/>\n";
      }

      //printf("AH 3\n");      
    }

    // maybe
    if (edge12->is_primary()
        && edge12->vertex1()
        && belongs(site1.source_category(), GEOMETRY_CATEGORY_POINT)
        && (edge12->vertex1()->x() == site1.point0().x() && edge12->vertex1()->y() == site1.point0().y())
        ) {
      edge12->twin()->foot(edge12->vertex1()->x(), edge12->vertex1()->y());
      if (showfeet) {
      event_log_ += "<circle id=\"footcornerm_" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"\" ";
      event_log_ += "cx=\""+to_str(edge12->vertex1()->x())+"\" cy=\""+to_str(edge12->vertex1()->y())+"\" r=\"110000\" class=\"evtfootc\"/>\n";
      }
    }
    // maybe
    if (edge23->is_primary()
        && edge23->vertex1()
        && belongs(site3.source_category(), GEOMETRY_CATEGORY_POINT)
        && (edge23->vertex1()->x() == site3.point0().x() && edge23->vertex1()->y() == site3.point0().y())
        ) {
      edge23->twin()->foot(edge23->vertex1()->x(), edge23->vertex1()->y());
      if (showfeet) {
      event_log_ += "<circle id=\"footcornerm_" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"\" ";
      event_log_ += "cx=\""+to_str(edge23->vertex1()->x())+"\" cy=\""+to_str(edge23->vertex1()->y())+"\" r=\"110000\" class=\"evtfootc\"/>\n";
      }
    }

    // special case derived from visual debug
    if (   belongs(site1.source_category(), GEOMETRY_CATEGORY_POINT)
        && belongs(site3.source_category(), GEOMETRY_CATEGORY_SEGMENT)
       ) {
      if (
         (    (site3.point0().x() == site1.point0().x() && site3.point0().y() == site1.point0().y())
           || (site3.point1().x() == site1.point0().x() && site3.point1().y() == site1.point0().y())
         )
         && edge23->is_primary()
         && (edge23->vertex0()->x() == site1.point0().x() && edge23->vertex0()->y() == site1.point0().y())
         ) {
        
        // visually, looked like this was already there
        // but maybe that was wrongly associated with another edge?
        // so set/reset to see...
        // yeah this wasn't set, even though there appeared to be a foot there
        // so that was probably a wrong foot from something else
        // and this is maybe right
        edge23->foot(site1.point0().x(), site1.point0().y());


        double x0 = site3.point0().x();
        double y0 = site3.point0().y();
        double x1 = site3.point1().x();
        double y1 = site3.point1().y();
        double x = edge23->vertex1()->x();
        double y = edge23->vertex1()->y();
        makefoot(x, y, x0, y0, x1, y1);
        edge23->twin()->foot(x, y);
        if (showfeet) {
        event_log_ += "<circle id=\"footcornerm_" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"\" ";
        event_log_ += "cx=\""+to_str(x)+"\" cy=\""+to_str(y)+"\" r=\"110000\" class=\"evtfootc\"/>\n";
        }
        reflect(x, y, edge23->vertex0()->x(), edge23->vertex0()->y(),
                      edge23->vertex1()->x(), edge23->vertex1()->y());
        edge23->next()->foot(x, y);
        if (showfeet) {
        event_log_ += "<circle id=\"footcornerm_" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"\" ";
        event_log_ += "cx=\""+to_str(x)+"\" cy=\""+to_str(y)+"\" r=\"110000\" class=\"evtfootc\"/>\n";
        }
        //printf("\nMAYBEMAYBE 1\n");
      }
    }
    // similar special case derived from above, but maybe doesn't happen?
    if (   belongs(site3.source_category(), GEOMETRY_CATEGORY_POINT)
        && belongs(site1.source_category(), GEOMETRY_CATEGORY_SEGMENT)
       ) {
      if (
         (site1.point0().x() == site3.point0().x() && site1.point0().y() == site3.point0().y())
         ||
         (site1.point1().x() == site3.point0().x() && site1.point1().y() == site3.point0().y())
         ) {
        //printf("\nMAYBEMAYBE 2\n");
      }
    }


    // for straight edges
    if (edge23->vertex1() && !edge23->twin()->foot() 
        && edge23->is_linear() 
        //&& edge23->twin()->is_primary() 
        && belongs(site3.source_category(), GEOMETRY_CATEGORY_SEGMENT)) {
      double x0 = site3.point0().x();
      double y0 = site3.point0().y();
      double x1 = site3.point1().x();
      double y1 = site3.point1().y();
      double x = edge23->vertex1()->x();
      double y = edge23->vertex1()->y();
      makefoot(x, y, x0, y0, x1, y1);
      edge23->twin()->foot(x, y);
      if (showfeet) {
      event_log_ += "<circle id=\"footcornerm_" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"\" ";
      event_log_ += "cx=\""+to_str(x)+"\" cy=\""+to_str(y)+"\" r=\"110000\" class=\"evtfootc\"/>\n";
      }
    }

    

    if (edge23->vertex1() && !new_edge2.foot() 
        && new_edge2.is_linear() 
        //&& new_edge2.is_primary() 
        && belongs(site3.source_category(), GEOMETRY_CATEGORY_SEGMENT)) {
      double x0 = site3.point0().x();
      double y0 = site3.point0().y();
      double x1 = site3.point1().x();
      double y1 = site3.point1().y();
      double x = new_vertex.x();
      double y = new_vertex.y();
      makefoot(x, y, x0, y0, x1, y1);
      new_edge2.foot(x, y);
      if (showfeet) {
      event_log_ += "<circle id=\"footcornerm_" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"\" ";
      event_log_ += "cx=\""+to_str(x)+"\" cy=\""+to_str(y)+"\" r=\"110000\" class=\"evtfootc\"/>\n";
      }
    }

    if (!edge12->foot() 
        && edge12->is_linear() 
        //&& edge12->is_primary() 
        && belongs(site1.source_category(), GEOMETRY_CATEGORY_SEGMENT)
        && edge12->vertex1()
        ) {
      double x0 = site1.point0().x();
      double y0 = site1.point0().y();
      double x1 = site1.point1().x();
      double y1 = site1.point1().y();
      double x = new_vertex.x();
      double y = new_vertex.y();
      makefoot(x, y, x0, y0, x1, y1);
      edge12->foot(x, y);
      if (showfeet) {
      event_log_ += "<circle id=\"footcornerm_" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"\" ";
      event_log_ += "cx=\""+to_str(x)+"\" cy=\""+to_str(y)+"\" r=\"330000\" class=\"evtfootc\"/>\n";
      }
      if (new_vertex.r() == 0 && edge12->vertex1()) {
        edge12->twin()->foot(edge12->vertex1()->x(), edge12->vertex1()->y());
        if (showfeet) {
        event_log_ += "<circle id=\"footcornerm_" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"\" ";
        event_log_ += "cx=\""+to_str(edge12->vertex1()->x())+"\" cy=\""+to_str(edge12->vertex1()->y())+"\" r=\"220000\" class=\"evtfootc\"/>\n";
        }
      }
    }

    // didn't see change with this
    // might be redundant - or not right
    // thinking not right, though it is picking up a case that the first
    // foot finding conditional should get
    // ... yeah doesn't seem right
    if (false
        && edge12->vertex1()
        && !edge12->next()->foot() 
        && edge12->is_linear() 
        //&& edge12->is_primary() 
        && belongs(site1.source_category(), GEOMETRY_CATEGORY_SEGMENT)) {
      double x0 = site1.point0().x();
      double y0 = site1.point0().y();
      double x1 = site1.point1().x();
      double y1 = site1.point1().y();
      double x = edge12->vertex1()->x();
      double y = edge12->vertex1()->y();
      makefoot(x, y, x0, y0, x1, y1);
      edge12->next()->foot(x, y);
      if (showfeet) {
      event_log_ += "<circle id=\"footcornerm_" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"\" ";
      event_log_ += "cx=\""+to_str(x)+"\" cy=\""+to_str(y)+"\" r=\"220000\" class=\"evtfootc\"/>\n";
      }
    }
    
    // trying to fill in corner feet
    // seems to work without creating feet where they don't belong
    if (!edge12->is_primary()
        &&  edge12->vertex1()
        && !edge12->twin()->prev()->is_primary()
        &&  edge12->next()->is_primary()
        && !edge12->next()->foot()
       ) {
      edge12->next()->foot(edge12->vertex1()->x(), edge12->vertex1()->y());
      if (showfeet) {
      event_log_ += "<circle id=\"footcornerm_" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"\" ";
      event_log_ += "cx=\""+to_str(edge12->vertex1()->x())+"\" cy=\""+to_str(edge12->vertex1()->y())+"\" r=\"30000\" class=\"evtfootc\"/>\n";
      }
    }
    if (   !edge23->is_primary()
        &&  edge23->vertex1()
        && !edge23->twin()->prev()->is_primary()
        &&  edge23->next()->is_primary()
        && !edge23->next()->foot()
       ) {
      edge23->next()->foot(edge23->vertex1()->x(), edge23->vertex1()->y());
      if (showfeet) {
      event_log_ += "<circle id=\"footcornerm_" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"\" ";
      event_log_ += "cx=\""+to_str(edge23->vertex1()->x())+"\" cy=\""+to_str(edge23->vertex1()->y())+"\" r=\"40000\" class=\"evtfootc\"/>\n";
      }
    }

    
    // another based on special case, that hopefully has general utility
    // this worked for that one special case
    if (  !edge23->foot()
        && edge23->vertex1()
        && edge23->is_linear()
        && edge23->twin()->next()->foot()
       ) {
        double x = edge23->twin()->next()->foot()->x();
        double y = edge23->twin()->next()->foot()->y();
        reflect(x, y, edge23->vertex0()->x(), edge23->vertex0()->y(),
                      edge23->vertex1()->x(), edge23->vertex1()->y());
        edge23->foot(x, y);
        if (showfeet) {
        event_log_ += "<circle id=\"footcornerm_" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"\" ";
        event_log_ += "cx=\""+to_str(x)+"\" cy=\""+to_str(y)+"\" r=\"40000\" class=\"evtfootc\"/>\n";
        }
    }
    // same as above but for edge12 - no test case to demonstrate it yet though
    if (  !edge12->foot()
        && edge12->vertex1()
        && edge12->is_linear()
        && edge12->twin()->next()->foot()
       ) {
        double x = edge12->twin()->next()->foot()->x();
        double y = edge12->twin()->next()->foot()->y();
        reflect(x, y, edge12->vertex0()->x(), edge12->vertex0()->y(),
                      edge12->vertex1()->x(), edge12->vertex1()->y());
        edge12->foot(x, y);
        if (showfeet) {
        event_log_ += "<circle id=\"footcornerm_" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"\" ";
        event_log_ += "cx=\""+to_str(x)+"\" cy=\""+to_str(y)+"\" r=\"40000\" class=\"evtfootc\"/>\n";
        }
    }

    // another special case
    // got it - that was last of first round of many missing feet
    if (   !edge23->foot()
        && !edge12->is_primary()
        //&& !edge23->is_curved()
        &&  edge23->is_primary()
        &&  edge12->vertex1()
       ) {
      edge23->foot(edge12->vertex1()->x(), edge12->vertex1()->y());
    }
    // might need similar but not same as above for similar edge12 case
    // if that's possible, but should demonstrate or illustrate need for that


    // curved edges
    
    // on t/test this might not have significant effect with or without
    // ... this fixes missing feet on full hex grid test
    if (
        //!edge12->foot() && 
        edge12->is_curved()) {
      if (belongs(site1.source_category(), GEOMETRY_CATEGORY_SEGMENT)) {
        double x0 = site1.point0().x();
        double y0 = site1.point0().y();
        double x1 = site1.point1().x();
        double y1 = site1.point1().y();
        double x = new_vertex.x();
        double y = new_vertex.y();
        makefoot(x, y, x0, y0, x1, y1);
        edge12->foot(x, y);
        //printf("pre3_1\n");
        if (showfeet) {
        event_log_ += "<circle id=\"footpre3_1_" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"\" ";
        event_log_ += "cx=\""+to_str(x)+"\" cy=\""+to_str(y)+"\" r=\"110000\" class=\"evtfoot2\"/>\n";
        }
      } else if (belongs(site1.source_category(), GEOMETRY_CATEGORY_POINT)) {
        edge12->foot(site1.point0().x(), site1.point0().y());
        //printf("pre3_2\n");
        if (showfeet) {
        event_log_ += "<circle id=\"footpre3_2_" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"\" ";
        event_log_ += "cx=\""+to_str(site1.point0().x())+"\" cy=\""+to_str(site1.point0().y())+"\" r=\"110000\" class=\"evtfoot2\"/>\n";
        }
      }
    }

    // on t/test this is not good
    if (false && 
        //!edge23->twin()->foot() && 
        edge23->is_curved()) {
      if (belongs(site3.source_category(), GEOMETRY_CATEGORY_SEGMENT)) {
        double x0 = site3.point0().x();
        double y0 = site3.point0().y();
        double x1 = site3.point1().x();
        double y1 = site3.point1().y();
        double x = new_vertex.x();
        double y = new_vertex.y();
        makefoot(x, y, x0, y0, x1, y1);
        //edge23->foot(x, y);
        edge23->twin()->foot(x, y);
        //printf("NEWpre3_1\n");
        if (showfeet) {
        event_log_ += "<circle id=\"footNEWpre3_1_" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"\" ";
        event_log_ += "cx=\""+to_str(x)+"\" cy=\""+to_str(y)+"\" r=\"110000\" class=\"evtfoot2\"/>\n";
        }
      }
      if (belongs(site3.source_category(), GEOMETRY_CATEGORY_POINT)) {
        //edge12->foot(site3.point0().x(), site3.point0().y());
        edge23->twin()->foot(site3.point0().x(), site3.point0().y());
        //printf("NEWpre3_2\n");
        if (showfeet) {
        event_log_ += "<circle id=\"footNEWpre3_2_" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"\" ";
        event_log_ += "cx=\""+to_str(site3.point0().x())+"\" cy=\""+to_str(site3.point0().y())+"\" r=\"110000\" class=\"evtfoot2\"/>\n";
        }
      }
    }
    
    if (edge12->twin()->foot() && edge12->twin()->cell()->contains_point()
        && !edge23->foot()) {
        edge23->foot(edge12->twin()->foot()->x(), edge12->twin()->foot()->y());
        //printf("around point\n");
        if (showfeet) {
        event_log_ += "<circle id=\"footNEWpre3_2_" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"\" ";
        event_log_ += "cx=\""+to_str(edge12->twin()->foot()->x())+"\" cy=\""+to_str(edge12->twin()->foot()->y())+"\" r=\"110000\" class=\"evtfoot2\"/>\n";
        }
    }

    if (edge23->foot() && edge23->cell()->contains_point()
        && edge23->next() && !edge23->next()->foot() 
        && edge23->next()->is_primary()) {
        // rare
        edge23->next()->foot(edge23->foot()->x(), edge23->foot()->y());
        //printf("around point SOME MORE\n");
        if (showfeet) {
        event_log_ += "<circle id=\"footNEWpre3_2_" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"\" ";
        event_log_ += "cx=\""+to_str(edge23->foot()->x())+"\" cy=\""+to_str(edge23->foot()->y())+"\" r=\"110000\" class=\"evtfoot2\"/>\n";
        }
    }

    if (edge23->twin()->foot() && edge23->twin()->cell()->contains_point()
        && !new_edge2.foot()) {
        new_edge2.foot(edge23->twin()->foot()->x(), edge23->twin()->foot()->y());
        //printf("around point TO NEW 2\n");
        if (showfeet) {
        event_log_ += "<circle id=\"footNEWpre3_2_" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"\" ";
        event_log_ += "cx=\""+to_str(edge23->twin()->foot()->x())+"\" cy=\""+to_str(edge23->twin()->foot()->y())+"\" r=\"110000\" class=\"evtfoot2\"/>\n";
        }
    }
    
    if (edge12->foot() && edge12->cell()->contains_point()
        && !new_edge1.foot()) {
        new_edge1.foot(edge12->foot()->x(), edge12->foot()->y());
        //printf("around point TO NEW 1\n");
        if (showfeet) {
        event_log_ += "<circle id=\"footNEWpre3_2_" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"\" ";
        event_log_ += "cx=\""+to_str(edge12->foot()->x())+"\" cy=\""+to_str(edge12->foot()->y())+"\" r=\"110000\" class=\"evtfoot2\"/>\n";
        }
    }

    // derived from lingering special case on hex fan housing (way above hex mesh)
    // ... but also handles at least three feet in t/test, so this is legit
    // ... yes this is happening in several cases that are also handled by other
    // cases above, and then a few those miss
    if (
       !edge23->foot()
       && edge23->is_curved() && new_edge2.is_curved()
       && edge23->twin()->cell()->contains_point()
       && new_edge2.cell()->contains_point()
       && edge12->foot()
       ) {
      double x = edge12->foot()->x();
      double y = edge12->foot()->y();
      reflect(x, y, edge12->vertex0()->x(), edge12->vertex0()->y(),
                   edge12->vertex1()->x(), edge12->vertex1()->y());
      edge23->foot(x, y);
      event_log_ += "<circle id=\"footNEWpre3_2_" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"\" ";
      event_log_ += "cx=\""+to_str(x)+"\" cy=\""+to_str(y)+"\" r=\"110000\" class=\"evtfootworking\"/>\n";
    }

    if (!is_linear) {
      if (belongs(site1.source_category(), GEOMETRY_CATEGORY_POINT)
          && belongs(site3.source_category(), GEOMETRY_CATEGORY_SEGMENT)
         ) {
        if (new_edge2.vertex0()) {
          double x0 = site3.point0().x();
          double y0 = site3.point0().y();
          double x1 = site3.point1().x();
          double y1 = site3.point1().y();
          double x = new_edge2.vertex0()->x();
          double y = new_edge2.vertex0()->y();
          //printf("mf3");
          makefoot(x, y, x0, y0, x1, y1);
          //printf("\n");
          new_edge2.foot(x, y);
         
          if (showfeet) {
          event_log_ += "<circle id=\"foot3_" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"\" ";
          event_log_ += "cx=\""+to_str(x)+"\" cy=\""+to_str(y)+"\" r=\"110000\" class=\"evtfoot2\"/>\n";
          }

        }
      }
      
      if (belongs(site3.source_category(), GEOMETRY_CATEGORY_POINT)
          && belongs(site1.source_category(), GEOMETRY_CATEGORY_SEGMENT)) {
        
        new_edge2.foot(site3.point0().x(),  site3.point0().y());
      
        if (showfeet) {
        event_log_ += "<circle id=\"foot_after_3_" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"\" ";
        event_log_ += "cx=\""+to_str(site3.point0().x())+"\" cy=\""+to_str(site3.point0().y())+"\" r=\"110000\" class=\"evtfoot2\"/>\n";
        }

        //printf("premf4\n");
        if (
          false
          //new_edge2.vertex0() 
          //&& new_edge2.vertex1()
          
        ) {
          double x0 = site1.point0().x();
          double y0 = site1.point0().y();
          double x1 = site1.point1().x();
          double y1 = site1.point1().y();
          double x = new_edge1.vertex0()->x();
          double y = new_edge1.vertex0()->y();
          //double x = new_vertex.x();
          //double y = new_vertex.y();
          //printf("mf4");
          makefoot(x, y, x0, y0, x1, y1);
          //printf("\n");
          // We were getting in here often, but this foot is probably 
          // getting reset most of the time by one of the other conditionals.
          // Only noticed one foot affected by perturbing this foot.
          new_edge1.foot(x, y);
          if (showfeet) {
          event_log_ += "<circle id=\"foot4_" + to_str((UV) &site1) + "_"+to_str((UV) &site3)+"\" ";
          event_log_ += "cx=\""+to_str(x)+"\" cy=\""+to_str(y)+"\" r=\"110000\" class=\"evtfoot2\"/>\n";
          }
        }
      }
    }
    
    event_log_ += "</g>\n";

    // Return a pointer to the new half-edge.
    return std::make_pair(&new_edge1, &new_edge2);
  }

  void _build() {
    // Remove degenerate edges.
    edge_iterator last_edge = edges_.begin();
    for (edge_iterator it = edges_.begin(); it != edges_.end(); it += 2) {
      const vertex_type* v1 = it->vertex0();
      const vertex_type* v2 = it->vertex1();
      if (v1 && v2 && vertex_equality_predicate_(*v1, *v2)) {
        remove_edge(&(*it));
      }
      else {
        if (it != last_edge) {
          edge_type* e1 = &(*last_edge = *it);
          edge_type* e2 = &(*(last_edge + 1) = *(it + 1));

          e1->twin(e2);
          e2->twin(e1);
          if (e1->prev()) {
            e1->prev()->next(e1);
            e2->next()->prev(e2);
          }
          if (e2->prev()) {
            e1->next()->prev(e1);
            e2->prev()->next(e2);
          }
        }
        last_edge += 2;
      }
    }
    edges_.erase(last_edge, edges_.end());

    // Set up incident edge pointers for cells and vertices.
    for (edge_iterator it = edges_.begin(); it != edges_.end(); ++it) {
      it->cell()->incident_edge(&(*it));
      if (it->vertex0()) {
        it->vertex0()->incident_edge(&(*it));
      }
    }

    // Remove degenerate vertices.
    vertex_iterator last_vertex = vertices_.begin();
    for (vertex_iterator it = vertices_.begin(); it != vertices_.end(); ++it) {
      if (it->incident_edge()) {
        if (it != last_vertex) {
          *last_vertex = *it;
          vertex_type* v = &(*last_vertex);
          edge_type* e = v->incident_edge();
          do {
            e->vertex0(v);
            e = e->rot_next();
          } while (e != v->incident_edge());
        }
        ++last_vertex;
      }
    }
    vertices_.erase(last_vertex, vertices_.end());

    // Set up next/prev pointers for infinite edges.
    if (vertices_.empty()) {
      if (!edges_.empty()) {
        // Update prev/next pointers for the line edges.
        edge_iterator edge_it = edges_.begin();
        edge_type* edge1 = &(*edge_it);
        edge1->next(edge1);
        edge1->prev(edge1);
        ++edge_it;
        edge1 = &(*edge_it);
        ++edge_it;

        while (edge_it != edges_.end()) {
          edge_type* edge2 = &(*edge_it);
          ++edge_it;

          edge1->next(edge2);
          edge1->prev(edge2);
          edge2->next(edge1);
          edge2->prev(edge1);

          edge1 = &(*edge_it);
          ++edge_it;
        }

        edge1->next(edge1);
        edge1->prev(edge1);
      }
    } else {
      // Update prev/next pointers for the ray edges.
      for (cell_iterator cell_it = cells_.begin();
         cell_it != cells_.end(); ++cell_it) {
        if (cell_it->is_degenerate())
          continue;
        // Move to the previous edge while
        // it is possible in the CW direction.
        edge_type* left_edge = cell_it->incident_edge();
        while (left_edge->prev() != NULL) {
          left_edge = left_edge->prev();
          // Terminate if this is not a boundary cell.
          if (left_edge == cell_it->incident_edge())
            break;
        }

        if (left_edge->prev() != NULL)
          continue;

        edge_type* right_edge = cell_it->incident_edge();
        while (right_edge->next() != NULL)
          right_edge = right_edge->next();
        left_edge->prev(right_edge);
        right_edge->next(left_edge);
      }
    }

    // The above gets us the complete Voronoi diagram.
    // Now we'll narrow that down to just the medial axis for the polygon.
    
    if (0) { // handy data dump to copy-paste between stages while debugging
    for (edge_iterator it = edges_.begin(); it != edges_.end(); ++it) {
      printf("edge %lld: %lld, %lld, %lld, %ld, %ld, %s, %s, %s, %s, %s\n",
          (long long unsigned int) &(*it),
          (long long unsigned int) it->twin(),
          (long long unsigned int) it->next(),
          (long long unsigned int) it->prev(),
          it->color(),
          it->cell()->source_index(),
          it->is_curved()?"curved":"      ",
          it->is_finite()?"finite":"      ",
          it->is_primary()?"primary":"      ",
          it->twin() == it->next() ? "next=twin":"twok",
          &(*it) == it->next() ? "next=itself":"nxtok"
      );
    }
    }

    // Mark edges exterior to the polygon by setting color attribute to 1.
    // (Adjacent vertices and cells are also marked.)
    for (edge_iterator it = edges_.begin(); it != edges_.end(); ++it) {
      if (!it->is_finite()) { mark_exterior(&(*it)); }
    }

    // Now all the cells associated with the polygon's outer contour segments 
    // have color() == 1, while all cells associated with holes still have 
    // color() == 0. This isn't always enough information to label all edges
    // inside holes correctly. We'll go ahead and label edges not associated
    // with the outer cells as edges in holes, and then later correct
    // mislabeled edges.

    for (edge_iterator it = edges_.begin(); it != edges_.end(); ++it) {
      if (  it->cell()->color() == 0
         // cells with color 0 at this point are either holes
         // or regions within the polygon associated with concavites
         && it->twin()->cell()->color() == 0
         // this avoids labeling edges coming directly off
         // of the inner side of the medial axis that surrounds a hole
         && (  it->next()->twin()->cell()->color() != 1
            && it->prev()->twin()->cell()->color() != 1)
         ) {
        it->color(1);
      }
    }

    // Now we find cells with a mix of primary edges labeled as inside and 
    // outside. Adjacent primary edges can't have different inside-outside 
    // status. We're sure about the edges we've labled as within the polygon 
    // so far. So we recursively label adjacent primary edges as within if they
    // don't have that label already, and non-primary edges associated with
    // curved edges get their labels fixed too.

    for (cell_iterator it = cells_.begin(); it != cells_.end(); ++it) {
      //printf("    cell source_index %ld\n",it->source_index());
      edge_type* e = it->incident_edge();
      do {
        if (e->is_primary() && e->next()->is_primary()) {
          if (e->color() == 0 && e->next()->color() != 0) {
            //printf("    start first recurse\n");
            mark_interior(e->next());
            } 
          if (e->color() != 0 && e->next()->color() == 0) {
            //printf("    start second recurse\n");
            mark_interior(e, true);
            }
        }
        e = e->next();
      } while (e != it->incident_edge());
    }

    // Deeper edges still escape recursive stuff above
    
    
    
    
    
    
    
// Adjust this or copy modify to capture a few more cases
// prob bring back req that at least one is primary
// and consider looking at consective primary by doing a next->twin->next
// jump over any secondary () and then if you do a 1 to 0 color conversion
// probably change color of the hopped-over secondary too.



    for (edge_iterator it = edges_.begin(); it != edges_.end(); ++it) {
      if (
          //it->is_primary() &&
           (  
           //it->next()->is_primary() && 
           it->next()->color() == 0
           && 
           //it->prev()->is_primary() && 
           it->prev()->color() == 0
           )
         ) {
        it->color(0);
        it->twin()->color(0);
      }
    }


    // yes, this fixed some missed cases
    for (edge_iterator it = edges_.begin(); it != edges_.end(); ++it) {
      if (
          it->is_primary() && it->color() == 0
          && !it->next()->is_primary() 
          && it->next()->twin()->next()->is_primary()
          && it->next()->twin()->next()->color() == 1
         ) {
        it->next()->twin()->next()->color(0);
        it->next()->twin()->next()->twin()->color(0);
        // the hopped-over non primaries should be changed too
        it->next()->color(0);
        it->next()->twin()->color(0);
      }
    }
    
    
    
    
    // Some non-primary edges of more complex cells don't get corrected
    // by the steps above. But now they're easy to identify and fix.
    
    for (edge_iterator it = edges_.begin(); it != edges_.end(); ++it) {
      if (!it->is_primary() &&
           (  it->rot_next()->is_primary() && it->rot_next()->color() == 0
           && it->rot_prev()->is_primary() && it->rot_prev()->color() == 0
           )
         ) {
        it->color(0);
        it->twin()->color(0);
      }
    }
//event_log_ += "<!-- magentas -->\n";

    // Still missing some - think it's secondaries.
    for (edge_iterator it = edges_.begin(); it != edges_.end(); ++it) {
      if (
          !it->is_primary() && it->color() == 1
          && it->vertex0() && it->vertex1()
           &&
           (
           (
            it->twin()->prev()->is_curved() &&
           //&& (it->twin()->prev()->is_curved() 
           //    || (it->twin()->prev()->is_finite() && it->twin()->prev()->is_primary()) )
           it->prev()->color() == 1 && it->prev()->is_primary()
           && it->twin()->prev()->color() == 0
           && it->twin()->prev()->prev()->color() == 0
           )
           //||
           //(
           //&& (it->prev()->is_curved() 
           //    || (it->prev()->is_finite() && it->prev()->is_primary()) )
           //it->next()->color() == 1 && it->next()->is_primary()
           //it->prev()->color() == 0
           //&& it->prev()->prev()->color() == 0
           //)
           )
         ) {
        printf("\n\nYES\n\n");
        //it->color(0);
        //it->twin()->color(0);
        if (it->vertex0() && it->vertex1()) {
event_log_ += "<line style=\"stroke-width:200000;stroke:magenta;\" ";
event_log_ += "x1=\""+to_str(it->vertex0()->x())+"\" y1=\""+to_str(it->vertex0()->y())+"\" ";
event_log_ += "x2=\""+to_str(it->vertex1()->x())+"\" y2=\""+to_str(it->vertex1()->y())+"\" />\n";
}
      }
    }


bool dbge = false;
if (dbge) {
for (edge_iterator it = edges_.begin(); it != edges_.end(); ++it) {
edge_type* edge = &(*it);
if (edge->color() == 1 && edge->vertex0() && edge->vertex1()) {
event_log_ += "<line style=\"stroke-width:300000;stroke:yellow;\" ";
event_log_ += "x1=\""+to_str(edge->vertex0()->x())+"\" y1=\""+to_str(edge->vertex0()->y())+"\" ";
event_log_ += "x2=\""+to_str(edge->vertex1()->x())+"\" y2=\""+to_str(edge->vertex1()->y())+"\" />\n";
if (edge->is_primary()) {
event_log_ += "<line style=\"stroke-width:200000;stroke:orange;\" ";
event_log_ += "x1=\""+to_str(edge->vertex0()->x())+"\" y1=\""+to_str(edge->vertex0()->y())+"\" ";
event_log_ += "x2=\""+to_str(edge->vertex1()->x())+"\" y2=\""+to_str(edge->vertex1()->y())+"\" />\n";
} else if (!edge->is_primary()) {
event_log_ += "<line style=\"stroke-width:200000;stroke:aqua;\" ";
event_log_ += "x1=\""+to_str(edge->vertex0()->x())+"\" y1=\""+to_str(edge->vertex0()->y())+"\" ";
event_log_ += "x2=\""+to_str(edge->vertex1()->x())+"\" y2=\""+to_str(edge->vertex1()->y())+"\" />\n";
}
}
}
}

    // Now all edges within the polygon have color() == 0 and all edges
    // outside of the polygon or inside holes in the polygon have 
    // color() == 1.

    /////////////
    // At this point we modify the half edge graph to better represent the 
    // the medial axis.
    // The main thing to do is update next() and prev() pointers to follow
    // along the primary edges that represent the medial axis, instead
    // of having them point just to the next/prev within each Voronoi cell.

    // Get the edge corresponding to the first polygon input segment
    // so the first loop we traverse is the outer polygon loop.
    // Currently it doesn't matter that we process that first, but we
    // may want to enhance the output data structure later to reflect 
    // inside vs outside/in-concavity/in-hole medial axis edge collections.

    edge_type * start_edge = NULL;
    for (edge_iterator it = edges_.begin(); it != edges_.end(); ++it) {
      if (it->cell()->source_index() == 0
          && it->color() == 0
          && it->is_primary()
         ) {
        start_edge = &(*it);
        break;
      }
    }

    // Walk the edge references and modify to represent medial axis.
    while (start_edge != NULL) {
      edge_type * edge = start_edge;
      //start_edge->color(2);
      do {
        // mark visited internal edges (will restore to 0 afterward)
        edge->color(2);

        // if next edge is within polygon
        if (edge->next()->color() == 0 || edge->next()->color() == 2) {
          if (edge->next()->is_primary()) { 
            // go to next edge within same cell 
            edge = edge->next();
          } else { 
            // skip over a non-primary edge to the primary edge that follows it
            edge_type* prev = edge;
            edge = edge->next()->twin()->next();
            // get foot from non-primary endpoint and
            // mirror foot info from the non-primary to the twin
            if (prev->twin()->vertex0() && prev->twin()->vertex1() && prev->next()->vertex1()) {
              // The reflect about line case is simple:
              double x = prev->next()->vertex1()->x();
              double y = prev->next()->vertex1()->y();
              
              if (!edge->foot()) {
                edge->foot(x, y);
              }
              
              double x0 = prev->twin()->vertex0()->x();
              double y0 = prev->twin()->vertex0()->y();
              double x1, y1;
              // would like to already have foot in place
              // but not quite there yet, and this performs well
              if (true || !prev->twin()->foot()) {
                if (!prev->twin()->is_curved()) {
                  double x1 = prev->twin()->vertex1()->x() + 0;
                  double y1 = prev->twin()->vertex1()->y() + 0;
                  reflect(x, y, x0, y0, x1, y1);
                  prev->twin()->foot(x, y);
                  //printf("reflect foot to line\n");
                } else {
                  // The case for a curved edge isn't as simple, but 
                  // it seems most feet in this case are already properly
                  // calculated by the event-processing code.
                  // It may be that we never get to, or should never need to 
                  // get into this else{}. Maybe eliminate this after foot-finding
                  // in the event-processing has been fully understood and
                  // implemented.
                  // ... "reflect foot for parabola" still happens a lot -
                  //     still depending on it
                  if (!prev->twin()->prev()->is_primary()) {
                    //printf("from prev twin prev non-primary\n");
                    prev->twin()->foot(prev->twin()->prev()->vertex0()->x(),
                                       prev->twin()->prev()->vertex0()->y()
                                      );       
                  }
                  else {
                    double x = prev->next()->vertex1()->x();
                    double y = prev->next()->vertex1()->y();
                    if (!edge->is_curved()) {
                      double x0 = edge->vertex0()->x();
                      double y0 = edge->vertex0()->y();
                      double x1 = edge->vertex1()->x();
                      double y1 = edge->vertex1()->y();
                      reflect(x, y, x0, y0, x1, y1);
                      prev->twin()->foot(x, y);
                      //printf("reflect foot for parabola\n");
                    }                  
                  }
                }
              }
            }
            // first make the clipped-out edges ahead link to themselves
            prev->next()->twin()->next(prev->next());
            prev->next()->prev(prev->next()->twin());
            // now link this to new next
            prev->next(edge);
            edge->prev(prev);
          }
        } else {
          // corner - end touches polygon, so turn around
          edge_type * prev = edge;
          edge = edge->twin();
if (dbge) {
event_log_ += "<line style=\"stroke-width:15000;stroke:red;\" ";
event_log_ += "x1=\""+to_str(prev->vertex0()->x())+"\" y1=\""+to_str(prev->vertex0()->y())+"\" ";
event_log_ += "x2=\""+to_str(prev->vertex1()->x())+"\" y2=\""+to_str(prev->vertex1()->y())+"\" />\n";
event_log_ += "<line style=\"stroke-width:15000;stroke:red;\" ";
event_log_ += "x1=\""+to_str(edge->vertex0()->x())+"\" y1=\""+to_str(edge->vertex0()->y())+"\" ";
event_log_ += "x2=\""+to_str(edge->vertex1()->x())+"\" y2=\""+to_str(edge->vertex1()->y())+"\" />\n";
}

// This may be obsolete now that we're doing corner foot processing in site events
// and it actually gives some bad results only in combination with the new
// processing there.
// It also may be that missing foot code later picks up some of these ? when we
// don't do this here.
if (false) {
          // figure feet
          double theta = atan2(edge->vertex1()->y() - edge->vertex0()->y(),
                               edge->vertex1()->x() - edge->vertex0()->x());
          double footx = prev->vertex0()->x() + prev->vertex0()->r();
          double footy = prev->vertex0()->y();
          // This should always come out <= 1 and > 0.
          // Sometimes it spills over just beyond 1 in the case of 
          // a corner so shallow it's practically flat.
          // So snap to 1 if > 1.
          double for_acos = prev->vertex0()->r() 
            / sqrt(pow(prev->vertex1()->x() - prev->vertex0()->x(),2)
                 + pow(prev->vertex1()->y() - prev->vertex0()->y(),2)
                  );
          if (for_acos > 1) { for_acos = 1; }
          double phi = acos( for_acos );

          rotate_2d(footx, footy, theta + phi, prev->vertex0()->x(), 
                                               prev->vertex0()->y()
          );
          
          if (!prev->foot()) {
          prev->foot(footx, footy);
          }

          rotate_2d(footx, footy, -2*phi, prev->vertex0()->x(), 
                                          prev->vertex0()->y()
          );
          
          if (!edge->next()->foot()) {
          edge->next()->foot(footx, footy);
          }

          if (!edge->foot()) {
          edge->foot(edge->vertex0()->x(), edge->vertex0()->y());
          }
}
          // first connect edges ahead to eachother
          prev->next()->prev(edge->prev());
          edge->prev()->next(prev->next());
if (dbge) {
if (prev->next()->vertex0() && prev->next()->vertex1()) {
event_log_ += "<line style=\"stroke-width:15000;stroke:blue;\" ";
event_log_ += "x1=\""+to_str(prev->next()->vertex0()->x())+"\" y1=\""+to_str(prev->next()->vertex0()->y())+"\" ";
event_log_ += "x2=\""+to_str(prev->next()->vertex1()->x())+"\" y2=\""+to_str(prev->next()->vertex1()->y())+"\" />\n";
}
if (edge->prev()->vertex0() && edge->prev()->vertex1()) {
event_log_ += "<line style=\"stroke-width:15000;stroke:blue;\" ";
event_log_ += "x1=\""+to_str(edge->prev()->vertex0()->x())+"\" y1=\""+to_str(edge->prev()->vertex0()->y())+"\" ";
event_log_ += "x2=\""+to_str(edge->prev()->vertex1()->x())+"\" y2=\""+to_str(edge->prev()->vertex1()->y())+"\" />\n";
}

if (prev->next()->prev()->vertex1() && prev->next()->prev()->vertex0()) {
event_log_ += "<line style=\"stroke-width:10000;stroke:gray;\" ";
event_log_ += "x1=\""+to_str(prev->next()->prev()->vertex0()->x())+"\" y1=\""+to_str(prev->next()->prev()->vertex0()->y())+"\" ";
event_log_ += "x2=\""+to_str(prev->next()->prev()->vertex1()->x())+"\" y2=\""+to_str(prev->next()->prev()->vertex1()->y())+"\" />\n";
}
if (edge->prev()->next()->vertex0() && edge->prev()->next()->vertex1()) {
event_log_ += "<line style=\"stroke-width:10000;stroke:gray;\" ";
event_log_ += "x1=\""+to_str(edge->prev()->next()->vertex0()->x())+"\" y1=\""+to_str(edge->prev()->next()->vertex0()->y())+"\" ";
event_log_ += "x2=\""+to_str(edge->prev()->next()->vertex1()->x())+"\" y2=\""+to_str(edge->prev()->next()->vertex1()->y())+"\" />\n";
}
}
          // now link the corner edges together
          prev->next(edge);
          edge->prev(prev);
if (dbge) {
event_log_ += "<line style=\"stroke-width:10000;stroke:pink;\" ";
event_log_ += "x1=\""+to_str(prev->vertex0()->x())+"\" y1=\""+to_str(prev->vertex0()->y())+"\" ";
event_log_ += "x2=\""+to_str(prev->vertex1()->x())+"\" y2=\""+to_str(prev->vertex1()->y())+"\" />\n";
event_log_ += "<line style=\"stroke-width:10000;stroke:pink;\" ";
event_log_ += "x1=\""+to_str(edge->vertex0()->x())+"\" y1=\""+to_str(edge->vertex0()->y())+"\" ";
event_log_ += "x2=\""+to_str(edge->vertex1()->x())+"\" y2=\""+to_str(edge->vertex1()->y())+"\" />\n";

event_log_ += "<line style=\"stroke-width:5000;stroke:purple;\" ";
event_log_ += "x1=\""+to_str(prev->next()->vertex0()->x())+"\" y1=\""+to_str(prev->next()->vertex0()->y())+"\" ";
event_log_ += "x2=\""+to_str(prev->next()->vertex1()->x())+"\" y2=\""+to_str(prev->next()->vertex1()->y())+"\" />\n";
event_log_ += "<line style=\"stroke-width:5000;stroke:purple;\" ";
event_log_ += "x1=\""+to_str(edge->prev()->vertex0()->x())+"\" y1=\""+to_str(edge->prev()->vertex0()->y())+"\" ";
event_log_ += "x2=\""+to_str(edge->prev()->vertex1()->x())+"\" y2=\""+to_str(edge->prev()->vertex1()->y())+"\" />\n";
}

        }
      } while (edge != start_edge && edge->color() != 2);

      // After the first run, any further runs are following internal hole 
      // loops. Find the first edge of the first/next hole.
      start_edge = NULL;
      for (edge_iterator it = edges_.begin(); it != edges_.end(); ++it) {
        if (it->color() == 0
            && it->is_primary()
           ) {
          start_edge = &(*it);
          break;
        }
      }
    }

    // Restore color() == 0 for internal edges.
    for (edge_iterator it = edges_.begin(); it != edges_.end(); ++it) {
      if (it->color() == 2) {
        it->color(0);
      }
    }

    // add some missing feet

    event_log_ += "<g id=\"missingfeet\">\n";

    start_edge = NULL;
    for (edge_iterator it = edges_.begin(); it != edges_.end(); ++it) {
      if (it->color() == 0
          && it->is_primary()
         ) {
        start_edge = &(*it);
        break;
      }
    }
    while (start_edge != NULL) {
      edge_type * edge = start_edge;
      do {
        if (!edge->foot()
           && edge->next()->foot()
           && edge->color() == 0
           && edge->is_primary()
           ) {
           
          if (edge->cell()->contains_point()) {
          edge->foot(edge->next()->foot()->x(), edge->next()->foot()->y());
          event_log_ += "<circle cx=\""+to_str(edge->next()->foot()->x())+"\" cy=\""+to_str(edge->next()->foot()->y())+"\" r=\"110000\" class=\"evtfoot3\"/>\n";
          }
          else {
            double x  = edge->vertex0()->x();
            double y  = edge->vertex0()->y();
            double x0 = edge->next()->foot()->x();
            double y0 = edge->next()->foot()->y();
            double x1 = edge->next()->vertex0()->x();
            double y1 = edge->next()->vertex0()->y();
            rotate_2d(x1, y1, 3.14159/2, x0, y0);
            makefoot(x, y, x0, y0, x1, y1);
            edge->foot(x, y);
            event_log_ += "<circle cx=\""+to_str(x)+"\" cy=\""+to_str(y)+"\" r=\"110000\" class=\"evtfoot3\"/>\n";
          }
          //printf("fixed missing foot for consecutive linear edges\n");
          if (! edge->prev()->foot()
              && ( !edge->prev()->is_curved()
                 || edge->prev()->cell()->contains_point() )
              && (edge->prev()->color() == 0 || edge->prev()->color() == 2)
              && edge->prev()->is_primary()) {
            edge = edge->prev();
            edge->color(0);
            edge = edge->prev();
            edge->color(0);
          }
        }
        // mark visited internal edges (will restore to 0 afterward)
        edge->color(2);
        edge = edge->next();
      } while (edge != start_edge && edge->color() != 2);
      
      // After the first run, any further runs are following internal hole 
      // loops. Find the first edge of the first/next hole.
      start_edge = NULL;
      for (edge_iterator it = edges_.begin(); it != edges_.end(); ++it) {
        if (it->color() == 0
            && it->is_primary()
           ) {
          start_edge = &(*it);
          break;
        }
      }
    }
    
    event_log_ += "</g>\n";

    // Restore color() == 0 for internal edges.
    for (edge_iterator it = edges_.begin(); it != edges_.end(); ++it) {
      if (it->color() == 2) {
        it->color(0);
      }
    }

    // check for any missing feet
    for (edge_iterator it = edges_.begin(); it != edges_.end(); ++it) {
      if (it->color() == 0
          && it->is_primary()
         ) {
        if (!it->foot()) {
          //croak("\nNO FOOT\n\n");          
          printf("NO FOOT\n");
          event_log_ += "<line class=\"edge_missing_foot\" ";
          event_log_ += " x1=\""+to_str(it->vertex0()->x())+"\" y1=\""+to_str(it->vertex0()->y())+"\"";
          event_log_ += " x2=\""+to_str(it->vertex1()->x())+"\" y2=\""+to_str(it->vertex1()->y())+"\"/>\n";

          // For debugging, put in placeholders for missing feet 
          // in some rediculous location, so we can still get some kind of 
          // output without a segfault.
          it->foot(0, 0);
          }
        if (it->next() &&  !it->next()->next()) {printf("NO NEXT\n\n");}
      }
    } 

    
    // Debug reporting
    /*
    if (0) {
      printf("original edges\n");
      printf("srcInd isInf curved   color");
      printf("  this     twin       next       prev        point\n");
      for (edge_iterator it = edges_.begin(); it != edges_.end(); ++it) {
        printf("%3d   %5s  %7s  %2d  ",
          it->cell()->source_index(),
          (it->is_finite() ? "     ":" INF "),
          (it->is_curved() ? " curve ":" line  "),
          it->color()
        );
        printf("%llu, %llu , %llu, %llu ",
          (unsigned long long int) &(*it),
          (unsigned long long int) it->twin(),
          (unsigned long long int) it->next(),
          (unsigned long long int) it->prev()
        );
      if (it->vertex0()) {
        printf("[%f , %f , %f]",
          it->vertex0()->x(),
          it->vertex0()->y(),
          it->vertex0()->r()
        );
      }
      else {printf("no vert0");}
      printf("\n");
      }
    }
    */
  }

 private:
  typedef typename cell_container_type::iterator cell_iterator;
  typedef typename vertex_container_type::iterator vertex_iterator;
  typedef typename edge_container_type::iterator edge_iterator;
  typedef typename TRAITS::vertex_equality_predicate_type
    vertex_equality_predicate_type;

  template <typename SEvent>
  bool is_primary_edge(const SEvent& site1, const SEvent& site2) const {
    bool flag1 = site1.is_segment();
    bool flag2 = site2.is_segment();
    if (flag1 && !flag2) {
      return (site1.point0() != site2.point0()) &&
             (site1.point1() != site2.point0());
    }
    if (!flag1 && flag2) {
      return (site2.point0() != site1.point0()) &&
             (site2.point1() != site1.point0());
    }
    return true;
  }

  template <typename SEvent>
  bool is_linear_edge(const SEvent& site1, const SEvent& site2) const {
    if (!is_primary_edge(site1, site2)) {
      return true;
    }
    return !(site1.is_segment() ^ site2.is_segment());
  }

  // Remove degenerate edge.
  void remove_edge(edge_type* edge) {
    
    // Are these two ifs necessary?
    // Put these in for debugging, where the problem was something else,
    // but these do fill in/transfer some missing feet.
    // After revising the foot-finding (trying to do it all in the sweepline
    // event processing), see if these are still needed.
    if (edge->foot() && !edge->next()->foot()) {
      edge->next()->foot(edge->foot()->x(), edge->foot()->y());
    }
    if (edge->twin()->foot() && !edge->twin()->next()->foot()) {
      edge->twin()->next()->foot(edge->twin()->foot()->x(), edge->twin()->foot()->y());
    }

    // Update the endpoints of the incident edges to the second vertex.
    vertex_type* vertex = edge->vertex0();
    edge_type* updated_edge = edge->twin()->rot_next();
    while (updated_edge != edge->twin()) {
      updated_edge->vertex0(vertex);
      updated_edge = updated_edge->rot_next();
    }
    
    edge_type* edge1 = edge;
    edge_type* edge2 = edge->twin();

    edge_type* edge1_rot_prev = edge1->rot_prev();
    edge_type* edge1_rot_next = edge1->rot_next();

    edge_type* edge2_rot_prev = edge2->rot_prev();
    edge_type* edge2_rot_next = edge2->rot_next();

    // Update prev/next pointers for the incident edges.
    edge1_rot_next->twin()->next(edge2_rot_prev);
    edge2_rot_prev->prev(edge1_rot_next->twin());
    edge1_rot_prev->prev(edge2_rot_next->twin());
    edge2_rot_next->twin()->next(edge1_rot_prev);

  }

  void mark_exterior(edge_type* edge) {
    if (edge->color() == 1) {
      return;
    }
    edge->color(1);
    edge->twin()->color(1);
    edge->cell()->color(1);
    edge->twin()->cell()->color(1);
    vertex_type* v = edge->vertex1();
    if (!v) {v = edge->vertex0();}
    if (v == NULL || !edge->is_primary()) {
      return;
    }
    v->color(1);
    edge_type* e = v->incident_edge();
    do {
      mark_exterior(e);
      e = e->rot_next();
    } while (e != v->incident_edge());
  }

  void mark_interior(edge_type* edge, bool backward = false) {
    // This function seems to work as intended, though it's still
    // on probation. The conditionals in the do {} while(); might not all be
    // correct or necessary (or they might be). 
    edge->color(0);
    edge->twin()->color(0);
    vertex_type* v = edge->vertex0();
    edge_type* e;
    if (edge->is_curved()) {
      edge_type* start_e = (edge->cell()->contains_point()) ? edge : edge->twin();
      e = start_e;
      do {
        if (!e->is_primary()) { 
          e->color(0);
          e->twin()->color(0);
        }
        e = e->next();
      } while (e != start_e);
    }

    if (!backward) {
      v = edge->vertex1();
    } 

    if (!v) {
      return;
    }
    e = v->incident_edge();
    v->color(0);
    e = v->incident_edge();
    do {
      if (e->is_primary() && e->next()->is_primary()) {
        if (e->color() == 0 && e->next()->color() != 0) {
          mark_interior(e->next());
        }
        if (e->color() != 0 && e->next()->color() == 0) {
          mark_interior(e, true);
        }
      }
      if (e->is_primary() && e->prev()->is_primary()) {
        if (e->color() == 0 && e->prev()->color() != 0) {
          mark_interior(e->prev(), true);
        }
        if (e->color() != 0 && e->prev()->color() == 0) {
          mark_interior(e);
        }
      }
      e = e->rot_next();
    } while (e != v->incident_edge());
  }

  bool is_exterior (const edge_type& e) { return (e->color() != 0); }

  void rotate_2d(double &x, double &y, const double theta, const double xo = 0, const double yo = 0) {
    double xp;
    x -= xo;
    y -= yo;
    xp = (x * cos(theta) - y * sin(theta)) + xo;
    y  = (y * cos(theta) + x * sin(theta)) + yo;
    x  = xp;
  }
  template <typename CT>
  void reflect(CT &x, CT &y, const CT x0, const CT y0, const CT x1, const CT y1) {
    double dy = (double) (y1 - y0);
    double dx = (double) (x1 - x0);
    if (dy == 0 && dx == 0) {return;}
    double theta = atan2(dy, dx);
    rotate_2d(x, y, -theta, x0, y0);
    y -= y0;
    y *= -1.0;
    y += y0;
    rotate_2d(x, y, theta, x0, y0);
  }

void makefoot(double & x, double & y, const double x0, const double y0,
                                      const double x1, const double y1) {
    // infinite slope case first
    if (x1 - x0 == 0) {
        x = x0;
    } else {
      double m  = (y1 - y0)/(x1 - x0);
      if (m == 0) {
          y = y0;
      }
      else {
        double intersect_x = ((m * x0) - y0 + ((1 / m) * x) + (y)) / (m + (1 / m));
          double intersect_y = -(x0 - intersect_x) * m + y0;
            x = intersect_x;
            y = intersect_y;
      }
    }
  }

  cell_container_type cells_;
  vertex_container_type vertices_;
  edge_container_type edges_;
  std::string event_log_;
  vertex_equality_predicate_type vertex_equality_predicate_;

  // Disallow copy constructor and operator=
  medial_axis(const medial_axis&);
  void operator=(const medial_axis&);
};
}  // polygon
}  // boost

#endif  // BOOST_POLYGON_MEDIAL_AXIS