/* See delvor.h */
/* Last edited on 2023-02-25 16:04:08 by stolfi */

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

#include <quad.h>
#include <delaunay.h>
#include <pswr.h>
#include <jsrandom.h>
#include <affirm.h>
#include <rn.h>

#include <drawfigs.h>

#include <delvor.h>

/* GENERIC PROTOTYPES */

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

void make_random_sites(int n, bool_t normal, delaunay_site_t sites[]);
void make_12_sites(int n, delaunay_site_t sites[]);

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

r2_t r2_from_hi2(hi2_point_t *ap);
hi2_point_t hi2_from_r2(r2_t *ap, int32_t W);
r2_t circumcenter(hi2_point_t *ap, hi2_point_t *bp, hi2_point_t *cp);

/* GENERIC PROCEDURES */

r2_t r2_from_hi2(hi2_point_t *ap)
  {
    int32_t W = ap->c[0], X = ap->c[1], Y = ap->c[2];
    r2_t rp;
    rp.c[0] = ((double)X)/((double)W);
    rp.c[1] = ((double)Y)/((double)W);
    return rp;
  }

hi2_point_t hi2_from_r2(r2_t *ap, int32_t W);
  {
    int X = (int)floor(W*ap->c[0] + 0.5);
    int Y = (int)floor(W*ap->c[1] + 0.5);
    return (hi2_point_t){{ W, X, Y }};
  }

r2_t circumcenter(hi2_point_t *ap, hi2_point_t *bp, hi2_point_t *cp)
  { 
    r2_t A = r2_from_hi2(ap);
    r2_t B = r2_from_hi2(bp);
    r2_t C = r2_from_hi2(cp);
    return r2_circumcenter(&A, &B, &C);
  }

void make_random_sites(int n, int normal, int32_t W, delaunay_site_t sites[])
  { int i;
        
    srandom(4615 + 4634);
    for (i=0; i < n; i++) 
      { sites[i].index = i;
        r2_t p;
        if (normal)
          { double sigma = 0.5/3;
            p.c[0] = fmax(-0.5, fmin(+0.5, sigma*dgaussrand()));
            p.c[0] = fmax(-0.5, fmin(+0.5, sigma*dgaussrand()));
          }
        else
          { p.c[1] = drandom() - 0.5; 
            p.c[1] = drandom() - 0.5;
          }
        sites[i].pt = hi2_from_r2(&p, W);
      }
  }
  
void make_12_sites (int n, delaunay_site_t sites[])
  {
    affirm(n == 12, "make_12_sites: bad n");
    int k = 0;
    sites[k] = (delaunay_site_t){ .index = k, .pt = (hi2_point_t){{1000, -750, +575}} }; k++;
    sites[k] = (delaunay_site_t){ .index = k, .pt = (hi2_point_t){{1000,  -50, +525}} }; k++;
    sites[k] = (delaunay_site_t){ .index = k, .pt = (hi2_point_t){{1000, +700, +750}} }; k++;
    sites[k] = (delaunay_site_t){ .index = k, .pt = (hi2_point_t){{1000, -650,  -75}} }; k++;
    sites[k] = (delaunay_site_t){ .index = k, .pt = (hi2_point_t){{1000, +325,    0}} }; k++;
    sites[k] = (delaunay_site_t){ .index = k, .pt = (hi2_point_t){{1000, +800, +175}} }; k++;
    sites[k] = (delaunay_site_t){ .index = k, .pt = (hi2_point_t){{1000, -400, -450}} }; k++;
    sites[k] = (delaunay_site_t){ .index = k, .pt = (hi2_point_t){{1000, +575, -425}} }; k++;
    sites[k] = (delaunay_site_t){ .index = k, .pt = (hi2_point_t){{1000, +225, -725}} }; k++;
    sites[k] = (delaunay_site_t){ .index = k, .pt = (hi2_point_t){{1000, -325, +350}} }; k++;
    sites[k] = (delaunay_site_t){ .index = k, .pt = (hi2_point_t){{1000, +500, +510}} }; k++;
    sites[k] = (delaunay_site_t){ .index = k, .pt = (hi2_point_t){{1000, -250,  -75}} }; k++;
    affirm(k == n, "make_12_sites: something amiss");
  }

/* Collecting lists of edges: */

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

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

