#include "project.h"

#ifndef PI
#define PI 3.141592654
#endif


typedef struct {
  int             x, y;
}               Point;

typedef struct contour_struct {
  float          *data;
  float         **dcol;
  unsigned char  *flags;
  unsigned char **fcol;
  int             x, y;
  int             x1, y1;
  Jwgline         current_line;
  int             current_line_size;
  int             lp;
  Handle          h;
}              *Contour;


#define INRANGE(l1,l2,c) ( ((c)<(l2)) && ((c)>=(l1)) )
#define RANGEX(c,x)	( ((x)>0) ? ( ((x)<((c)->x1)) ? (x):((c)->x1)) : 1)
#define RANGEY(c,y)	(((y)>0) ?  ( ((y)<((c)->y1)) ? (y):((c)->y1)) : 1)

#define FLAGS(c,p) ((c->fcol[p.y])[p.x])
#define DATA(c,p) ((c->dcol[p.y])[p.x])

#define DEF_LINE_SIZE 8192





static void 
cont_reset_point(Contour c)
{
  c->current_line.npts = 0;
}

static void 
cont_add_point(Contour c, Jwgpos pos)
{
  Jwgpos         *ptr;
  int             n;

  if (c->current_line.npts == c->current_line_size) {
    /* uh-oh we're about to over flow */

    ptr = c->current_line.data;


    n = c->current_line_size;
    printf("JWG: CONTOUR: Increasing maximum line size from %d to %d\n",
	   n, n * 2);

    c->current_line_size = n * 2;

    c->current_line.data = malloc(sizeof(Jwgpos) * n * 2);
    bcopy(ptr, c->current_line.data, sizeof(Jwgpos) * n);

    free(ptr);

  }
  jwg_xform(c->h, &pos);

  c->current_line.data[c->current_line.npts] = pos;
  c->current_line.npts++;

}

static unsigned char
find_direction(Point p1, Point p2)
{
  p2.x -= p1.x;
  p2.y -= p1.y;

  switch (p2.y) {
  case -1:
    switch (p2.x) {
    case -1:
      return (1);
    case 0:
      return (2);
    case 1:
      return (4);
    }
  case 0:
    switch (p2.x) {
    case -1:
      return (128);
    case 1:
      return (8);
    }
  case 1:
    switch (p2.x) {
    case -1:
      return (64);
    case 0:
      return (32);
    case 1:
      return (16);
    }
  }


  fprintf(stderr, "Contouring cockup type 3 ... (bye bye)\n");
  abort();

}

static Point
add_direction(Point p, int d)
{
  Point           dirs[8] = {{-1, -1}, {0, -1}, {1, -1}, {1, 0}, {1, 1}, {0, 1}, {-1, 1}, {-1, 0}};

  p.x = dirs[d].x;
  p.y = dirs[d].y;

  return (p);
}


static int
points_marked(Contour c, Point p1, Point p2)
{
  return (c->fcol[p1.y])[p1.x] & find_direction(p1, p2);
}

static void
mark_points(Contour c, Point p1, Point p2)
{
  (c->fcol[p1.y])[p1.x] |= find_direction(p1, p2);
  (c->fcol[p2.y])[p2.x] |= find_direction(p2, p1);
}



static void trace(Contour c, Point p1, Point p2, Point pl, float level)
{
  Point           pn;
  float           v1, v2, vn;
  Jwgpos          pos;
  float           r;

printf(".");
fflush(stdout);

  if (!c->lp)
  {
  pos.x = 1;
  pos.y = 1;

  cont_add_point(c, pos);
  }

  v1 = DATA(c, p1);
  v2 = DATA(c, p2);

  while (1) {

    /*
     * printf("(%d,%d,%f) (%d,%d,%f) (%d,%d,%f)\n", p1.x, p1.y, DATA(c, p1),
     * p2.x, p2.y, DATA(c, p2), pl.x, pl.y, DATA(c, pl));
     */

    r = (v2 - v1);
    if (r == 0.0) {
      r = 0.5;
    } else {
      r = (level - v1) / r;
    }

    pos.x = r * (float) RANGEX(c, p2.x);
    pos.y = r * (float) RANGEY(c, p2.y);

    r = 1.0 - r;

    pos.x += r * (float) RANGEX(c, p1.x);
    pos.y += r * (float) RANGEY(c, p1.y);

    cont_add_point(c, pos);







    if (DATA(c, p1) > level) {
      printf("Contouring cockup type 1 -- bye bye\n");
      abort();
    }
    if (DATA(c, p2) <= level) {
      printf("Contouring cockup type 2 -- bye bye\n");
      abort();
    }
    if (points_marked(c, p1, p2)) {
      if (c->lp) {
	jwg_write_line(c->h, &c->current_line);
	c->current_line.npts = 0;
      } else {
	pos.x = 1;
	pos.y = 1;

	cont_add_point(c, pos);
      }


      return;
    }
    mark_points(c, p1, p2);

    pn.x = p1.x + p2.x - pl.x;
    pn.y = p1.y + p2.y - pl.y;

    vn = DATA(c, pn);

    if (vn <= level) {
      pl = p1;
      p1 = pn;
      v1 = vn;
    } else {
      pl = p2;
      p2 = pn;
      v2 = vn;
    }




  }

}





