/*
 * voiro.c
 *
 * compute the voronoi diagram and the convex hull of a set of random points.
 * Uses ideas from an O(n log n) algorithm, not sure if it is actually really
 * O (n log n) though.
 *
 * has numerical robustness issues.
 *
 * (c) 2006 Simon Budig <simon@budig.de>, released under the GPL v2.0.
 *
 *
 * compile with the following command:
 *  gcc -Wall -g -o voiro `pkg-config --libs --cflags cairo gtk+-2.0` voiro.c
 *
 * you can optionally give the number of points on the commandline,
 * a number < 2 gives a testcase, where the numerical issues become apparent.
 *
 * The following hotkeys are there:
 *
 *   d        - dump the coordinates of the points
 *   f, F11   - toggle fullscreen
 *   p        - switch between three ways to draw the dots
 *   t        - toggle text labels on the points
 *   c        - toggle drawing of convex hull
 *   v        - toggle drawing of voronoi diagram
 *   <Right>  - next step in the algorithm
 *   <left>   - previous step in the algorithm
 *   <Return> - finish the algorithm
 *   n        - restart the algorithm, repeated pressing generates new pointsets
 *   q, <Esc> - quit.
 */

#include <gtk/gtk.h>
#include <gdk/gdkkeysyms.h>
#include <cairo.h>
#include <math.h>
#include <stdlib.h>
#include <stdlib.h>

typedef struct _xy_pair         Point;
typedef struct _voiro_diagram   VoiroDiagram;
typedef struct _voiro_convex    VoiroConvex;
typedef struct _voiro_boundary  VoiroBoundary;
typedef struct _voiro_context   VoiroContext;
typedef enum   _voiro_runmode   VoiroRunmode;
typedef enum   _voiro_boundtype VoiroBoundtype;

enum _voiro_runmode
{
  VOIRO_STEPPING,
  VOIRO_FINISH,
  VOIRO_RESTART,
  VOIRO_NEW,
  VOIRO_QUIT,
};

enum _voiro_boundtype
{
  VOIRO_INVALID,
  VOIRO_SEGMENT,
  VOIRO_RAY,
  VOIRO_LINE,
};


struct _xy_pair
{
  double x, y;
};

struct _voiro_diagram
{
  gint            max_level;
  gint            n_sites;
  Point          *sites;
  gint           *sorted_order;
  VoiroConvex    *convex;
  VoiroBoundary  *mergeline;
  VoiroBoundary **bounds;
};

struct _voiro_convex
{
  VoiroDiagram *diagram;

  VoiroConvex  *left;
  VoiroConvex  *right;

  gint          level;

  gint          n_upper;
  gint         *upper;
  gint          n_lower;
  gint         *lower;
};

struct _voiro_boundary
{
  gint           site;
  gint           neighbor;

  gboolean       valid;

  Point          p0;
  gboolean       end0;
  Point          p1;
  gboolean       end1;

  VoiroBoundary *next;
  VoiroBoundary *nnext;
};

struct _voiro_context
{
  VoiroDiagram *diagram;

  gdouble       min_x;
  gdouble       max_x;
  gdouble       min_y;
  gdouble       max_y;
  gdouble       pixelwidth;

  gint          show_level;
  gboolean      fullscreen;
  gboolean      dot_mode;
  gboolean      show_text;
  gboolean      show_convex;
  gboolean      show_voronoi;

  gboolean      in_algorithm;
  VoiroRunmode  runmode;
  gint          showstep;
  gint          step;

  GtkWidget    *window;
  GtkWidget    *darea;
  cairo_t      *cr;
};


void
voiro_algorithm_step (VoiroContext *new_context)
{
  static VoiroContext *context = NULL;

  if (new_context)
    {
      context = new_context;
      return;
    }

  if (context->runmode == VOIRO_STEPPING &&
      context->step >= context->showstep)
    gtk_main ();

  context->step++;
}


gint
ori (Point p, Point q, Point r)
{
  gdouble dx1, dy1, dx2, dy2;

  dx1 = q.x - p.x;
  dy1 = q.y - p.y;
  dx2 = r.x - p.x;
  dy2 = r.y - p.y;

  if (dx1*dy2 > dy1*dx2)
    return 1;
  else if (dx1*dy2 < dy1*dx2)
    return -1;
  else if ((dx1*dx2 < 0) || (dy1*dy2 < 0))
    return -1;
  else if ((dx1*dx1+dy1*dy1) >= (dx2*dx2 + dy2*dy2))
    return 0;
  else
    return 1;
}


void
voiro_convex_free (VoiroConvex *convex)
{
  if (convex->left)
    voiro_convex_free (convex->left);
  if (convex->right)
    voiro_convex_free (convex->right);

  g_free (convex->upper);
  g_free (convex->lower);
  g_free (convex);
}