/* NEAREST POINT */

void nearest_draw(PSStream *ps);
void nearest_make_sites(int n, delaunay_site_t sites[]);

#define NEAREST_NSITES 12

void nearest_make_sites(int n, delaunay_site_t sites[])
  { make_12_sites(n, sites); }
 
void nearest_draw(PSStream *ps)
  {
    delaunay_site_t sites[NEAREST_NSITES];
    int i, ibest = -1;
    double dbest;

    nearest_make_sites(NEAREST_NSITES, sites);

    plot_sites(ps, 0.5, NEAREST_NSITES, sites);

    hi2_point_t q = (hi2_point_t) {{100, +20, +25}};
    
    urat64_t d2best = urat64_INF;
    for (i=0; i<NEAREST_NSITES; i++)
      { urat64_t d2 = hi2_dist_sqr(&q, &(sites[i].pt));
        if (urat64_compare(d2, d2best) < 0) { d2best = d2; ibest = i;}
      }
      
    pswr_set_pen(ps, 0,0,0, 0.1, 1.0, 1.0);
    pswr_segment(ps, p.c[0], p.c[1], sites[ibest].p.c[0], sites[ibest].p.c[1]);
    pswr_set_pen(ps, 0,0,0, 0.1, 0.0, 0.0);
    

    pswr_set_fill_color(ps, 0.5,0.5,0.5);
    pswr_dot(ps, p.c[0], p.c[1], 0.5, TRUE, TRUE);
  }

/* VORONOI DIAGRAM */

void voronoi_draw(PSStream *ps);
void voronoi_make_sites(int n, delaunay_site_t sites[]);
void voronoi_plot_diagram(PSStream *ps, quad_arc_t e);
void voronoi_draw_edge(quad_arc_t e, void *closure);

#define VORONOI_NSITES 12

void voronoi_make_sites(int n, delaunay_site_t sites[])
  { make_12_sites(n, sites); }
 
void voronoi_draw_edge(quad_arc_t e, void *closure)
  {
    PSStream *ps = (PSStream *)closure;
    quad_arc_t d = quad_tor(e);
    delaunay_site_t *os = quad_org(d);
    delaunay_site_t *ds = quad_dst(d);
    delaunay_site_t *ls = quad_dst(quad_onext(d));
    delaunay_site_t *rs = quad_dst(quad_onext(quad_sym(d)));
    hi2_point_t *op = &(os->pt);
    hi2_point_t *dp = &(ds->pt);
    hi2_point_t *lp = &(ls->pt);
    hi2_point_t *rp = &(rs->pt);
    r2_point_t ap, bp;
    
    if (hi2_orient(op, dp, lp) > 0)
      { ap = circumcenter(op, dp, lp); }
    else
      { ap = perp_dir(op, dp, 10.0); }
      
    if (r2_orient(dp, op, rp) > 0)
      { bp = circumcenter(dp, op, rp); }
    else
      { bp = perp_dir(dp, op, 10.0); }
      
    pswr_segment(ps, ap.c[0], ap.c[1], bp.c[0], bp.c[1]);
  }

void voronoi_plot_diagram(PSStream *ps, quad_arc_t e)
  { 
    quad_enum(e, voronoi_draw_edge, (void *)ps);
  }

void voronoi_draw(PSStream *ps)
  {
    delaunay_site_t sites[VORONOI_NSITES];
    quad_arc_t e;

    voronoi_make_sites(VORONOI_NSITES, sites);

    e = delaunay_build(sites, VORONOI_NSITES);

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

    pswr_set_pen(ps, 0,0,0, 0.0, 0.0, 0.0);
    voronoi_plot_diagram(ps, quad_rot(e));

    plot_sites(ps, 0.5, VORONOI_NSITES, sites);
  }

/* DELAUNAY DIAGRAM */

void delaunay_draw(PSStream *ps);
void delaunay_make_sites(int n, delaunay_site_t sites[]);
void delaunay_plot_diagram(PSStream *ps, quad_arc_t e);
void delaunay_draw_edge(quad_arc_t e, void *closure);
void delaunay_draw_edge_circle(quad_arc_t e, void *closure);
void delaunay_draw_face_circle(quad_arc_t e, void *closure);
void delaunay_trace_bubble(quad_arc_t e, void *closure);
quad_arc_t delaunay_pick_nice_edge(int ne, quad_arc_t e[]);

