Geografic-Information-System/lib/gis/plot.c

/*!
 * \file lib/gis/plot.c
 *
 * \brief GIS Library - Plotting functions.
 *
 * Plot lines and filled polygons. Input space is current
 * window. Output space and output functions are user
 * defined. Converts input east,north lines and polygons to output x,y
 * and calls user supplied line drawing routines to do the plotting.
 *
 * Handles global wrap-around for lat-lon locations.
 *
 * Does not perform window clipping.
 * Clipping must be done by the line draw routines supplied by the user.
 *
 * Note:
 *  Hopefully, cartographic style projection plotting will be added later.
 *
 * (C) 2001-2008, 2013 by the GRASS Development Team
 *
 * This program is free software under the GNU General Public License
 * (>=v2). Read the file COPYING that comes with GRASS for details.
 *
 * \author Original author CERL
 */

#include <stdlib.h>
#include <math.h>
#include <grass/gis.h>

static void fastline(double, double, double, double);
static void slowline(double, double, double, double);
static void plot_line(double, double, double, double, void (*)());
static double wrap_east(double, double);
static int edge(double, double, double, double);
static int edge_point(double, int);

static int edge_order(const void *, const void *);
static void row_solid_fill(int, double, double);
static void row_dotted_fill(int, double, double);
static int ifloor(double);
static int iceil(double);

struct point {
    double x;
    int y;
};
#define POINT struct point

static struct state {
    struct Cell_head window;
    double xconv, yconv;
    double left, right, top, bottom;
    int ymin, ymax;
    int dotted_fill_gap;

    POINT *P;
    int np;
    int npalloc;

    void (*row_fill)(int, double, double);
    int (*move)(int, int);
    int (*cont)(int, int);
} state;

static struct state *st = &state;

#define OK            0
#define TOO_FEW_EDGES 2
#define NO_MEMORY     1
#define OUT_OF_SYNC  -1

/*!
 * \brief Initialize plotting routines
 *
 * Initializes the plotting capability. This routine must be called
 * once before calling the G_plot_*() routines described below.  The
 * parameters <i>t, b, l, r</i> are the top, bottom, left, and right
 * of the output x,y coordinate space. They are not integers, but
 * doubles to allow for subpixel registration of the input and output
 * coordinate spaces. The input coordinate space is assumed to be the
 * current GRASS region, and the routines supports both planimetric
 * and latitude-longitude coordinate systems.

 * <b>Move</b> and <b>Cont</b> are subroutines that will draw lines in x,y
 * space. They will be called as follows:
 * - Move(x, y) move to x,y (no draw)
 * - Cont(x, y) draw from previous position to x,y. Cont(~) is responsible for clipping
 *
 * \param t,b,l,r top, bottom, left, right
 * \param move Move function
 * \param Cont Cont function
 */
void G_setup_plot(double t, double b, double l, double r,
		  int (*Move) (int, int), int (*Cont) (int, int))
{
    G_get_set_window(&st->window);

    st->left = l;
    st->right = r;
    st->top = t;
    st->bottom = b;

    st->xconv = (st->right - st->left) / (st->window.east - st->window.west);
    st->yconv = (st->bottom - st->top) / (st->window.north - st->window.south);

    if (st->top < st->bottom) {
	st->ymin = iceil(st->top);
	st->ymax = ifloor(st->bottom);
    }
    else {
	st->ymin = iceil(st->bottom);
	st->ymax = ifloor(st->top);
    }

    st->move = Move;
    st->cont = Cont;
}

/*!
 * \brief Set row_fill routine to row_solid_fill or row_dotted_fill
 *
 * After calling this function, G_plot_polygon() and G_plot_area()
 * fill shapes with solid or dotted lines. If gap is greater than
 * zero, this value will be used for row_dotted_fill.  Otherwise,
 * row_solid_fill is used.
 *
 * \param gap
 */
void G_setup_fill(int gap)
{
    if (gap > 0) {
	st->row_fill = row_dotted_fill;
	st->dotted_fill_gap = gap + 1;
    }
    else
	st->row_fill = row_solid_fill;
}

#define X(e) (st->left + st->xconv * ((e) - st->window.west))
#define Y(n) (st->top + st->yconv * (st->window.north - (n)))

#define EAST(x) (st->window.west + ((x)-st->left)/st->xconv)
#define NORTH(y) (st->window.north - ((y)-st->top)/st->yconv)


