plot.c

Go to the documentation of this file.
00001 /*****************************************************************
00002  * Plot lines and filled polygons. Input space is database window.
00003  * Output space and output functions are user defined.
00004  * Converts input east,north lines and polygons to output x,y
00005  * and calls user supplied line drawing routines to do the plotting.
00006  *
00007  * Handles global wrap-around for lat-lon databases.
00008  *
00009  * Does not perform window clipping.
00010  * Clipping must be done by the line draw routines supplied by the user.
00011  *
00012  * Note:
00013  *  Hopefully, cartographic style projection plotting will be added later.
00014  *******************************************************************/
00015 #include <grass/gis.h>
00016 #include <stdlib.h>
00017 
00018 static double xconv, yconv;
00019 static double left, right, top, bottom;
00020 static int ymin, ymax;
00021 static struct Cell_head window;
00022 static int fastline(double,double,double,double);
00023 static int slowline(double,double,double,double);
00024 static int plot_line(double,double,double,double,int (*)());
00025 static double nearest(double,double);
00026 static int edge (double,double,double,double);
00027 static int edge_point(double,int);
00028 
00029 #define POINT struct point
00030 POINT
00031 {
00032     double x;
00033     int y;
00034 };
00035 static int edge_order(const void *, const void *);
00036 static int row_solid_fill(int,double,double);
00037 static int row_dotted_fill(int,double,double);
00038 static int dotted_fill_gap = 2;
00039 static int ifloor(double);
00040 static int iceil(double);
00041 static int (*row_fill)() = row_solid_fill;
00042 static int (*move)() = NULL;
00043 static int (*cont)() = NULL;
00044 
00045 static double fabs(double x)
00046 {
00047     return x>0?x:-x;
00048 }
00049 
00064 extern double G_adjust_easting();
00065 
00066 /*
00067  * G_setup_plot (t, b, l, r, Move, Cont)
00068  *     double t, b, l, r;
00069  *     int (*Move)(), (*Cont)();
00070  *
00071  * initialize the plotting capability.
00072  *    t,b,l,r:   top, bottom, left, right of the output x,y coordinate space.
00073  *    Move,Cont: subroutines that will draw lines in x,y space.
00074  *       Move(x,y)   move to x,y (no draw)
00075  *       Cont(x,y)   draw from previous position to x,y
00076  * Notes:
00077  *   Cont() is responsible for clipping.
00078  *   The t,b,l,r are only used to compute coordinate transformations.
00079  *   The input space is assumed to be the current GRASS window.
00080  */
00081 
00103 int G_setup_plot (
00104     double t,double b,double l,double r,
00105     int (*Move)(),
00106     int (*Cont)())
00107 {
00108     G_get_set_window (&window);
00109 
00110     left = l;
00111     right = r;
00112     top = t;
00113     bottom = b;
00114 
00115     xconv = (right-left)/(window.east-window.west);
00116     yconv = (bottom-top)/(window.north-window.south);
00117 
00118     if (top < bottom)
00119     {
00120         ymin = iceil(top);
00121         ymax = ifloor(bottom);
00122     }
00123     else
00124     {
00125         ymin = iceil(bottom);
00126         ymax = ifloor(top);
00127     }
00128 
00129     move = Move;
00130     cont = Cont;
00131 
00132     return 0;
00133 }
00134 
00146 int G_setup_fill(int gap)
00147 {
00148     if (gap > 0) {
00149         row_fill = row_dotted_fill;
00150         dotted_fill_gap = gap + 1;
00151     } else
00152         row_fill = row_solid_fill;
00153 
00154     return 0;
00155 }
00156 
00157 #define X(e) (left + xconv * ((e) - window.west))
00158 #define Y(n) (top + yconv * (window.north - (n)))
00159 
00160 #define EAST(x) (window.west + ((x)-left)/xconv)
00161 #define NORTH(y) (window.north - ((y)-top)/yconv)
00162 
00163 
00177 int G_plot_where_xy  (double east, double north, int *x, int *y)
00178 
00179 {
00180     *x = ifloor(X(G_adjust_easting(east,&window))+0.5);
00181     *y = ifloor(Y(north)+0.5);
00182 
00183     return 0;
00184 }
00185 
00186 
00200 int G_plot_where_en  (int x, int  y, double *east, double *north)
00201 
00202 {
00203     *east = G_adjust_easting(EAST(x),&window);
00204     *north = NORTH(y);
00205 
00206     return 0;
00207 }
00208 
00209 int G_plot_point  (double east, double north)
00210 
00211 {
00212     int x,y;
00213 
00214     G_plot_where_xy(east,north,&x,&y);
00215     move (x,y);
00216     cont (x,y);
00217 
00218     return 0;
00219 }
00220 
00221 /*
00222  * Line in map coordinates is plotted in output x,y coordinates
00223  * This routine handles global wrap-around for lat-long databses.
00224  *
00225  */
00226 
00242 int G_plot_line  (double east1, double north1, double east2, double north2)
00243 
00244 {
00245     int fastline();
00246     plot_line (east1, north1, east2, north2, fastline);
00247 
00248     return 0;
00249 }
00250 
00251 int G_plot_line2  (double east1, double north1, double east2, double north2)
00252 
00253 {
00254     int slowline();
00255     plot_line (east1, north1, east2, north2, slowline);
00256 
00257     return 0;
00258 }
00259 
00260 /* fastline converts double rows/cols to ints then plots
00261  * this is ok for graphics, but not the best for vector to raster
00262  */
00263 
00264 static int fastline(double x1,double y1,double x2,double y2)
00265 {
00266     move (ifloor(x1+0.5),ifloor(y1+0.5));
00267     cont (ifloor(x2+0.5),ifloor(y2+0.5));
00268 
00269     return 0;
00270 }
00271 
00272 /* NOTE (shapiro): 
00273  *   I think the adding of 0.5 in slowline is not correct
00274  *   the output window (left, right, top, bottom) should already
00275  *   be adjusted for this: left=-0.5; right = window.cols-0.5;
00276  */
00277 
00278 static int slowline(double x1,double y1,double x2,double y2)
00279 {
00280     double dx, dy;
00281     double m,b;
00282     int xstart, xstop, ystart, ystop;
00283 
00284     dx = x2-x1;
00285     dy = y2-y1;
00286 
00287     if (fabs(dx) > fabs(dy))
00288     {
00289         m = dy/dx;
00290         b = y1 - m * x1;
00291 
00292         if (x1 > x2)
00293         {
00294             xstart = iceil (x2-0.5);
00295             xstop  = ifloor (x1+0.5);
00296         }
00297         else
00298         {
00299             xstart = iceil (x1-0.5);
00300             xstop = ifloor (x2+0.5);
00301         }
00302         if (xstart <= xstop)
00303         {
00304             ystart = ifloor(m * xstart + b + 0.5);
00305             move (xstart, ystart);
00306             while (xstart <= xstop)
00307             {
00308                 cont (xstart++, ystart);
00309                 ystart = ifloor(m * xstart + b + 0.5);
00310             }
00311         }
00312     }
00313     else
00314     {
00315         if(dx==dy) /* they both might be 0 */
00316            m = 1.;
00317         else 
00318            m = dx/dy;
00319         b = x1 - m * y1;
00320 
00321         if (y1 > y2)
00322         {
00323             ystart = iceil (y2-0.5);
00324             ystop  = ifloor (y1+0.5);
00325         }
00326         else
00327         {
00328             ystart = iceil (y1-0.5);
00329             ystop = ifloor (y2+0.5);
00330         }
00331         if (ystart <= ystop)
00332         {
00333             xstart = ifloor(m * ystart + b + 0.5);
00334             move (xstart, ystart);
00335             while (ystart <= ystop)
00336             {
00337                 cont (xstart, ystart++);
00338                 xstart = ifloor(m * ystart + b + 0.5);
00339             }
00340         }
00341     }
00342 
00343     return 0;
00344 }
00345 
00346 static int plot_line(double east1,double north1,double east2,double north2,
00347     int (*line)())
00348 {
00349     double x1,x2,y1,y2;
00350 
00351     y1 = Y(north1);
00352     y2 = Y(north2);
00353 
00354     if (window.proj == PROJECTION_LL)
00355     {
00356         if (east1 > east2)
00357             while ((east1-east2) > 180)
00358                 east2 += 360;
00359         else if (east2 > east1)
00360             while ((east2-east1) > 180)
00361                 east1 += 360;
00362         while (east1 > window.east)
00363         {
00364             east1 -= 360.0;
00365             east2 -= 360.0;
00366         }
00367         while (east1 < window.west)
00368         {
00369             east1 += 360.0;
00370             east2 += 360.0;
00371         }
00372         x1 = X(east1);
00373         x2 = X(east2);
00374 
00375         line (x1,y1,x2,y2);
00376 
00377         if (east2 > window.east || east2 < window.west)
00378         {
00379             while (east2 > window.east)
00380             {
00381                 east1 -= 360.0;
00382                 east2 -= 360.0;
00383             }
00384             while (east2 < window.west)
00385             {
00386                 east1 += 360.0;
00387                 east2 += 360.0;
00388             }
00389             x1 = X(east1);
00390             x2 = X(east2);
00391             line (x1,y1,x2,y2);
00392         }
00393     }
00394     else
00395     {
00396         x1 = X(east1);
00397         x2 = X(east2);
00398         line (x1,y1,x2,y2);
00399     }
00400 
00401     return 0;
00402 }
00403 
00404 /*
00405  * G_plot_polygon (x, y, n)
00406  * 
00407  *    double *x       x coordinates of vertices
00408  *    double *y       y coordinates of vertices
00409  *    int n           number of verticies
00410  *
00411  * polygon fill from map coordinate space to plot x,y space.
00412  * for lat-lon, handles global wrap-around as well as polar polygons.
00413  *
00414  * returns 0 ok, 2 n<3, -1 weird internal error, 1 no memory
00415  */
00416 
00417 static POINT *P;
00418 static int np;
00419 static int npalloc = 0;
00420 
00421 #define OK 0
00422 #define TOO_FEW_EDGES 2
00423 #define NO_MEMORY 1
00424 #define OUT_OF_SYNC -1
00425 
00426 static double nearest(double e0,double e1)
00427 {
00428     while (e0 - e1 > 180)
00429         e1 += 360.0;
00430     while (e1 - e0 > 180)
00431         e1 -= 360.0;
00432     
00433     return e1;
00434 }
00435 
00436 
00449 int G_plot_polygon (
00450     double *x,double *y,
00451     int n)
00452 {
00453     int i;
00454     int pole;
00455     double x0,x1;
00456     double y0,y1;
00457     double shift,E,W=0L;
00458     double e0,e1;
00459     int shift1, shift2;
00460 
00461     if (n < 3)
00462         return TOO_FEW_EDGES;
00463 
00464 /* traverse the perimeter */
00465 
00466     np = 0;
00467     shift1 = 0;
00468 
00469 /* global wrap-around for lat-lon, part1 */
00470     if (window.proj == PROJECTION_LL)
00471     {
00472         /*
00473         pole = G_pole_in_polygon(x,y,n);
00474         */
00475         pole = 0;
00476 
00477         e0 = x[n-1];
00478         E = W = e0;
00479 
00480         x0 = X(e0);
00481         y0 = Y(y[n-1]);
00482 
00483         if (pole &&!edge (x0, y0, x0, Y(90.0*pole)))
00484                 return NO_MEMORY;
00485 
00486         for (i = 0; i < n; i++)
00487         {
00488             e1 = nearest (e0, x[i]);
00489             if (e1 > E) E = e1;
00490             if (e1 < W) W = e1;
00491 
00492             x1 = X(e1);
00493             y1 = Y(y[i]);
00494 
00495             if(!edge (x0, y0, x1, y1))
00496                 return NO_MEMORY;
00497 
00498             x0 = x1;
00499             y0 = y1;
00500             e0 = e1;
00501         }
00502         if (pole &&!edge (x0, y0, x0, Y(90.0*pole)))
00503                 return NO_MEMORY;
00504 
00505         shift = 0;        /* shift into window */
00506         while (E+shift > window.east)
00507             shift -= 360.0;
00508         while (E+shift < window.west)
00509             shift += 360.0;
00510         shift1 = X(x[n-1]+shift) - X(x[n-1]);
00511     }
00512     else
00513     {
00514         x0 = X(x[n-1]);
00515         y0 = Y(y[n-1]);
00516 
00517         for (i = 0; i < n; i++)
00518         {
00519             x1 = X(x[i]);
00520             y1 = Y(y[i]);
00521             if(!edge (x0, y0, x1, y1))
00522                 return NO_MEMORY;
00523             x0 = x1;
00524             y0 = y1;
00525         }
00526     }
00527 
00528 /* check if perimeter has odd number of points */
00529     if (np%2)
00530         return OUT_OF_SYNC;
00531 
00532 /* sort the edge points by col(x) and then by row(y) */
00533     qsort (P, np, sizeof(POINT), &edge_order);
00534 
00535 /* plot */
00536     for (i = 1; i < np; i += 2)
00537     {
00538         if (P[i].y != P[i-1].y)
00539             return OUT_OF_SYNC;
00540         row_fill (P[i].y, P[i-1].x+shift1, P[i].x+shift1);
00541     }
00542     if (window.proj == PROJECTION_LL)   /* now do wrap-around, part 2 */
00543     {
00544         shift = 0;
00545         while (W+shift < window.west)
00546             shift += 360.0;
00547         while (W+shift > window.east)
00548             shift -= 360.0;
00549         shift2 = X(x[n-1]+shift) - X(x[n-1]);
00550         if (shift2 != shift1)
00551         {
00552             for (i = 1; i < np; i += 2)
00553             {
00554                 row_fill (P[i].y, P[i-1].x+shift2, P[i].x+shift2);
00555             }
00556         }
00557     }
00558     return OK;
00559 }
00560 
00561 /*
00562  * G_plot_area (xs, ys, rpnts, rings)
00563  *      double **xs;  -- pointer to pointer for X's
00564  *      double **ys;  -- pointer to pointer for Y's
00565  *      int *rpnts;   -- array of ints w/ num points per ring
00566  *      int rings;    -- number of rings
00567  *
00568  * Essentially a copy of G_plot_polygon, with minor mods to
00569  * handle a set of polygons.  return values are the same.
00570  */
00571 
00587 int G_plot_area (double **xs, double **ys, int *rpnts, int rings)
00588 {
00589     int i, j, n;
00590     int pole;
00591     double x0,x1, *x;
00592     double y0,y1, *y;
00593     double shift,E,W=0L;
00594     double e0,e1;
00595     int *shift1 = NULL, shift2;
00596 
00597 /* traverse the perimeter */
00598 
00599     np = 0;
00600     shift1 = (int *) G_calloc (sizeof(int), rings);
00601     
00602     for (j = 0; j < rings; j++) 
00603     {
00604         n = rpnts[j];
00605         
00606         if (n < 3)
00607             return TOO_FEW_EDGES;
00608         
00609         x = xs[j];
00610         y = ys[j];
00611         
00612     /* global wrap-around for lat-lon, part1 */
00613         if (window.proj == PROJECTION_LL)
00614         {
00615             /*
00616             pole = G_pole_in_polygon(x,y,n);
00617             */
00618             pole = 0;
00619 
00620             e0 = x[n-1];
00621             E = W = e0;
00622 
00623             x0 = X(e0);
00624             y0 = Y(y[n-1]);
00625 
00626             if (pole &&!edge (x0, y0, x0, Y(90.0*pole)))
00627                     return NO_MEMORY;
00628 
00629             for (i = 0; i < n; i++)
00630             {
00631                 e1 = nearest (e0, x[i]);
00632                 if (e1 > E) E = e1;
00633                 if (e1 < W) W = e1;
00634 
00635                 x1 = X(e1);
00636                 y1 = Y(y[i]);
00637 
00638                 if(!edge (x0, y0, x1, y1))
00639                     return NO_MEMORY;
00640 
00641                 x0 = x1;
00642                 y0 = y1;
00643                 e0 = e1;
00644             }
00645             if (pole &&!edge (x0, y0, x0, Y(90.0*pole)))
00646                     return NO_MEMORY;
00647 
00648             shift = 0;        /* shift into window */
00649             while (E+shift > window.east)
00650                 shift -= 360.0;
00651             while (E+shift < window.west)
00652                 shift += 360.0;
00653             shift1[j] = X(x[n-1]+shift) - X(x[n-1]);
00654         }
00655         else
00656         {
00657             x0 = X(x[n-1]);
00658             y0 = Y(y[n-1]);
00659 
00660             for (i = 0; i < n; i++)
00661             {
00662                 x1 = X(x[i]);
00663                 y1 = Y(y[i]);
00664                 if(!edge (x0, y0, x1, y1))
00665                     return NO_MEMORY;
00666                 x0 = x1;
00667                 y0 = y1;
00668             }
00669         }
00670     } /* for() */
00671     
00672 /* check if perimeter has odd number of points */
00673     if (np%2)
00674         return OUT_OF_SYNC;
00675 
00676 /* sort the edge points by col(x) and then by row(y) */
00677     qsort (P, np, sizeof(POINT), &edge_order);
00678 
00679 /* plot */
00680     for (j = 0; j < rings; j++)
00681     {
00682         for (i = 1; i < np; i += 2)
00683         {
00684             if (P[i].y != P[i-1].y)
00685                 return OUT_OF_SYNC;
00686             row_fill (P[i].y, P[i-1].x+shift1[j], P[i].x+shift1[j]);
00687         }
00688         if (window.proj == PROJECTION_LL)       /* now do wrap-around, part 2 */
00689         {
00690                 n = rpnts[j];
00691                 x = xs[j];
00692                 y = ys[j];
00693 
00694                 shift = 0;
00695                 while (W+shift < window.west)
00696                     shift += 360.0;
00697                 while (W+shift > window.east)
00698                     shift -= 360.0;
00699                 shift2 = X(x[n-1]+shift) - X(x[n-1]);
00700                 if (shift2 != shift1[j])
00701                 {
00702                     for (i = 1; i < np; i += 2)
00703                     {
00704                         row_fill (P[i].y, P[i-1].x+shift2, P[i].x+shift2);
00705                     }
00706                 }
00707         }
00708     }
00709     G_free (shift1);
00710     return OK;
00711 
00712 }
00713                           
00714 static int edge (double x0,double y0,double x1,double y1)
00715 {
00716     register double m;
00717     double dy, x;
00718     int ystart, ystop;
00719 
00720 
00721 /* tolerance to avoid FPE */
00722     dy = y0 - y1;
00723     if (fabs(dy) < 1e-10)
00724         return 1;
00725 
00726     m = (x0 - x1) / dy;
00727 
00728     if (y0 < y1)
00729     {
00730         ystart = iceil  (y0);
00731         ystop  = ifloor (y1);
00732         if (ystop == y1) ystop--; /* if line stops at row center, don't include point */
00733     }
00734     else
00735     {
00736         ystart = iceil  (y1);
00737         ystop  = ifloor (y0);
00738         if (ystop == y0) ystop--; /* if line stops at row center, don't include point */
00739     }
00740     if (ystart > ystop)
00741         return 1;       /* does not cross center line of row */
00742 
00743     x = m * (ystart - y0) + x0;
00744     while (ystart <= ystop)
00745     {
00746         if(!edge_point (x, ystart++))
00747             return 0;
00748         x += m;
00749     }
00750     return 1;
00751 }
00752 
00753 static int edge_point( double x, register int y)
00754 {
00755 
00756     if (y < ymin || y > ymax)
00757         return 1;
00758     if (np >= npalloc)
00759     {
00760         if (npalloc > 0)
00761         {
00762             npalloc *= 2;
00763             P = (POINT *) G_realloc (P, npalloc * sizeof (POINT));
00764         }
00765         else
00766         {
00767             npalloc = 32;
00768             P = (POINT *) G_malloc (npalloc * sizeof (POINT));
00769         }
00770         if (P == NULL)
00771         {
00772             npalloc = 0;
00773             return 0;
00774         }
00775     }
00776     P[np].x   = x;
00777     P[np++].y = y;
00778     return 1;
00779 }
00780 
00781 static int edge_order(const void *aa, const void *bb)
00782 {
00783     const struct point *a = aa, *b = bb;
00784     if (a->y < b->y) return (-1);
00785     if (a->y > b->y) return (1);
00786 
00787     if (a->x < b->x) return (-1);
00788     if (a->x > b->x) return (1);
00789 
00790     return (0);
00791 }
00792 
00793 static int row_solid_fill(int y,double x1,double x2)
00794 {
00795     int i1,i2;
00796 
00797     i1 = iceil(x1);
00798     i2 = ifloor(x2);
00799     if (i1 <= i2)
00800     {
00801         move (i1, y);
00802         cont (i2, y);
00803     }
00804 
00805     return 0;
00806 }
00807 
00808 static int row_dotted_fill(int y,double x1,double x2)
00809 {
00810     int i1,i2,i;
00811 
00812     if(y != iceil(y / dotted_fill_gap) * dotted_fill_gap)
00813             return 0;
00814 
00815     i1 = iceil(x1 / dotted_fill_gap) * dotted_fill_gap;
00816     i2 = ifloor(x2);
00817     if (i1 <= i2)
00818     {
00819         for(i = i1; i <= i2; i += dotted_fill_gap)
00820         {
00821             move (i, y);
00822             cont (i, y);
00823         }
00824     }
00825 
00826     return 0;
00827 }
00828 
00829 static int ifloor(double x)
00830 {
00831     int i;
00832     i = (int) x;
00833     if (i > x) i--;
00834     return i;
00835 }
00836 
00837 static int iceil(double x)
00838 {
00839     int i;
00840 
00841     i = (int) x;
00842     if (i < x) i++;
00843     return i;
00844 }
00845 
00846 /*
00847  * G_plot_fx(e1,e2)
00848  *
00849  * plot f(x) from x=e1 to x=e2
00850  */
00851 
00852 
00864 int G_plot_fx (
00865     double (*f)(),
00866     double east1,double east2)
00867 {
00868     double east,north,north1;
00869     double incr;
00870 
00871 
00872     incr = fabs(1.0 / xconv) ;
00873 
00874     east  = east1;
00875     north = f(east1);
00876 
00877     if (east1 > east2)
00878     {
00879         while ((east1 -= incr) > east2)
00880         {
00881             north1 = f(east1);
00882             G_plot_line (east, north, east1, north1);
00883             north = north1;
00884             east  = east1;
00885         }
00886     }
00887     else
00888     {
00889         while ((east1 += incr) < east2)
00890         {
00891             north1 = f(east1);
00892             G_plot_line (east, north, east1, north1);
00893             north = north1;
00894             east  = east1;
00895         }
00896     }
00897     G_plot_line (east, north, east2, f(east2));
00898 
00899     return 0;
00900 }

Generated on Wed Dec 19 14:59:06 2007 for GRASS by  doxygen 1.5.4