/*******************************************************************
 * Line drawing in the current window.
 *
 * Clip window:
 *   D_set_clip_window (top, bottom ,left, right)
 *      establish clipping region for subseqent line drawing.
 *   D_set_clip_window_to_map_window ()
 *     set clipping to pixels corresponding to the current map region
 *     (default)
 *   D_set_clip_window_to_screen_window ()
 *     set clipping to full extent of the window (ie disables clipping on screen)
 *
 * Moves.
 *   D_move_abs(x,y)   move to x,y.
 *   D_move_rel(x,y)   move to +x,+y.
 *      Set current position. Position is not clipped.
 *
 * Draw line 
 *   D_cont_abs(x,y)   draw to x,y.
 *   D_cont_rel(x,y)   draw to +x,+y.
 *      Line draw from current position. New postion is not clipped.
 *      The lines drawn are clipped however.
 *      Return values indicate the nature of the clipping:
 *        0 no clipping
 *        1 part of the line is drawn
 *       -1 none of the line is drawn
 *   
 *
 */

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

#include <grass/gis.h>
#include <grass/raster.h>
#include <grass/display.h>

struct rectangle
{
    double left;
    double rite;
    double bot;
    double top;
};

struct vector
{
    double x, y;
};

struct plane
{
    double x, y, k;
};

static struct vector cur;

static struct rectangle clip;

static struct plane pl_left = { -1, 0, 0 };
static struct plane pl_rite = { 1, 0, 0 };
static struct plane pl_bot = { 0, -1, 0 };
static struct plane pl_top = { 0, 1, 0 };

static int window_set;

#define min(x,y) ((x) < (y) ? (x) : (y))
#define max(x,y) ((x) > (y) ? (x) : (y))

#define round(x) ((int) floor(0.5 + (x)))

static int *xi, *yi;
static int nalloc_i;

static double *xf, *yf;
static int nalloc_f;

static void alloc_int(int n)
{

    if (nalloc_i >= n)
	return;

    nalloc_i = n;
    xi = G_realloc(xi, nalloc_i * sizeof(int));
    yi = G_realloc(yi, nalloc_i * sizeof(int));
}

static void alloc_float(int n)
{

    if (nalloc_f >= n)
	return;

    nalloc_f = n + 10;
    xf = G_realloc(xf, nalloc_f * sizeof(double));
    yf = G_realloc(yf, nalloc_f * sizeof(double));
}

static void dealloc_float(const double **x, const double **y, int release)
{
    if (release) {
	G_free(*(double **)x);
	G_free(*(double **)y);
    }

    *x = xf;
    *y = yf;

    nalloc_f = 0;

    xf = NULL;
    yf = NULL;
}

static int do_filter(int *x, int *y, int n)
{
    int i, j;

    for (i = 0, j = 1; j < n; j++) {
	if (x[j] == x[i] && y[j] == y[i])
	    continue;
	i++;
	if (i == j)
	    continue;
	x[i] = x[j];
	y[i] = y[j];
    }

    return i + 1;
}

static void do_round(const double *x, const double *y, int n)
{
    int i;

    alloc_int(n);

    for (i = 0; i < n; i++) {
	xi[i] = round(D_u_to_d_col(x[i]));
	yi[i] = round(D_u_to_d_row(y[i]));
    }
}

static void do_floor(const double *x, const double *y, int n)
{
    int i;

    alloc_int(n);

    for (i = 0; i < n; i++) {
	xi[i] = floor(D_u_to_d_col(x[i]));
	yi[i] = floor(D_u_to_d_row(y[i]));
    }
}

static double dist_plane(double x, double y, const struct plane *p)
{
    return x * p->x + y * p->y + p->k;
}

static double interpolate(double a, double b, double ka, double kb)
{
    return (a * kb - b * ka) / (kb - ka);
}