/*!
 * \brief Converts east,north to x,y
 *
 * The map coordinates <i>east,north</i> are converted
 * to pixel coordinates <i>x,y</i>.
 *
 * \param east easting
 * \param north nothing
 * \param x x coordinate
 * \param y y coordinate
 */
void G_plot_where_xy(double east, double north, int *x, int *y)
{
    *x = ifloor(X(G_adjust_easting(east, &st->window)) + 0.5);
    *y = ifloor(Y(north) + 0.5);
}

/*!
 * \brief Converts x,y to east,north
 *
 * The pixel coordinates <i>x,y</i> are converted to map
 * coordinates <i>east,north</i>.
 *
 * \param x x coordinate
 * \param y y coordinate
 * \param east easting
 * \param north northing
 */

void G_plot_where_en(int x, int y, double *east, double *north)
{
    *east = G_adjust_easting(EAST(x), &st->window);
    *north = NORTH(y);
}

/*!
  \brief Plot point

  \param east easting
  \param north northing
*/
void G_plot_point(double east, double north)
{
    int x, y;

    G_plot_where_xy(east, north, &x, &y);
    st->move(x, y);
    st->cont(x, y);
}

/*!
 * \brief Plot line between latlon coordinates (fastline)
 *
 * A line from <i>east1,north1</i> to <i>east2,north2</i> is plotted
 * in output x,y coordinates (e.g. pixels for graphics.) This routine
 * handles global wrap-around for latitude-longitude databases.
 *
 * \param east1, north1 first point (start line node)
 * \param east2, north2 second point (end line node)
 */
void G_plot_line(double east1, double north1, double east2, double north2)
{
    plot_line(east1, north1, east2, north2, fastline);
}

/*!
 * \brief Plot line between latlon coordinates (slowline)
 *
 * A line from <i>east1,north1</i> to <i>east2,north2</i> is plotted
 * in output x,y coordinates (e.g. pixels for graphics.) This routine
 * handles global wrap-around for latitude-longitude databases.
 *
 * \param east1, north1 first point (start line node)
 * \param east2, north2 second point (end line node)
 */
void G_plot_line2(double east1, double north1, double east2, double north2)
{
    plot_line(east1, north1, east2, north2, slowline);
}

/* fastline converts double rows/cols to ints then plots
 * this is ok for graphics, but not the best for vector to raster
 */
static void fastline(double x1, double y1, double x2, double y2)
{
    st->move(ifloor(x1 + 0.5), ifloor(y1 + 0.5));
    st->cont(ifloor(x2 + 0.5), ifloor(y2 + 0.5));
}

/* NOTE (shapiro): 
 *   I think the adding of 0.5 in slowline is not correct
 *   the output window (left, right, top, bottom) should already
 *   be adjusted for this: left=-0.5; right = window.cols-0.5;
 */

static void slowline(double x1, double y1, double x2, double y2)
{
    double dx, dy;
    double m, b;
    int xstart, xstop, ystart, ystop;

    dx = x2 - x1;
    dy = y2 - y1;

    if (fabs(dx) > fabs(dy)) {
	m = dy / dx;
	b = y1 - m * x1;

	if (x1 > x2) {
	    xstart = iceil(x2 - 0.5);
	    xstop = ifloor(x1 + 0.5);
	}
	else {
	    xstart = iceil(x1 - 0.5);
	    xstop = ifloor(x2 + 0.5);
	}
	if (xstart <= xstop) {
	    ystart = ifloor(m * xstart + b + 0.5);
	    st->move(xstart, ystart);
	    while (xstart <= xstop) {
		st->cont(xstart++, ystart);
		ystart = ifloor(m * xstart + b + 0.5);
	    }
	}
    }
    else {
	if (dx == dy)		/* they both might be 0 */
	    m = 1.;
	else
	    m = dx / dy;
	b = x1 - m * y1;

	if (y1 > y2) {
	    ystart = iceil(y2 - 0.5);
	    ystop = ifloor(y1 + 0.5);
	}
	else {
	    ystart = iceil(y1 - 0.5);
	    ystop = ifloor(y2 + 0.5);
	}
	if (ystart <= ystop) {
	    xstart = ifloor(m * ystart + b + 0.5);
	    st->move(xstart, ystart);
	    while (ystart <= ystop) {
		st->cont(xstart, ystart++);
		xstart = ifloor(m * ystart + b + 0.5);
	    }
	}
    }
}

