/*
 * Routines for bitmap manipulation
 *
  gsumi version 0.5

  Copyright 1997 Owen Taylor <owt1@cornell.edu>

  Heavily based on

  xink version 0.02

  Copyright 1997 Raph Levien <raph@acm.org>

  This code is free for commercial and non-commercial use or
  redistribution, as long as the source code release, startup screen,
  or product packaging includes this copyright notice.
*/

#include "gsumi.h"
#include "rect.h"

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

#define ROUND(A) floor(A+0.5)

static Blob *
blob_new (int y, int height)
{
  Blob *result;

  result = malloc (sizeof (Blob) +  sizeof(BlobSpan) * height);
  result->y = y;
  result->height = height;

  return result;
}

typedef enum {
  NONE = 0,
  LEFT = 1 << 0,
  RIGHT = 1 << 1,
} EdgeType;

Blob *
blob_convex_union (Blob *b1, Blob *b2)
{
  Blob *result;
  int y, x1, x2, y1, y2, i1, i2;
  int i, j;
  int start;
  EdgeType *present;

  /* Create the storage for the result */

  y = MIN(b1->y,b2->y);
  result = blob_new (y, MAX(b1->y+b1->height,b2->y+b2->height)-y);

  if (result->height == 0)
    return result;

  present = malloc (sizeof (EdgeType) * result->height);
  memset (present, 0, result->height * sizeof(EdgeType));

  /* Initialize spans from original objects */

  for (i=0, j=b1->y-y; i<b1->height; i++,j++)
    {
      if (b1->data[i].right >= b1->data[i].left)
	{
	  present[j] = LEFT | RIGHT;
	  result->data[j].left = b1->data[i].left;
	  result->data[j].right = b1->data[i].right;
	}
    }

  for (i=0, j=b2->y-y; i<b2->height; i++,j++)
    {
      if (b2->data[i].right >= b2->data[i].left)
	{
	  if (present[j])
	    {
	      if (result->data[j].left > b2->data[i].left)
		result->data[j].left = b2->data[i].left;
	      if (result->data[j].right < b2->data[i].right)
		result->data[j].right = b2->data[i].right;
	    }
	  else
	    {
	      present[j] = LEFT | RIGHT;
	      result->data[j].left = b2->data[i].left;
	      result->data[j].right = b2->data[i].right;
	    }
	}
    }

  /* Now walk through edges, deleting points that aren't on convex hull */

  start = 0;
  while (!(present[start])) start++;

  /*    left edge */

  i1 = start-1; 
  i2 = start;
  x1 = result->data[start].left - result->data[start].right;
  y1 = 0;
  
  for (i=start+1;i<result->height;i++)
    {
      if (!(present[i] & LEFT))
	continue;

      x2 = result->data[i].left - result->data[i2].left;
      y2 = i-i2;
      
      while (x2*y1 - x1*y2 < 0) /* clockwise rotation */
	{
	  present[i2] &= ~LEFT;
	  i2 = i1;
	  while (!(present[--i1] & LEFT) && i1>=start);

	  if (i1<start)
	    {
	      x1 = result->data[start].left - result->data[start].right;
	      y1 = 0;
	    }
	  else
	    {
	      x1 = result->data[i2].left - result->data[i1].left;
	      y1 = i2 - i1;
	    }
	  x2 = result->data[i].left - result->data[i2].left;
	  y2 = i - i2;
	}
      x1 = x2;
      y1 = y2;
      i1 = i2;
      i2 = i;
    }

  /*     Right edge */

  i1 = start -1; 
  i2 = start;
  x1 = result->data[start].right - result->data[start].left;
  y1 = 0;
  
  for (i=start+1;i<result->height;i++)
    {
      if (!(present[i] & RIGHT))
	continue;

      x2 = result->data[i].right - result->data[i2].right;
      y2 = i-i2;
      
      while (x2*y1 - x1*y2 > 0) /* counter-clockwise rotation */
	{
	  present[i2] &= ~RIGHT;
	  i2 = i1;
	  while (!(present[--i1] & RIGHT) && i1>=start);

	  if (i1<start)
	    {
	      x1 = result->data[start].right - result->data[start].left;
	      y1 = 0;
	    }
	  else
	    {
	      x1 = result->data[i2].right - result->data[i1].right;
	      y1 = i2 - i1;
	    }
	  x2 = result->data[i].right - result->data[i2].right;
	  y2 = i - i2;
	}
      x1 = x2;
      y1 = y2;
      i1 = i2;
      i2 = i;
    }

  /* Restore edges of spans that were deleted in last step or never present */

  /* We fill only interior regions of convex hull, as if we were filling
     polygons. But since we draw ellipses with nearest points, not interior
     points, maybe it would look better if we did the same here. Probably
     not a big deal either way after anti-aliasing */

  /*     left edge */
  for (i1=start; i1<result->height-1; i1++)
    {
      /* Find empty gaps */
      if (!(present[i1+1] & LEFT))
	{
	  int increment;	/* fractional part */
	  int denom;		/* denominator of fraction */
	  int step;		/* integral step */
	  int frac;		/* fractional step */
	  int reverse;

	  /* find bottom of gap */
	  i2 = i1+2;
	  while (!(present[i2] & LEFT) && i2 < result->height) i2++;
	  
	  if (i2 < result->height)
	    {
	      denom = i2-i1;
	      x1 = result->data[i1].left;
	      x2 = result->data[i2].left;
	      step = (x2-x1)/denom;
	      frac = x2-x1 - step*denom;
	      if (frac < 0)
		{
		  frac = -frac;
		  reverse = 1;
		}
	      else
		reverse = 0;
	      
	      increment = 0;
	      for (i=i1+1; i<i2; i++)
		{
		  x1 += step;
		  increment += frac;
		  if (increment >= denom)
		    {
		      increment -= denom;
		      x1 += reverse ? -1 : 1;
		    }
		  if (increment == 0 || reverse)
		    result->data[i].left = x1;
		  else
		    result->data[i].left = x1 + 1;
		}
	    }
	  i1 = i2-1;		/* advance to next possibility */
	}
    }

  /*     right edge */
  for (i1=start; i1<result->height-1; i1++)
    {
      /* Find empty gaps */
      if (!(present[i1+1] & RIGHT))
	{
	  int increment;	/* fractional part */
	  int denom;		/* denominator of fraction */
	  int step;		/* integral step */
	  int frac;		/* fractional step */
	  int reverse;

	  /* find bottom of gap */
	  i2 = i1+2;
	  while (!(present[i2] & RIGHT) && i2 < result->height) i2++;
	  
	  if (i2 < result->height)
	    {
	      denom = i2-i1;
	      x1 = result->data[i1].right;
	      x2 = result->data[i2].right;
	      step = (x2-x1)/denom;
	      frac = x2-x1 - step*denom;
	      if (frac < 0)
		{
		  frac = -frac;
		  reverse = 1;
		}
	      else
		reverse = 0;
	      
	      increment = 0;
	      for (i=i1+1; i<i2; i++)
		{
		  x1 += step;
		  increment += frac;
		  if (increment >= denom)
		    {
		      increment -= denom;
		      x1 += reverse ? -1 : 1;
		    }
		  if (reverse && increment != 0)
		    result->data[i].right = x1 - 1;
		  else
		    result->data[i].right = x1;
		}
	    }
	  i1 = i2-1;		/* advance to next possibility */
	}
    }
  
  /* Mark empty lines at top and bottom as unused */
  for (i=0;i<start;i++)
    {
      result->data[i].left = 0;
      result->data[i].right = -1;
    }
  for (i=result->height-1;!present[i];i--)
    {
      result->data[i].left = 0;
      result->data[i].right = -1;
    }

  free (present);
  return result;
}