#define DELAUNAY_NSITES 12

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

void delaunay_draw_edge(quad_arc_t e, void *closure)
  {
    PSStream *ps = (PSStream *)closure;
    delaunay_site_t *os = quad_org(e);
    delaunay_site_t *ds = quad_dst(e);
    Point *op = &(os->p);
    Point *dp = &(ds->p);
    pswr_segment(ps, op->c[0], op->c[1], dp->c[0], dp->c[1]);
  }

void delaunay_draw_edge_circle(quad_arc_t e, void *closure)
  {
    PSStream *ps = (PSStream *)closure;
    delaunay_site_t *ss = quad_org(e);
    delaunay_site_t *ts = quad_dst(e);
    delaunay_site_t *ls = quad_dst(quad_onext(e));
    delaunay_site_t *rs = quad_dst(quad_onext(quad_sym(e)));
    Point *sp = &(ss->p);
    Point *tp = &(ts->p);
    Point *lp = &(ls->p);
    Point *rp = &(rs->p);
    Point a, b, c1, c2, m, c;
    double rad;
    
    if (r2_orient(sp, tp, lp) > 0)
      { a = circumcenter(sp, tp, lp); }
    else
      { a = perp_dir(sp, tp, 1.0); }
      
    if (r2_orient(tp, sp, rp) > 0)
      { b = circumcenter(tp, sp, rp); }
    else
      { b = perp_dir(tp, sp, 1.0); }
      
    c1 = (Point) {{0.66*a.c[0] + 0.34*b.c[0], 0.66*a.c[1] + 0.34*b.c[1]}};
    c2 = (Point) {{0.34*a.c[0] + 0.66*b.c[0], 0.34*a.c[1] + 0.66*b.c[1]}};
    m =  (Point) {{0.5 * sp->c[0] + 0.5 * tp->c[0], 0.5 * sp->c[1] + 0.5 * tp->c[1]}};
    
    c = (r2_dist(&c1, &m) > r2_dist(&c2, &m) ? c1 : c2);
      
    rad = r2_dist(sp, &c);
    pswr_set_pen(ps, 0,0,0, 0.0, 1.5, 1.5);
    pswr_set_fill_color(ps, 0.75,0.75,0.75);
    pswr_circle(ps, c.c[0], c.c[1], rad, TRUE, TRUE);
    pswr_set_fill_color(ps, 0,0,0);
    pswr_dot(ps, c.c[0], c.c[1], 0.25, TRUE, FALSE);
  }

void delaunay_trace_bubble(quad_arc_t e, void *closure)
  {
    PSStream *ps = (PSStream *)closure;
    delaunay_site_t *ss = quad_org(e);
    delaunay_site_t *ts = quad_dst(e);
    delaunay_site_t *ls = quad_dst(quad_onext(e));
    delaunay_site_t *rs = quad_dst(quad_onext(quad_sym(e)));
    Point *sp = &(ss->p);
    Point *tp = &(ts->p);
    Point *lp = &(ls->p);
    Point *rp = &(rs->p);
    Point a, b, c;
    double rad;
    double dstep = 0.05;
    int i;
    int nsteps;
    
    if (r2_orient(sp, tp, lp) > 0)
      { a = circumcenter(sp, tp, lp); }
    else
      { a = perp_dir(sp, tp, 0.2); }
      
    if (r2_orient(tp, sp, rp) > 0)
      { b = circumcenter(tp, sp, rp); }
    else
      { b = perp_dir(tp, sp, 0.2); }
      
    { double d = r2_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.c[0] + beta*b.c[0], alfa*a.c[1] + beta*b.c[1]}};
        rad = r2_dist(sp, &c);
        pswr_set_pen(ps, 0,0,0, 0.0, 1.0, 2.0);
        pswr_circle(ps, c.c[0], c.c[1], rad, TRUE, TRUE);
      }
  }