static void plot_line(double east1, double north1, double east2, double north2,
		      void (*line)(double, double, double, double))
{
    double x1, x2, y1, y2;

    y1 = Y(north1);
    y2 = Y(north2);

    if (st->window.proj == PROJECTION_LL) {
	if (east1 > east2)
	    while ((east1 - east2) > 180)
		east2 += 360;
	else if (east2 > east1)
	    while ((east2 - east1) > 180)
		east1 += 360;
	while (east1 > st->window.east) {
	    east1 -= 360.0;
	    east2 -= 360.0;
	}
	while (east1 < st->window.west) {
	    east1 += 360.0;
	    east2 += 360.0;
	}
	x1 = X(east1);
	x2 = X(east2);

	line(x1, y1, x2, y2);

	if (east2 > st->window.east || east2 < st->window.west) {
	    while (east2 > st->window.east) {
		east1 -= 360.0;
		east2 -= 360.0;
	    }
	    while (east2 < st->window.west) {
		east1 += 360.0;
		east2 += 360.0;
	    }
	    x1 = X(east1);
	    x2 = X(east2);
	    line(x1, y1, x2, y2);
	}
    }
    else {
	x1 = X(east1);
	x2 = X(east2);
	line(x1, y1, x2, y2);
    }
}

static double wrap_east(double e0, double e1)
{
    while (e0 - e1 > 180)
	e1 += 360.0;
    while (e1 - e0 > 180)
	e1 -= 360.0;

    return e1;
}


/*!
 * \brief Plot filled polygon with n vertices
 *
 * The polygon, described by the <i>n</i> vertices
 * <i>east,north</i>, is plotted in the output x,y space as a filled polygon.
 *
 * \param x coordinates of vertices
 * \param y coordinates of vertices
 * \param n number of vertices
 *
 * \return 0 on success
 * \return 2 n < 3
 * \return -1 weird internal error
 * \return 1 no memory
 */
int G_plot_polygon(const double *x, const double *y, int n)
{
    int i;
    int pole;
    double x0, x1;
    double y0, y1;
    double shift, E, W = 0L;
    double e0, e1;
    int shift1, shift2;

    if (!st->row_fill)
	st->row_fill = row_solid_fill;

    if (n < 3)
	return TOO_FEW_EDGES;

    /* traverse the perimeter */

    st->np = 0;
    shift1 = 0;

    /* global wrap-around for lat-lon, part1 */
    if (st->window.proj == PROJECTION_LL) {
	/*
	   pole = G_pole_in_polygon(x,y,n);
	 */
	pole = 0;

	e0 = x[n - 1];
	E = W = e0;

	x0 = X(e0);
	y0 = Y(y[n - 1]);

	if (pole && !edge(x0, y0, x0, Y(90.0 * pole)))
	    return NO_MEMORY;

	for (i = 0; i < n; i++) {
	    e1 = wrap_east(e0, x[i]);
	    if (e1 > E)
		E = e1;
	    if (e1 < W)
		W = e1;

	    x1 = X(e1);
	    y1 = Y(y[i]);

	    if (!edge(x0, y0, x1, y1))
		return NO_MEMORY;

	    x0 = x1;
	    y0 = y1;
	    e0 = e1;
	}
	if (pole && !edge(x0, y0, x0, Y(90.0 * pole)))
	    return NO_MEMORY;

	shift = 0;		/* shift into window */
	while (E + shift > st->window.east)
	    shift -= 360.0;
	while (E + shift < st->window.west)
	    shift += 360.0;
	shift1 = X(x[n - 1] + shift) - X(x[n - 1]);
    }
    else {
	x0 = X(x[n - 1]);
	y0 = Y(y[n - 1]);

	for (i = 0; i < n; i++) {
	    x1 = X(x[i]);
	    y1 = Y(y[i]);
	    if (!edge(x0, y0, x1, y1))
		return NO_MEMORY;
	    x0 = x1;
	    y0 = y1;
	}
    }

    /* check if perimeter has odd number of points */
    if (st->np & 1) {
	G_warning("Weird internal error: perimeter has odd number of points");
	return OUT_OF_SYNC;
    }

    /* sort the edge points by col(x) and then by row(y) */
    qsort(st->P, st->np, sizeof(POINT), edge_order);

    /* plot */
    for (i = 1; i < st->np; i += 2) {
	if (st->P[i].y != st->P[i - 1].y) {
	    G_warning("Weird internal error: edge leaves row");
	    return OUT_OF_SYNC;
	}
	st->row_fill(st->P[i].y, st->P[i - 1].x + shift1, st->P[i].x + shift1);
    }
    if (st->window.proj == PROJECTION_LL) {	/* now do wrap-around, part 2 */
	shift = 0;
	while (W + shift < st->window.west)
	    shift += 360.0;
	while (W + shift > st->window.east)
	    shift -= 360.0;
	shift2 = X(x[n - 1] + shift) - X(x[n - 1]);
	if (shift2 != shift1) {
	    for (i = 1; i < st->np; i += 2) {
		st->row_fill(st->P[i].y, st->P[i - 1].x + shift2, st->P[i].x + shift2);
	    }
	}
    }
    return OK;
}

