/* See cirpack.h */
/* Last edited on 2005-02-04 07:11:54 by stolfi */ 

#include <i2.h>
#include <r2.h>
#include <interval.h>
#include <hexpaint.h>
#include <cpkbasic.h>
#include <pswr.h>
#include <pswr_color.h>
#include <pswr_iso.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>

#define MAX_CANVAS_SIZE 2048

hxp_canvas_t hxp_canvas_new(r2_t org, r2_t step, i2_t size, hxp_pixel_t v)
  { 
    int nx = size.c[0], ny = size.c[1];
    demand(nx <= MAX_CANVAS_SIZE, "canvas too big");
    demand(ny <= MAX_CANVAS_SIZE, "canvas too big");
    int nt = nx*ny;
    hxp_pixel_t *pix = (hxp_pixel_t *)notnull(malloc(nt*sizeof(hxp_pixel_t)), "no mem");
    int i;
    for (i = 0; i < nt; i++) { pix[i] = v; }
    return (hxp_canvas_t){ pix, org, step, size };
  }

r2_t hxp_pixel_pos(hxp_canvas_t *s, i2_t p)
  {
    int jx = p.c[0], iy = p.c[1]; 
    double dx = s->step.c[0], dy = s->step.c[1];
    double x = s->org.c[0] + (jx + (double)(iy & 1)/2.0)*dx;
    double y = s->org.c[1] + iy*dy;
    return (r2_t){{x,y}};
  }

int cmpy (r2_t *a, r2_t *b);
  /* Compares points {a} and {b} in Y-coordinate, then X: */

int cmpy (r2_t *a, r2_t *b)
  { 
    if (a->c[1] < b->c[1]) 
      { return -1; }
    else if (a->c[1] > b->c[1])
      { return +1; }
    else if (a->c[0] < b->c[0]) 
      { return -1; }
    else if (a->c[0] > b->c[0])
      { return +1; }
    else
      { return 0; }
  }

