#include "drawfigs.h"
#include <ioprotos.h>
#include <quad.h>
#include <delaunay.h>
#include <js.h>
#include <ps.h>

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>

/* GENERIC PROTOTYPES */

typedef site_struct Point;

typedef struct {
    int n_edges;
    int max_edges;
    edge_ref *edges;
  } Edge_list;

void make_sites (int n, int normal, float fraction, Point sites[]);
void make_12_sites (int n, Point sites[]);
void plot_sites(FILE* f, float radius, int n, Point sites[]);
float ccw(Point *a, Point *b, Point *c);
Point circ_center(Point *a, Point *b, Point *c);
Point perp_dir(Point *a, Point *b, float rbig);
double dist(Point *p, Point *q);

Edge_list *collect_edges(edge_ref e, int max_edges);
void collect_edge(edge_ref e, void *cls);

/* GENERIC PROCEDURES */

void make_sites (int n, int normal, float fraction, Point sites[])
  { int i;
    float x, y, r, rmax = 0.0;

    srandom(4615 + 4634);
    if (normal)
      { for (i=0; i < n; i++) 
          { sites[i].x = nrandom();
            sites[i].y = nrandom();
          }
      }
    else
      { for (i=0 ; i < n; i++)
          { sites[i].x = frandom() - 0.5; 
            sites[i].y = frandom() - 0.5;
          }
      }

    for (i=0; i < n; i++) 
      { x = sites[i].x;
        y = sites[i].y;
        fprintf(stderr, "s[%d] = (%f %f)\n", i, x, y);
        if (x < 0.0) x = -x;
        if (y < 0.0) y = -y;
        r = (x > y ? x : y);
        if (r > rmax) rmax = r;
      }

    rmax /= fraction;

    for (i=0; i < n; i++) 
      { sites[i].x /= rmax;
        sites[i].y /= rmax;
      }
  }
  
void make_12_sites (int n, Point sites[])
  {
    sites[ 0] = (Point){-0.750, +0.575};
    sites[ 1] = (Point){-0.050, +0.525};
    sites[ 2] = (Point){+0.700, +0.750};
    sites[ 3] = (Point){-0.650, -0.075};
    sites[ 4] = (Point){+0.325, +0.000};
    sites[ 5] = (Point){+0.800, +0.175};
    sites[ 6] = (Point){-0.400, -0.450};
    sites[ 7] = (Point){+0.575, -0.425};
    sites[ 8] = (Point){+0.225, -0.725};
    sites[ 9] = (Point){-0.325, +0.350};
    sites[10] = (Point){+0.500, +0.510};
    sites[11] = (Point){-0.250, -0.075};
    assert(n == 12, "make_sites: bad n");
  }

void plot_sites(FILE* f, float radius, int n, Point sites[])
  { int i;
    for (i=0; i<n; i++)
      { ps_fill_circle(f, sites[i].x, sites[i].y, radius, 0.0); }
  }

double dist(Point *p, Point *q)
 { double dx = ((double) p->x) - q->x;
   double dy = ((double) p->y) - q->y;
   double d = sqrt(dx*dx + dy*dy);
   return (d);
 }

/* Counterclockwise triangle predicate: */

/* Defined in delaunay.c

float ccw(Point *a, Point *b, Point *c)
  {
    double x1 = a->x, y1 = a->y;
    double x2 = b->x, y2 = b->y;
    double x3 = c->x, y3 = c->y;

    return (x2*y3-y2*x3) - (x1*y3-y1*x3) + (x1*y2-y1*x2);
  }

*/

/* Circumcenter: */

Point circ_center(Point *a, Point *b, Point *c)
  {
    Point ctr;

    double px[3], py[3];
    double dx[3], dy[3], mx[3], my[3], q[3];
    double det[3], rx[3], ry[3];
    double swx, swy, sw;
    int i, j;
    
    px[0] = ((double) a->x); py[0] = ((double) a->y); 
    px[1] = ((double) b->x); py[1] = ((double) b->y); 
    px[2] = ((double) c->x); py[2] = ((double) c->y); 
    
    for (i=0; i<3; i++)
      { j = (i==2 ? 0 : i+1);
        dx[i] = px[j] - px[i];
        dy[i] = py[j] - py[i];
        mx[i] = (px[j] + px[i])/2.0;
        my[i] = (py[j] + py[i])/2.0;
        q[i] = dx[i]*mx[i] + dy[i]*my[i];
      }
      
    swx = 0.0; swy = 0.0; sw = 0.0;

    for (i=0; i<3; i++)
      { j = (i==2 ? 0 : i+1);
        det[i] = dx[i]*dy[j] - dx[j]*dy[i];
        rx[i] = q[i]*dy[j] - q[j]*dy[i];
        ry[i] = dx[i]*q[j] - dx[j]*q[i];
        swx += det[i] * rx[i];
        swy += det[i] * ry[i];
        sw += det[i] * det[i];
        /*
          fprintf(stderr, 
            "circ_center: xc = %f yc = %f det = %f\n", 
            rx[i]/det[i], ry[i]/det[i], det[i]
          );
        */
      }
      
    fprintf(stderr, "\n");
        
    ctr.x = swx/sw;
    ctr.y = swy/sw;
    return (ctr);
  }