void delaunay_draw_face_circle(quad_arc_t e, void *closure)
  {
    PSStream *ps = (PSStream *)closure;
    delaunay_site_t *ss = quad_org(e);
    delaunay_site_t *ts = quad_dst(e);
    delaunay_site_t *ls = quad_dst(quad_onext(e));
    
    Point *sp = &(ss->p);
    Point *tp = &(ts->p);
    Point *lp = &(ls->p);
    Point c;
    double rad;
    
    if (r2_orient(sp, tp, lp) > 0)
      { c = circumcenter(sp, tp, lp); }
    else
      { affirm(0, "draw_face_circle: bad face"); }
      
    rad = r2_dist(sp, &c);
    pswr_set_pen(ps, 0,0,0, 0.0, 1.5, 1.5);
    pswr_set_fill_color(ps, 0.75,0.75,0.75); 
    pswr_circle(ps, c.c[0], c.c[1], rad, TRUE, TRUE);
    pswr_set_fill_color(ps, 0,0,0); 
    pswr_dot(ps, c.c[0], c.c[1], 0.1, TRUE, FALSE);
  }

void delaunay_plot_diagram(PSStream *ps, quad_arc_t e)
  { 
    quad_enum(e, delaunay_draw_edge, (void *) ps);
  }

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

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

    for (i=0; i<ne; i++)
      { quad_arc_t g = e[i];
        delaunay_site_t *ss = quad_org(g);
        delaunay_site_t *ts = quad_dst(g);
        Point *sp = &(ss->p);
        Point *tp = &(ts->p);
        Point m = {{(sp->c[0] + tp->c[0])/2.0, (sp->c[1] + tp->c[1])/2.0}};
        double d = r2_dist(&m, &o);
        double h = r2_dist(sp, tp);
        if ((d < dbest) && (h > hmin)) { dbest = d; ebest = g; }
      }

    return (ebest);
  }

void delaunay_draw(PSStream *ps)
  {
    delaunay_site_t sites[DELAUNAY_NSITES];
    int max_edges = 3*(DELAUNAY_NSITES-2);
    quad_arc_t 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 */
        pswr_set_pen(ps, 0,0,0, 0.0, 1.0, 2.0);
        voronoi_plot_diagram(ps, quad_rot(e));
      }
    else if (delaunay_variant == 1)
      { /* Delaunay with an edge circle */
        Edge_list *el = collect_edges(e, max_edges);
        quad_arc_t g = delaunay_pick_nice_edge(el->n_edges, el->edges);
        delaunay_draw_edge_circle(g, (void *)ps);
        pswr_set_pen(ps, 0,0,0, 0.0, 1.0, 2.0);
        voronoi_plot_diagram(ps, quad_rot(e));
      }
    else if (delaunay_variant == 2)
      { /* Delaunay with a face circle */
        Edge_list *el = collect_edges(e, max_edges);
        quad_arc_t g = delaunay_pick_nice_edge(el->n_edges, el->edges);
        g = quad_sym(quad_oprev(quad_sym(g)));
        delaunay_draw_face_circle(g, (void *)ps); 
        pswr_set_pen(ps, 0,0,0, 0.0, 1.0, 2.0);
        voronoi_plot_diagram(ps, quad_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;
        
        /* Delete edges that cross xmedian and are above ycut: */
        pswr_set_pen(ps, 0,0,0, 0.5, 0.0, 0.0);
        for (i=0; i<el->n_edges; i++)
          { quad_arc_t g = el->edges[i];
            delaunay_site_t *ss = quad_org(g);
            delaunay_site_t *ts = quad_dst(g);
            if (ss->p.c[0] > ts->p.c[0]) { delaunay_site_t *aux = ss; ss = ts; ts = aux;}
            Point *sp = &(ss->p);
            Point *tp = &(ts->p);
            if ((sp->c[0] < xmedian) && (tp->c[0] >= xmedian))
              { if (r2_orient(sp, tp, &cutpoint) < 0)
                  { quad_destroy_edge(g); }
                else
                  { delaunay_draw_edge(g, (void *)ps);
                    delaunay_trace_bubble(g, (void*)ps); }
              }
          }
        /* Draw cut line: */
        pswr_set_pen(ps, 0,0,0, 0.0, 2.0, 1.0);
        pswr_segment(ps, xmedian, -1.0, xmedian, +1.0);
      }
    else
      { affirm(0, "delaunay_variant is bogus"); }

    pswr_set_pen(ps, 0,0,0, 0.0, 0.0, 0.0);
    delaunay_plot_diagram(ps, e);

    plot_sites(ps, 0.5, DELAUNAY_NSITES, sites);

  }

