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

/*****************************************************************
 * Plot lines and filled polygons. Input space is database 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 databases.
 *
 * 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.
 *******************************************************************/
#include <stdlib.h>
#include <math.h>
#include <grass/gis.h>

static double xconv, yconv;
static double left, right, top, bottom;
static int ymin, ymax;
static struct Cell_head window;
static int fastline(double, double, double, double);
static int slowline(double, double, double, double);
static int plot_line(double, double, double, double, int (*)());
static double nearest(double, double);
static int edge(double, double, double, double);
static int edge_point(double, int);

#define POINT struct point
POINT {
    double x;
    int y;
};
static int edge_order(const void *, const void *);
static int row_solid_fill(int, double, double);
static int row_dotted_fill(int, double, double);
static int dotted_fill_gap = 2;
static int ifloor(double);
static int iceil(double);
static int (*row_fill) () = row_solid_fill;
static int (*move) (int, int);
static int (*cont) (int, int);

/*!
 * \brief returns east larger than west
 *
 * If the region projection is
 * PROJECTION_LL, then this routine returns an equivalent <b>east</b> that is
 * larger, but no more than 360 degrees larger, than the coordinate for the
 * western edge of the region. Otherwise no adjustment is made and the original
 * <b>east</b> is returned.
 *
 *  \param east
 *  \param region
 *  \return double
 */

/*
 * G_setup_plot (t, b, l, r, Move, Cont)
 *     double t, b, l, r;
 *     int (*Move)(), (*Cont)();
 *
 * initialize the plotting capability.
 *    t,b,l,r:   top, bottom, left, right of the output x,y coordinate space.
 *    Move,Cont: subroutines that will draw lines in x,y space.
 *       Move(x,y)   move to x,y (no draw)
 *       Cont(x,y)   draw from previous position to x,y
 * Notes:
 *   Cont() is responsible for clipping.
 *   The t,b,l,r are only used to compute coordinate transformations.
 *   The input space is assumed to be the current GRASS window.
 */

/*!
 * \brief initialize plotting routines
 *
 * Initializes the plotting
 * capability. This routine must be called once before calling the
 * <b>G_plot_*(~)</b> routines described below.
 * The parameters <b>t, b, l, r</b> 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 ~
 *  \return int
 */

int G_setup_plot(double t, double b, double l, double r,
		 int (*Move) (int, int), int (*Cont) (int, int))
{
    G_get_set_window(&window);

    left = l;
    right = r;
    top = t;
    bottom = b;

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

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

    move = Move;
    cont = Cont;

    return 0;
}

/*!
 * \brief set row_fill routine to row_solid_fill or row_dotted_fill
 *
 * After calling this function, <b>G_plot_polygon()</b> and
 * <b>G_plot_area()</b> 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 int
 *  \return int
 */
int G_setup_fill(int gap)
{
    if (gap > 0) {
	row_fill = row_dotted_fill;
	dotted_fill_gap = gap + 1;
    }
    else
	row_fill = row_solid_fill;

    return 0;
}

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

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


/*!
 * \brief east,north to x,y
 *
 * The map coordinates <b>east,north</b> are converted
 * to pixel coordinates <b>x,y.</b>
 *
 *  \param east
 *  \param north
 *  \param x
 *  \param y
 *  \return int
 */

int G_plot_where_xy(double east, double north, int *x, int *y)
{
    *x = ifloor(X(G_adjust_easting(east, &window)) + 0.5);
    *y = ifloor(Y(north) + 0.5);

    return 0;
}


/*!
 * \brief x,y to east,north
 *
 * The pixel coordinates <b>x,y</b> are converted to map
 * coordinates <b>east,north.</b>
 *
 *  \param x
 *  \param y
 *  \param east
 *  \param north
 *  \return int
 */

int G_plot_where_en(int x, int y, double *east, double *north)
{
    *east = G_adjust_easting(EAST(x), &window);
    *north = NORTH(y);

    return 0;
}

int G_plot_point(double east, double north)
{
    int x, y;

    G_plot_where_xy(east, north, &x, &y);
    move(x, y);
    cont(x, y);

    return 0;
}

/*
 * Line in map coordinates is plotted in output x,y coordinates
 * This routine handles global wrap-around for lat-long databses.
 *
 */

/*!
 * \brief plot line between latlon coordinates
 *
 * A line from <b>east1,north1</b>
 * to <b>east2,north2</b> is plotted in output x,y coordinates (e.g. pixels for
 * graphics.) This routine handles global wrap-around for latitude-longitude
 * databases.
 *
 *  \param east1
 *  \param north1
 *  \param east2
 *  \param north2
 *  \return int
 */

int G_plot_line(double east1, double north1, double east2, double north2)
{
    return plot_line(east1, north1, east2, north2, fastline);
}