/* Perpendicular direction (point at infinity) */

Point perp_dir(Point *a, Point *b, float rbig)
  {
    Point prp;
    
    double ax = (double) a->x;
    double ay = (double) a->y;
    double bx = (double) b->x;
    double by = (double) b->y;
    
    double dx = bx - ax;
    double dy = by - ay;
    
    double d = sqrt(dx*dx + dy*dy);
    
    double mx = (ax + bx)/2.0;
    double my = (ay + by)/2.0;
    
    prp.x = mx - (dy/d)*rbig;
    prp.y = my + (dx/d)*rbig;

    return (prp);
  }

/* Collecting lists of edges: */

Edge_list *collect_edges(edge_ref e, int max_edges)
  { Edge_list *el = (Edge_list *) malloc(sizeof(Edge_list));
    
    el->n_edges = 0;
    el->max_edges = max_edges;
    el->edges = (edge_ref *) malloc(max_edges * sizeof(edge_ref));
    
    quad_enum(e, collect_edge, (void *) el);
    return(el);
    }

void collect_edge(edge_ref e, void *cls)
  { Edge_list *el = (Edge_list *) cls;
    assert ((el->n_edges < el->max_edges), "collect_edges: array overflow");
    (el->edges)[el->n_edges] = e;
    (el->n_edges)++;
  }

/* NEAREST POINT */

void nearest_draw (FILE *f);
void nearest_make_sites (int n, Point sites[]);

#define NEAREST_NSITES 12

void nearest_make_sites (int n, Point sites[])
  { make_12_sites(n, sites); }
 
void nearest_draw (FILE *f)
  {
    Point sites[NEAREST_NSITES];
    Point p;
    int i, ibest;
    double dbest;

    nearest_make_sites (NEAREST_NSITES, sites);

    plot_sites(f, 0.02, NEAREST_NSITES, sites);

    p = (Point) {+0.20, +0.25};
    
    dbest = 1.0e20;
    for (i=0; i<NEAREST_NSITES; i++)
      { double d = dist(&p, &(sites[i]));
        if (d < dbest) { dbest = d; ibest = i;}
      }
      
    ps_set_pen(f, 0.0, 0.0, 1.0, 1.0);
    ps_draw_segment(f, p.x, p.y, sites[ibest].x, sites[ibest].y);
    ps_set_pen(f, 0.0, 0.0, 0.0, 0.0);

    ps_fill_and_draw_circle(f, p.x, p.y, 0.02, 0.5);
  }

/* VORONOI DIAGRAM */

void voronoi_draw (FILE *f);
void voronoi_make_sites (int n, Point sites[]);
void voronoi_plot_diagram (FILE *f, edge_ref e);
void voronoi_draw_edge (edge_ref e, void *closure);

#define VORONOI_NSITES 12

void voronoi_make_sites (int n, Point sites[])
  { make_12_sites(n, sites); }
 
void voronoi_draw_edge (edge_ref e, void *closure)
  {
    FILE *f = (FILE *)closure;
    edge_ref d = TOR(e);
    Point *op = ORG(d);
    Point *dp = DEST(d);
    Point *lp = DEST(ONEXT(d));
    Point *rp = DEST(ONEXT(SYM(d)));
    Point ap, bp;
    
    if (ccw(op, dp, lp) > 0.0)
      { ap = circ_center(op, dp, lp); }
    else
      { ap = perp_dir(op, dp, 10.0); }
      
    if (ccw(dp, op, rp) > 0.0)
      { bp = circ_center(dp, op, rp); }
    else
      { bp = perp_dir(dp, op, 10.0); }
      
    ps_draw_segment(f, ap.x, ap.y, bp.x, bp.y);
  }