/*!
 * \brief Plot multiple polygons
 *
 * Like G_plot_polygon(), except it takes a set of polygons, each with
 * npts[<i>i</i>] vertices, where the number of polygons is specified
 * with the <i>rings</i> argument. It is especially useful for
 * plotting vector areas with interior islands.
 *
 * \param xs pointer to pointer for X's
 * \param ys pointer to pointer for Y's
 * \param rpnts array of ints w/ num points per ring
 * \param rings number of rings
 *
 * \return 0 on success
 * \return 2 n < 3
 * \return -1 weird internal error
 * \return 1 no memory
 */
int G_plot_area(double *const *xs, double *const *ys, int *rpnts, int rings)
{
    int i, j, n;
    int pole;
    double x0, x1, *x;
    double y0, y1, *y;
    double shift, E, W = 0L;
    double e0, e1;
    int *shift1 = NULL, shift2;

    if (!st->row_fill)
	st->row_fill = row_solid_fill;

    /* traverse the perimeter */

    st->np = 0;
    shift1 = (int *)G_calloc(sizeof(int), rings);

    for (j = 0; j < rings; j++) {
	n = rpnts[j];

	if (n < 3)
	    return TOO_FEW_EDGES;

	x = xs[j];
	y = ys[j];

	/* global wrap-around for lat-lon, part1 */
	if (st->window.proj == PROJECTION_LL) {
	    /*
	       pole = G_pole_in_polygon(x,y,n);
	     */
	    pole = 0;

	    e0 = x[n - 1];
	    E = W = e0;

	    x0 = X(e0);
	    y0 = Y(y[n - 1]);

	    if (pole && !edge(x0, y0, x0, Y(90.0 * pole)))
		return NO_MEMORY;

	    for (i = 0; i < n; i++) {
		e1 = wrap_east(e0, x[i]);
		if (e1 > E)
		    E = e1;
		if (e1 < W)
		    W = e1;

		x1 = X(e1);
		y1 = Y(y[i]);

		if (!edge(x0, y0, x1, y1))
		    return NO_MEMORY;

		x0 = x1;
		y0 = y1;
		e0 = e1;
	    }
	    if (pole && !edge(x0, y0, x0, Y(90.0 * pole)))
		return NO_MEMORY;

	    shift = 0;		/* shift into window */
	    while (E + shift > st->window.east)
		shift -= 360.0;
	    while (E + shift < st->window.west)
		shift += 360.0;
	    shift1[j] = X(x[n - 1] + shift) - X(x[n - 1]);
	}
	else {
	    x0 = X(x[n - 1]);
	    y0 = Y(y[n - 1]);

	    for (i = 0; i < n; i++) {
		x1 = X(x[i]);
		y1 = Y(y[i]);
		if (!edge(x0, y0, x1, y1))
		    return NO_MEMORY;
		x0 = x1;
		y0 = y1;
	    }
	}
    }				/* for() */

    /* check if perimeter has odd number of points */
    if (st->np & 1) {
	G_warning("Weird internal error: perimeter has odd number of points");
	return OUT_OF_SYNC;
    }

    /* sort the edge points by col(x) and then by row(y) */
    qsort(st->P, st->np, sizeof(POINT), &edge_order);

    /* plot */
    for (j = 0; j < rings; j++) {
	for (i = 1; i < st->np; i += 2) {
	    if (st->P[i].y != st->P[i - 1].y) {
		G_warning("Weird internal error: edge leaves row");
		return OUT_OF_SYNC;
	    }
	    st->row_fill(st->P[i].y, st->P[i - 1].x + shift1[j], st->P[i].x + shift1[j]);
	}
	if (st->window.proj == PROJECTION_LL) {	/* now do wrap-around, part 2 */
	    n = rpnts[j];
	    x = xs[j];
	    y = ys[j];

	    shift = 0;
	    while (W + shift < st->window.west)
		shift += 360.0;
	    while (W + shift > st->window.east)
		shift -= 360.0;
	    shift2 = X(x[n - 1] + shift) - X(x[n - 1]);
	    if (shift2 != shift1[j]) {
		for (i = 1; i < st->np; i += 2) {
		    st->row_fill(st->P[i].y, st->P[i - 1].x + shift2, st->P[i].x + shift2);
		}
	    }
	}
    }
    G_free(shift1);
    return OK;

}