/* GABRIEL DIAGRAM */

void gabriel_draw(PSStream *ps);
void gabriel_make_sites(int n, delaunay_site_t sites[]);
quad_arc_t gabriel_from_delaunay(quad_arc_t e, int n, delaunay_site_t sites[]);
void gabriel_plot_diagram(PSStream *ps, quad_arc_t e);
void gabriel_analyze_edge(quad_arc_t e, void *closure);
void gabriel_draw_edge(quad_arc_t e, void *closure);
void gabriel_draw_edge_circle(quad_arc_t e, void *closure);

#define GABRIEL_NSITES 12

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

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

void gabriel_analyze_edge(quad_arc_t e, void *closure)
  {
    gabriel_data_t *cls = (gabriel_data_t *) closure;
    delaunay_site_t *ss = quad_org(e);
    delaunay_site_t *ts = quad_dst(e);
    Point *sp = &(ss->p);
    Point *tp = &(ts->p);
    Point m = {{(sp->c[0] + tp->c[0])/2.0, (sp->c[1] + tp->c[1])/2.0}};
    double r = r2_dist(sp, &m);
    
    int i;
    int nin = 0;
    
    for (i=0; i < (cls->n); i++)
      { Point *u = &(cls->sites[i].p);
        double d = r2_dist(u, &m);
        if ((d < r) && (u != sp) && (u != tp)) 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)++;
      }
  }

quad_arc_t gabriel_from_delaunay(quad_arc_t e, int n, delaunay_site_t sites[])
  {
    int i;

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

    affirm(cls.n_gut_edges > 0, "gabriel_from_delaunay: no edges found");
    
    /* Delete non-edges: */
    for (i=0; i < cls.n_bad_edges; i++)
      { quad_arc_t g = cls.bad_edges[i];
        quad_destroy_edge(g);
      }
      
    /* Return some edge: */
    { int ne = cls.n_gut_edges;
      double dbest = 1.0e20;
      quad_arc_t ebest = cls.gut_edges[0]; /* Just in case */
      Point o = {{0.0, 0.0}};
      double hmin = 0.20;
      
      for (i=0; i<ne; i++)
        { quad_arc_t g = cls.gut_edges[i];
          delaunay_site_t *ss = quad_org(g);
          delaunay_site_t *ts = quad_dst(g);
          Point *sp = &(ss->p);
          Point *tp = &(ts->p);
          Point m = {{(sp->c[0] + tp->c[0])/2.0, (sp->c[1] + tp->c[1])/2.0}};
          double d = r2_dist(&m, &o);
          double h = r2_dist(sp, tp);
          if ((d < dbest) && (h > hmin)) { dbest = d; ebest = g; }
        }
      
      return (ebest);
    }
  }

void gabriel_draw_edge(quad_arc_t e, void *closure)
  {
    PSStream *ps = (PSStream *)closure;
    delaunay_site_t *op = quad_org(e);
    delaunay_site_t *dp = quad_dst(e);
    pswr_segment(ps, op->p.c[0], op->p.c[1], dp->p.c[0], dp->p.c[1]);
  }

void gabriel_draw_edge_circle(quad_arc_t e, void *closure)
  {
    PSStream *ps = (PSStream *)closure;
    delaunay_site_t *ss = quad_org(e);
    delaunay_site_t *ts = quad_dst(e);
    Point *sp = &(ss->p);
    Point *tp = &(ts->p);
    Point m = {{(sp->c[0] + tp->c[0])/2.0, (sp->c[1] + tp->c[1])/2.0}};
    double r = r2_dist(sp, &m);
    pswr_set_pen(ps, 0,0,0, 0.0, 1.5, 1.5);
    pswr_set_fill_color(ps, 0.75,0.75,0.75);
    pswr_circle(ps, m.c[0], m.c[1], r, TRUE, TRUE);
  }

void gabriel_plot_diagram(PSStream *ps, quad_arc_t e)
  { 
    quad_enum(e, gabriel_draw_edge, (void *)ps);
  }