/****************************************************************
 * Code to scan convert an arbitrary ellipse into a Blob. Based
 * on Van Aken's conic algorithm in Foley and Van Damn 
 ****************************************************************/

/* Return octant from gradient */
static int
blob_get_octant (int D, int E)
{
  if (D>=0)
    {
    if (E<0)
      return (D<-E) ? 1 : 2;
    else
      return (D>E) ? 3 : 4;
    }
  else
    if (E>0)
      return (-D<E) ? 5 : 6;
    else
      return (-D>-E) ? 7 : 8;
}

static void
blob_conic_add_pixel (Blob *b, EdgeType *present, int x, int y, int octant)
{
  /*  printf ("%d %d\n",x,y); */
  if (y<b->y || y>=b->y+b->height)
    {
      /*      g_warning("Out of bounds!\n"); */
    }
  else
    {
      if (octant <= 4)
	{
	  if (present[y-b->y] & RIGHT)
	    b->data[y-b->y].right = MAX(b->data[y-b->y].right,x);
	  else
	    {
	      b->data[y-b->y].right = x;
	      present[y-b->y] |= RIGHT;
	    }
	}
      else
	{
	  if (present[y-b->y] & LEFT)
	    b->data[y-b->y].left = MIN(b->data[y-b->y].left,x);
	  else
	    {
	      b->data[y-b->y].left = x;
	      present[y-b->y] |= LEFT;
	    }
	}
    }
}