static int clip_plane(struct vector *a, struct vector *b,
		      const struct plane *p, int *clipped)
{
    double ka = dist_plane(a->x, a->y, p);
    double kb = dist_plane(b->x, b->y, p);
    double kab;

    /* both outside */
    if (ka > 0 && kb > 0)
	return 1;

    /* both inside */
    if (ka <= 0 && kb <= 0)
	return 0;

    *clipped = 1;

    /* a outside - swap a and b */
    if (ka >= 0) {
	struct vector *t;
	double kt;

	t = a;
	a = b;
	b = t;
	kt = ka;
	ka = kb;
	kb = kt;
    }

    kab = kb - ka;

    b->x = interpolate(a->x, b->x, ka, kb);
    b->y = interpolate(a->y, b->y, ka, kb);

    return 0;
}

static int do_clip(struct vector *a, struct vector *b)
{
    int clipped = 0;

    if (a->x < clip.left && b->x < clip.left)
	return -1;
    if (a->x > clip.rite && b->x > clip.rite)
	return -1;
    if (a->y < clip.bot && b->y < clip.bot)
	return -1;
    if (a->y > clip.top && b->y > clip.top)
	return -1;

    if (clip_plane(a, b, &pl_left, &clipped))
	return -1;
    if (clip_plane(a, b, &pl_rite, &clipped))
	return -1;
    if (clip_plane(a, b, &pl_bot, &clipped))
	return -1;
    if (clip_plane(a, b, &pl_top, &clipped))
	return -1;

    return clipped;
}

static int shift_count(double dx)
{
    return (int)floor(dx / 360);
}

static double shift_angle(double dx)
{
    return shift_count(dx) * 360;
}

static double coerce(double x)
{
    x += 180;
    x -= shift_angle(x);
    x -= 180;
    return x;
}

static int euclidify(double *x, const double *y, int n, int no_pole)
{
    double ux0 = clip.left;
    double ux1 = clip.rite;
    double x0, x1;
    int lo, hi, count;
    int i;

    x0 = x1 = x[0];

    for (i = 1; i < n; i++) {
	if (fabs(y[i]) < 89.9)
	    x[i] = x[i - 1] + coerce(x[i] - x[i - 1]);

	x0 = min(x0, x[i]);
	x1 = max(x1, x[i]);
    }

    if (no_pole && fabs(x[n - 1] - x[0]) > 180)
	return 0;

    lo = -shift_count(ux1 - x0);
    hi = shift_count(x1 - ux0);
    count = hi - lo + 1;

    for (i = 0; i < n; i++)
	x[i] -= lo * 360;

    return count;
}

static void do_ll_wrap(const double *x, const double *y, int n,
		       void (*func) (const double *, const double *, int))
{
    double *xx = G_malloc(n * sizeof(double));
    int count, i;

    memcpy(xx, x, n * sizeof(double));
    count = euclidify(xx, y, n, 0);

    for (i = 0; i < count; i++) {
	int j;

	(*func) (xx, y, n);

	for (j = 0; j < n; j++)
	    xx[j] -= 360;
    }

    G_free(xx);
}

/*!
 * \brief set clipping window
 *
 * Sets the clipping window to the pixel window that corresponds
 * to the current database region. This is the default.
 *
 *  \param top
 *  \param bottom
 *  \param left
 *  \param right
 */

void D_set_clip(double t, double b, double l, double r)
{
    clip.left = min(l, r);
    clip.rite = max(l, r);
    clip.bot = min(b, t);
    clip.top = max(b, t);

    pl_left.k = clip.left;
    pl_rite.k = -clip.rite;
    pl_bot.k = clip.bot;
    pl_top.k = -clip.top;

    window_set = 1;
}

/*!
 * \brief set clipping window to map window
 *
 * Sets the clipping window to the pixel window that corresponds to the
 * current database region. This is the default.
 *
 *  \param ~
 */

void D_clip_to_map(void)
{
    D_set_clip(D_get_u_north(),
	       D_get_u_south(), D_get_u_west(), D_get_u_east());
}

/*!
 * \brief move to pixel
 *
 * Move without drawing to
 * pixel location <b>x,y</b>, even if it falls outside the clipping window.
 *
 *  \param x
 *  \param y
 */

void D_move_clip(double x, double y)
{
    cur.x = x;
    cur.y = y;
}