void gabriel_draw(PSStream *ps)
  {
    delaunay_site_t sites[GABRIEL_NSITES];
    quad_arc_t 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 *)ps);
    
    pswr_set_pen(ps, 0,0,0, 0.0, 0.0, 0.0);
    gabriel_plot_diagram(ps, e);

    plot_sites(ps, 0.5, GABRIEL_NSITES, sites);

  }

/* RELATIVE NEIGHBORHOOD GRAPH */

void relnei_draw(PSStream *ps);
void relnei_make_sites(int n, delaunay_site_t sites[]);
quad_arc_t relnei_from_delaunay(quad_arc_t e, int n, delaunay_site_t sites[]);
void relnei_plot_diagram(PSStream *ps, quad_arc_t e);
void relnei_analyze_edge(quad_arc_t e, void *closure);
void relnei_draw_edge(quad_arc_t e, void *closure);
void relnei_draw_edge_lune(quad_arc_t e, void *closure);

#define RELNEI_NSITES 12

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

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

void relnei_analyze_edge(quad_arc_t e, void *closure)
  {
    relnei_data_t *cls = (relnei_data_t *) closure;
    delaunay_site_t *ss = quad_org(e);
    delaunay_site_t *ts = quad_dst(e);
    Point *sp = &(ss->p);
    Point *tp = &(ts->p);
    double r = r2_dist(sp, tp);
    
    int i;
    int nin = 0;
    
    for (i=0; i < (cls->n); i++)
      { Point *u = &(cls->sites[i].p);
        double ds = r2_dist(u, sp);
        double dt = r2_dist(u, tp);
        if ((u != sp) && (u != tp) && (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)++;
      }
  }

quad_arc_t relnei_from_delaunay(quad_arc_t e, int n, delaunay_site_t sites[])
  {
    int i;

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

    affirm(cls.n_gut_edges > 0, "relnei_from_delaunay: no edges found");
    
    /* Delete non-edges: */
    for (i=0; i < cls.n_bad_edges; i++)
      { quad_arc_t g = cls.bad_edges[i];
        quad_destroy_edge(g);
      }
      
    /* Return some edge: */
    { int ne = cls.n_gut_edges;
      double dbest = 1.0e20;
      quad_arc_t ebest = cls.gut_edges[0]; /* Just in case */
      Point o = {{0.0, 0.0}};
      double hmin = 0.20;
      
      for (i=0; i<ne; i++)
        { quad_arc_t g = cls.gut_edges[i];
          delaunay_site_t *ss = quad_org(g);
          delaunay_site_t *ts = quad_dst(g);
          Point *sp = &(ss->p);
          Point *tp = &(ts->p);
          Point m = {{(sp->c[0] + tp->c[0])/2.0, (sp->c[1] + tp->c[1])/2.0}};
          double d = r2_dist(&m, &o);
          double h = r2_dist(sp, tp);
          if ((d < dbest) && (h > hmin)) { dbest = d; ebest = g; }
        }
      
      return (ebest);
    }
  }

void relnei_draw_edge(quad_arc_t e, void *closure)
  {
    PSStream *ps = (PSStream *)closure;
    delaunay_site_t *op = quad_org(e);
    delaunay_site_t *dp = quad_dst(e);
    pswr_segment(ps, op->p.c[0], op->p.c[1], dp->p.c[0], dp->p.c[1]);
  }

void relnei_draw_edge_lune(quad_arc_t e, void *closure)
  {
    PSStream *ps = (PSStream *)closure;
    delaunay_site_t *ss = quad_org(e);
    delaunay_site_t *ts = quad_dst(e);
    Point *sp = &(ss->p);
    Point *tp = &(ts->p);
    Point m = {{(sp->c[0] + tp->c[0])/2.0, (sp->c[1] + tp->c[1])/2.0}};
    double r = r2_dist(sp, &m);
    double tilt = atan2(sp->c[1] - m.c[1], sp->c[0] - m.c[0]);

    pswr_set_pen(ps, 0,0,0, 0.0, 1.5, 1.5);
    pswr_set_fill_color(ps, 0.75,0.75,0.75);
    pswr_lune(ps, m.c[0], m.c[1], r, tilt, TRUE, TRUE);
  }

void relnei_plot_diagram(PSStream *ps, quad_arc_t e)
  { 
    quad_enum(e, relnei_draw_edge, (void *)ps);
  }

void relnei_draw(PSStream *ps)
  {
    delaunay_site_t sites[RELNEI_NSITES];
    quad_arc_t 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 *)ps);
    
    pswr_set_pen(ps, 0,0,0, 0.0, 0.0, 0.0);
    relnei_plot_diagram(ps, e);

    plot_sites(ps, 0.5, RELNEI_NSITES, sites);

  }