void hxp_paint_polygon
  ( r2_vec_t p,      /* Vertices of polygon, ccw around interior. */
    hxp_canvas_t *s, /* Where to paint. */
    hxp_pixel_t val  /* Value to paint. */
  )
  {
    /* Get hold of params: */
    int np = p.nel;      /* Number of vertices in polygon. */
    int nx = s->size.c[0]; /* Num of points per canvas row. */
    int ny = s->size.c[1]; /* Num of rows in canvas. */
    double ox = s->org.c[0], oy = s->org.c[1];    /* Lower left corner of grid. */
    double dx = s->step.c[0], dy = s->step.c[1];  /* Steps of canvas's grid. */

    /* Min and max scanline intercepted by polygon: */
    int iyminp = ny, iymaxp = -1;

    /* Make a list of all edge crossings for each scanline: */
    int ncross[ny];          /* Number of crossings in each scanline. */
    double_vec_t cross[ny];  /* Abscissas of crossings in each scanline. */
    int iy;
    for (iy = 0; iy < ny; iy++) { ncross[iy] = 0; }

    int ip;
    for (ip = 0; ip < np; ip++)
      { 
        /* Get an edge {u,v} of the polygon: */
        r2_t u = p.el[ip]; 
        r2_t v = p.el[(ip+1)%np]; 

        /* Make sure {u} is the lower endpoint: */
        if (u.c[1] > v.c[1]) { r2_t t = u; u = v; v = t; }

        /* Compute range {iymin..iymax} of scanlines crossed by edge. */
        /* We assume that the polygon is shifted down by epsilon. */
        int iymin = (int)ceil((u.c[1] - oy)/dy);
        int iymax = (int)ceil((v.c[1] - oy)/dy)-1;
        /* Intersect {iymin..iymax} with {0..ny-1}: */
        if (iymin < 0) { iymin = 0; }
        if (iymax > ny-1) { iymax = ny-1; }
        if (iymin <= iymax)
          { 
            /* Remember the polygon's Y extent: */
            if (iymin < iyminp) { iyminp = iymin; }
            if (iymax > iymaxp) { iymaxp = iymax; }
          
            /* Get edge's slope {Dx/Dy}: */
            double slope = (v.c[0] - u.c[0])/(v.c[1] - u.c[1]); 
            /* Add crossings to the crossings lists: */
            for (iy = iymin; iy <= iymax; iy++)
              { double y = oy + iy*dy;
                double x = u.c[0] + slope*(y - u.c[1]);
                int *nc = &(ncross[iy]);
                double_vec_t *cr = &(cross[iy]);
                if (*nc == 0) 
                  { *cr =  double_vec_new(6); }
                else
                  { double_vec_expand(cr, *nc); }
                cr->el[*nc] = x;
                (*nc)++;
              }
          }
      }

    /* Now fill each scanline according to its crossings: */
    for (iy = iyminp; iy <= iymaxp; iy++)
      { 
        /* Get the crossings of this scanline: */
        int nc = ncross[iy];
        assert(nc > 0); /* Polygon is connected. */
        assert(nc % 2 == 0); /* Jordan's curve theorem. */
        double *cx = cross[iy].el;
        /* Get pixels of this scanline: */
        hxp_pixel_t *pix = &(s->pix[iy*nx]); 
        /* Sort crossings (insertion sort should be fast enough): */
        int ic;
        for (ic = 1; ic < nc; ic++)
          { double cxi = cx[ic];
            int jc = ic; 
            while ((jc > 0) && (cx[jc-1] > cxi))
              { cx[jc] = cx[jc-1]; jc--; }
            cx[jc] = cxi;
          }
        /* Process crossing by even-odd rule: */
        for (ic = 0; ic < nc; ic += 2)
          { 
            /* Get next "inside" interval: */
            double xmin = cx[ic], xmax = cx[ic+1]; 
            /* Get range {ixmin..ixmax} of pixels contained in interval: */
            /* We assume that the polygon is shifted left by delta: */
            int ixmin = (int)ceil((xmin - ox)/dx - (iy & 1)/2.0);
            int ixmax = (int)ceil((xmax - ox)/dx - (iy & 1)/2.0) - 1;
            /* Intersect {ixmin..ixmax} with {0..nx-1}: */
            if (ixmin < 0) { ixmin = 0; }
            if (ixmax > nx-1) { ixmax = nx-1; }
            /* Set pixels: */
            int ix; 
            for (ix = ixmin; ix <= ixmax; ix++) { pix[ix] = val; }
          }
        /* Free crossings list: */
        free(cross[iy].el);
      }
  }       
  
void hxp_paint_circle
  ( r2_t ctr,        /* Center of circle. */
    double rad,      /* Radius of circle. */
    hxp_canvas_t *s, /* Where to paint. */
    hxp_pixel_t val  /* Value to paint. */
  )
  {
    /* Get hold of params: */
    int nx = s->size.c[0]; /* Num of points per canvas row. */
    int ny = s->size.c[1]; /* Num of rows in canvas. */
    double ox = s->org.c[0], oy = s->org.c[1];    /* Lower left corner of grid. */
    double dx = s->step.c[0], dy = s->step.c[1];  /* Steps of canvas's grid. */
    double rad2 = rad*rad;

    /* Compute range {iymin..iymax} of scanlines crossed by circle. */
    double ymin = ctr.c[1] - rad;
    double ymax = ctr.c[1] + rad;
    
    /* We assume that the circle is fattened by epsilon. */
    int iymin = (int)ceil((ymin - oy)/dy);
    int iymax = (int)floor((ymax - oy)/dy);
    /* Intersect {imin..imax} with {0..ny-1}: */
    if (iymin < 0) { iymin = 0; }
    if (iymax > ny-1) { iymax = ny-1; }
    if (iymin > iymax) { return; }
    /* Enumerate relevant scanlines: */
    int iy;
    for (iy = iymin; iy <= iymax; iy++)
      { /* Get pixels of this scanline: */
        hxp_pixel_t *pix = &(s->pix[iy*nx]); 
        /* Find half-extend {hx} of circle on scanling {iy}: */
        double y = oy + iy*dy, y2 = y*y;
        double hx = (y2 > rad2 ? 0.0 : sqrt(rad2 - y2));
        double xmin = ctr.c[0] - hx, xmax = ctr.c[0] + hx; 
        /* Get range {ixmin..(ixmax+1)-1} of pixels contained in interval: */
        /* We assume that the polygon is shifted left by delta: */
        int ixmin = (int)ceil((xmin - ox)/dx - (iy & 1)/2.0);
        int ixmax = (int)floor((xmax - ox)/dx - (iy & 1)/2.0);
        /* Intersect {ixmin..ixmax} with {0..nx-1}: */
        if (ixmin < 0) { ixmin = 0; }
        if (ixmax > nx-1) { ixmax = nx-1; }
        int ix; 
        for (ix = ixmin; ix < ixmax; ix++) { pix[ix] = val; }
      }
  }