void voronoi_plot_diagram (FILE *f, edge_ref e)
  { 
    quad_enum (e, voronoi_draw_edge, (void *) f);
  }

void voronoi_draw (FILE *f)
  {
    Point sites[VORONOI_NSITES];
    edge_ref e;

    voronoi_make_sites (VORONOI_NSITES, sites);

    e = delaunay_build (sites, VORONOI_NSITES);

    printf("Delaunay returned %x:%d \n", (e & 0xfffffffcu), (e&3));

    ps_set_pen (f, 0.0, 0.0, 0.0, 0.0);
    voronoi_plot_diagram (f, ROT(e));

    plot_sites(f, 0.02, VORONOI_NSITES, sites);
  }

/* DELAUNAY DIAGRAM */

void delaunay_draw (FILE *f);
void delaunay_make_sites (int n, Point sites[]);
void delaunay_plot_diagram (FILE *f, edge_ref e);
void delaunay_draw_edge (edge_ref e, void *closure);
void delaunay_draw_edge_circle (edge_ref e, void *closure);
void delaunay_draw_face_circle (edge_ref e, void *closure);
void delaunay_trace_bubble (edge_ref e, void *closure);
edge_ref delaunay_pick_nice_edge(int ne, edge_ref e[]);

#define DELAUNAY_NSITES 12

void delaunay_make_sites (int n, Point sites[])
  { make_12_sites(n, sites); }

void delaunay_draw_edge (edge_ref e, void *closure)
  {
    FILE *f = (FILE *)closure;
    Point *op = ORG(e);
    Point *dp = DEST(e);
    ps_draw_segment(f, op->x, op->y, dp->x, dp->y);
  }

void delaunay_draw_edge_circle (edge_ref e, void *closure)
  {
    FILE *f = (FILE *)closure;
    Point *s = ORG(e);
    Point *t = DEST(e);
    Point *l = DEST(ONEXT(e));
    Point *r = DEST(ONEXT(SYM(e)));
    Point a, b, c1, c2, m, c;
    double rad;
    
    if (ccw(s, t, l) > 0.0)
      { a = circ_center(s, t, l); }
    else
      { a = perp_dir(s, t, 1.0); }
      
    if (ccw(t, s, r) > 0.0)
      { b = circ_center(t, s, r); }
    else
      { b = perp_dir(t, s, 1.0); }
      
    c1 = (Point) {0.66*a.x + 0.34*b.x, 0.66*a.y + 0.34*b.y};
    c2 = (Point) {0.34*a.x + 0.66*b.x, 0.34*a.y + 0.66*b.y};
    m =  (Point) {0.5 * s->x + 0.5 * t->x, 0.5 * s->y + 0.5 * t->y};
    
    c = (dist(&c1, &m) > dist(&c2, &m) ? c1 : c2);
      
    rad = dist(s, &c);
    ps_set_pen (f, 0.0, 0.0, 1.5, 1.5);
    ps_fill_and_draw_circle (f, c.x, c.y, rad, 0.75);
    ps_fill_circle(f, c.x, c.y, 0.01, 0.0);
  }

void delaunay_trace_bubble (edge_ref e, void *closure)
  {
    FILE *f = (FILE *)closure;
    Point *s = ORG(e);
    Point *t = DEST(e);
    Point *l = DEST(ONEXT(e));
    Point *r = DEST(ONEXT(SYM(e)));
    Point a, b, c;
    double rad;
    double dstep = 0.05;
    int i;
    int nsteps;
    
    if (ccw(s, t, l) > 0.0)
      { a = circ_center(s, t, l); }
    else
      { a = perp_dir(s, t, 0.2); }
      
    if (ccw(t, s, r) > 0.0)
      { b = circ_center(t, s, r); }
    else
      { b = perp_dir(t, s, 0.2); }
      
    { double d = dist(&a, &b);
      nsteps  = (int) (d/dstep);
      if (nsteps < 1) nsteps = 1;
    }
      
    for (i=0; i<nsteps; i++)
      { double alfa = (((double) i)  + 0.5)/((double) nsteps);
        double beta = 1.0 - alfa;
        c = (Point) {alfa*a.x + beta*b.x, alfa*a.y + beta*b.y};
        rad = dist(s, &c);
        ps_set_pen (f, 0.0, 0.0, 1.0, 2.0);
        ps_draw_circle (f, c.x, c.y, rad);
      }
  }