static void
blob_conic (Blob *b, int xs, int ys,
	    int A, int B, int C, int D, int E, int F)
{
  int x,y;			/* current point */
  int octant;			/* current octant */
  int dxsquare, dysquare;	/* change in (x,y) for square moves */
  int dxdiag, dydiag;		/* change in (x,y) for diagonal moves */
  int d,u,v,k1,k2,k3;		/* decision variables and increments */
  int octantCount;		/* number of octants to be drawn */
  int count;			/* number of steps for last octant */
  int fx,dfx;			/* check x-boundaries */
  int fy,dfy;			/* check y-boundaries */
  int ddf;			/* second difference for boundary-checking */
  int tmp, i;

  EdgeType *present;
  
  present = malloc (sizeof (EdgeType) * b->height);
  memset (present, 0, b->height * sizeof(EdgeType));

  octant = blob_get_octant (D,E);

  switch (octant)
    {
    case 1:
      d = ROUND (A+B/2.+C/4.+D+E/2.+F);
      u = ROUND (A+B/2.+D);
      v = ROUND (A+B/2.+D+E);
      k1 = 2*A;
      k2 = 2*A + B;
      k3 = k2 + B + 2*C;
      dxsquare = 1;
      dysquare = 0;
      dxdiag = 1;
      dydiag = 1;
      break;
    case 2:
      d = ROUND (A/4.+B/2.+C+D/2.+E+F);
      u = ROUND (B/2.+C+E);
      v = ROUND (B/2.+C+D+E);
      k1 = 2*C;
      k2 = B + 2*C;
      k3 = 2*A + 2*B + 2*C;
      dxsquare = 0;
      dysquare = 1;
      dxdiag = 1;
      dydiag = 1;
      break;
    case 3:
      d = ROUND (A/4.-B/2.+C-D/2.+E+F);
      u = ROUND (-B/2.+C+E);
      v = ROUND (-B/2.+C-D+E);
      k1 = 2*C;
      k2 = 2*C - B;
      k3 = 2*A - 2*B + 2*C;
      dxsquare = 0;
      dysquare = 1;
      dxdiag = -1;
      dydiag = 1;
      break;
    case 4:
      d = ROUND (A-B/2.+C/4.-D+E/2.+F);
      u = ROUND (A-B/2.-D);
      v = ROUND (A-B/2.-D+E);
      k1 = 2*A;
      k2 = 2*A - B;
      k3 = k2 - B + 2*C;
      dxsquare = -1;
      dysquare = 0;
      dxdiag = -1;
      dydiag = 1;
      break;
    case 5:
      d = ROUND (A+B/2.+C/4.-D-E/2.+F);
      u = ROUND (A+B/2.-D);
      v = ROUND (A+B/2.-D-E);
      k1 = 2*A;
      k2 = 2*A + B;
      k3 = k2 + B + 2*C;
      dxsquare = -1;
      dysquare = 0;
      dxdiag = -1;
      dydiag = -1;
      break;
    case 6:
      d = ROUND (A/4.+B/2.+C-D/2.-E+F);
      u = ROUND (B/2.+C-E);
      v = ROUND (B/2.+C-D-E);
      k1 = 2*C;
      k2 = B + 2*C;
      k3 = 2*A + 2*B + 2*C;
      dxsquare = 0;
      dysquare = -1;
      dxdiag = -1;
      dydiag = -1;
      break;
    case 7:
      d = ROUND (A/4.-B/2.+C+D/2.-E+F);
      u = ROUND (-B/2.+C-E);
      v = ROUND (-B/2.+C+D-E);
      k1 = 2*C;
      k2 = 2*C - B;
      k3 = 2*A - 2*B + 2*C;
      dxsquare = 0;
      dysquare = -1;
      dxdiag = 1;
      dydiag = -1;
      break;
    default:			/* case 8: */
      d = ROUND (A-B/2.+C/4.+D-E/2.+F);
      u = ROUND (A-B/2.+D);
      v = ROUND (A-B/2.+D-E);
      k1 = 2*A;
      k2 = 2*A - B;
      k3 = k2 - B + 2*C;
      dxsquare = 1;
      dysquare = 0;
      dxdiag = 1;
      dydiag = -1;
      break;
    }

  octantCount = 8;
  x = xs;
  y = ys;
  count = 0;			/* ignore until last octant */

  /* Initialize boundary checking - we keep track of the discriminants
     for the conic as quadratics in x and y, and when they go negative
     we know we are beyond the boundaries of the conic. */

  ddf = 2*(B*B - 4*A*C);
  fx = E*E - 4*C*F;
  dfx = ddf/2 + dxdiag * (2*B*E - 4*C*D);
  fy = D*D - 4*A*F;
  dfy = ddf/2 + dydiag * (2*B*D - 4*A*E);

  while (1)
    {
      if (octantCount == 0)
	{
	  /* figure out remaining steps in square direction */
	  switch (octant)
	    {
	    case 1:
	    case 8:
	      count = xs - x;
	      break;
	    case 2:
	    case 3:
	      count = ys - y;
	      break;
	    case 4:
	    case 5:
	      count = x - xs;
	      break;
	    case 6:
	    case 7:
	      count = y - ys;
	      break;
	    }
	  /*	  if (count < 0)
	    g_warning("Negative count (%d) in octant %d\n",count,octant); */
	  if (count <= 0)
	    goto done;
	     
	}
      if (octant %2)	/* odd octants */
	{
	  while (v < k2/2)
	    {
	      blob_conic_add_pixel (b, present, x, y, octant);
	      if (d<0)
		{
		  x += dxsquare; y += dysquare;
		  u += k1;
		  v += k2;
		  d += u;
		  if (dxsquare)
		    {
		      fx += dfx;
		      dfx += ddf;
		    }
		  if (dysquare)
		    {
		      fy += dfy;
		      dfy += ddf;
		    }
		}
	      else
		{
		  x += dxdiag; y += dydiag;
		  u += k2;
		  v += k3;
		  d += v;
		  fx += dfx;
		  dfx += ddf;
		  fy += dfy;
		  dfy += ddf;
#ifdef NOHOP
		  /* Check for octant hopping */
		  if (u < v && fx >= 0 && fy >= 0)
		    {
		      /* hopped diag across ellipse, take square instead */
		      x += dxsquare - dxdiag;
		      y += dysquare - dydiag;
		      d -= v;
		      u += k1 - k2;
		      v += k2 - k3;
		      d += u;
		      dfx -= ddf;
		      fx -= dfx;
		      dfy -= ddf;
		      fy -= dfy;
		      if (dxsquare)
			{
			  fx += dfx;
			  dfx += ddf;
			}
		      if (dysquare)
			{
			  fy += dfy;
			  dfy += ddf;
			}
		    }
#endif
		}
	      if (count && --count == 0)
		goto done;
	    }
	  /* We now cross diagonal octant boundary */
	  d = ROUND (d - u + v/2. - k2/2. + 3*k3/8.);
	  u = ROUND (-u + v - k2/2. + k3/2.);
	  v = ROUND (v - k2 + k3/2.); /* could be v + A - C */
	  k1 = k1 - 2*k2 + k3;
	  k2 = k3 - k2;
	  tmp = dxsquare;
	  dxsquare = -dysquare;
	  dysquare = tmp;
	}
      else			/* Even octants */
	{
	  while (u < k2 /2)
	    {
	      blob_conic_add_pixel (b, present, x, y, octant);
	      if (d<0)
		{
		  x += dxdiag; y += dydiag;
		  u += k2;
		  v += k3;
		  d += v;
		  fx += dfx;
		  dfx += ddf;
		  fy += dfy;
		  dfy += ddf;
		}
	      else
		{
		  x += dxsquare; y += dysquare;
		  u += k1;
		  v += k2;
		  d += u;
		  if (dxsquare)
		    {
		      fx += dfx;
		      dfx += ddf;
		    }
		  if (dysquare)
		    {
		      fy += dfy;
		      dfy += ddf;
		    }
#ifdef NOHOP
		  /* Check for octant hopping */
		   if ((fx >= 0) && (fy >= 0) &&
		       ((((octant == 2) || (octant == 6)) &&
			 2*u > ROUND(v + B/2 + C)) ||
			(((octant == 4) || (octant == 8)) &&
			 2*u > ROUND(v - B/2 + A))))
		       {
			 /* hopped square across ellipse, take diag instead */
			 x += dxdiag - dxsquare;
			 y += dydiag - dysquare;
			 d -= u;
			 u += k2 - k1;
			 v += k3 - k2;
			 d += v;
			 if (dxsquare)
			   {
			     dfx -= ddf;
			     fx -= dfx;
			   }
			 if (dysquare)
			   {
			     dfy -= ddf;
			     fy -= dfy;
			   }
			 fx += dfx;
			 dfx += ddf;
			 fy += dfy;
			 dfy += ddf;
		       }
#endif
		}
	      if (count && --count <= 0)
		goto done;
	    }
	  /* We now cross square octant boundary */
	  d = d + u - v + k1 - k2;
	  v = 2*u - v + k1 - k2;
	  /* Do v first; it depends on u */
	  u = u + k1 - k2;
	  k3 = 4 * (k1 - k2) + k3;
	  k2 = 2 * k1 - k2;
	  tmp = dxdiag;
	  dxdiag = -dydiag;
	  dydiag = tmp;
	  /* Change appropriate boundary checking increment */
	  if (octant == 2 || octant == 6)
	    dfx = -dfx + ddf;
	  else
	    dfy = -dfy + ddf;
	}
      octant++;
      if (octant > 8) octant -= 8;
      octantCount--;
    }

done:				/* jump out of two levels */
  
  for (i=0; i<b->height; i++)
    if (present[i] != (LEFT | RIGHT))
      {
	b->data[i].left = 0;
	b->data[i].right = -1;
      }

  free (present);
}