void hxp_paint_stroke
  ( r2_t a,          /* Vertex of triangle. */
    r2_t b,          /* Vertex of triangle. */
    double r,        /* Half-width of stroke. */
    hxp_canvas_t *s, /* Where to paint. */
    hxp_pixel_t val  /* Value to paint. */
  )
  {
    /* Zero-area figures contain no pixels: */
    if (r == 0) { return; }
    /* Get the displacement vector from {a} to {b}: */
    r2_t d; r2_sub(&b, &a, &d);
    /* Zero-area figures contain no pixels: */
    if ((d.c[0] == 0) && (d.c[1] == 0)) { return; }
    /* Get its direction: */
    r2_dir(&d, &d);
    /* Get a vector {v} with length {r} perpendicular to {d}: */
    r2_t v = (r2_t){{-r*d.c[1], r*d.c[0]}};
    /* Build the polygon: */
    r2_vec_t p = r2_vec_new(4);
    r2_add(&a, &v, &(p.el[0]));
    r2_add(&b, &v, &(p.el[1]));
    r2_sub(&b, &v, &(p.el[2]));
    r2_sub(&a, &v, &(p.el[3]));
    /* Plot the polygon: */
    hxp_paint_polygon(p, s, val);
    free(p.el);
  }

void hxp_plot_canvas
  ( PSStream *ps,     /* Where to plot. */
    hxp_canvas_t *s,  /* The canvas to plot. */
    double rad(hxp_pixel_t v) /* Radius function. */
  )
  { 
    /* Get hold of params: */
    int nx = s->size.c[0]; /* Num of points per canvas row. */
    int ny = s->size.c[1]; /* Num of rows in canvas. */

    /* Draw pixels: */
    pswr_set_pen(ps, 0.000,0.000,0.000, 0.10, 0.0, 0.0);
    pswr_set_fill_color(ps, 0.000,0.250,0.000);
    int ix, iy;
    for (iy = 0; iy < ny; iy++)
      { for (ix = 0; ix < nx; ix++)
          { r2_t ctr = hxp_pixel_pos(s, (i2_t){{ix,iy}});
            double r = rad(s->pix[iy*nx + ix]);
            if (r > 0) 
              { pswr_circle(ps, ctr.c[0], ctr.c[1], r, TRUE, TRUE); }
          }
      }

  }

i2_vec_t hxp_list_pixels (hxp_canvas_t *s, hxp_pixel_t v)
  {
    /* Count pixels that are equal to {v}: */
    int nx = s->size.c[0], ny = s->size.c[1]; 
    int nv = 0;
    int iy;
    for (iy = 0; iy < ny; iy++)
      { hxp_pixel_t *pix = &(s->pix[iy*nx]); 
        int ix;
        for (ix = 0; ix < nx; ix++)
          { if (pix[ix] == v) { nv++; } }
      }

    /* List their positions: */
    i2_vec_t K = i2_vec_new(nv);
    nv = 0; 
    for (iy = 0; iy < ny; iy++)
      { hxp_pixel_t *pix = &(s->pix[iy*nx]); 
        int ix;
        for (ix = 0; ix < nx; ix++)
          { if (pix[ix] == v) 
              { K.el[nv] = (i2_t){{ix,iy}}; nv++; }
          }
      }
    return K;
  }