/*!
 * \brief line to x,y
 *
 * Draws a line from the
 * current position to pixel location <b>x,y.</b> Any part of the line that
 * falls outside the clipping window is not drawn.
 * <b>Note.</b> The new position is <b>x,y</b>, even if it falls outside the
 * clipping window. Returns 0 if the line was contained entirely in the clipping
 * window, 1 if the line had to be clipped to draw it.
 *
 *  \param x
 *  \param y
 *  \return int
 */

static int line_clip(double x1, double y1, double x2, double y2)
{
    struct vector a, b;
    int clipped;

    a.x = x1;
    a.y = y1;

    b.x = x2;
    b.y = y2;

    clipped = do_clip(&a, &b);

    if (clipped >= 0) {
	int x1 = round(D_u_to_d_col(a.x));
	int y1 = round(D_u_to_d_row(a.y));
	int x2 = round(D_u_to_d_col(b.x));
	int y2 = round(D_u_to_d_row(b.y));

	R_move_abs(x1, y1);
	R_cont_abs(x2, y2);
    }

    return clipped;
}

static int line_clip_ll(double ax, double ay, double bx, double by)
{
    double ux0 = clip.left;
    double ux1 = clip.rite;
    double x0, x1;
    int lo, hi, i;
    int ret;

    bx = ax + coerce(bx - ax);

    x0 = min(ax, bx);
    x1 = max(ax, bx);

    lo = -shift_count(ux1 - x0);
    hi = shift_count(x1 - ux0);

    ret = 0;

    for (i = lo; i <= hi; i++)
	ret |= line_clip(ax + i * 360, ay, bx + i * 360, by);

    return ret;
}

int D_cont_clip(double x, double y)
{
    int ret;

    if (!window_set)
	D_clip_to_map();

    if (D_is_lat_lon())
	ret = line_clip_ll(cur.x, cur.y, x, y);
    else
	ret = line_clip(cur.x, cur.y, x, y);

    cur.x = x;
    cur.y = y;

    return ret;
}

void D_polydots_clip(const double *x, const double *y, int n)
{
    double ux0 = clip.left;
    int i, j;

    if (!window_set)
	D_clip_to_map();

    alloc_float(n);

    for (i = j = 0; i < n; i++) {
	double xx = x[i];
	double yy = y[i];

	if (D_is_lat_lon())
	    xx -= shift_angle(x[i] - ux0);

	if (xx < clip.left || xx > clip.rite)
	    continue;
	if (yy < clip.bot || yy > clip.top)
	    continue;

	xf[j] = xx;
	yf[j] = yy;
	j++;
    }

    do_floor(xf, yf, n);
    n = do_filter(xi, yi, n);

    R_polydots_abs(xi, yi, j);
}

static int cull_polyline_plane(int *pn, const double *x, const double *y,
			       const struct plane *p)
{
    int n = *pn;
    int last = -1;
    int prev = 0;
    double x0 = x[prev];
    double y0 = y[prev];
    double d0 = dist_plane(x0, y0, p);
    int i, j;

    for (i = 0, j = 0; i < n; i++) {
	double x1 = x[i];
	double y1 = y[i];
	double d1 = dist_plane(x1, y1, p);
	int in0 = d0 <= 0;
	int in1 = d1 <= 0;

	if (!in0 && in1 && last != prev) {	/* entering */
	    alloc_float(j + 1);
	    xf[j] = x0;
	    yf[j] = y0;
	    j++;
	    last = prev;
	}

	if (in1 || in0) {	/* inside or leaving */
	    alloc_float(j + 1);
	    xf[j] = x1;
	    yf[j] = y1;
	    j++;
	    last = i;
	}

	x0 = x1;
	y0 = y1;
	d0 = d1;
	prev = i;
    }

    *pn = j;

    return (j == 0);
}