/* Scan convert an ellipse specified by _offsets_ of major and
   minor axes, and by center into a blob */
Blob *
blob_ellipse (double xc, double yc, double xp, double yp, double xq, double yq)
{
  Blob *r;
  double A,B,C,D,E,F;		/* coefficients of conic */
  double xprod,tmp;
  double height;
  double dpx, dpy;
  int y;

  xprod = xp*yq - xq*yp;

  if (xprod == 0)		/* colinear points */
    {
      g_print("Colinear points!\n");
      g_debug("gsumi");
    }
  
  if (xprod < 0)
    {
      tmp = xp; xp = xq; xq = tmp;
      tmp = yp; yp = yq; yq = tmp;
      xprod = -xprod;
    }
  
  A = yp*yp + yq*yq;
  B = -2 * (xp*yp + xq*yq);
  C = xp*xp + xq*xq;
  D = 2*yq*xprod;
  E = -2*xq*xprod;
  F = 0;
  /* Now offset the ellipse so that the center is exact, but the
     starting point is no longer exactly on the ellipse */

  dpx = ROUND(xp+xc)-xp-xc;
  dpy = ROUND(yp+yc)-yp-yc;

  F += dpx*(A*dpx+D+B*dpy) + dpy*(C*dpy+E);
  D += 2*A*dpx + B*dpy;
  E += B*dpx + 2*C*dpy;

  height = sqrt(A);
  y = floor(yc-height-0.5);
  /* We allow an extra pixel of slop on top and bottom to deal with
     round-off error */
  r = blob_new (y-1,ceil(yc+height+0.5)-y+3);

  /* Although it seems that multiplying A-F by a constant would improve
     things, this seems not to be the case in practice. Several test
     cases showed about equal improvement and degradation */
  blob_conic (r,ROUND(xp+xc),ROUND(yp+yc),ROUND(16*A),ROUND(16*B),ROUND(16*C),
	      ROUND(16*D),ROUND(16*E),ROUND(16*F));

  return r;
}

void
blob_bounds(Blob *b, rect *r)
{
  int i;

  i = 0;
  while (i<b->height && b->data[i].left > b->data[i].right)
    i++;

  if (i<b->height)
    {
      r->y0 = b->y + i;
      r->x0 = b->data[i].left;
      r->x1 = b->data[i].right + 1;
      while (i<b->height && b->data[i].left <= b->data[i].right)
	{
	  r->x0 = MIN(b->data[i].left, r->x0);
	  r->x1 = MAX(b->data[i].right+1, r->x1);
	  i++;
	}
      r->y1 = b->y + i;
    }
  else
    {
      r->x0 = r->y0 = 0;
      r->x1 = r->y1 = 0;
    }
}

void
blob_dump(Blob *b) {
  int i,j;
  for (i=0; i<b->height; i++)
    {
      for (j=0;j<b->data[i].left;j++)
	putchar(' ');
      for (j=b->data[i].left;j<=b->data[i].right;j++)
	putchar('*');
      putchar('\n');
    }
}