GRASS Programmer's Manual 6.4.1(2011)
|
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 }