static void polyline_cull(const double *x, const double *y, int n)
{
    alloc_float(n + 10);

    if (cull_polyline_plane(&n, x, y, &pl_left))
	return;

    dealloc_float(&x, &y, 0);

    if (cull_polyline_plane(&n, x, y, &pl_rite))
	return;

    dealloc_float(&x, &y, 1);

    if (cull_polyline_plane(&n, x, y, &pl_bot))
	return;

    dealloc_float(&x, &y, 1);

    if (cull_polyline_plane(&n, x, y, &pl_top))
	return;

    dealloc_float(&x, &y, 1);

    do_floor(x, y, n);
    n = do_filter(xi, yi, n);

    R_polyline_abs(xi, yi, n);
}

void D_polyline_cull(const double *x, const double *y, int n)
{
    if (n < 2)
	return;

    if (!window_set)
	D_clip_to_map();

    if (D_is_lat_lon())
	do_ll_wrap(x, y, n, polyline_cull);
    else
	polyline_cull(x, y, n);
}

static void polyline_clip(const double *x, const double *y, int n)
{
    int i;

    for (i = 1; i < n; i++)
	line_clip(x[i - 1], y[i - 1], x[i], y[i]);
}

void D_polyline_clip(const double *x, const double *y, int n)
{
    if (n < 2)
	return;

    if (!window_set)
	D_clip_to_map();

    if (D_is_lat_lon())
	do_ll_wrap(x, y, n, polyline_clip);
    else
	polyline_clip(x, y, n);
}

static int cull_polygon_plane(int *pn, const double *x, const double *y,
			      const struct plane *p)
{
    int n = *pn;
    int last = -1;
    int prev = n - 1;
    double x0 = x[prev];
    double y0 = y[prev];
    double d0 = dist_plane(x0, y0, p);
    int i, j;

    for (i = j = 0; i < n; i++) {
	double x1 = x[i];
	double y1 = y[i];
	double d1 = dist_plane(x1, y1, p);
	int in0 = d0 <= 0;
	int in1 = d1 <= 0;

	if (!in0 && in1 && last != prev) {	/* entering */
	    alloc_float(j + 1);
	    xf[j] = x0;
	    yf[j] = y0;
	    j++;
	    last = prev;
	}

	if (in1 || in0) {	/* inside or leaving */
	    alloc_float(j + 1);
	    xf[j] = x1;
	    yf[j] = y1;
	    j++;
	    last = i;
	}

	x0 = x1;
	y0 = y1;
	d0 = d1;
	prev = i;
    }

    *pn = j;

    return (j == 0);
}

static void polygon_cull(const double *x, const double *y, int n)
{
    alloc_float(n + 10);

    if (cull_polygon_plane(&n, x, y, &pl_left))
	return;

    dealloc_float(&x, &y, 0);

    if (cull_polygon_plane(&n, x, y, &pl_rite))
	return;

    dealloc_float(&x, &y, 1);

    if (cull_polygon_plane(&n, x, y, &pl_bot))
	return;

    dealloc_float(&x, &y, 1);

    if (cull_polygon_plane(&n, x, y, &pl_top))
	return;

    dealloc_float(&x, &y, 1);

    do_round(x, y, n);
    n = do_filter(xi, yi, n);

    R_polygon_abs(xi, yi, n);
}

void D_polygon_cull(const double *x, const double *y, int n)
{
    if (!window_set)
	D_clip_to_map();

    if (D_is_lat_lon())
	do_ll_wrap(x, y, n, polygon_cull);
    else
	polygon_cull(x, y, n);
}

static int clip_polygon_plane(int *pn, const double *x, const double *y,
			      const struct plane *p)
{
    int n = *pn;
    double x0 = x[n - 1];
    double y0 = y[n - 1];
    double d0 = dist_plane(x0, y0, p);
    int i, j;

    for (i = j = 0; i < n; i++) {
	double x1 = x[i];
	double y1 = y[i];
	double d1 = dist_plane(x1, y1, p);
	int in0 = d0 <= 0;
	int in1 = d1 <= 0;

	if (in0 != in1) {	/* edge crossing */
	    alloc_float(j + 1);
	    xf[j] = interpolate(x0, x1, d0, d1);
	    yf[j] = interpolate(y0, y1, d0, d1);
	    j++;
	}

	if (in1) {		/* point inside */
	    alloc_float(j + 1);
	    xf[j] = x[i];
	    yf[j] = y[i];
	    j++;
	}

	x0 = x1;
	y0 = y1;
	d0 = d1;
    }

    *pn = j;

    return (j == 0);
}

