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 "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_fill(int,double,double);
00037 static int ifloor(double);
00038 static int iceil(double);
00039 static int (*move)() = NULL;
00040 static int (*cont)() = NULL;
00041 
00042 static double fabs(double x)
00043 {
00044     return x>0?x:-x;
00045 }
00046 
00061 extern double G_adjust_easting();
00062 
00063 /*
00064  * G_setup_plot (t, b, l, r, Move, Cont)
00065  *     double t, b, l, r;
00066  *     int (*Move)(), (*Cont)();
00067  *
00068  * initialize the plotting capability.
00069  *    t,b,l,r:   top, bottom, left, right of the output x,y coordinate space.
00070  *    Move,Cont: subroutines that will draw lines in x,y space.
00071  *       Move(x,y)   move to x,y (no draw)
00072  *       Cont(x,y)   draw from previous position to x,y
00073  * Notes:
00074  *   Cont() is responsible for clipping.
00075  *   The t,b,l,r are only used to compute coordinate transformations.
00076  *   The input space is assumed to be the current GRASS window.
00077  */
00078 
00100 int G_setup_plot (
00101     double t,double b,double l,double r,
00102     int (*Move)(),
00103     int (*Cont)())
00104 {
00105     G_get_set_window (&window);
00106 
00107     left = l;
00108     right = r;
00109     top = t;
00110     bottom = b;
00111 
00112     xconv = (right-left)/(window.east-window.west);
00113     yconv = (bottom-top)/(window.north-window.south);
00114 
00115     if (top < bottom)
00116     {
00117         ymin = iceil(top);
00118         ymax = ifloor(bottom);
00119     }
00120     else
00121     {
00122         ymin = iceil(bottom);
00123         ymax = ifloor(top);
00124     }
00125 
00126     move = Move;
00127     cont = Cont;
00128 
00129     return 0;
00130 }
00131 
00132 #define X(e) (left + xconv * ((e) - window.west))
00133 #define Y(n) (top + yconv * (window.north - (n)))
00134 
00135 #define EAST(x) (window.west + ((x)-left)/xconv)
00136 #define NORTH(y) (window.north - ((y)-top)/yconv)
00137 
00138 
00152 int G_plot_where_xy (east, north, x, y)
00153     double east, north;
00154     int *x, *y;
00155 {
00156     *x = ifloor(X(G_adjust_easting(east,&window))+0.5);
00157     *y = ifloor(Y(north)+0.5);
00158 
00159     return 0;
00160 }
00161 
00162 
00176 int G_plot_where_en (x, y, east, north)
00177     double *east, *north;
00178 {
00179     *east = G_adjust_easting(EAST(x),&window);
00180     *north = NORTH(y);
00181 
00182     return 0;
00183 }
00184 
00185 int G_plot_point (east, north)
00186     double east, north;
00187 {
00188     int x,y;
00189 
00190     G_plot_where_xy(east,north,&x,&y);
00191     move (x,y);
00192     cont (x,y);
00193 
00194     return 0;
00195 }
00196 
00197 /*
00198  * Line in map coordinates is plotted in output x,y coordinates
00199  * This routine handles global wrap-around for lat-long databses.
00200  *
00201  */
00202 
00218 int G_plot_line (east1, north1, east2, north2)
00219     double east1, north1, east2, north2;
00220 {
00221     int fastline();
00222     plot_line (east1, north1, east2, north2, fastline);
00223 
00224     return 0;
00225 }
00226 
00227 int G_plot_line2 (east1, north1, east2, north2)
00228     double east1, north1, east2, north2;
00229 {
00230     int slowline();
00231     plot_line (east1, north1, east2, north2, slowline);
00232 
00233     return 0;
00234 }
00235 
00236 /* fastline converts double rows/cols to ints then plots
00237  * this is ok for graphics, but not the best for vector to raster
00238  */
00239 
00240 static int fastline(double x1,double y1,double x2,double y2)
00241 {
00242     move (ifloor(x1+0.5),ifloor(y1+0.5));
00243     cont (ifloor(x2+0.5),ifloor(y2+0.5));
00244 
00245     return 0;
00246 }
00247 
00248 /* NOTE (shapiro): 
00249  *   I think the adding of 0.5 in slowline is not correct
00250  *   the output window (left, right, top, bottom) should already
00251  *   be adjusted for this: left=-0.5; right = window.cols-0.5;
00252  */
00253 
00254 static int slowline(double x1,double y1,double x2,double y2)
00255 {
00256     double dx, dy;
00257     double m,b;
00258     int xstart, xstop, ystart, ystop;
00259 
00260     dx = x2-x1;
00261     dy = y2-y1;
00262 
00263     if (fabs(dx) > fabs(dy))
00264     {
00265         m = dy/dx;
00266         b = y1 - m * x1;
00267 
00268         if (x1 > x2)
00269         {
00270             xstart = iceil (x2-0.5);
00271             xstop  = ifloor (x1+0.5);
00272         }
00273         else
00274         {
00275             xstart = iceil (x1-0.5);
00276             xstop = ifloor (x2+0.5);
00277         }
00278         if (xstart <= xstop)
00279         {
00280             ystart = ifloor(m * xstart + b + 0.5);
00281             move (xstart, ystart);
00282             while (xstart <= xstop)
00283             {
00284                 cont (xstart++, ystart);
00285                 ystart = ifloor(m * xstart + b + 0.5);
00286             }
00287         }
00288     }
00289     else
00290     {
00291         if(dx==dy) /* they both might be 0 */
00292            m = 1.;
00293         else 
00294            m = dx/dy;
00295         b = x1 - m * y1;
00296 
00297         if (y1 > y2)
00298         {
00299             ystart = iceil (y2-0.5);
00300             ystop  = ifloor (y1+0.5);
00301         }
00302         else
00303         {
00304             ystart = iceil (y1-0.5);
00305             ystop = ifloor (y2+0.5);
00306         }
00307         if (ystart <= ystop)
00308         {
00309             xstart = ifloor(m * ystart + b + 0.5);
00310             move (xstart, ystart);
00311             while (ystart <= ystop)
00312             {
00313                 cont (xstart, ystart++);
00314                 xstart = ifloor(m * ystart + b + 0.5);
00315             }
00316         }
00317     }
00318 
00319     return 0;
00320 }
00321 
00322 static int plot_line(double east1,double north1,double east2,double north2,
00323     int (*line)())
00324 {
00325     double x1,x2,y1,y2;
00326 
00327     y1 = Y(north1);
00328     y2 = Y(north2);
00329 
00330     if (window.proj == PROJECTION_LL)
00331     {
00332         if (east1 > east2)
00333             while ((east1-east2) > 180)
00334                 east2 += 360;
00335         else if (east2 > east1)
00336             while ((east2-east1) > 180)
00337                 east1 += 360;
00338         while (east1 > window.east)
00339         {
00340             east1 -= 360.0;
00341             east2 -= 360.0;
00342         }
00343         while (east1 < window.west)
00344         {
00345             east1 += 360.0;
00346             east2 += 360.0;
00347         }
00348         x1 = X(east1);
00349         x2 = X(east2);
00350 
00351         line (x1,y1,x2,y2);
00352 
00353         if (east2 > window.east || east2 < window.west)
00354         {
00355             while (east2 > window.east)
00356             {
00357                 east1 -= 360.0;
00358                 east2 -= 360.0;
00359             }
00360             while (east2 < window.west)
00361             {
00362                 east1 += 360.0;
00363                 east2 += 360.0;
00364             }
00365             x1 = X(east1);
00366             x2 = X(east2);
00367             line (x1,y1,x2,y2);
00368         }
00369     }
00370     else
00371     {
00372         x1 = X(east1);
00373         x2 = X(east2);
00374         line (x1,y1,x2,y2);
00375     }
00376 
00377     return 0;
00378 }
00379 
00380 /*
00381  * G_plot_polygon (x, y, n)
00382  * 
00383  *    double *x       x coordinates of vertices
00384  *    double *y       y coordinates of vertices
00385  *    int n           number of verticies
00386  *
00387  * polygon fill from map coordinate space to plot x,y space.
00388  * for lat-lon, handles global wrap-around as well as polar polygons.
00389  *
00390  * returns 0 ok, 2 n<3, -1 weird internal error, 1 no memory
00391  */
00392 
00393 static POINT *P;
00394 static int np;
00395 static int npalloc = 0;
00396 
00397 #define OK 0
00398 #define TOO_FEW_EDGES 2
00399 #define NO_MEMORY 1
00400 #define OUT_OF_SYNC -1
00401 
00402 static double nearest(double e0,double e1)
00403 {
00404     while (e0 - e1 > 180)
00405         e1 += 360.0;
00406     while (e1 - e0 > 180)
00407         e1 -= 360.0;
00408     
00409     return e1;
00410 }
00411 
00412 
00425 int G_plot_polygon (
00426     double *x,double *y,
00427     int n)
00428 {
00429     int i;
00430     int pole;
00431     double x0,x1;
00432     double y0,y1;
00433     double shift,E,W=0L;
00434     double e0,e1;
00435     int shift1, shift2;
00436 
00437     if (n < 3)
00438         return TOO_FEW_EDGES;
00439 
00440 /* traverse the perimeter */
00441 
00442     np = 0;
00443     shift1 = 0;
00444 
00445 /* global wrap-around for lat-lon, part1 */
00446     if (window.proj == PROJECTION_LL)
00447     {
00448         /*
00449         pole = G_pole_in_polygon(x,y,n);
00450         */
00451         pole = 0;
00452 
00453         e0 = x[n-1];
00454         E = W = e0;
00455 
00456         x0 = X(e0);
00457         y0 = Y(y[n-1]);
00458 
00459         if (pole &&!edge (x0, y0, x0, Y(90.0*pole)))
00460                 return NO_MEMORY;
00461 
00462         for (i = 0; i < n; i++)
00463         {
00464             e1 = nearest (e0, x[i]);
00465             if (e1 > E) E = e1;
00466             if (e1 < W) W = e1;
00467 
00468             x1 = X(e1);
00469             y1 = Y(y[i]);
00470 
00471             if(!edge (x0, y0, x1, y1))
00472                 return NO_MEMORY;
00473 
00474             x0 = x1;
00475             y0 = y1;
00476             e0 = e1;
00477         }
00478         if (pole &&!edge (x0, y0, x0, Y(90.0*pole)))
00479                 return NO_MEMORY;
00480 
00481         shift = 0;        /* shift into window */
00482         while (E+shift > window.east)
00483             shift -= 360.0;
00484         while (E+shift < window.west)
00485             shift += 360.0;
00486         shift1 = X(x[n-1]+shift) - X(x[n-1]);
00487     }
00488     else
00489     {
00490         x0 = X(x[n-1]);
00491         y0 = Y(y[n-1]);
00492 
00493         for (i = 0; i < n; i++)
00494         {
00495             x1 = X(x[i]);
00496             y1 = Y(y[i]);
00497             if(!edge (x0, y0, x1, y1))
00498                 return NO_MEMORY;
00499             x0 = x1;
00500             y0 = y1;
00501         }
00502     }
00503 
00504 /* check if perimeter has odd number of points */
00505     if (np%2)
00506         return OUT_OF_SYNC;
00507 
00508 /* sort the edge points by col(x) and then by row(y) */
00509     qsort (P, np, sizeof(POINT), &edge_order);
00510 
00511 /* plot */
00512     for (i = 1; i < np; i += 2)
00513     {
00514         if (P[i].y != P[i-1].y)
00515             return OUT_OF_SYNC;
00516         row_fill (P[i].y, P[i-1].x+shift1, P[i].x+shift1);
00517     }
00518     if (window.proj == PROJECTION_LL)   /* now do wrap-around, part 2 */
00519     {
00520         shift = 0;
00521         while (W+shift < window.west)
00522             shift += 360.0;
00523         while (W+shift > window.east)
00524             shift -= 360.0;
00525         shift2 = X(x[n-1]+shift) - X(x[n-1]);
00526         if (shift2 != shift1)
00527         {
00528             for (i = 1; i < np; i += 2)
00529             {
00530                 row_fill (P[i].y, P[i-1].x+shift2, P[i].x+shift2);
00531             }
00532         }
00533     }
00534     return OK;
00535 }
00536 
00537 /*
00538  * G_plot_area (xs, ys, rpnts, rings)
00539  *      double **xs;  -- pointer to pointer for X's
00540  *      double **ys;  -- pointer to pointer for Y's
00541  *      int *rpnts;   -- array of ints w/ num points per ring
00542  *      int rings;    -- number of rings
00543  *
00544  * Essentially a copy of G_plot_polygon, with minor mods to
00545  * handle a set of polygons.  return values are the same.
00546  */
00547 
00563 int G_plot_area (double **xs, double **ys, int *rpnts, int rings)
00564 {
00565     int i, j, n;
00566     int pole;
00567     double x0,x1, *x;
00568     double y0,y1, *y;
00569     double shift,E,W=0L;
00570     double e0,e1;
00571     int *shift1 = NULL, shift2;
00572 
00573 /* traverse the perimeter */
00574 
00575     np = 0;
00576     shift1 = (int *) G_calloc (sizeof(int), rings);
00577     
00578     for (j = 0; j < rings; j++) 
00579     {
00580         n = rpnts[j];
00581         
00582         if (n < 3)
00583             return TOO_FEW_EDGES;
00584         
00585         x = xs[j];
00586         y = ys[j];
00587         
00588     /* global wrap-around for lat-lon, part1 */
00589         if (window.proj == PROJECTION_LL)
00590         {
00591             /*
00592             pole = G_pole_in_polygon(x,y,n);
00593             */
00594             pole = 0;
00595 
00596             e0 = x[n-1];
00597             E = W = e0;
00598 
00599             x0 = X(e0);
00600             y0 = Y(y[n-1]);
00601 
00602             if (pole &&!edge (x0, y0, x0, Y(90.0*pole)))
00603                     return NO_MEMORY;
00604 
00605             for (i = 0; i < n; i++)
00606             {
00607                 e1 = nearest (e0, x[i]);
00608                 if (e1 > E) E = e1;
00609                 if (e1 < W) W = e1;
00610 
00611                 x1 = X(e1);
00612                 y1 = Y(y[i]);
00613 
00614                 if(!edge (x0, y0, x1, y1))
00615                     return NO_MEMORY;
00616 
00617                 x0 = x1;
00618                 y0 = y1;
00619                 e0 = e1;
00620             }
00621             if (pole &&!edge (x0, y0, x0, Y(90.0*pole)))
00622                     return NO_MEMORY;
00623 
00624             shift = 0;        /* shift into window */
00625             while (E+shift > window.east)
00626                 shift -= 360.0;
00627             while (E+shift < window.west)
00628                 shift += 360.0;
00629             shift1[j] = X(x[n-1]+shift) - X(x[n-1]);
00630         }
00631         else
00632         {
00633             x0 = X(x[n-1]);
00634             y0 = Y(y[n-1]);
00635 
00636             for (i = 0; i < n; i++)
00637             {
00638                 x1 = X(x[i]);
00639                 y1 = Y(y[i]);
00640                 if(!edge (x0, y0, x1, y1))
00641                     return NO_MEMORY;
00642                 x0 = x1;
00643                 y0 = y1;
00644             }
00645         }
00646     } /* for() */
00647     
00648 /* check if perimeter has odd number of points */
00649     if (np%2)
00650         return OUT_OF_SYNC;
00651 
00652 /* sort the edge points by col(x) and then by row(y) */
00653     qsort (P, np, sizeof(POINT), &edge_order);
00654 
00655 /* plot */
00656     for (j = 0; j < rings; j++)
00657     {
00658         for (i = 1; i < np; i += 2)
00659         {
00660             if (P[i].y != P[i-1].y)
00661                 return OUT_OF_SYNC;
00662             row_fill (P[i].y, P[i-1].x+shift1[j], P[i].x+shift1[j]);
00663         }
00664         if (window.proj == PROJECTION_LL)       /* now do wrap-around, part 2 */
00665         {
00666                 n = rpnts[j];
00667                 x = xs[j];
00668                 y = ys[j];
00669 
00670                 shift = 0;
00671                 while (W+shift < window.west)
00672                     shift += 360.0;
00673                 while (W+shift > window.east)
00674                     shift -= 360.0;
00675                 shift2 = X(x[n-1]+shift) - X(x[n-1]);
00676                 if (shift2 != shift1[j])
00677                 {
00678                     for (i = 1; i < np; i += 2)
00679                     {
00680                         row_fill (P[i].y, P[i-1].x+shift2, P[i].x+shift2);
00681                     }
00682                 }
00683         }
00684     }
00685     free (shift1);
00686     return OK;
00687 
00688 }
00689                           
00690 static int edge (double x0,double y0,double x1,double y1)
00691 {
00692     register double m;
00693     double dy, x;
00694     int ystart, ystop;
00695 
00696 
00697 /* tolerance to avoid FPE */
00698     dy = y0 - y1;
00699     if (fabs(dy) < 1e-10)
00700         return 1;
00701 
00702     m = (x0 - x1) / dy;
00703 
00704     if (y0 < y1)
00705     {
00706         ystart = iceil  (y0);
00707         ystop  = ifloor (y1);
00708         if (ystop == y1) ystop--; /* if line stops at row center, don't include point */
00709     }
00710     else
00711     {
00712         ystart = iceil  (y1);
00713         ystop  = ifloor (y0);
00714         if (ystop == y0) ystop--; /* if line stops at row center, don't include point */
00715     }
00716     if (ystart > ystop)
00717         return 1;       /* does not cross center line of row */
00718 
00719     x = m * (ystart - y0) + x0;
00720     while (ystart <= ystop)
00721     {
00722         if(!edge_point (x, ystart++))
00723             return 0;
00724         x += m;
00725     }
00726     return 1;
00727 }
00728 
00729 static int edge_point( double x, register int y)
00730 {
00731 
00732     if (y < ymin || y > ymax)
00733         return 1;
00734     if (np >= npalloc)
00735     {
00736         if (npalloc > 0)
00737         {
00738             npalloc *= 2;
00739             P = (POINT *) realloc (P, npalloc * sizeof (POINT));
00740         }
00741         else
00742         {
00743             npalloc = 32;
00744             P = (POINT *) malloc (npalloc * sizeof (POINT));
00745         }
00746         if (P == NULL)
00747         {
00748             npalloc = 0;
00749             return 0;
00750         }
00751     }
00752     P[np].x   = x;
00753     P[np++].y = y;
00754     return 1;
00755 }
00756 
00757 static int edge_order(const void *aa, const void *bb)
00758 {
00759     const struct point *a = aa, *b = bb;
00760     if (a->y < b->y) return (-1);
00761     if (a->y > b->y) return (1);
00762 
00763     if (a->x < b->x) return (-1);
00764     if (a->x > b->x) return (1);
00765 
00766     return (0);
00767 }
00768 
00769 static int row_fill(int y,double x1,double x2)
00770 {
00771     int i1,i2;
00772 
00773     i1 = iceil(x1);
00774     i2 = ifloor(x2);
00775     if (i1 <= i2)
00776     {
00777         move (i1, y);
00778         cont (i2, y);
00779     }
00780 
00781     return 0;
00782 }
00783 
00784 static int ifloor(double x)
00785 {
00786     int i;
00787     i = (int) x;
00788     if (i > x) i--;
00789     return i;
00790 }
00791 
00792 static int iceil(double x)
00793 {
00794     int i;
00795 
00796     i = (int) x;
00797     if (i < x) i++;
00798     return i;
00799 }
00800 
00801 /*
00802  * G_plot_fx(e1,e2)
00803  *
00804  * plot f(x) from x=e1 to x=e2
00805  */
00806 
00807 
00819 int G_plot_fx (
00820     double (*f)(),
00821     double east1,double east2)
00822 {
00823     double east,north,north1;
00824     double incr;
00825 
00826 
00827     incr = fabs(1.0 / xconv) ;
00828 
00829     east  = east1;
00830     north = f(east1);
00831 
00832     if (east1 > east2)
00833     {
00834         while ((east1 -= incr) > east2)
00835         {
00836             north1 = f(east1);
00837             G_plot_line (east, north, east1, north1);
00838             north = north1;
00839             east  = east1;
00840         }
00841     }
00842     else
00843     {
00844         while ((east1 += incr) < east2)
00845         {
00846             north1 = f(east1);
00847             G_plot_line (east, north, east1, north1);
00848             north = north1;
00849             east  = east1;
00850         }
00851     }
00852     G_plot_line (east, north, east2, f(east2));
00853 
00854     return 0;
00855 }

Generated on Mon Jan 1 19:49:25 2007 for GRASS by  doxygen 1.5.1