void delaunay_draw_face_circle (edge_ref e, void *closure)
  {
    FILE *f = (FILE *)closure;
    Point *s = ORG(e);
    Point *t = DEST(e);
    Point *l = DEST(ONEXT(e));
    Point c;
    double rad;
    
    if (ccw(s, t, l) > 0.0)
      { c = circ_center(s, t, l); }
    else
      { assert(0, "draw_face_circle: bad face"); }
      
    rad = dist(s, &c);
    ps_set_pen (f, 0.0, 0.0, 1.5, 1.5);
    ps_fill_and_draw_circle (f, c.x, c.y, rad, 0.75);
    ps_fill_circle(f, c.x, c.y, 0.01, 0.0);
  }

void delaunay_plot_diagram (FILE *f, edge_ref e)
  { 
    quad_enum (e, delaunay_draw_edge, (void *) f);
  }

int delaunay_variant; /* Ugly hack: to get three pics out of delaunay_draw */

edge_ref delaunay_pick_nice_edge(int ne, edge_ref e[])
  /* Return an edge suitable for showing the empty circle property: */
  {
    double dbest = 1.0e20;
    edge_ref ebest = e[0]; /* Just in case */
    Point o = {0.0, 0.0};
    double hmin = 0.20;
    int i;

    for (i=0; i<ne; i++)
      { edge_ref g = e[i];
        Point *s = ORG(g);
        Point *t = DEST(g);
        Point m = {(s->x + t->x)/2.0, (s->y + t->y)/2.0};
        double d = dist(&m, &o);
        double h = dist(s, t);
        if ((d < dbest) && (h > hmin)) { dbest = d; ebest = g; }
      }

    return (ebest);
  }

void delaunay_draw (FILE *f)
  {
    Point sites[DELAUNAY_NSITES];
    int max_edges = 3*(DELAUNAY_NSITES-2);
    edge_ref e;

    delaunay_make_sites (DELAUNAY_NSITES, sites);

    e = delaunay_build (sites, DELAUNAY_NSITES);

    printf("Delaunay returned %x:%d \n", (e & 0xfffffffcu), (e&3));

    if (delaunay_variant == 0)
      { /* Delaunay with Voronoi */
        ps_set_pen(f, 0.0, 0.0, 1.0, 2.0);
        voronoi_plot_diagram (f, ROT(e));
      }
    else if (delaunay_variant == 1)
      { /* Delaunay with an edge circle */
        Edge_list *el = collect_edges(e, max_edges);
        edge_ref g = delaunay_pick_nice_edge(el->n_edges, el->edges);
        delaunay_draw_edge_circle (g, (void *) f);
        ps_set_pen(f, 0.0, 0.0, 1.0, 2.0);
        voronoi_plot_diagram (f, ROT(e));
      }
    else if (delaunay_variant == 2)
      { /* Delaunay with a face circle */
        Edge_list *el = collect_edges(e, max_edges);
        edge_ref g = delaunay_pick_nice_edge(el->n_edges, el->edges);
        g = SYM(OPREV(SYM(g)));
        delaunay_draw_face_circle (g, (void *) f); 
        ps_set_pen(f, 0.0, 0.0, 1.0, 2.0);
        voronoi_plot_diagram (f, ROT(e));
      }
    else if (delaunay_variant == 3)
      { /* Pure delaunay diagram */ }
    else if (delaunay_variant == 4)
      { /* Recursive Delaunay with bubble */ 
        Edge_list *el = collect_edges(e, max_edges);
        double xmedian = 0.0; /* Should compute it... */
        double ycut = 0.0;
        Point cutpoint = (Point){xmedian, ycut};
        int i;
        
        /* Draw cut line: */
        ps_set_pen(f, 0.0, 0.0, 2.0, 1.0);
        ps_draw_segment(f, xmedian, -1.0, xmedian, +1.0);
        
        /* Delete edges that cross xmedian and are above ycut: */
        ps_set_pen(f, 0.0, 0.5, 0.0, 0.0);
        for (i=0; i<el->n_edges; i++)
          { edge_ref g = el->edges[i];
            Point *s = ORG(g);
            Point *t = DEST(g);
            if (s->x > t->x) { Point *aux = s; s = t; t = aux;}
            if ((s->x < xmedian) && (t->x >= xmedian))
              { if (ccw(s, t, &cutpoint) < 0.0)
                  { destroy_edge(g); }
                else
                  { delaunay_draw_edge(g, (void *) f);
                    delaunay_trace_bubble(g, (void*) f); }
              }
          }
      }
    else
      { assert (0, "delaunay_variant is bogus"); }

    ps_set_pen (f, 0.0, 0.0, 0.0, 0.0);
    delaunay_plot_diagram (f, e);

    plot_sites(f, 0.02, DELAUNAY_NSITES, sites);

  }