void
voiro_sort_points (Point *points, gint npoints)
{
  Point *part1, *part2;
  gint size1, size2;
  gint i1, i2;

  if (npoints < 2)
    return;

  size1 = npoints/2;
  part1 = g_memdup (points, size1 * sizeof (Point));
  if (size1 > 1)
    voiro_sort_points (part1, size1);

  size2 = npoints - size1;
  part2 = g_memdup (&points[size1], size2 * sizeof (Point));
  if (size2 > 1)
    voiro_sort_points (part2, size2);

  i1 = i2 = 0;

  while (i1 < size1 && i2 < size2)
    {
      if (part1[i1].x < part2[i2].x)
        {
          points[i1 + i2] = part1[i1];
          i1 ++;
        }
      else if (part1[i1].x > part2[i2].x)
        {
          points[i1 + i2] = part2[i2];
          i2 ++;
        }
      else
        {
          if (part1[i1].y <= part2[i2].y)
            {
              points[i1 + i2] = part1[i1];
              i1++;
            }
          else
            {
              points[i1 + i2] = part2[i2];
              i2++;
            }
        }
    }

  for (; i1 < size1; i1++)
    points[i1 + i2] = part1[i1];

  for (; i2 < size2; i2++)
    points[i1 + i2] = part2[i2];

  g_free (part1);
  g_free (part2);
}


gboolean
voiro_bounds_intersect (VoiroDiagram  *diagram,
                        VoiroBoundary *bound1,
                        VoiroBoundary *bound2,
                        gdouble       *ret_t1,
                        gdouble       *ret_t2,
                        Point         *intersect)
{
  gdouble det, t1, t2;

  if (!bound1->valid)
    return FALSE;

  det = ((bound1->p1.y - bound1->p0.y) * (bound2->p1.x - bound2->p0.x) -
         (bound1->p1.x - bound1->p0.x) * (bound2->p1.y - bound2->p0.y));

  if (det == 0)
    return FALSE;

  t1 = (bound1->p0.x - bound2->p0.x) * (bound2->p1.y - bound2->p0.y) -
       (bound1->p0.y - bound2->p0.y) * (bound2->p1.x - bound2->p0.x);

  t2 = (bound2->p0.y - bound1->p0.y) * (bound1->p1.x - bound1->p0.x) -
       (bound2->p0.x - bound1->p0.x) * (bound1->p1.y - bound1->p0.y);

  t1 /= det;
  t2 /= det;

  if (ret_t1)
    *ret_t1 = t1;

  if (ret_t2)
    *ret_t2 = t2;

  if (intersect)
    {
      intersect->x = bound1->p0.x + t1 * (bound1->p1.x - bound1->p0.x);
      intersect->y = bound1->p0.y + t1 * (bound1->p1.y - bound1->p0.y);
    }

  if ((bound1->end0 && t1 < 0.0) ||
      (bound2->end0 && t2 < 0.0))
    return FALSE;

  if ((bound1->end1 && t1 > 1.0) ||
      (bound2->end1 && t2 > 1.0))
    return FALSE;

  return TRUE;
}

void
voiro_bounds_cut_cw (VoiroBoundary *first,
                     gdouble        ft,
                     VoiroBoundary *second,
                     gdouble        st,
                     Point          intersection)
{
  Point dir;
  gboolean tmp;

  if (ft < 0)
    {
      dir.x = first->p1.x - first->p0.x;
      dir.y = first->p1.y - first->p0.y;

      first->p0.x = intersection.x - dir.x;
      first->p0.y = intersection.y - dir.y;
    }

  first->p1 = intersection;
  first->end1 = TRUE;

  dir.x = second->p1.x - second->p0.x + first->p1.x;
  dir.y = second->p1.y - second->p0.y + first->p1.y;

  if (ori (first->p0, first->p1, dir) < 0)
    {
      /* flip direction of second boundary */

      dir = second->p0;
      second->p0 = second->p1;
      second->p1 = dir;

      tmp = second->end0;
      second->end0 = second->end1;
      second->end1 = tmp;

      st = 1.0 - st;
    }

  if (st > 1.0)
    {
      dir.x = second->p1.x - second->p0.x;
      dir.y = second->p1.y - second->p0.y;

      second->p1.x = intersection.x + dir.x;
      second->p1.y = intersection.y + dir.y;
    }

  second->p0 = intersection;
  second->end0 = TRUE;
}