static void possible(Contour c, int x1, int y1, int x2, int y2, int xl, int yl, float level)
{
  Point           p1, p2, pl;

  p1.x = x1;
  p1.y = y1;
  p2.x = x2;
  p2.y = y2;
  pl.x = xl;
  pl.y = yl;

  if (points_marked(c, p1, p2))
    return;

  trace(c, p1, p2, pl, level);


}


static void contour(Contour c, float level)
{
  int             x, y;
  float         **c0, **c1;
  float          *r00, *r10, *r11, *r01;

  c1 = c0 = c->dcol;
  c1--;

  for (y = 0; y < c->y; ++y) {
    r00 = r01 = *(c0);
    r01--;
    if (y) {
      r10 = r11 = *(c1);
      r11--;
    }
    for (x = 0; x < c->x; ++x) {

      if (x) {
	if (INRANGE(*r00, *r01, level))
	  possible(c, x, y, x - 1, y, x, y - 1, level);
	if (INRANGE(*r01, *r00, level))
	  possible(c, x - 1, y, x, y, x - 1, y + 1, level);
      }
      if (y) {
	if (INRANGE(*r00, *r10, level))
	  possible(c, x, y - 1, x, y, x - 1, y, level);
	if (INRANGE(*r10, *r00, level))
	  possible(c, x, y, x, y - 1, x + 1, y - 1, level);
      }
      r00++;
      r10++;
      r11++;
      r01++;
    }
    c0++;
    c1++;
  }

}



INTERNAL void
do_contour(Handle h, float level, int side, int lp)
{
  Contour         c = (Contour) malloc(sizeof(struct contour_struct));
  int             i;
  float           hug = 1.0e10 * (float) side;

 
  c->data = h->data->data;
  c->x1 = c->x = h->data->w;
  c->y1 = c->y = h->data->h;

 
  c->dcol = (float **) malloc(sizeof(float *) * c->y);



  c->dcol[0] = c->data;
  for (i = 1; i < c->y; ++i)
    c->dcol[i] = c->dcol[i - 1] + c->x;

  c->flags = (unsigned char *) malloc(c->x * c->y);
  bzero(c->flags, c->x * c->y);

  c->fcol = (unsigned char **) malloc(c->y * sizeof(unsigned char *));

  c->fcol[0] = c->flags;
  for (i = 1; i < c->y; ++i)
    c->fcol[i] = c->fcol[i - 1] + c->x;

  c->x1--;
  c->y1--;

  for (i = 0; i < c->y; ++i) {
    (c->dcol[i])[0] = hug;
    (c->dcol[i])[c->x1] = hug;
  }
  for (i = 0; i < c->x; ++i) {
    (c->dcol[0])[i] = hug;
    (c->dcol[c->y1])[i] = hug;
  }

  c->x1--;
  c->y1--;

  /*
   * { Point p;
   * 
   * for (p.y=0;p.y<c->y;++p.y) { for (p.x=0;p.x<c->x;++p.x) {
   * 
   * printf("%.1f ",DATA(c,p));
   * 
   * } printf("\n");
   * 
   * }
   * } */

  c->current_line.npts = 0;
  c->current_line_size = DEF_LINE_SIZE;
  c->current_line.data = (Jwgpos *) malloc(sizeof(Jwgpos) * DEF_LINE_SIZE);

  c->h = h;
  c->lp = lp;

  contour(c, level);
  printf("\n");


  if (!c->lp) 
     jwg_write_polygon(c->h, &c->current_line);

  free(c->flags);
  free(c->fcol);
  free(c->dcol);
  free(c->current_line.data);
  
}