/* GABRIEL DIAGRAM */

void gabriel_draw (FILE *f);
void gabriel_make_sites (int n, Point sites[]);
edge_ref gabriel_from_delaunay (edge_ref e, int n, Point sites[]);
void gabriel_plot_diagram (FILE *f, edge_ref e);
void gabriel_analyze_edge (edge_ref e, void *closure);
void gabriel_draw_edge (edge_ref e, void *closure);
void gabriel_draw_edge_circle (edge_ref e, void *closure);

#define GABRIEL_NSITES 12

void gabriel_make_sites (int n, Point sites[])
  { make_12_sites(n, sites); }

typedef struct{
    int n_gut_edges; /* Gabriel edge counter */
    int n_bad_edges; /* Delaunay edges that are not gabriel */
    int n; /* Number of sites */
    Point *sites;
    edge_ref *gut_edges;
    edge_ref *bad_edges;
  } gabriel_data;

void gabriel_analyze_edge (edge_ref e, void *closure)
  {
    gabriel_data *cls = (gabriel_data *) closure;
    Point *s = ORG(e);
    Point *t = DEST(e);
    Point m = {(s->x + t->x)/2.0, (s->y + t->y)/2.0};
    double r = dist(s, &m);
    
    int i;
    int nin = 0;
    
    for (i=0; i < (cls->n); i++)
      { Point *u = &((cls->sites)[i]);
        double d = dist(u, &m);
        if ((d < r) && (u != s) && (u != t)) nin++;
      }
    
    if (nin == 0)
      { cls->gut_edges[cls->n_gut_edges] = e;
        (cls->n_gut_edges)++;
      }
    else
      { cls->bad_edges[cls->n_bad_edges] = e;
        (cls->n_bad_edges)++;
      }
  }

edge_ref gabriel_from_delaunay (edge_ref e, int n, Point sites[])
  {
    int i;

    int max_edges = 3*(n-2);
    gabriel_data cls;
    cls.n = n;
    cls.sites = sites;
    cls.n_gut_edges = 0;
    cls.n_bad_edges = 0;
    cls.gut_edges = (edge_ref *) malloc(max_edges * sizeof(edge_ref));
    cls.bad_edges = (edge_ref *) malloc(max_edges * sizeof(edge_ref));
    
    quad_enum (e, gabriel_analyze_edge, (void *) &cls);

    assert(cls.n_gut_edges > 0, "gabriel_from_delaunay: no edges found");
    
    /* Delete non-edges: */
    for (i=0; i < cls.n_bad_edges; i++)
      { edge_ref g = cls.bad_edges[i];
        destroy_edge(g);
      }
      
    /* Return some edge: */
    { int ne = cls.n_gut_edges;
      double dbest = 1.0e20;
      edge_ref ebest = cls.gut_edges[0]; /* Just in case */
      Point o = {0.0, 0.0};
      double hmin = 0.20;
      
      for (i=0; i<ne; i++)
        { edge_ref g = cls.gut_edges[i];
          Point *s = ORG(g);
          Point *t = DEST(g);
          Point m = {(s->x + t->x)/2.0, (s->y + t->y)/2.0};
          double d = dist(&m, &o);
          double h = dist(s, t);
          if ((d < dbest) && (h > hmin)) { dbest = d; ebest = g; }
        }
      
      return (ebest);
    }
  }

void gabriel_draw_edge (edge_ref e, void *closure)
  {
    FILE *f = (FILE *)closure;
    Point *op = ORG(e);
    Point *dp = DEST(e);
    ps_draw_segment(f, op->x, op->y, dp->x, dp->y);
  }

void gabriel_draw_edge_circle (edge_ref e, void *closure)
  {
    FILE *f = (FILE *)closure;
    Point *s = ORG(e);
    Point *t = DEST(e);
    Point m = {(s->x + t->x)/2.0, (s->y + t->y)/2.0};
    double r = dist(s, &m);
    ps_set_pen (f, 0.0, 0.0, 1.5, 1.5);
    ps_fill_and_draw_circle (f, m.x, m.y, r, 0.75);
  }