void
voiro_bounds_cut_ccw (VoiroBoundary *first,
                      gdouble        ft,
                      VoiroBoundary *second,
                      gdouble        st,
                      Point          intersection)
{
  Point dir;
  gboolean tmp;

  if (ft < 0)
    {
      dir.x = first->p1.x - first->p0.x;
      dir.y = first->p1.y - first->p0.y;

      first->p0.x = intersection.x - dir.x;
      first->p0.y = intersection.y - dir.y;
    }

  first->p1 = intersection;
  first->end1 = TRUE;

  dir.x = second->p1.x - second->p0.x + first->p1.x;
  dir.y = second->p1.y - second->p0.y + first->p1.y;

  if (ori (first->p0, first->p1, dir) > 0)
    {
      /* flip direction of second boundary */

      dir = second->p0;
      second->p0 = second->p1;
      second->p1 = dir;

      tmp = second->end0;
      second->end0 = second->end1;
      second->end1 = tmp;

      st = 1.0 - st;
    }

  if (st > 1.0)
    {
      dir.x = second->p1.x - second->p0.x;
      dir.y = second->p1.y - second->p0.y;

      second->p1.x = intersection.x + dir.x;
      second->p1.y = intersection.y + dir.y;
    }

  second->p0 = intersection;
  second->end0 = TRUE;
}


void
voiro_bounds_check_ori (VoiroDiagram  *diagram,
                        gint           site,
                        VoiroBoundary *boundary)
{
  VoiroBoundary *b;
  gint o;

  o = ori (boundary->p0, boundary->p1, diagram->sites[site]);

  for (b = diagram->bounds[site]; b != NULL;
       b = (b->site == site) ? b->next : b->nnext)
    {
      if (ori (boundary->p0, boundary->p1, b->p0) == o ||
          ori (boundary->p0, boundary->p1, b->p1) == o)
        continue;

      if (voiro_bounds_intersect (diagram, boundary, b, NULL, NULL, NULL))
        continue;

      b->valid = FALSE;
    }
}


void
voiro_insert_line (VoiroDiagram  *diagram,
                   VoiroBoundary *boundary)
{
  boundary->next = diagram->bounds[boundary->site];
  diagram->bounds[boundary->site] = boundary;

  boundary->nnext = diagram->bounds[boundary->neighbor];
  diagram->bounds[boundary->neighbor] = boundary;
}


VoiroBoundary *
voiro_build_middle (VoiroDiagram *diagram,
                    gint          site,
                    gint          neighbor,
                    Point        *p0)
{
  VoiroBoundary *boundary;
  gdouble dx, dy;

  g_return_val_if_fail (site != neighbor, NULL);

  boundary = g_new0 (VoiroBoundary, 1);
  boundary->site = site;
  boundary->neighbor = neighbor;
  boundary->next = NULL;
  boundary->nnext = NULL;

  dx = diagram->sites[neighbor].x - diagram->sites[site].x;
  dy = diagram->sites[neighbor].y - diagram->sites[site].y;

  if (p0)
    {
      boundary->p0.x = p0->x;
      boundary->p0.y = p0->y;
      boundary->end0 = TRUE;
    }
  else
    {
      boundary->p0.x = diagram->sites[site].x + dx / 2;
      boundary->p0.y = diagram->sites[site].y + dy / 2;
      boundary->end0 = FALSE;
    }

  boundary->p1.x = boundary->p0.x - dy;
  boundary->p1.y = boundary->p0.y + dx;

  boundary->end1 = FALSE;
  boundary->valid = TRUE;

  return boundary;
}


void
voiro_merge_support_ccw (VoiroDiagram  *diagram,
                         gint          *start,
                         gint           n_start,
                         gint          *end,
                         gint           n_end,
                         gint          *ret_startpt,
                         gint          *ret_endpt,
                         gint         **ret_merged,
                         gint          *ret_n_merged)
{
  gint i, si, ei;
  gboolean busy;

  /* merge two ccw boundaries into one */
  si = n_start - 1;
  ei = 0;

  busy = TRUE;
  while (busy)
    {
      busy = FALSE;
      while (ei+1 < n_end &&
             ori (diagram->sites[start[si]],
                  diagram->sites[end[ei]],
                  diagram->sites[end[ei+1]]) >= 0)
        {
          ei ++;
          busy = TRUE;
        }

      while (si > 0 &&
             ori (diagram->sites[start[si-1]],
                  diagram->sites[start[si]],
                  diagram->sites[end[ei]]) >= 0)
        {
          si --;
          busy = TRUE;
        }
    }

  *ret_n_merged = si + 1 + (n_end - ei);
  *ret_merged = g_new0 (gint, *ret_n_merged);
  for (i = 0; i <= si; i++)
    (*ret_merged)[i] = start[i];
  for (i = ei; i < n_end; i++)
    (*ret_merged)[si+1+i-ei] = end[i];

  *ret_startpt = si;
  *ret_endpt   = ei;
}