static void polygon_clip(const double *x, const double *y, int n)
{
    alloc_float(n + 10);

    if (clip_polygon_plane(&n, x, y, &pl_left))
	return;

    dealloc_float(&x, &y, 0);

    if (clip_polygon_plane(&n, x, y, &pl_rite))
	return;

    dealloc_float(&x, &y, 1);

    if (clip_polygon_plane(&n, x, y, &pl_bot))
	return;

    dealloc_float(&x, &y, 1);

    if (clip_polygon_plane(&n, x, y, &pl_top))
	return;

    dealloc_float(&x, &y, 1);

    do_round(x, y, n);
    n = do_filter(xi, yi, n);

    R_polygon_abs(xi, yi, n);
}

void D_polygon_clip(const double *x, const double *y, int n)
{
    if (!window_set)
	D_clip_to_map();

    if (D_is_lat_lon())
	do_ll_wrap(x, y, n, polygon_clip);
    else
	polygon_clip(x, y, n);
}

static void box_clip(double x1, double y1, double x2, double y2)
{
    double t, b, l, r;
    int ti, bi, li, ri;

    l = max(clip.left, min(x1, x2));
    r = min(clip.rite, max(x1, x2));
    b = max(clip.bot, min(y1, y2));
    t = min(clip.top, max(y1, y2));

    li = round(D_u_to_d_col(l));
    ri = round(D_u_to_d_col(r));
    bi = round(D_u_to_d_row(b));
    ti = round(D_u_to_d_row(t));

    R_box_abs(li, ti, ri, bi);
}

static void box_clip_ll(double x1, double y1, double x2, double y2)
{
    double ux0 = clip.left;
    double ux1 = clip.rite;
    int lo, hi, i;

    x2 = x1 + coerce(x2 - x1);

    lo = -shift_count(ux1 - x1);
    hi = shift_count(x2 - ux0);

    for (i = lo; i <= hi; i++)
	box_clip(x1 + i * 360, y1, x2 + i * 360, y2);
}

void D_box_clip(double x1, double y1, double x2, double y2)
{
    if (!window_set)
	D_clip_to_map();

    if (D_is_lat_lon())
	box_clip_ll(x1, y1, x2, y2);
    else
	box_clip(x1, y1, x2, y2);
}

void D_move(double x, double y)
{
    int xi = round(D_u_to_d_col(x));
    int yi = round(D_u_to_d_row(y));

    R_move_abs(xi, yi);
}

void D_cont(double x, double y)
{
    int xi = round(D_u_to_d_col(x));
    int yi = round(D_u_to_d_row(y));

    R_cont_abs(xi, yi);
}

void D_polydots(const double *x, const double *y, int n)
{
    do_floor(x, y, n);
    n = do_filter(xi, yi, n);
    R_polydots_abs(xi, yi, n);
}

void D_polyline(const double *x, const double *y, int n)
{
    do_floor(x, y, n);
    n = do_filter(xi, yi, n);
    R_polyline_abs(xi, yi, n);
}

void D_polygon(const double *x, const double *y, int n)
{
    do_round(x, y, n);
    n = do_filter(xi, yi, n);
    R_polygon_abs(xi, yi, n);
}

void D_box(double x1, double y1, double x2, double y2)
{
    double l = min(x1, x2);
    double r = max(x1, x2);
    double b = min(y1, y2);
    double t = max(y1, y2);
    int li = round(D_u_to_d_col(l));
    int ri = round(D_u_to_d_col(r));
    int bi = round(D_u_to_d_row(b));
    int ti = round(D_u_to_d_row(t));

    R_box_abs(li, ti, ri, bi);
}

void D_line_width(double d)
{
    int w = round(d);

    if (w < 0)
	w = 0;

    R_line_width(w);
}