static int edge(double x0, double y0, double x1, double y1)
{
    double m, d;
    double x;
    int ystart, ystop;
    int exp;

    /* tolerance to avoid FPE */
    d = GRASS_EPSILON;
    if (y0 != y1) {
	if (fabs(y0) > fabs(y1))
	    d = fabs(y0);
	else
	    d = fabs(y1);

	d = frexp(d, &exp);
	exp -= 53;
	d = ldexp(d, exp);
    }

    if (fabs(y0 - y1) < d)
	return 1;

    if (y0 < y1) {
	ystart = iceil(y0);
	ystop = ifloor(y1);
	if (ystop == y1)
	    ystop--;		/* if line stops at row center, don't include point */
    }
    else {
	ystart = iceil(y1);
	ystop = ifloor(y0);
	if (ystop == y0)
	    ystop--;		/* if line stops at row center, don't include point */
    }

    if (ystart > ystop)
	return 1;		/* does not cross center line of row */

    m = (x0 - x1) / (y0 - y1);
    x = m * (ystart - y0) + x0;
    while (ystart <= ystop) {
	if (!edge_point(x, ystart++))
	    return 0;
	x += m;
    }

    return 1;
}

static int edge_point(double x, int y)
{

    if (y < st->ymin || y > st->ymax)
	return 1;
    if (st->np >= st->npalloc) {
	if (st->npalloc > 0) {
	    st->npalloc *= 2;
	    st->P = (POINT *) G_realloc(st->P, st->npalloc * sizeof(POINT));
	}
	else {
	    st->npalloc = 32;
	    st->P = (POINT *) G_malloc(st->npalloc * sizeof(POINT));
	}
	if (st->P == NULL) {
	    st->npalloc = 0;
	    return 0;
	}
    }
    st->P[st->np].x = x;
    st->P[st->np++].y = y;
    return 1;
}

static int edge_order(const void *aa, const void *bb)
{
    const struct point *a = aa, *b = bb;

    if (a->y < b->y)
	return (-1);
    if (a->y > b->y)
	return (1);

    if (a->x < b->x)
	return (-1);
    if (a->x > b->x)
	return (1);

    return (0);
}

static void row_solid_fill(int y, double x1, double x2)
{
    int i1, i2;

    i1 = iceil(x1);
    i2 = ifloor(x2);
    if (i1 <= i2) {
	st->move(i1, y);
	st->cont(i2, y);
    }
}

static void row_dotted_fill(int y, double x1, double x2)
{
    int i1, i2, i;

    if (y != iceil(y / st->dotted_fill_gap) * st->dotted_fill_gap)
	return;

    i1 = iceil(x1 / st->dotted_fill_gap) * st->dotted_fill_gap;
    i2 = ifloor(x2);
    if (i1 <= i2) {
	for (i = i1; i <= i2; i += st->dotted_fill_gap) {
	    st->move(i, y);
	    st->cont(i, y);
	}
    }
}

static int ifloor(double x)
{
    int i;

    i = (int)x;
    if (i > x)
	i--;
    return i;
}

static int iceil(double x)
{
    int i;

    i = (int)x;
    if (i < x)
	i++;
    return i;
}

/*!
 * \brief Plot f(east1) to f(east2)
 *
 * The function <i>f(east)</i> is plotted from <i>east1</i> to
 * <i>east2</i>. The function <i>f(east)</i> must return the map
 * northing coordinate associated with east.
 *
 * \param f plotting function
 * \param east1 easting (first point)
 * \param east2 easting (second point)
 */
void G_plot_fx(double (*f) (double), double east1, double east2)
{
    double east, north, north1;
    double incr;


    incr = fabs(1.0 / st->xconv);

    east = east1;
    north = f(east1);

    if (east1 > east2) {
	while ((east1 -= incr) > east2) {
	    north1 = f(east1);
	    G_plot_line(east, north, east1, north1);
	    north = north1;
	    east = east1;
	}
    }
    else {
	while ((east1 += incr) < east2) {
	    north1 = f(east1);
	    G_plot_line(east, north, east1, north1);
	    north = north1;
	    east = east1;
	}
    }

    G_plot_line(east, north, east2, f(east2));
}