/* this needs sorted points */
VoiroConvex *
voiro_convex_hull (VoiroDiagram  *diagram,
                   gint          *points,
                   gint           n_points,
                   VoiroConvex  **place)
{
  VoiroConvex *convex, *l_convex, *r_convex;
  gint size1, size2;
  gint  i, li, ri, nb, lower_li, lower_ri;
  VoiroBoundary *middle, *b, *min_bound, *forbidden;
  gdouble        tm, tb, min_tm, min_tb;
  Point          inter, min_inter;

  g_return_val_if_fail (n_points > 0, NULL);

  convex = g_new0 (VoiroConvex, 1);
  *place = convex;

  if (n_points <= 2)
    {
      convex->n_upper = convex->n_lower = n_points;
      convex->upper = g_new (gint, n_points);
      convex->lower = g_new (gint, n_points);
      for (i = 0; i < n_points; i++)
        convex->upper[i] = convex->lower[n_points - 1 - i] = points[i];

      convex->left = convex->right = NULL;
      convex->level = 1;
      diagram->max_level = MAX (convex->level, diagram->max_level);

      if (n_points == 2)
        {
          middle = voiro_build_middle (diagram, points[0], points[1], NULL);
          voiro_insert_line (diagram, middle);
        }
      voiro_algorithm_step (NULL);
      return convex;
    }

  /* Divide and Conquer */

  size1 = n_points / 2;
  size2 = n_points - size1;

  l_convex = voiro_convex_hull (diagram, points, size1, &convex->left);
  r_convex = voiro_convex_hull (diagram, &points[size1], size2, &convex->right);

  voiro_merge_support_ccw (diagram,
                           r_convex->lower, r_convex->n_lower,
                           l_convex->lower, l_convex->n_lower,
                           &lower_ri, &lower_li,
                           &convex->lower, &convex->n_lower);

  voiro_merge_support_ccw (diagram,
                           l_convex->upper, l_convex->n_upper,
                           r_convex->upper, r_convex->n_upper,
                           &li, &ri,
                           &convex->upper, &convex->n_upper);


  diagram->mergeline = NULL;

  lower_li = l_convex->lower[lower_li];
  lower_ri = r_convex->lower[lower_ri];

  li = l_convex->upper[li];
  ri = r_convex->upper[ri];

  diagram->mergeline = middle = voiro_build_middle (diagram, ri, li, NULL);
  middle->next = NULL;

  forbidden = NULL;

  while (TRUE)
    {
      min_bound = NULL;

      /* look for the topmost boundary intersecting middle */
      for (b = diagram->bounds[li]; b != NULL;
           b = (b->site == li) ? b->next : b->nnext)
        {
          nb = (b->site == li) ? b->neighbor : b->site;

          if (b != forbidden &&
              nb != ri &&
              nb != li &&
              voiro_bounds_intersect (diagram,
                                      middle, b,
                                      &tm, &tb, &inter))
            {
              if (min_bound == NULL || tm < min_tm)
                {
                  min_bound = b;
                  min_tm = tm;
                  min_tb = tb;
                  min_inter = inter;
                }
            }
        }

      for (b = diagram->bounds[ri]; b != NULL;
           b = (b->site == ri) ? b->next : b->nnext)
        {
          nb = (b->site == ri) ? b->neighbor : b->site;

          if (b != forbidden &&
              nb != li &&
              nb != ri &&
              voiro_bounds_intersect (diagram,
                                      middle, b,
                                      &tm, &tb, &inter))
            {
              if (min_bound == NULL || tm < min_tm)
                {
                  min_bound = b;
                  min_tm = tm;
                  min_tb = tb;
                  min_inter = inter;
                }
            }
        }

      if (min_bound)
        {
          if (min_bound->site == ri ||
              min_bound->neighbor == ri)
            {
              voiro_bounds_cut_cw (middle, min_tm,
                                   min_bound, min_tb, min_inter);
              voiro_bounds_check_ori (diagram, ri, middle);
              voiro_algorithm_step (NULL);

              forbidden = min_bound;

              ri = (min_bound->site == ri) ? min_bound->neighbor : min_bound->site;
            }
          else if (min_bound->site == li ||
                   min_bound->neighbor == li)
            {
              voiro_bounds_cut_ccw (middle, min_tm,
                                    min_bound, min_tb, min_inter);
              voiro_bounds_check_ori (diagram, li, middle);
              voiro_algorithm_step (NULL);

              forbidden = min_bound;

              li = (min_bound->site == li) ? min_bound->neighbor : min_bound->site;
            }
          else
            {
              g_printerr ("Found a weird line: (%d) %d, %d...\n",
                          min_bound->site, ri, li);
            }

          middle = voiro_build_middle (diagram, ri, li, &min_inter);
          middle->next = diagram->mergeline;
          diagram->mergeline = middle;
        }
      else
        {
          break;
        }

      voiro_algorithm_step (NULL);

      if (li == lower_li && ri == lower_ri)
        break;
    }

  for (middle = diagram->mergeline; middle != NULL; )
    {
      b = middle->next;
      voiro_insert_line (diagram, middle);
      middle = b;
    }

  diagram->mergeline = NULL;

  convex->level = MAX (l_convex->level, r_convex->level) + 1;
  diagram->max_level = MAX (convex->level, diagram->max_level);

  voiro_algorithm_step (NULL);
  return convex;
}


