00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
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
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
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
00223
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
00261
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
00273
00274
00275
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)
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
00406
00407
00408
00409
00410
00411
00412
00413
00414
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
00465
00466 np = 0;
00467 shift1 = 0;
00468
00469
00470 if (window.proj == PROJECTION_LL)
00471 {
00472
00473
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;
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
00529 if (np%2)
00530 return OUT_OF_SYNC;
00531
00532
00533 qsort (P, np, sizeof(POINT), &edge_order);
00534
00535
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)
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
00563
00564
00565
00566
00567
00568
00569
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
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
00613 if (window.proj == PROJECTION_LL)
00614 {
00615
00616
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;
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 }
00671
00672
00673 if (np%2)
00674 return OUT_OF_SYNC;
00675
00676
00677 qsort (P, np, sizeof(POINT), &edge_order);
00678
00679
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)
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
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--;
00733 }
00734 else
00735 {
00736 ystart = iceil (y1);
00737 ystop = ifloor (y0);
00738 if (ystop == y0) ystop--;
00739 }
00740 if (ystart > ystop)
00741 return 1;
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
00848
00849
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 }