int G_plot_line2(double east1, double north1, double east2, double north2)
{
    return 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 int fastline(double x1, double y1, double x2, double y2)
{
    move(ifloor(x1 + 0.5), ifloor(y1 + 0.5));
    cont(ifloor(x2 + 0.5), ifloor(y2 + 0.5));

    return 0;
}

/* 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 int 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);
	    move(xstart, ystart);
	    while (xstart <= xstop) {
		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);
	    move(xstart, ystart);
	    while (ystart <= ystop) {
		cont(xstart, ystart++);
		xstart = ifloor(m * ystart + b + 0.5);
	    }
	}
    }

    return 0;
}

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

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

    if (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 > window.east) {
	    east1 -= 360.0;
	    east2 -= 360.0;
	}
	while (east1 < window.west) {
	    east1 += 360.0;
	    east2 += 360.0;
	}
	x1 = X(east1);
	x2 = X(east2);

	line(x1, y1, x2, y2);

	if (east2 > window.east || east2 < window.west) {
	    while (east2 > window.east) {
		east1 -= 360.0;
		east2 -= 360.0;
	    }
	    while (east2 < 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);
    }

    return 0;
}

/*
 * G_plot_polygon (x, y, n)
 * 
 *    double *x       x coordinates of vertices
 *    double *y       y coordinates of vertices
 *    int n           number of verticies
 *
 * polygon fill from map coordinate space to plot x,y space.
 * for lat-lon, handles global wrap-around as well as polar polygons.
 *
 * returns 0 ok, 2 n<3, -1 weird internal error, 1 no memory
 */

static POINT *P;
static int np;
static int npalloc = 0;

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

static double nearest(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 <b>n</b> vertices
 * <b>east,north</b>, is plotted in the output x,y space as a filled polygon.
 *
 *  \param east
 *  \param north
 *  \param n
 *  \return int
 */

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 (n < 3)
	return TOO_FEW_EDGES;

    /* traverse the perimeter */

    np = 0;
    shift1 = 0;

    /* global wrap-around for lat-lon, part1 */
    if (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 = nearest(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 > window.east)
	    shift -= 360.0;
	while (E + shift < 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 (np % 2)
	return OUT_OF_SYNC;

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

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

/*
 * G_plot_area (xs, ys, rpnts, rings)
 *      double **xs;  -- pointer to pointer for X's
 *      double **ys;  -- pointer to pointer for Y's
 *      int *rpnts;   -- array of ints w/ num points per ring
 *      int rings;    -- number of rings
 *
 * Essentially a copy of G_plot_polygon, with minor mods to
 * handle a set of polygons.  return values are the same.
 */

/*!
 * \brief plot multiple polygons
 *
 * Like G_plot_polygon, except it takes a set of polygons,
 * each with \textbf{npts[<i>i</i>]} vertices, where the number of polygons 
 * is specified with the <b>rings</b> argument.  It is especially useful for 
 * plotting vector areas with interior islands.
 *
 *  \param xs
 *  \param ys
 *  \param npts
 *  \param rings
 *  \return int
 */

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;

    /* traverse the perimeter */

    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 (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 = nearest(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 > window.east)
		shift -= 360.0;
	    while (E + shift < 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 (np % 2)
	return OUT_OF_SYNC;

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

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

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

}

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


    /* tolerance to avoid FPE */
    dy = y0 - y1;
    if (fabs(dy) < 1e-10)
	return 1;

    m = (x0 - x1) / dy;

    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 */

    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 < ymin || y > ymax)
	return 1;
    if (np >= npalloc) {
	if (npalloc > 0) {
	    npalloc *= 2;
	    P = (POINT *) G_realloc(P, npalloc * sizeof(POINT));
	}
	else {
	    npalloc = 32;
	    P = (POINT *) G_malloc(npalloc * sizeof(POINT));
	}
	if (P == NULL) {
	    npalloc = 0;
	    return 0;
	}
    }
    P[np].x = x;
    P[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 int row_solid_fill(int y, double x1, double x2)
{
    int i1, i2;

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

    return 0;
}

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

    if (y != iceil(y / dotted_fill_gap) * dotted_fill_gap)
	return 0;

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

    return 0;
}

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;
}

/*
 * G_plot_fx(e1,e2)
 *
 * plot f(x) from x=e1 to x=e2
 */


/*!
 * \brief plot f(east1) to f(east2)
 *
 * The function <b>f(east)</b> is plotted from
 * <b>east1</b> to <b>east2.</b> The function <b>f(east)</b> must return
 * the map northing coordinate associated with east.
 *
 *  \param ~
 *  \return int
 */

int G_plot_fx(double (*f) (double), double east1, double east2)
{
    double east, north, north1;
    double incr;


    incr = fabs(1.0 / 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));

    return 0;
}