/* blob.c: routines for manipulating scan converted convex
 *         polygons.
 *
 * Copyright 1998-1999, Owen Taylor <otaylor@gtk.org>
 *
 * > Please contact the above author before modifying the copy <
 * > of this file in the GIMP distribution. Thanks.            <
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
*/

#include "config.h"

#include <glib-object.h>

#include "libgimpmath/gimpmath.h"

#include "paint-types.h"

#include "gimpink-blob.h"


typedef enum
{
  EDGE_NONE  = 0,
  EDGE_LEFT  = 1 << 0,
  EDGE_RIGHT = 1 << 1
} EdgeType;


/*  local function prototypes  */

static Blob * blob_new            (gint      y,
                                   gint      height);
static void   blob_fill           (Blob     *b,
                                   EdgeType *present);
static void   blob_make_convex    (Blob     *b,
                                   EdgeType *present);

#if 0
static void   blob_line_add_pixel (Blob     *b,
                                   gint      x,
                                   gint      y);
static void   blob_line           (Blob     *b,
                                   gint      x0,
                                   gint      y0,
                                   gint      x1,
                                   gint      y1);
#endif


/*  public functions  */

/* Return blob for the given (convex) polygon
 */
Blob *
blob_polygon (BlobPoint *points,
              gint       npoints)
{
  Blob     *result;
  EdgeType *present;
  gint      i;
  gint      im1;
  gint      ip1;
  gint      ymin, ymax;

  ymax = points[0].y;
  ymin = points[0].y;

  for (i = 1; i < npoints; i++)
    {
      if (points[i].y > ymax)
        ymax = points[i].y;
      if (points[i].y < ymin)
        ymin = points[i].y;
    }

  result = blob_new (ymin, ymax - ymin + 1);
  present = g_new0 (EdgeType, result->height);

  im1 = npoints - 1;
  i = 0;
  ip1 = 1;

  for (; i < npoints ; i++)
    {
      gint sides = 0;
      gint j     = points[i].y - ymin;

      if (points[i].y < points[im1].y)
        sides |= EDGE_RIGHT;
      else if (points[i].y > points[im1].y)
        sides |= EDGE_LEFT;

      if (points[ip1].y < points[i].y)
        sides |= EDGE_RIGHT;
      else if (points[ip1].y > points[i].y)
        sides |= EDGE_LEFT;

      if (sides & EDGE_RIGHT)
        {
          if (present[j] & EDGE_RIGHT)
            {
              result->data[j].right = MAX (result->data[j].right, points[i].x);
            }
          else
            {
              present[j] |= EDGE_RIGHT;
              result->data[j].right = points[i].x;
            }
        }

      if (sides & EDGE_LEFT)
        {
          if (present[j] & EDGE_LEFT)
            {
              result->data[j].left = MIN (result->data[j].left, points[i].x);
            }
          else
            {
              present[j] |= EDGE_LEFT;
              result->data[j].left = points[i].x;
            }
        }

      im1 = i;
      ip1++;
      if (ip1 == npoints)
        ip1 = 0;
    }

  blob_fill (result, present);
  g_free (present);

  return result;
}

/* Scan convert a square specified by _offsets_ of major and minor
 * axes, and by center into a blob
 */
Blob *
blob_square (gdouble xc,
             gdouble yc,
             gdouble xp,
             gdouble yp,
             gdouble xq,
             gdouble yq)
{
  BlobPoint points[4];

  /* Make sure we order points ccw */

  if (xp * yq - yq * xp < 0)
    {
      xq = -xq;
      yq = -yq;
    }

  points[0].x = xc + xp + xq;
  points[0].y = yc + yp + yq;
  points[1].x = xc + xp - xq;
  points[1].y = yc + yp - yq;
  points[2].x = xc - xp - xq;
  points[2].y = yc - yp - yq;
  points[3].x = xc - xp + xq;
  points[3].y = yc - yp + yq;

  return blob_polygon (points, 4);
}

/* Scan convert a diamond specified by _offsets_ of major and minor
 * axes, and by center into a blob
 */