void
voiro_cleanup_model (VoiroDiagram *diagram)
{
  gint i;

  if (diagram->bounds)
    {
      for (i = 0; i < diagram->n_sites; i++)
        {
          VoiroBoundary *b, *next;

          b = diagram->bounds[i];
          while (b)
            {
              next = b->site == i ? b->next : b->nnext;
              if (b->site <= i && b->neighbor <= i)
                g_free (b);
              b = next;
            }

          diagram->bounds[i] = NULL;
        }
    }

  if (diagram->convex)
    {
      voiro_convex_free (diagram->convex);
      diagram->convex = NULL;
    }

  diagram->convex = NULL;
}


void
voiro_init_model (VoiroDiagram *diagram,
                  gint          n_sites)
{
  gdouble r, phi;
  gint i;
  gboolean test = FALSE;

  voiro_cleanup_model (diagram);

  if (diagram->sites)
    g_free (diagram->sites);

  if (diagram->sorted_order)
    g_free (diagram->sorted_order);

  if (diagram->bounds)
    g_free (diagram->bounds);
  if (n_sites < 2)
    {
      test = TRUE;
      n_sites = 8;
    }

  diagram->n_sites = n_sites;
  diagram->sites  = g_new (Point, diagram->n_sites);

  if (test)
    {
      for (i = 0; i < n_sites; i++)
        {
          diagram->sites[i].x = cos ((gdouble) i / n_sites * G_PI * 2);
          diagram->sites[i].y = sin ((gdouble) i / n_sites * G_PI * 2);
        }
    }
  else
    {
      for (i = 0; i < diagram->n_sites; i++)
        {
          r = sqrt (g_random_double_range (0, 1));
          phi = g_random_double_range (0, 2*G_PI);
          diagram->sites[i].x = r * cos (phi);
          diagram->sites[i].y = r * sin (phi);
        }
    }

  diagram->sorted_order = g_new0 (gint, diagram->n_sites);
  for (i = 0; i < diagram->n_sites; i++)
    diagram->sorted_order[i] = i;

  diagram->bounds = g_new0 (VoiroBoundary *, diagram->n_sites);

  voiro_sort_points (diagram->sites, diagram->n_sites);
}


void
voiro_draw_begin (VoiroContext *context)
{
  GtkAllocation *alloc;
  gdouble scaling;
  gdouble x0, x1, y0, y1;

  g_return_if_fail (context->cr == NULL);

  context->cr = gdk_cairo_create (context->darea->window);

  alloc = &context->darea->allocation;
  scaling = MIN (alloc->width, alloc->height) / 2.2;
  cairo_new_path (context->cr);
  cairo_rectangle (context->cr, 0, 0, alloc->width, alloc->height);
  cairo_clip (context->cr);
  cairo_stroke (context->cr);

  cairo_translate (context->cr, alloc->width / 2, alloc->height / 2);
  cairo_scale (context->cr, scaling, -scaling);

  x0 = y0 = 1.0;
  cairo_device_to_user_distance (context->cr, &x0, &y0);
  context->pixelwidth = MAX (x0, y0);

  cairo_set_font_size (context->cr, context->pixelwidth * 14);
  x0 = y0 = 0.0;
  x1 = alloc->width;
  y1 = alloc->height;
  cairo_device_to_user (context->cr, &x0, &y0);
  cairo_device_to_user (context->cr, &x1, &y1);
  context->min_x = MIN (x0, x1);
  context->min_y = MIN (y0, y1);
  context->max_x = MAX (x0, x1);
  context->max_y = MAX (y0, y1);
}


void
voiro_draw_end (VoiroContext *context)
{
  cairo_destroy (context->cr);
  context->cr = NULL;
}


void
voiro_draw_point (VoiroContext *context,
                  Point         point,
                  gint          nr,
                  gint          mode)
{
  cairo_t *cr = context->cr;
  char str[6];

  switch (mode)
    {
      case 0:
        break;
      case 1:
        cairo_set_source_rgba (cr, 0.0, 0.0, 0.0, 1.0);
        cairo_new_path (cr);
        cairo_arc (cr, point.x, point.y, 1*context->pixelwidth, 0, 2*G_PI);
        cairo_fill (cr);
        break;
      case 2:
        cairo_set_source_rgba (cr, 1.0, 1.0, 1.0, 1.0);
        cairo_new_path (cr);
        cairo_arc (cr, point.x, point.y, 3*context->pixelwidth, 0, 2*G_PI);
        cairo_fill_preserve (cr);
        cairo_set_source_rgba (cr, 0.0, 0.0, 0.0, 1.0);
        cairo_set_line_width (cr, 1.2*context->pixelwidth);
        cairo_stroke (cr);
        break;
    }

  if (nr < 0)
    return;

  snprintf (str, 6, "%d", nr);
  cairo_save (context->cr);
  cairo_move_to (context->cr,
                 point.x + 4 * context->pixelwidth,
                 point.y - 14 * context->pixelwidth);
  cairo_scale (context->cr, 1.0, -1.0);
  cairo_show_text (context->cr, str);
  cairo_restore (context->cr);
}