void gabriel_plot_diagram (FILE *f, edge_ref e)
  { 
    quad_enum (e, gabriel_draw_edge, (void *) f);
  }

void gabriel_draw (FILE *f)
  {
    Point sites[GABRIEL_NSITES];
    edge_ref e;
    
    gabriel_make_sites (GABRIEL_NSITES, sites);

    e = delaunay_build (sites, GABRIEL_NSITES);

    printf("Delaunay returned %x:%d \n", (e & 0xfffffffcu), (e&3));

    e = gabriel_from_delaunay (e, GABRIEL_NSITES, sites);

    printf("Gabriel returned %x:%d \n", (e & 0xfffffffcu), (e&3));

    gabriel_draw_edge_circle (e, (void *) f);
    
    ps_set_pen (f, 0.0, 0.0, 0.0, 0.0);
    gabriel_plot_diagram (f, e);

    plot_sites(f, 0.02, GABRIEL_NSITES, sites);

  }

/* RELATIVE NEIGHBORHOOD GRAPH */

void relnei_draw (FILE *f);
void relnei_make_sites (int n, Point sites[]);
edge_ref relnei_from_delaunay (edge_ref e, int n, Point sites[]);
void relnei_plot_diagram (FILE *f, edge_ref e);
void relnei_analyze_edge (edge_ref e, void *closure);
void relnei_draw_edge (edge_ref e, void *closure);
void relnei_draw_edge_lune (edge_ref e, void *closure);

#define RELNEI_NSITES 12

void relnei_make_sites (int n, Point sites[])
  { make_12_sites(n, sites); }

typedef struct{
    int n_gut_edges; /* Relnei edge counter */
    int n_bad_edges; /* Delaunay edges that are not relnei */
    int n; /* Number of sites */
    Point *sites;
    edge_ref *gut_edges;
    edge_ref *bad_edges;
  } relnei_data;

void relnei_analyze_edge (edge_ref e, void *closure)
  {
    relnei_data *cls = (relnei_data *) closure;
    Point *s = ORG(e);
    Point *t = DEST(e);
    double r = dist(s, t);
    
    int i;
    int nin = 0;
    
    for (i=0; i < (cls->n); i++)
      { Point *u = &((cls->sites)[i]);
        double ds = dist(u, s);
        double dt = dist(u, t);
        if ((u != s) && (u != t) && (ds < r) && (dt < r)) nin++;
      }
    
    if (nin == 0)
      { cls->gut_edges[cls->n_gut_edges] = e;
        (cls->n_gut_edges)++;
      }
    else
      { cls->bad_edges[cls->n_bad_edges] = e;
        (cls->n_bad_edges)++;
      }
  }

edge_ref relnei_from_delaunay (edge_ref e, int n, Point sites[])
  {
    int i;

    int max_edges = 3*(n-2);
    relnei_data cls;
    cls.n = n;
    cls.sites = sites;
    cls.n_gut_edges = 0;
    cls.n_bad_edges = 0;
    cls.gut_edges = (edge_ref *) malloc(max_edges * sizeof(edge_ref));
    cls.bad_edges = (edge_ref *) malloc(max_edges * sizeof(edge_ref));
    
    quad_enum (e, relnei_analyze_edge, (void *) &cls);

    assert(cls.n_gut_edges > 0, "relnei_from_delaunay: no edges found");
    
    /* Delete non-edges: */
    for (i=0; i < cls.n_bad_edges; i++)
      { edge_ref g = cls.bad_edges[i];
        destroy_edge(g);
      }
      
    /* Return some edge: */
    { int ne = cls.n_gut_edges;
      double dbest = 1.0e20;
      edge_ref ebest = cls.gut_edges[0]; /* Just in case */
      Point o = {0.0, 0.0};
      double hmin = 0.20;
      
      for (i=0; i<ne; i++)
        { edge_ref g = cls.gut_edges[i];
          Point *s = ORG(g);
          Point *t = DEST(g);
          Point m = {(s->x + t->x)/2.0, (s->y + t->y)/2.0};
          double d = dist(&m, &o);
          double h = dist(s, t);
          if ((d < dbest) && (h > hmin)) { dbest = d; ebest = g; }
        }
      
      return (ebest);
    }
  }

void relnei_draw_edge (edge_ref e, void *closure)
  {
    FILE *f = (FILE *)closure;
    Point *op = ORG(e);
    Point *dp = DEST(e);
    ps_draw_segment(f, op->x, op->y, dp->x, dp->y);
  }