Blob *
blob_diamond (gdouble xc,
              gdouble yc,
              gdouble xp,
              gdouble yp,
              gdouble xq,
              gdouble yq)
{
  BlobPoint points[4];

  /* Make sure we order points ccw */

  if (xp * yq - yq * xp < 0)
    {
      xq = -xq;
      yq = -yq;
    }

  points[0].x = xc + xp;
  points[0].y = yc + yp;
  points[1].x = xc - xq;
  points[1].y = yc - yq;
  points[2].x = xc - xp;
  points[2].y = yc - yp;
  points[3].x = xc + xq;
  points[3].y = yc + yq;

  return blob_polygon (points, 4);
}


#define TABLE_SIZE 256

#define ELLIPSE_SHIFT 2
#define TABLE_SHIFT   12
#define TOTAL_SHIFT   (ELLIPSE_SHIFT + TABLE_SHIFT)

/*
 * The choose of this values limits the maximal image_size to
 * 16384 x 16384 pixels. The values will overflow as soon as
 * x or y > INT_MAX / (1 << (ELLIPSE_SHIFT + TABLE_SHIFT)) / SUBSAMPLE
 *
 * Alternatively the code could be change the code as follows:
 *
 *   xc_base = floor (xc)
 *   xc_shift = 0.5 + (xc - xc_base) * (1 << TOTAL_SHIFT);
 *
 *    gint x = xc_base + (xc_shift + c * xp_shift + s * xq_shift +
 *             (1 << (TOTAL_SHIFT - 1))) >> TOTAL_SHIFT;
 *
 * which would change the limit from the image to the ellipse size
 */

static gboolean trig_initialized = FALSE;
static gint     trig_table[TABLE_SIZE];

/* Scan convert an ellipse specified by _offsets_ of major and
 * minor axes, and by center into a blob
 */
Blob *
blob_ellipse (gdouble xc,
              gdouble yc,
              gdouble xp,
              gdouble yp,
              gdouble xq,
              gdouble yq)
{
  Blob     *result;
  EdgeType *present;
  gint      i;
  gdouble   r1, r2;
  gint      maxy, miny;
  gint      step;
  gdouble   max_radius;

  gint      xc_shift, yc_shift;
  gint      xp_shift, yp_shift;
  gint      xq_shift, yq_shift;

  if (! trig_initialized)
    {
      trig_initialized = TRUE;

      for (i = 0; i < 256; i++)
        trig_table[i] = 0.5 + sin (i * (G_PI / 128.0)) * (1 << TABLE_SHIFT);
    }

  /* Make sure we traverse ellipse in ccw direction */

  if (xp * yq - yq * xp < 0)
    {
      xq = -xq;
      yq = -yq;
    }

  /* Compute bounds as if we were drawing a rectangle */

  maxy = ceil  (yc + fabs (yp) + fabs (yq));
  miny = floor (yc - fabs (yp) - fabs (yq));

  result = blob_new (miny, maxy - miny + 1);
  present = g_new0 (EdgeType, result->height);

  /* Figure out a step that will draw most of the points */

  r1 = sqrt (xp * xp + yp * yp);
  r2 = sqrt (xq * xq + yq * yq);
  max_radius = MAX (r1, r2);
  step = TABLE_SIZE;

  while (step > 1 && (TABLE_SIZE / step < 4 * max_radius))
    step >>= 1;

  /* Fill in the edge points */

  xc_shift = 0.5 + xc * (1 << TOTAL_SHIFT);
  yc_shift = 0.5 + yc * (1 << TOTAL_SHIFT);
  xp_shift = 0.5 + xp * (1 << ELLIPSE_SHIFT);
  yp_shift = 0.5 + yp * (1 << ELLIPSE_SHIFT);
  xq_shift = 0.5 + xq * (1 << ELLIPSE_SHIFT);
  yq_shift = 0.5 + yq * (1 << ELLIPSE_SHIFT);

  for (i = 0 ; i < TABLE_SIZE ; i += step)
    {
      gint s = trig_table[i];
      gint c = trig_table[(TABLE_SIZE + TABLE_SIZE / 4 - i) % TABLE_SIZE];

      gint x = (xc_shift + c * xp_shift + s * xq_shift +
                (1 << (TOTAL_SHIFT - 1))) >> TOTAL_SHIFT;
      gint y = ((yc_shift + c * yp_shift + s * yq_shift +
                (1 << (TOTAL_SHIFT - 1))) >> TOTAL_SHIFT) - result->y;

      gint dydi = c * yq_shift  - s * yp_shift;

      if (dydi <= 0) /* left edge */
        {
          if (present[y] & EDGE_LEFT)
            {
              result->data[y].left = MIN (result->data[y].left, x);
            }
          else
            {
              present[y] |= EDGE_LEFT;
              result->data[y].left = x;
            }
        }

      if (dydi >= 0) /* right edge */
        {
          if (present[y] & EDGE_RIGHT)
            {
              result->data[y].right = MAX (result->data[y].right, x);
            }
          else
            {
              present[y] |= EDGE_RIGHT;
              result->data[y].right = x;
            }
        }
    }

  /* Now fill in missing points */

  blob_fill (result, present);
  g_free (present);

  return result;
}

