00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
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
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
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
00199
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
00237
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
00249
00250
00251
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)
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
00382
00383
00384
00385
00386
00387
00388
00389
00390
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
00441
00442 np = 0;
00443 shift1 = 0;
00444
00445
00446 if (window.proj == PROJECTION_LL)
00447 {
00448
00449
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;
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
00505 if (np%2)
00506 return OUT_OF_SYNC;
00507
00508
00509 qsort (P, np, sizeof(POINT), &edge_order);
00510
00511
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)
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
00539
00540
00541
00542
00543
00544
00545
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
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
00589 if (window.proj == PROJECTION_LL)
00590 {
00591
00592
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;
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 }
00647
00648
00649 if (np%2)
00650 return OUT_OF_SYNC;
00651
00652
00653 qsort (P, np, sizeof(POINT), &edge_order);
00654
00655
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)
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
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--;
00709 }
00710 else
00711 {
00712 ystart = iceil (y1);
00713 ystop = ifloor (y0);
00714 if (ystop == y0) ystop--;
00715 }
00716 if (ystart > ystop)
00717 return 1;
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
00803
00804
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 }