void
voiro_draw_polyline (VoiroContext *context,
                     gint         *points,
                     gint          n_points)
{
  cairo_t *cr = context->cr;
  gint i;
  g_return_if_fail (n_points > 0);

  cairo_set_line_width (cr, context->pixelwidth);

  for (i = 0; i < n_points; i++)
    {
      cairo_arc (cr, context->diagram->sites[points[i]].x,
                 context->diagram->sites[points[i]].y,
                 5*context->pixelwidth, 0, 2*G_PI);
      cairo_stroke (cr);
    }

  cairo_move_to (cr,
                 context->diagram->sites[points[0]].x,
                 context->diagram->sites[points[0]].y);

  for (i = 1; i < n_points; i++)
    cairo_line_to (cr,
                   context->diagram->sites[points[i]].x,
                   context->diagram->sites[points[i]].y);

  cairo_stroke (cr);
}


void
voiro_draw_convex (VoiroContext *context,
                   VoiroConvex  *convex,
                   gint          level)
{
  if ((level == 0 && convex->level) ||
      (level > 0 && convex->level == level))
    {
      cairo_set_source_rgba (context->cr, 0.0, 0.0, 1.0, 0.5);

      if (convex->n_upper)
        voiro_draw_polyline (context, convex->upper, convex->n_upper);
      if (convex->n_lower && convex->level > 1)
        voiro_draw_polyline (context, convex->lower, convex->n_lower);
    }
  else
    {
      if (convex->left)
        voiro_draw_convex (context, convex->left, level);
      if (convex->right)
        voiro_draw_convex (context, convex->right, level);
    }

}


gboolean
voiro_clip_boundary (VoiroContext  *context,
                     VoiroBoundary *boundary,
                     Point         *point0,
                     Point         *point1)
{
  Point p0, p1, dir, ptmp;
  gint code0, code1, ctmp;
  gdouble t;

#define CLIP_TOP    1
#define CLIP_BOTTOM 2
#define CLIP_RIGHT  4
#define CLIP_LEFT   8

  /* Cohen Sutherland, adjusted for unbounded lines */

  if (!boundary->valid)
    return FALSE;

  code0 = code1 = 0;

  p0 = boundary->p0;
  p1 = boundary->p1;
  dir.x = p1.x - p0.x;
  dir.y = p1.y - p0.y;

  /* tweak code for unbounded lines */
  if (boundary->end0 == FALSE)
    {
      if (dir.x > 0)
        code0 |= CLIP_LEFT;
      else if (dir.x < 0)
        code0 |= CLIP_RIGHT;
      if (dir.y > 0)
        code0 |= CLIP_BOTTOM;
      else if (dir.y < 0)
        code0 |= CLIP_TOP;
    }

  if (boundary->end1 == FALSE)
    {
      if (dir.x > 0)
        code1 |= CLIP_RIGHT;
      else if (dir.x < 0)
        code1 |= CLIP_LEFT;
      if (dir.y > 0)
        code1 |= CLIP_TOP;
      else if (dir.y < 0)
        code1 |= CLIP_BOTTOM;
    }

  while (1)
    {
      if (p0.x < context->min_x)
        code0 |= CLIP_LEFT;
      else if (p0.x > context->max_x)
        code0 |= CLIP_RIGHT;
      if (p0.y < context->min_y)
        code0 |= CLIP_BOTTOM;
      else if (p0.y > context->max_x)
        code0 |= CLIP_TOP;

      if (p1.x < context->min_x)
        code1 |= CLIP_LEFT;
      if (p1.x > context->max_x)
        code1 |= CLIP_RIGHT;
      if (p1.y < context->min_y)
        code1 |= CLIP_BOTTOM;
      if (p1.y > context->max_x)
        code1 |= CLIP_TOP;

      if ((code0 & code1) != 0)
        {
          return FALSE;
        }

      if (code0 == 0 && code1 == 0)
        {
          *point0 = p0;
          *point1 = p1;
          return TRUE;
        }

      if (code0 == 0)
        {
          /* nothing to clip on first endpoint, swap */
          ptmp = p0; p0 = p1; p1 = ptmp;
          ctmp = code0; code0 = code1; code1 = ctmp;
        }

      if (code0 & CLIP_LEFT)
        {
          t = (context->min_x - p0.x) / (p1.x - p0.x);
          p0.y += t * (p1.y - p0.y);
          p0.x = context->min_x;

          code0 &= ~CLIP_LEFT;
        }
      else if (code0 & CLIP_RIGHT)
        {
          t = (context->max_x - p0.x) / (p1.x - p0.x);
          p0.y += t * (p1.y - p0.y);
          p0.x = context->max_x;

          code0 &= ~CLIP_RIGHT;
        }
      else if (code0 & CLIP_BOTTOM)
        {
          t = (context->min_y - p0.y) / (p1.y - p0.y);
          p0.x += t * (p1.x - p0.x);
          p0.y = context->min_y;

          code0 &= ~CLIP_BOTTOM;
        }
      else if (code0 & CLIP_TOP)
        {
          t = (context->max_y - p0.y) / (p1.y - p0.y);
          p0.x += t * (p1.x - p0.x);
          p0.y = context->max_y;

          code0 &= ~CLIP_TOP;
        }
    }

#undef CLIP_LEFT
#undef CLIP_RIGHT
#undef CLIP_TOP
#undef CLIP_BOTTOM
}