void
blob_bounds (Blob *b,
             gint *x,
             gint *y,
             gint *width,
             gint *height)
{
  gint i;
  gint x0, x1, y0, y1;

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

  if (i < b->height)
    {
      y0 = b->y + i;
      x0 = b->data[i].left;
      x1 = b->data[i].right + 1;

      while (i < b->height && b->data[i].left <= b->data[i].right)
        {
          x0 = MIN (b->data[i].left,      x0);
          x1 = MAX (b->data[i].right + 1, x1);
          i++;
        }

      y1 = b->y + i;
    }
  else
    {
      x0 = y0 = 0;
      x1 = y1 = 0;
    }

  *x = x0;
  *y = y0;
  *width = x1 - x0;
  *height = y1 - y0;
}

Blob *
blob_convex_union (Blob *b1,
                   Blob *b2)
{
  Blob     *result;
  gint      y;
  gint      i, j;
  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 = g_new0 (EdgeType, result->height);

  /* 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] = EDGE_LEFT | EDGE_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] = EDGE_LEFT | EDGE_RIGHT;
              result->data[j].left = b2->data[i].left;
              result->data[j].right = b2->data[i].right;
            }
        }
    }

  blob_make_convex (result, present);

  g_free (present);

  return result;
}

Blob *
blob_duplicate (Blob *b)
{
  g_return_val_if_fail (b != NULL, NULL);

  return g_memdup (b, sizeof (Blob) +  sizeof (BlobSpan) * (b->height - 1));
}

#if 0
void
blob_dump (Blob *b)
{
  gint 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');
    }
}
#endif


/*  private functions  */

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

  result = g_malloc (sizeof (Blob) +  sizeof (BlobSpan) * (height - 1));

  result->y      = y;
  result->height = height;

  return result;
}

static void
blob_fill (Blob     *b,
           EdgeType *present)
{
  gint start;
  gint x1, x2, i1, i2;
  gint i;

  /* Mark empty lines at top and bottom as unused */

  start = 0;
  while (! present[start])
    {
      b->data[start].left  = 0;
      b->data[start].right = -1;
      start++;
    }

  if (present[start] != (EDGE_RIGHT | EDGE_LEFT))
    {
      if (present[start] == EDGE_RIGHT)
        b->data[start].left = b->data[start].right;
      else
        b->data[start].right = b->data[start].left;

      present[start] = EDGE_RIGHT | EDGE_LEFT;
    }

  for (i = b->height - 1; ! present[i]; i--)
    {
      b->data[i].left  = 0;
      b->data[i].right = -1;
    }

  if (present[i] != (EDGE_RIGHT | EDGE_LEFT))
    {
      if (present[i] == EDGE_RIGHT)
        b->data[i].left = b->data[i].right;
      else
        b->data[i].right = b->data[i].left;

      present[i] = EDGE_RIGHT | EDGE_LEFT;
    }


  /* Restore missing edges  */

  /* 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 < b->height - 2; i1++)
    {
      /* Find empty gaps */
      if (! (present[i1 + 1] & EDGE_LEFT))
        {
          gint increment; /* fractional part */
          gint denom;     /* denominator of fraction */
          gint step;      /* integral step */
          gint frac;      /* fractional step */
          gint reverse;

          /* find bottom of gap */
          i2 = i1 + 2;
          while (i2 < b->height && ! (present[i2] & EDGE_LEFT))
            i2++;

          if (i2 < b->height)
            {
              denom = i2 - i1;
              x1 = b->data[i1].left;
              x2 = b->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)
                    b->data[i].left = x1;
                  else
                    b->data[i].left = x1 + 1;
                }
            }

          i1 = i2 - 1;        /* advance to next possibility */
        }
    }

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

          /* find bottom of gap */
          i2 = i1 + 2;
          while (i2 < b->height && ! (present[i2] & EDGE_RIGHT))
            i2++;

          if (i2 < b->height)
            {
              denom = i2 - i1;
              x1 = b->data[i1].right;
              x2 = b->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)
                    b->data[i].right = x1 - 1;
                  else
                    b->data[i].right = x1;
                }
            }

          i1 = i2 - 1;        /* advance to next possibility */
        }
    }

}

