[GRASSLIST:10387] Re: G_plot_polygon() suggestion for vector areas
Huidae Cho
grass4u at gmail.com
Fri Feb 17 20:17:22 EST 2006
I've added G_setup_fill() function to fill a polygon with dots. I chose dotted
fill over dashed fill because of kind of the moire effect. You can see a
screenshot of d.vect with "fill=2 width=1" option. "fill=2 width=1" means that
it fills areas with 1 pixel dots and 2 pixel skip.
There is a limitation: if two vectors are overlaid with the same fill, only the
latter will be shown.
Please find attached the figure.
Huidae Cho
On Fri, Feb 17, 2006 at 11:45:25PM +0100, Markus Neteler wrote:
> Hi Huidae,
>
> you may remember the hatching of vector maps:
> do you consider to integrate this addition somehow in CVS?
>
> Markus
>
> On Sun, Aug 21, 2005 at 03:14:20PM -0500, Huidae Cho wrote:
> > In your screenshot, each symbol has a different color, so some symbols
> > are nearly invisible. I think this is a kind of side effect of dashed
> > fill because symbols are drawn in black color after filling polygons.
> >
> > d.vect width=2 draws BLACK thick symbols.
> >
> > Huidae Cho
> >
> > On Sun, Aug 21, 2005 at 05:38:15PM +0200, Markus Neteler wrote:
> > > Sorry, I hit the "send button" too quickly.
> > > Attached a screenshot. The basic/x symbols come nicely,
> > > it probably depends on the screen size.
> > >
> > > Markus
> > >
> > >
> > > On Sun, Aug 21, 2005 at 05:35:11PM +0200, Markus Neteler wrote:
> > > > Thanks, I have tried it with spearfish
> > > >
> > > >
> > > > d.rast elevation.dem
> > > > d.vect fields -c
> > > > d.vect archsites icon=basic/x
> > > >
> > > > Looks very good to me!
> > > >
> > > > On Thu, Aug 18, 2005 at 12:02:10PM -0500, Huidae Cho wrote:
> > > > > Find attached lib/gis/plot.c. You will find some problems with thin
> > > > > lines drawn on dash filled area. For example, some basic/x symbols are
> > > > > almost invisible when drawn on filled area.
> > > > >
> > > > > The row_fill function pointer is set to row_dashed_fill by default. Give
> > > > > it a try.
> > > > >
> > > > > Huidae Cho
> > > > >
> > > > >
> > > > > On Thu, Aug 18, 2005 at 10:08:01AM +0200, Markus Neteler wrote:
> > > > > > Hi Huidae,
> > > > > >
> > > > > > On Wed, Aug 17, 2005 at 05:26:21PM -0500, Huidae Cho wrote:
> > > > > > > Hi Markus,
> > > > > > >
> > > > > > > My first name is Huidae and the last name is Cho, so "Huidae Cho" is
> > > > > > > already western style. Just call me Huidae /heedae/.
> > > > > >
> > > > > > I see - thanks for the update!
> > > > > >
> > > > > > > Slanted dashed line is difficult to me, instead, I applied horizontal
> > > > > > > dashed lines to row_fill() and got the attached image.
> > > > > > >
> > > > > > > What do you think?
> > > > > >
> > > > > > Cool - I like it :-) Exactly doing the job...
> > > > > >
> > > > > > Markus
> > > > > >
> > > > > > > Huidae Cho
> > > > > > >
> > > > > > > On Wed, Aug 17, 2005 at 08:52:45PM +0200, Markus Neteler wrote:
> > > > > > > > Hi Huidae Cho
> > > > > > > >
> > > > > > > > (first of all, I would like to ask you if I should call
> > > > > > > > you huidae or Cho or Huidae Cho? I know that Chinese
> > > > > > > > people have the problem that the Western people always
> > > > > > > > confuse their name. However, I don't really know the order
> > > > > > > > of yours - sorry, I am asking a bit late :-)
> > > > > > > >
> > > > > > > > I was thinking of vector area drawing in d.vect.
> > > > > > > > I got the impression that you are fairly familiar with
> > > > > > > > the R_*() functions due to the latest great changes.
> > > > > > > > Dreaming of semi-transparent vector areas, I thought if
> > > > > > > > it would be possible to clone G_plot_polygon() in
> > > > > > > > lib/gis/plot.c to something like G_plot_polygon_stripes()
> > > > > > > > which would render an area such as
> > > > > > > >
> > > > > > > > +---------+
> > > > > > > > |/////////|
> > > > > > > > |/////////|
> > > > > > > > |/////////|
> > > > > > > > +---------+
> > > > > > > >
> > > > > > > > so that we could still see a bit of the underlying map
> > > > > > > > such as DEM or sat image.
> > > > > > > >
> > > > > > > > Do you think that's feasible?
> > > > > > > >
> > > > > > > > Markus
> > > > > > > >
> > > > > >
> > > > > >
> > > > > >
> > > > > > --
> > > > > > Markus Neteler <neteler itc it> http://mpa.itc.it
> > > > > > ITC-irst - Centro per la Ricerca Scientifica e Tecnologica
> > > > > > MPBA - Predictive Models for Biol. & Environ. Data Analysis
> > > > > > Via Sommarive, 18 - 38050 Povo (Trento), Italy
> > > >
> > > > > /*****************************************************************
> > > > > * 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 "gis.h"
> > > > > #include <stdlib.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_dashed_fill(int,double,double);
> > > > > static int ifloor(double);
> > > > > static int iceil(double);
> > > > > static int (*row_fill)() = row_dashed_fill;
> > > > > static int (*move)() = NULL;
> > > > > static int (*cont)() = NULL;
> > > > >
> > > > > static double fabs(double x)
> > > > > {
> > > > > return x>0?x:-x;
> > > > > }
> > > > >
> > > > > /*!
> > > > > * \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
> > > > > */
> > > > >
> > > > > extern double G_adjust_easting();
> > > > >
> > > > > /*
> > > > > * 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 (*Cont)())
> > > > > {
> > > > > 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
> > > > > *
> > > > > * After calling this function, <b>G_plot_polygon()</b> and
> > > > > * <b>G_plot_area()</b> fill shapes with solid (opaque) lines.
> > > > > *
> > > > > * \param void
> > > > > * \return int
> > > > > */
> > > > > int G_solid_fill(void)
> > > > > {
> > > > > row_fill = row_solid_fill;
> > > > >
> > > > > return 0;
> > > > > }
> > > > >
> > > > > /*!
> > > > > * \brief set row_fill routine to row_dashed_fill
> > > > > *
> > > > > * After calling this function, <b>G_plot_polygon()</b> and
> > > > > * <b>G_plot_area()</b> fill shapes with dashed (kind of transparent) lines.
> > > > > *
> > > > > * \param void
> > > > > * \return int
> > > > > */
> > > > > int G_dashed_fill(void)
> > > > > {
> > > > > row_fill = row_dashed_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 (east, north, x, y)
> > > > > double east, north;
> > > > > int *x, *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 (x, y, east, north)
> > > > > double *east, *north;
> > > > > {
> > > > > *east = G_adjust_easting(EAST(x),&window);
> > > > > *north = NORTH(y);
> > > > >
> > > > > return 0;
> > > > > }
> > > > >
> > > > > int G_plot_point (east, north)
> > > > > double east, 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 (east1, north1, east2, north2)
> > > > > double east1, north1, east2, north2;
> > > > > {
> > > > > int fastline();
> > > > > plot_line (east1, north1, east2, north2, fastline);
> > > > >
> > > > > return 0;
> > > > > }
> > > > >
> > > > > int G_plot_line2 (east1, north1, east2, north2)
> > > > > double east1, north1, east2, north2;
> > > > > {
> > > > > int slowline();
> > > > > plot_line (east1, north1, east2, north2, slowline);
> > > > >
> > > > > return 0;
> > > > > }
> > > > >
> > > > > /* 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 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 (
> > > > > double *x,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 **xs, double **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);
> > > > > }
> > > > > }
> > > > > }
> > > > > }
> > > > > free (shift1);
> > > > > return OK;
> > > > >
> > > > > }
> > > > >
> > > > > static int edge (double x0,double y0,double x1,double y1)
> > > > > {
> > > > > register 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, register int y)
> > > > > {
> > > > >
> > > > > if (y < ymin || y > ymax)
> > > > > return 1;
> > > > > if (np >= npalloc)
> > > > > {
> > > > > if (npalloc > 0)
> > > > > {
> > > > > npalloc *= 2;
> > > > > P = (POINT *) realloc (P, npalloc * sizeof (POINT));
> > > > > }
> > > > > else
> > > > > {
> > > > > npalloc = 32;
> > > > > P = (POINT *) 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_dashed_fill(int y,double x1,double x2)
> > > > > {
> > > > > int i1, i2, i;
> > > > >
> > > > > i1 = iceil(x1);
> > > > > i2 = ifloor(x2);
> > > > > if (i1 <= i2)
> > > > > {
> > > > > for(i = i1 + (i1 + y) % 2; i <= i2; i += 2)
> > > > > {
> > > > > 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 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;
> > > > > }
> > > >
> > > >
> > > > --
> > > > Markus Neteler <neteler itc it> http://mpa.itc.it
> > > > ITC-irst - Centro per la Ricerca Scientifica e Tecnologica
> > > > MPBA - Predictive Models for Biol. & Environ. Data Analysis
> > > > Via Sommarive, 18 - 38050 Povo (Trento), Italy
> > >
> > > --
> > > Markus Neteler <neteler itc it> http://mpa.itc.it
> > > ITC-irst - Centro per la Ricerca Scientifica e Tecnologica
> > > MPBA - Predictive Models for Biol. & Environ. Data Analysis
> > > Via Sommarive, 18 - 38050 Povo (Trento), Italy
> >
> >
>
> --
> Markus Neteler <neteler itc it> http://mpa.itc.it
> ITC-irst - Centro per la Ricerca Scientifica e Tecnologica
> MPBA - Predictive Models for Biol. & Environ. Data Analysis
> Via Sommarive, 18 - 38050 Povo (Trento), Italy
-------------- next part --------------
A non-text attachment was scrubbed...
Name: d.vect_fill.png
Type: image/png
Size: 117069 bytes
Desc: not available
Url : http://lists.osgeo.org/pipermail/grass-user/attachments/20060217/c224b487/d.vect_fill.png
More information about the grass-user
mailing list