Loading...
Searching...
No Matches
plot.c
Go to the documentation of this file.
/*!
* \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 (*)(double, double, double, double));
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));
}
G_free
void G_free(void *buf)
Free allocated memory.
Definition alloc.c:147
NULL
#define NULL
Definition ccmath.h:32
b
double b
Definition driver/set_window.c:5
l
double l
Definition driver/set_window.c:5
t
double t
Definition driver/set_window.c:5
r
double r
Definition driver/set_window.c:5
G_warning
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition gis/error.c:203
G_get_set_window
void G_get_set_window(struct Cell_head *window)
Get the current working window (region).
Definition gis/set_window.c:35
G_plot_area
int G_plot_area(double *const *xs, double *const *ys, int *rpnts, int rings)
Plot multiple polygons.
Definition plot.c:502
G_plot_line2
void G_plot_line2(double east1, double north1, double east2, double north2)
Plot line between latlon coordinates (slowline).
Definition plot.c:223
G_plot_polygon
int G_plot_polygon(const double *x, const double *y, int n)
Plot filled polygon with n vertices.
Definition plot.c:372
OUT_OF_SYNC
#define OUT_OF_SYNC
Definition plot.c:73
G_plot_where_en
void G_plot_where_en(int x, int y, double *east, double *north)
Converts x,y to east,north.
Definition plot.c:177
X
#define X(e)
Definition plot.c:143
POINT
#define POINT
Definition plot.c:50
G_plot_where_xy
void G_plot_where_xy(double east, double north, int *x, int *y)
Converts east,north to x,y.
Definition plot.c:160
G_setup_fill
void G_setup_fill(int gap)
Set row_fill routine to row_solid_fill or row_dotted_fill.
Definition plot.c:133
G_setup_plot
void G_setup_plot(double t, double b, double l, double r, int(*Move)(int, int), int(*Cont)(int, int))
Initialize plotting routines.
Definition plot.c:97
NORTH
#define NORTH(y)
Definition plot.c:147
G_plot_fx
void G_plot_fx(double(*f)(double), double east1, double east2)
Plot f(east1) to f(east2).
Definition plot.c:780
G_plot_line
void G_plot_line(double east1, double north1, double east2, double north2)
Plot line between latlon coordinates (fastline).
Definition plot.c:208
Y
#define Y(n)
Definition plot.c:144
NO_MEMORY
#define NO_MEMORY
Definition plot.c:72
OK
#define OK
Definition plot.c:70
G_plot_point
void G_plot_point(double east, double north)
Plot point.
Definition plot.c:189
EAST
#define EAST(x)
Definition plot.c:146
TOO_FEW_EDGES
#define TOO_FEW_EDGES
Definition plot.c:71
G_adjust_easting
double G_adjust_easting(double east, const struct Cell_head *window)
Returns east not smaller than west.
Definition window_map.c:58
x
#define x