static void
blob_make_convex (Blob     *b,
                  EdgeType *present)
{
  gint x1, x2, y1, y2, i1, i2;
  gint i;
  gint start;

  /* 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 = b->data[start].left - b->data[start].right;
  y1 = 0;

  for (i = start + 1; i < b->height; i++)
    {
      if (! (present[i] & EDGE_LEFT))
        continue;

      x2 = b->data[i].left - b->data[i2].left;
      y2 = i - i2;

      while (x2 * y1 - x1 * y2 < 0) /* clockwise rotation */
        {
          present[i2] &= ~EDGE_LEFT;
          i2 = i1;
          while ((--i1) >= start && (! (present[i1] & EDGE_LEFT)));

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

      x1 = x2;
      y1 = y2;
      i1 = i2;
      i2 = i;
    }

  /*     Right edge */

  i1 = start -1;
  i2 = start;
  x1 = b->data[start].right - b->data[start].left;
  y1 = 0;

  for (i = start + 1; i < b->height; i++)
    {
      if (! (present[i] & EDGE_RIGHT))
        continue;

      x2 = b->data[i].right - b->data[i2].right;
      y2 = i - i2;

      while (x2 * y1 - x1 * y2 > 0) /* counter-clockwise rotation */
        {
          present[i2] &= ~EDGE_RIGHT;
          i2 = i1;
          while ((--i1) >= start && (! (present[i1] & EDGE_RIGHT)));

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

          x2 = b->data[i].right - b->data[i2].right;
          y2 = i - i2;
        }

      x1 = x2;
      y1 = y2;
      i1 = i2;
      i2 = i;
    }

  blob_fill (b, present);
}


#if 0

static void
blob_line_add_pixel (Blob *b,
                     gint  x,
                     gint  y)
{
  if (b->data[y - b->y].left > b->data[y - b->y].right)
    {
      b->data[y - b->y].left = b->data[y - b->y].right = x;
    }
  else
    {
      b->data[y - b->y].left  = MIN (b->data[y - b->y].left,  x);
      b->data[y - b->y].right = MAX (b->data[y - b->y].right, x);
    }
}

static void
blob_line (Blob *b,
           gint  x0,
           gint  y0,
           gint  x1,
           gint  y1)
{
  gint dx, dy, d;
  gint incrE, incrNE;
  gint x, y;

  gint xstep = 1;
  gint ystep = 1;

  dx = x1 - x0;
  dy = y1 - y0;

  if (dx < 0)
    {
      dx = -dx;
      xstep = -1;
    }

  if (dy < 0)
    {
      dy = -dy;
      ystep = -1;
    }

  /*  for (y = y0; y != y1 + ystep ; y += ystep)
    {
      b->data[y-b->y].left = 0;
      b->data[y-b->y].right = -1;
      }*/

  x = x0;
  y = y0;

  if (dy < dx)
    {
      d = 2 * dy - dx;        /* initial value of d */
      incrE  = 2 * dy;        /* increment used for move to E */
      incrNE = 2 * (dy - dx); /* increment used for move to NE */

      blob_line_add_pixel (b, x, y);

      while (x != x1)
        {
          if (d <= 0)
            {
              d += incrE;
              x += xstep;
            }
          else
            {
              d += incrNE;
              x += xstep;
              y += ystep;
            }

          blob_line_add_pixel (b, x, y);
        }
    }
  else
    {
      d = 2 * dx - dy;        /* initial value of d */
      incrE  = 2 * dx;        /* increment used for move to E */
      incrNE = 2 * (dx - dy); /* increment used for move to NE */

      blob_line_add_pixel (b, x, y);

      while (y != y1)
        {
          if (d <= 0)
            {
              d += incrE;
              y += ystep;
            }
          else
            {
              d += incrNE;
              x += xstep;
              y += ystep;
            }

          blob_line_add_pixel (b, x, y);
        }
    }
}

#endif