/* MINIMUM SPANNING TREE */

void mintree_draw(PSStream *ps);
void mintree_make_sites(int n, delaunay_site_t sites[]);
quad_arc_t mintree_from_delaunay(quad_arc_t e, int n, delaunay_site_t sites[]);
void mintree_plot_diagram(PSStream *ps, quad_arc_t e);
void mintree_save_edge(quad_arc_t e, void *closure);
void mintree_sort_edges(int ne, quad_arc_t e[]);
int mintree_find_root(delaunay_site_t *s, int n, delaunay_site_t sites[], int vp[]);
void mintree_draw_edge(quad_arc_t e, void *closure);

#define MINTREE_NSITES 12

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

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

void mintree_save_edge(quad_arc_t 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, quad_arc_t e[])
  /* Sort e[] by increasing length */
  { int i, j;
    double *r = rn_alloc(ne);
    for (i=1; i<ne; i++)
      { quad_arc_t ei = e[i]; 
        delaunay_site_t *s = quad_org(ei);
        delaunay_site_t *t = quad_dst(ei);
        double ri = r2_dist(&(s->p), &(t->p));
        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(delaunay_site_t *s, int n, delaunay_site_t sites[], int vp[])
  { int i;
    for (i=0; (i<n) && (s != &(sites[i])); i++) { }
    affirm(i<n, "find_root: didn\'t find it");
    while (vp[i] != i) i = vp[i];
    return (i);
  }

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

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

    mintree_data cls;
    cls.n_del_edges = 0;
    cls.del_edges = (quad_arc_t *) malloc(max_edges * sizeof(quad_arc_t));
    
    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++)
        { quad_arc_t ei = (cls.del_edges)[ie];
          int is = mintree_find_root(quad_org(ei), n, sites, vp);
          int it = mintree_find_root(quad_dst(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++;
            }
        }
    }
    
    affirm(n_gut_edges > 0, "mintree_from_delaunay: no edges found");
    
    /* Delete bad edges: */    
    { int ie;
      for (ie=0; ie < n_bad_edges; ie++)
        { quad_arc_t g = bad_edges[ie];
          quad_destroy_edge(g);
        }
    }
      
    /* Return some edge: */
    { int ie;
      int ne = n_gut_edges;
      double dbest = 1.0e20;
      quad_arc_t ebest = gut_edges[0]; /* Just in case */
      Point o = {{0.0, 0.0}};
      double hmin = 0.20;
      
      for (ie=0; ie<ne; ie++)
        { quad_arc_t g = gut_edges[ie];
          delaunay_site_t *s = quad_org(g);
          delaunay_site_t *t = quad_dst(g);
          Point m = {{(s->p.c[0] + t->p.c[0])/2.0, (s->p.c[1] + t->p.c[1])/2.0}};
          double d = r2_dist(&m, &o);
          double h = r2_dist(&(s->p), &(t->p));
          if ((d < dbest) && (h > hmin)) { dbest = d; ebest = g; }
        }
      
      return (ebest);
    }
  }

void mintree_draw_edge(quad_arc_t e, void *closure)
  {
    PSStream *ps = (PSStream *)closure;
    delaunay_site_t *op = quad_org(e);
    delaunay_site_t *dp = quad_dst(e);
    pswr_segment(ps, op->p.c[0], op->p.c[1], dp->p.c[0], dp->p.c[1]);
  }

void mintree_plot_diagram(PSStream *ps, quad_arc_t e)
  { 
    quad_enum(e, mintree_draw_edge, (void *)ps);
  }

void mintree_draw(PSStream *ps)
  {
    delaunay_site_t sites[MINTREE_NSITES];
    quad_arc_t 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));

    pswr_set_pen(ps, 0,0,0, 0.0, 0.0, 0.0);
    mintree_plot_diagram(ps, e);

    plot_sites(ps, 0.5, MINTREE_NSITES, sites);

  }

/* MAIN PROGRAM */

void delvor_draw_all(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); 
      }
  }