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

#include <cirpack.h>
#include <hexpaint.h>
#include <interval.h>
#include <i2.h>
#include <r2.h>
#include <r2x2.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

r2_vec_t cpk_map_points(r2_vec_t P, r2_t *o, r2x2_t *M);
  /* Maps each point of {P} by subtracting {o} from it,
    and multiplying the result (as a row vector) by
    the matrix {M}.  Infinite points are returned unchanged. */
    
void cpk_get_bbox(r2_vec_t R, interval_t B[]);
  /* Stores in {B[0]} and {B[1]} the X and Y ranges (i.e. the bounding
   box) of all points {R[0..]}. Ignores points that have infinite
   coords. If there are no finite points, returns an empty box. */

r2_vec_t cpk_get_valid_points_canon
  ( double step,    /* Grid step. */
    r2_vec_t R,     /* Vertices of polygon, ccw around interior. */
    r2_vec_t P,     /* Forbidden points. */
    double dP,      /* Minimum distance from forbidden points. */
    r2_vec_t L,     /* Vertices of forbidden polygonal lines. */
    double dL       /* Minimum distance from forbidden lines. */
  );
  /* Like {cpk_get_valid_points}, but uses the canonical hex grid with
    the given {step} that includes the origin and has nearest
    neighbors in the X direction. */
  
r2_vec_t cpk_get_valid_points
  ( r2_t o,         /* Grid's reference point. */
    r2_t u,         /* Grid's X axis. */
    double step,    /* Grid step. */
    r2_vec_t R,     /* Vertices of polygon, ccw around interior. */
    r2_vec_t P,     /* Forbidden points. */
    double dP,      /* Minimum distance from forbidden points. */
    r2_vec_t L,     /* Vertices of forbidden polygonal lines. */
    double dL       /* Minimum distance from forbidden lines. */
  )
  {
    /* Normalize {u} to unit length: */
    r2_dir(&u, &u);
    /* Get the direction {v} orthogonal to {u}, ccw from it: */
    double ux = u.c[0], uy = u.c[1];
    double vx = -uy, vy = ux;
    
    /* Rewrite the polygon and the fixed sites in the {o,u,v} coord system: */
    r2x2_t M = (r2x2_t){{{ux,vx},{uy,vy}}};
    r2_vec_t RM = cpk_map_points(R, &o, &M);
    r2_vec_t PM = cpk_map_points(P, &o, &M);
    r2_vec_t LM = cpk_map_points(L, &o, &M);
    
    /* Find valid points in this system: */
    r2_vec_t VM = cpk_get_valid_points_canon(step, RM, PM, dP, LM, dL);
    
    /* Unmap the solution: */
    M = (r2x2_t){{{ux,uy},{vx,vy}}};
    r2_neg(&o, &o);
    r2_vec_t V = cpk_map_points(VM, &o, &M);
    
    /* Housecleaning: */
    free(RM.el); free(PM.el); free(LM.el); free(VM.el);
    
    return V;
  }

r2_vec_t cpk_map_points(r2_vec_t P, r2_t *o, r2x2_t *M) 
  { r2_vec_t Q = r2_vec_new(P.nel);
    int k;
    for (k = 0; k < P.nel; k++) 
      { r2_t *Pk = &(P.el[k]);
        r2_t *Qk = &(Q.el[k]);
        if ((fabs(Pk->c[0]) == INF) || (fabs(Pk->c[1]) == INF))
          { *Qk = *Pk; }
        else
          { r2_t d; r2_sub(Pk, o, &d); 
            r2x2_map_row(&d, M, Qk); 
          }
      }
    return Q;
  }    
  