void relnei_draw_edge_lune (edge_ref e, void *closure)
  {
    FILE *f = (FILE *)closure;
    Point *s = ORG(e);
    Point *t = DEST(e);
    Point m = {(s->x + t->x)/2.0, (s->y + t->y)/2.0};
    double r = dist(s, &m);
    double tilt = atan2(s->y - m.y, s->x - m.x);

    ps_set_pen (f, 0.0, 0.0, 1.5, 1.5);
    ps_fill_and_draw_lune (f, m.x, m.y, r, tilt, 0.75);
  }

void relnei_plot_diagram (FILE *f, edge_ref e)
  { 
    quad_enum (e, relnei_draw_edge, (void *) f);
  }

void relnei_draw (FILE *f)
  {
    Point sites[RELNEI_NSITES];
    edge_ref e;
    
    relnei_make_sites (RELNEI_NSITES, sites);

    e = delaunay_build (sites, RELNEI_NSITES);

    printf("Delaunay returned %x:%d \n", (e & 0xfffffffcu), (e&3));

    e = relnei_from_delaunay (e, RELNEI_NSITES, sites);

    printf("Relnei returned %x:%d \n", (e & 0xfffffffcu), (e&3));

    relnei_draw_edge_lune (e, (void *) f);
    
    ps_set_pen (f, 0.0, 0.0, 0.0, 0.0);
    relnei_plot_diagram (f, e);

    plot_sites(f, 0.02, RELNEI_NSITES, sites);

  }

/* MINIMUM SPANNING TREE */

void mintree_draw (FILE *f);
void mintree_make_sites (int n, Point sites[]);
edge_ref mintree_from_delaunay (edge_ref e, int n, Point sites[]);
void mintree_plot_diagram (FILE *f, edge_ref e);
void mintree_save_edge (edge_ref e, void *closure);
void mintree_sort_edges (int ne, edge_ref e[]);
int mintree_find_root(Point *s, int n, Point sites[], int vp[]);
void mintree_draw_edge (edge_ref e, void *closure);

#define MINTREE_NSITES 12

void mintree_make_sites (int n, Point sites[])
  { make_12_sites(n, sites); }

typedef struct{
    int n_del_edges; /* Delaunay edges */
    edge_ref *del_edges;  /* All Delaunay edges */
  } mintree_data;

void mintree_save_edge (edge_ref e, void *closure)
  {
    mintree_data *cls = (mintree_data *) closure;
    cls->del_edges[cls->n_del_edges] = e;
    (cls->n_del_edges)++;
  }

void mintree_sort_edges (int ne, edge_ref e[])
  /* Sort e[] by increasing length */
  { int i, j;
    double *r = (double *) malloc(ne * sizeof(double));
    for (i=1; i<ne; i++)
      { edge_ref ei = e[i]; 
        Point *s = ORG(ei);
        Point *t = DEST(ei);
        double ri = dist(s, t);
        for (j=i; (j>0) && (r[j-1] > ri); j--)
          { r[j] = r[j-1]; e[j] = e[j-1]; }
        r[j] = ri; e[j] = ei;
      }
  }

int mintree_find_root(Point *s, int n, Point sites[], int vp[])
  { int i;
    for (i=0; (i<n) && (s != &(sites[i])); i++) { }
    assert (i<n, "find_root: didn\'t find it");
    while (vp[i] != i) i = vp[i];
    return (i);
  }

