00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
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
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
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
00211
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
00241
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
00253
00254
00255
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)
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
00369
00370
00371
00372
00373
00374
00375
00376
00377
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
00426
00427 np = 0;
00428 shift1 = 0;
00429
00430
00431 if (window.proj == PROJECTION_LL) {
00432
00433
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;
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
00488 if (np % 2)
00489 return OUT_OF_SYNC;
00490
00491
00492 qsort(P, np, sizeof(POINT), &edge_order);
00493
00494
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) {
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
00518
00519
00520
00521
00522
00523
00524
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
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
00567 if (window.proj == PROJECTION_LL) {
00568
00569
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;
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 }
00623
00624
00625 if (np % 2)
00626 return OUT_OF_SYNC;
00627
00628
00629 qsort(P, np, sizeof(POINT), &edge_order);
00630
00631
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) {
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
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--;
00680 }
00681 else {
00682 ystart = iceil(y1);
00683 ystop = ifloor(y0);
00684 if (ystop == y0)
00685 ystop--;
00686 }
00687 if (ystart > ystop)
00688 return 1;
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
00795
00796
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 }