r2_vec_t cpk_get_valid_points_canon
  ( double step,    /* Grid step. */
    r2_vec_t R,     /* Vertices of polygon, ccw around interior. */
    r2_vec_t P,     /* Forbidden points. */
    double dP,      /* Minimum distance from forbidden points. */
    r2_vec_t L,     /* Vertices of forbidden polygonal lines. */
    double dL       /* Minimum distance from forbidden lines. */
  )
  {
    /* Grid steps in {x} and {y}: */
    double dy = step; 
    double dx = step*sqrt(3)/2;
    
    /* Find the polygon's axis-aligned bounding box and its sizes: */
    interval_t r[2]; /* Range of {p}'s {u} and {v} coords. */
    cpk_get_bbox(R, r);
    
    /* Round the corners to grid points on even rows: */
    int ixlo = (int)ceil(lo(r[0])/dx);
    int iylo = 2*(int)ceil(lo(r[1])/dy/2);
    
    int ixhi = (int)ceil(hi(r[0])/dx);
    int iyhi = 2*(int)ceil(hi(r[1])/dy/2);
    
    /* Compute the grid dimensions: */
    int nx = ixhi - ixlo;
    int ny = iyhi - iylo;
    
    /* Allocate the grid: */
    hxp_canvas_t s = hxp_canvas_new
      ( (r2_t){{ixlo*dx, iylo*dy}},
        (r2_t){{dx,dy}},
        (i2_t){{nx,ny}},
        0
      );
        
    /* Paint the polygon over the grid: */
    hxp_paint_polygon(R, &s, 1);
    
    /* Erase all pixels near the forbidden points: */
    int i;
    for (i = 0; i < P.nel; i++)
      { hxp_paint_circle(P.el[i], dP, &s, 0); }
    
    /* Erase all pixels near the forbidden lines: */
    r2_t *a = NULL;
    for (i = 0; i <= L.nel; i++)
      { 
        r2_t *b = &(L.el[i]); 
        if ((fabs(a->c[0]) != INF) && (fabs(a->c[1]) != INF)) 
          { hxp_paint_circle(*b, dL, &s, 0); 
            if (a != NULL) 
              { hxp_paint_stroke(*a, *b, dL, &s, 0); }
            a = b;
          }
        else
        { a = NULL; }
      }

    /* Get all pixels that are still {1}. */
    i2_vec_t K = hxp_list_pixels(&s, 1);
    
    /* Get their positions */
    r2_vec_t V = r2_vec_new(K.nel);
    for (i = 0; i < K.nel; i++)
      { V.el[i] = hxp_pixel_pos(&s, K.el[i]); }

    /* Housecleaning: */
    free(s.pix); free(K.el);
    
    return V;
  }

i2_vec_t cpk_incompat_edges (r2_vec_t V, double d_min)
  { /* Should use a bucket grid. */
    /* Slow and dumb for now: */
    i2_vec_t e = i2_vec_new(V.nel);
    int ne = 0;
    int j;
    for (j = 1; j < V.nel; j++)
      { r2_t *Vj = &(V.el[j]);
        int i;
        for (i = 0; i < j; i++)
          { r2_t *Vi = &(V.el[i]);
            if (r2_dist(Vi, Vj) < d_min)
              { i2_vec_expand(&e, ne);
                e.el[ne] = (i2_t){{i,j}};
                ne++;
              }
          }
      }
    i2_vec_trim(&e, ne);
    return e;
  }
    
void cpk_get_bbox(r2_vec_t R, interval_t B[])
  {
    B[0] = (interval_t){{INF,-INF}};
    B[1] = (interval_t){{INF,-INF}};
    int i;
    for (i = 0; i < R.nel; i++)
      {
        r2_t *Ri = &(R.el[i]);
        double x = Ri->c[0], y = Ri->c[1];
        if ((fabs(x) != INF) && (fabs(y) != INF)) 
          {   
            if (x < lo(B[0])) { lo(B[0]) = x; }
            if (x > hi(B[0])) { hi(B[0]) = x; }
            if (y < lo(B[1])) { lo(B[1]) = y; }
            if (y > hi(B[1])) { hi(B[1]) = y; }
          }
      }
  }