• Main Page
  • Related Pages
  • Data Structures
  • Files
  • File List
  • Globals

plot.c

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

Generated on Wed Oct 13 2010 12:09:30 for GRASS Programmer's Manual by  doxygen 1.7.1