void
voiro_draw_voronoi (VoiroContext *context)
{
  VoiroBoundary *boundary;
  Point p0, p1;
  gint i;

  cairo_set_line_width (context->cr, context->pixelwidth);
  cairo_set_source_rgba (context->cr, 0.0, 0.5, 0.0, 1.0);

  cairo_new_path (context->cr);

  for (i = 0; i < context->diagram->n_sites; i++)
    {
      boundary = context->diagram->bounds[i];
      while (boundary)
        {
          if (boundary->site == i && boundary->valid)
            {
              if (voiro_clip_boundary (context, boundary, &p0, &p1))
                {
                  cairo_move_to (context->cr, p0.x, p0.y);
                  cairo_line_to (context->cr, p1.x, p1.y);
                }
            }

          boundary = (i == boundary->site) ? boundary->next : boundary->nnext;
        }
    }

  cairo_stroke (context->cr);

  cairo_set_source_rgba (context->cr, 1.0, 0.0, 0.0, 1.0);

  cairo_new_path (context->cr);

  boundary = context->diagram->mergeline;
  while (boundary)
    {
      if (boundary->valid)
        {
          if (voiro_clip_boundary (context, boundary, &p0, &p1))
            {
              cairo_move_to (context->cr, p0.x, p0.y);
              cairo_line_to (context->cr, p1.x, p1.y);
            }
        }

      boundary = boundary->next;
    }

  cairo_stroke (context->cr);
}


gboolean
voiro_expose_event (GtkWidget      *widget,
                    GdkEventExpose *event,
                    gpointer       *data)
{
  VoiroContext *context = (VoiroContext *) data;
  VoiroDiagram *diagram = context->diagram;
  gint i;

  voiro_draw_begin (context);

  if (diagram->convex && context->show_convex)
    voiro_draw_convex (context, diagram->convex, context->show_level);

  if (context->show_voronoi)
    voiro_draw_voronoi (context);

  for (i = 0; i < diagram->n_sites; i++)
    {
      voiro_draw_point (context, diagram->sites[i],
                        context->show_text ? i : -1,
                        context->dot_mode);
    }

  voiro_draw_end (context);

  return FALSE;
}


void
voiro_button_event (GtkWidget      *widget,
                    GdkEventButton *event,
                    gpointer       *data)
{
  VoiroContext *context = (VoiroContext *) data;
  gdouble x, y;

  x = event->x;
  y = event->y;

  voiro_draw_begin (context);
  cairo_device_to_user (context->cr, &x, &y);
  voiro_draw_end (context);
}


void
voiro_motion_notify (GtkWidget      *widget,
                     GdkEventMotion *event,
                     gpointer       *data)
{
  return;
}