edge_ref mintree_from_delaunay (edge_ref e, int n, Point sites[])
  {
    int max_edges = 3*(n-2);

    int n_gut_edges = 0;
    edge_ref *gut_edges = (edge_ref *) malloc(max_edges * sizeof(edge_ref));
    int n_bad_edges = 0;
    edge_ref *bad_edges = (edge_ref *) malloc(max_edges * sizeof(edge_ref));

    mintree_data cls;
    cls.n_del_edges = 0;
    cls.del_edges = (edge_ref *) malloc(max_edges * sizeof(edge_ref));
    
    quad_enum (e, mintree_save_edge, (void *) &cls);
    
    /* Kruskal's algorithm (or was it Prim's?) */
    { 
      int *vp = (int *) malloc(n * sizeof(int));
      /* sites[vp[i]] is the parent of sites[i] in the union-find tree */
    
      int iv, ie;

      for (iv=0; iv<n; iv++) vp[iv] = iv;

      mintree_sort_edges (cls.n_del_edges, cls.del_edges);

      for (ie=0; ie<cls.n_del_edges; ie++)
        { edge_ref ei = (cls.del_edges)[ie];
          int is = mintree_find_root(ORG(ei), n, sites, vp);
          int it = mintree_find_root(DEST(ei), n, sites, vp);
          if (is != it)
            { /* Edge ei does not close a cycle */
              gut_edges[n_gut_edges] = ei;
              n_gut_edges++;
              if (is < it) vp[it] = is; else vp[is] = it;
            }
          else
            { /* Edge ei closes a cycle */
              bad_edges[n_bad_edges] = ei;
              n_bad_edges++;
            }
        }
    }
    
    assert(n_gut_edges > 0, "mintree_from_delaunay: no edges found");
    
    /* Delete bad edges: */    
    { int ie;
      for (ie=0; ie < n_bad_edges; ie++)
        { edge_ref g = bad_edges[ie];
          destroy_edge(g);
        }
    }
      
    /* Return some edge: */
    { int ie;
      int ne = n_gut_edges;
      double dbest = 1.0e20;
      edge_ref ebest = gut_edges[0]; /* Just in case */
      Point o = {0.0, 0.0};
      double hmin = 0.20;
      
      for (ie=0; ie<ne; ie++)
        { edge_ref g = gut_edges[ie];
          Point *s = ORG(g);
          Point *t = DEST(g);
          Point m = {(s->x + t->x)/2.0, (s->y + t->y)/2.0};
          double d = dist(&m, &o);
          double h = dist(s, t);
          if ((d < dbest) && (h > hmin)) { dbest = d; ebest = g; }
        }
      
      return (ebest);
    }
  }

void mintree_draw_edge (edge_ref e, void *closure)
  {
    FILE *f = (FILE *)closure;
    Point *op = ORG(e);
    Point *dp = DEST(e);
    ps_draw_segment(f, op->x, op->y, dp->x, dp->y);
  }

void mintree_plot_diagram (FILE *f, edge_ref e)
  { 
    quad_enum (e, mintree_draw_edge, (void *) f);
  }

void mintree_draw (FILE *f)
  {
    Point sites[MINTREE_NSITES];
    edge_ref e;
    
    mintree_make_sites (MINTREE_NSITES, sites);

    e = delaunay_build (sites, MINTREE_NSITES);

    printf("Delaunay returned %x:%d \n", (e & 0xfffffffcu), (e&3));

    e = mintree_from_delaunay (e, MINTREE_NSITES, sites);

    printf("Mintree returned %x:%d \n", (e & 0xfffffffcu), (e&3));

    ps_set_pen (f, 0.0, 0.0, 0.0, 0.0);
    mintree_plot_diagram (f, e);

    plot_sites(f, 0.02, MINTREE_NSITES, sites);

  }

/* MAIN PROGRAM */

int main(void)
  {
    int eps_format;
    for (eps_format=0; eps_format < 2; eps_format++)
      { gen_figure("del-nearest",  -1.0, +1.0, -1.0, +1.0, nearest_draw,  eps_format); 
        gen_figure("del-voronoi",  -1.0, +1.0, -1.0, +1.0, voronoi_draw,  eps_format); 
        delaunay_variant = 0;
        gen_figure("del-delaunay", -1.0, +1.0, -1.0, +1.0, delaunay_draw, eps_format); 
        delaunay_variant = 1;
        gen_figure("del-delcirc",  -1.0, +1.0, -1.0, +1.0, delaunay_draw, eps_format); 
        delaunay_variant = 2;
        gen_figure("del-delface",  -1.0, +1.0, -1.0, +1.0, delaunay_draw, eps_format); 
        delaunay_variant = 3;
        gen_figure("del-delpure",  -1.0, +1.0, -1.0, +1.0, delaunay_draw, eps_format); 
        delaunay_variant = 4;
        gen_figure("del-delbubb",  -1.0, +1.0, -1.0, +1.0, delaunay_draw, eps_format); 
        gen_figure("del-gabriel",  -1.0, +1.0, -1.0, +1.0, gabriel_draw,  eps_format); 
        gen_figure("del-relnei",   -1.0, +1.0, -1.0, +1.0, relnei_draw,   eps_format); 
        gen_figure("del-mintree",  -1.0, +1.0, -1.0, +1.0, mintree_draw,  eps_format); 
      }
    return(0);
  }