gboolean
voiro_key_press_event (GtkWidget   *widget,
                       GdkEventKey *event,
                       gpointer    *data)
{
  VoiroContext *context = (VoiroContext *) data;
  gint i;

  switch (event->keyval)
    {
      case GDK_BackSpace:
        if (context->show_level == 0)
          context->show_level = context->diagram->max_level;
        else
          context->show_level = MAX (1, context->show_level - 1);
        gtk_widget_queue_draw (context->darea);
        break;
      case GDK_space:
        if (context->show_level >= context->diagram->max_level ||
            context->show_level == 0)
          context->show_level = 0;
        else
          context->show_level = MIN (context->diagram->max_level,
                                     context->show_level + 1);
        gtk_widget_queue_draw (context->darea);
        break;

      case GDK_D:
      case GDK_d:
        for (i = 0; i < context->diagram->n_sites; i++)
          g_printerr ("      diagram->sites[%d].x = %6.3f;\n"
                      "      diagram->sites[%d].y = %6.3f;\n",
                      i, context->diagram->sites[i].x,
                      i, context->diagram->sites[i].y);
        break;

      case GDK_F:
      case GDK_f:
      case GDK_F11:
        if (context->fullscreen)
          gtk_window_unfullscreen (GTK_WINDOW (context->window));
        else
          gtk_window_fullscreen (GTK_WINDOW (context->window));
        context->fullscreen = !context->fullscreen;
        break;

      case GDK_P:
      case GDK_p:
        context->dot_mode = (context->dot_mode + 1) % 3;
        gtk_widget_queue_draw (context->darea);
        break;

      case GDK_T:
      case GDK_t:
        context->show_text = !context->show_text;
        gtk_widget_queue_draw (context->darea);
        break;

      case GDK_C:
      case GDK_c:
        context->show_convex = !context->show_convex;
        gtk_widget_queue_draw (context->darea);
        break;

      case GDK_V:
      case GDK_v:
        context->show_voronoi = !context->show_voronoi;
        gtk_widget_queue_draw (context->darea);
        break;

      case GDK_Right:
        context->runmode = VOIRO_RESTART;
        if (context->showstep == 0 ||
            context->in_algorithm)
          context->showstep ++;
        gtk_main_quit ();
        gtk_widget_queue_draw (context->darea);
        break;
      case GDK_Left:
        context->runmode = VOIRO_RESTART;
        if (context->showstep > 0)
          context->showstep --;
        gtk_main_quit ();
        gtk_widget_queue_draw (context->darea);
        break;

      case GDK_Return:
        context->runmode = VOIRO_FINISH;
        gtk_main_quit ();
        gtk_widget_queue_draw (context->darea);
        break;

      case GDK_N:
      case GDK_n:
        context->runmode = VOIRO_NEW;
        context->showstep = 0;
        gtk_main_quit ();
        gtk_widget_queue_draw (context->darea);
        break;

      case GDK_Q:
      case GDK_q:
      case GDK_Escape:
        context->runmode = VOIRO_QUIT;
        gtk_main_quit ();
        break;
      default:
        return FALSE;
    }

  return TRUE;
}


void
voiro_create_window (VoiroContext *context)
{
  GtkWidget *vbox;

  context->window = g_object_new (GTK_TYPE_WINDOW,
                                  "title", "Voiro",
                                  "default-width", 800,
                                  "default-height", 600,
                                  NULL);

  g_signal_connect (context->window,
                    "destroy", G_CALLBACK (gtk_main_quit), NULL);

  vbox = g_object_new (GTK_TYPE_VBOX,
                       "visible", TRUE,
                       "parent", context->window,
                       NULL);

  context->darea = g_object_new (GTK_TYPE_DRAWING_AREA,
                                 "width-request", 100,
                                 "height-request", 100,
                                 "visible", TRUE,
                                 "parent", vbox,
                                 NULL);
  g_object_connect (G_OBJECT (context->darea),
                    "signal::expose-event",
                        G_CALLBACK (voiro_expose_event), context,
                    "signal::motion_notify_event",
                        G_CALLBACK (voiro_motion_notify), context,
                    "signal::button_press_event",
                        G_CALLBACK (voiro_button_event), context,
                    NULL);

  g_object_connect (G_OBJECT (context->window),
                    "signal::key_press_event",
                        G_CALLBACK (voiro_key_press_event), context,
                    NULL);

  gtk_widget_set_events (context->darea,
                         gtk_widget_get_events (context->darea)
                         | GDK_LEAVE_NOTIFY_MASK
                         | GDK_BUTTON_PRESS_MASK
                         | GDK_BUTTON_RELEASE_MASK
                         | GDK_KEY_PRESS_MASK
                         | GDK_POINTER_MOTION_MASK
                         | GDK_POINTER_MOTION_HINT_MASK);
}


gint
main (gint argc, char *argv[])
{
  VoiroContext *context;
  gint n_sites = 15;

  gtk_init (&argc, &argv);

  context = g_new0 (VoiroContext, 1);
  context->diagram = g_new0 (VoiroDiagram, 1);

  context->dot_mode = 2;
  context->show_text = FALSE;
  context->show_convex = TRUE;
  context->show_voronoi = TRUE;

  if (argc > 1)
    n_sites = atoi (argv[1]);

  voiro_create_window (context);
  voiro_algorithm_step (context);
  gtk_widget_show (context->window);

  voiro_init_model (context->diagram, n_sites);
  context->showstep = 0;

  while (context->runmode != VOIRO_QUIT)
    {
      context->runmode = VOIRO_STEPPING;
      gtk_widget_queue_draw (context->darea);

      context->step = 0;
      voiro_cleanup_model (context->diagram);
      voiro_algorithm_step (NULL);
      context->in_algorithm = TRUE;
      voiro_convex_hull (context->diagram,
                         context->diagram->sorted_order,
                         context->diagram->n_sites,
                         &context->diagram->convex);
      context->in_algorithm = FALSE;

      switch (context->runmode)
        {
          case VOIRO_STEPPING:
            gtk_main ();
            break;
          case VOIRO_FINISH:
            context->showstep = context->step;
            break;
          case VOIRO_RESTART:
            break;
          case VOIRO_NEW:
            voiro_init_model (context->diagram, n_sites);
            break;
          case VOIRO_QUIT:
            break;
        }
    }

  return 0;
}