GRASS Programmer's Manual  6.4.2(2012)
draw2.c
Go to the documentation of this file.
00001 
00002 /*******************************************************************
00003  * Line drawing in the current window.
00004  *
00005  * Clip window:
00006  *   D_set_clip_window (top, bottom ,left, right)
00007  *      establish clipping region for subseqent line drawing.
00008  *   D_set_clip_window_to_map_window ()
00009  *     set clipping to pixels corresponding to the current map region
00010  *     (default)
00011  *   D_set_clip_window_to_screen_window ()
00012  *     set clipping to full extent of the window (ie disables clipping on screen)
00013  *
00014  * Moves.
00015  *   D_move_abs(x,y)   move to x,y.
00016  *   D_move_rel(x,y)   move to +x,+y.
00017  *      Set current position. Position is not clipped.
00018  *
00019  * Draw line 
00020  *   D_cont_abs(x,y)   draw to x,y.
00021  *   D_cont_rel(x,y)   draw to +x,+y.
00022  *      Line draw from current position. New postion is not clipped.
00023  *      The lines drawn are clipped however.
00024  *      Return values indicate the nature of the clipping:
00025  *        0 no clipping
00026  *        1 part of the line is drawn
00027  *       -1 none of the line is drawn
00028  *   
00029  *
00030  */
00031 
00032 #include <math.h>
00033 #include <string.h>
00034 
00035 #include <grass/gis.h>
00036 #include <grass/raster.h>
00037 #include <grass/display.h>
00038 
00039 struct rectangle
00040 {
00041     double left;
00042     double rite;
00043     double bot;
00044     double top;
00045 };
00046 
00047 struct vector
00048 {
00049     double x, y;
00050 };
00051 
00052 struct plane
00053 {
00054     double x, y, k;
00055 };
00056 
00057 static struct vector cur;
00058 
00059 static struct rectangle clip;
00060 
00061 static struct plane pl_left = { -1, 0, 0 };
00062 static struct plane pl_rite = { 1, 0, 0 };
00063 static struct plane pl_bot = { 0, -1, 0 };
00064 static struct plane pl_top = { 0, 1, 0 };
00065 
00066 static int window_set;
00067 
00068 #define min(x,y) ((x) < (y) ? (x) : (y))
00069 #define max(x,y) ((x) > (y) ? (x) : (y))
00070 
00071 #define round(x) ((int) floor(0.5 + (x)))
00072 
00073 static int *xi, *yi;
00074 static int nalloc_i;
00075 
00076 static double *xf, *yf;
00077 static int nalloc_f;
00078 
00079 static void alloc_int(int n)
00080 {
00081 
00082     if (nalloc_i >= n)
00083         return;
00084 
00085     nalloc_i = n;
00086     xi = G_realloc(xi, nalloc_i * sizeof(int));
00087     yi = G_realloc(yi, nalloc_i * sizeof(int));
00088 }
00089 
00090 static void alloc_float(int n)
00091 {
00092 
00093     if (nalloc_f >= n)
00094         return;
00095 
00096     nalloc_f = n + 10;
00097     xf = G_realloc(xf, nalloc_f * sizeof(double));
00098     yf = G_realloc(yf, nalloc_f * sizeof(double));
00099 }
00100 
00101 static void dealloc_float(const double **x, const double **y, int release)
00102 {
00103     if (release) {
00104         G_free(*(double **)x);
00105         G_free(*(double **)y);
00106     }
00107 
00108     *x = xf;
00109     *y = yf;
00110 
00111     nalloc_f = 0;
00112 
00113     xf = NULL;
00114     yf = NULL;
00115 }
00116 
00117 static int do_filter(int *x, int *y, int n)
00118 {
00119     int i, j;
00120 
00121     for (i = 0, j = 1; j < n; j++) {
00122         if (x[j] == x[i] && y[j] == y[i])
00123             continue;
00124         i++;
00125         if (i == j)
00126             continue;
00127         x[i] = x[j];
00128         y[i] = y[j];
00129     }
00130 
00131     return i + 1;
00132 }
00133 
00134 static void do_round(const double *x, const double *y, int n)
00135 {
00136     int i;
00137 
00138     alloc_int(n);
00139 
00140     for (i = 0; i < n; i++) {
00141         xi[i] = round(D_u_to_d_col(x[i]));
00142         yi[i] = round(D_u_to_d_row(y[i]));
00143     }
00144 }
00145 
00146 static void do_floor(const double *x, const double *y, int n)
00147 {
00148     int i;
00149 
00150     alloc_int(n);
00151 
00152     for (i = 0; i < n; i++) {
00153         xi[i] = floor(D_u_to_d_col(x[i]));
00154         yi[i] = floor(D_u_to_d_row(y[i]));
00155     }
00156 }
00157 
00158 static double dist_plane(double x, double y, const struct plane *p)
00159 {
00160     return x * p->x + y * p->y + p->k;
00161 }
00162 
00163 static double interpolate(double a, double b, double ka, double kb)
00164 {
00165     return (a * kb - b * ka) / (kb - ka);
00166 }
00167 
00168 static int clip_plane(struct vector *a, struct vector *b,
00169                       const struct plane *p, int *clipped)
00170 {
00171     double ka = dist_plane(a->x, a->y, p);
00172     double kb = dist_plane(b->x, b->y, p);
00173     double kab;
00174 
00175     /* both outside */
00176     if (ka > 0 && kb > 0)
00177         return 1;
00178 
00179     /* both inside */
00180     if (ka <= 0 && kb <= 0)
00181         return 0;
00182 
00183     *clipped = 1;
00184 
00185     /* a outside - swap a and b */
00186     if (ka >= 0) {
00187         struct vector *t;
00188         double kt;
00189 
00190         t = a;
00191         a = b;
00192         b = t;
00193         kt = ka;
00194         ka = kb;
00195         kb = kt;
00196     }
00197 
00198     kab = kb - ka;
00199 
00200     b->x = interpolate(a->x, b->x, ka, kb);
00201     b->y = interpolate(a->y, b->y, ka, kb);
00202 
00203     return 0;
00204 }
00205 
00206 static int do_clip(struct vector *a, struct vector *b)
00207 {
00208     int clipped = 0;
00209 
00210     if (a->x < clip.left && b->x < clip.left)
00211         return -1;
00212     if (a->x > clip.rite && b->x > clip.rite)
00213         return -1;
00214     if (a->y < clip.bot && b->y < clip.bot)
00215         return -1;
00216     if (a->y > clip.top && b->y > clip.top)
00217         return -1;
00218 
00219     if (clip_plane(a, b, &pl_left, &clipped))
00220         return -1;
00221     if (clip_plane(a, b, &pl_rite, &clipped))
00222         return -1;
00223     if (clip_plane(a, b, &pl_bot, &clipped))
00224         return -1;
00225     if (clip_plane(a, b, &pl_top, &clipped))
00226         return -1;
00227 
00228     return clipped;
00229 }
00230 
00231 static int shift_count(double dx)
00232 {
00233     return (int)floor(dx / 360);
00234 }
00235 
00236 static double shift_angle(double dx)
00237 {
00238     return shift_count(dx) * 360;
00239 }
00240 
00241 static double coerce(double x)
00242 {
00243     x += 180;
00244     x -= shift_angle(x);
00245     x -= 180;
00246     return x;
00247 }
00248 
00249 static int euclidify(double *x, const double *y, int n, int no_pole)
00250 {
00251     double ux0 = clip.left;
00252     double ux1 = clip.rite;
00253     double x0, x1;
00254     int lo, hi, count;
00255     int i;
00256 
00257     x0 = x1 = x[0];
00258 
00259     for (i = 1; i < n; i++) {
00260         if (fabs(y[i]) < 89.9)
00261             x[i] = x[i - 1] + coerce(x[i] - x[i - 1]);
00262 
00263         x0 = min(x0, x[i]);
00264         x1 = max(x1, x[i]);
00265     }
00266 
00267     if (no_pole && fabs(x[n - 1] - x[0]) > 180)
00268         return 0;
00269 
00270     lo = -shift_count(ux1 - x0);
00271     hi = shift_count(x1 - ux0);
00272     count = hi - lo + 1;
00273 
00274     for (i = 0; i < n; i++)
00275         x[i] -= lo * 360;
00276 
00277     return count;
00278 }
00279 
00280 static void do_ll_wrap(const double *x, const double *y, int n,
00281                        void (*func) (const double *, const double *, int))
00282 {
00283     double *xx = G_malloc(n * sizeof(double));
00284     int count, i;
00285 
00286     memcpy(xx, x, n * sizeof(double));
00287     count = euclidify(xx, y, n, 0);
00288 
00289     for (i = 0; i < count; i++) {
00290         int j;
00291 
00292         (*func) (xx, y, n);
00293 
00294         for (j = 0; j < n; j++)
00295             xx[j] -= 360;
00296     }
00297 
00298     G_free(xx);
00299 }
00300 
00313 void D_set_clip(double t, double b, double l, double r)
00314 {
00315     clip.left = min(l, r);
00316     clip.rite = max(l, r);
00317     clip.bot = min(b, t);
00318     clip.top = max(b, t);
00319 
00320     pl_left.k = clip.left;
00321     pl_rite.k = -clip.rite;
00322     pl_bot.k = clip.bot;
00323     pl_top.k = -clip.top;
00324 
00325     window_set = 1;
00326 }
00327 
00337 void D_clip_to_map(void)
00338 {
00339     D_set_clip(D_get_u_north(),
00340                D_get_u_south(), D_get_u_west(), D_get_u_east());
00341 }
00342 
00353 void D_move_clip(double x, double y)
00354 {
00355     cur.x = x;
00356     cur.y = y;
00357 }
00358 
00374 static int line_clip(double x1, double y1, double x2, double y2)
00375 {
00376     struct vector a, b;
00377     int clipped;
00378 
00379     a.x = x1;
00380     a.y = y1;
00381 
00382     b.x = x2;
00383     b.y = y2;
00384 
00385     clipped = do_clip(&a, &b);
00386 
00387     if (clipped >= 0) {
00388         int x1 = round(D_u_to_d_col(a.x));
00389         int y1 = round(D_u_to_d_row(a.y));
00390         int x2 = round(D_u_to_d_col(b.x));
00391         int y2 = round(D_u_to_d_row(b.y));
00392 
00393         R_move_abs(x1, y1);
00394         R_cont_abs(x2, y2);
00395     }
00396 
00397     return clipped;
00398 }
00399 
00400 static int line_clip_ll(double ax, double ay, double bx, double by)
00401 {
00402     double ux0 = clip.left;
00403     double ux1 = clip.rite;
00404     double x0, x1;
00405     int lo, hi, i;
00406     int ret;
00407 
00408     bx = ax + coerce(bx - ax);
00409 
00410     x0 = min(ax, bx);
00411     x1 = max(ax, bx);
00412 
00413     lo = -shift_count(ux1 - x0);
00414     hi = shift_count(x1 - ux0);
00415 
00416     ret = 0;
00417 
00418     for (i = lo; i <= hi; i++)
00419         ret |= line_clip(ax + i * 360, ay, bx + i * 360, by);
00420 
00421     return ret;
00422 }
00423 
00424 int D_cont_clip(double x, double y)
00425 {
00426     int ret;
00427 
00428     if (!window_set)
00429         D_clip_to_map();
00430 
00431     if (D_is_lat_lon())
00432         ret = line_clip_ll(cur.x, cur.y, x, y);
00433     else
00434         ret = line_clip(cur.x, cur.y, x, y);
00435 
00436     cur.x = x;
00437     cur.y = y;
00438 
00439     return ret;
00440 }
00441 
00442 void D_polydots_clip(const double *x, const double *y, int n)
00443 {
00444     double ux0 = clip.left;
00445     int i, j;
00446 
00447     if (!window_set)
00448         D_clip_to_map();
00449 
00450     alloc_float(n);
00451 
00452     for (i = j = 0; i < n; i++) {
00453         double xx = x[i];
00454         double yy = y[i];
00455 
00456         if (D_is_lat_lon())
00457             xx -= shift_angle(x[i] - ux0);
00458 
00459         if (xx < clip.left || xx > clip.rite)
00460             continue;
00461         if (yy < clip.bot || yy > clip.top)
00462             continue;
00463 
00464         xf[j] = xx;
00465         yf[j] = yy;
00466         j++;
00467     }
00468 
00469     do_floor(xf, yf, n);
00470     n = do_filter(xi, yi, n);
00471 
00472     R_polydots_abs(xi, yi, j);
00473 }
00474 
00475 static int cull_polyline_plane(int *pn, const double *x, const double *y,
00476                                const struct plane *p)
00477 {
00478     int n = *pn;
00479     int last = -1;
00480     int prev = 0;
00481     double x0 = x[prev];
00482     double y0 = y[prev];
00483     double d0 = dist_plane(x0, y0, p);
00484     int i, j;
00485 
00486     for (i = 0, j = 0; i < n; i++) {
00487         double x1 = x[i];
00488         double y1 = y[i];
00489         double d1 = dist_plane(x1, y1, p);
00490         int in0 = d0 <= 0;
00491         int in1 = d1 <= 0;
00492 
00493         if (!in0 && in1 && last != prev) {      /* entering */
00494             alloc_float(j + 1);
00495             xf[j] = x0;
00496             yf[j] = y0;
00497             j++;
00498             last = prev;
00499         }
00500 
00501         if (in1 || in0) {       /* inside or leaving */
00502             alloc_float(j + 1);
00503             xf[j] = x1;
00504             yf[j] = y1;
00505             j++;
00506             last = i;
00507         }
00508 
00509         x0 = x1;
00510         y0 = y1;
00511         d0 = d1;
00512         prev = i;
00513     }
00514 
00515     *pn = j;
00516 
00517     return (j == 0);
00518 }
00519 
00520 static void polyline_cull(const double *x, const double *y, int n)
00521 {
00522     alloc_float(n + 10);
00523 
00524     if (cull_polyline_plane(&n, x, y, &pl_left))
00525         return;
00526 
00527     dealloc_float(&x, &y, 0);
00528 
00529     if (cull_polyline_plane(&n, x, y, &pl_rite))
00530         return;
00531 
00532     dealloc_float(&x, &y, 1);
00533 
00534     if (cull_polyline_plane(&n, x, y, &pl_bot))
00535         return;
00536 
00537     dealloc_float(&x, &y, 1);
00538 
00539     if (cull_polyline_plane(&n, x, y, &pl_top))
00540         return;
00541 
00542     dealloc_float(&x, &y, 1);
00543 
00544     do_floor(x, y, n);
00545     n = do_filter(xi, yi, n);
00546 
00547     R_polyline_abs(xi, yi, n);
00548 }
00549 
00550 void D_polyline_cull(const double *x, const double *y, int n)
00551 {
00552     if (n < 2)
00553         return;
00554 
00555     if (!window_set)
00556         D_clip_to_map();
00557 
00558     if (D_is_lat_lon())
00559         do_ll_wrap(x, y, n, polyline_cull);
00560     else
00561         polyline_cull(x, y, n);
00562 }
00563 
00564 static void polyline_clip(const double *x, const double *y, int n)
00565 {
00566     int i;
00567 
00568     for (i = 1; i < n; i++)
00569         line_clip(x[i - 1], y[i - 1], x[i], y[i]);
00570 }
00571 
00572 void D_polyline_clip(const double *x, const double *y, int n)
00573 {
00574     if (n < 2)
00575         return;
00576 
00577     if (!window_set)
00578         D_clip_to_map();
00579 
00580     if (D_is_lat_lon())
00581         do_ll_wrap(x, y, n, polyline_clip);
00582     else
00583         polyline_clip(x, y, n);
00584 }
00585 
00586 static int cull_polygon_plane(int *pn, const double *x, const double *y,
00587                               const struct plane *p)
00588 {
00589     int n = *pn;
00590     int last = -1;
00591     int prev = n - 1;
00592     double x0 = x[prev];
00593     double y0 = y[prev];
00594     double d0 = dist_plane(x0, y0, p);
00595     int i, j;
00596 
00597     for (i = j = 0; i < n; i++) {
00598         double x1 = x[i];
00599         double y1 = y[i];
00600         double d1 = dist_plane(x1, y1, p);
00601         int in0 = d0 <= 0;
00602         int in1 = d1 <= 0;
00603 
00604         if (!in0 && in1 && last != prev) {      /* entering */
00605             alloc_float(j + 1);
00606             xf[j] = x0;
00607             yf[j] = y0;
00608             j++;
00609             last = prev;
00610         }
00611 
00612         if (in1 || in0) {       /* inside or leaving */
00613             alloc_float(j + 1);
00614             xf[j] = x1;
00615             yf[j] = y1;
00616             j++;
00617             last = i;
00618         }
00619 
00620         x0 = x1;
00621         y0 = y1;
00622         d0 = d1;
00623         prev = i;
00624     }
00625 
00626     *pn = j;
00627 
00628     return (j == 0);
00629 }
00630 
00631 static void polygon_cull(const double *x, const double *y, int n)
00632 {
00633     alloc_float(n + 10);
00634 
00635     if (cull_polygon_plane(&n, x, y, &pl_left))
00636         return;
00637 
00638     dealloc_float(&x, &y, 0);
00639 
00640     if (cull_polygon_plane(&n, x, y, &pl_rite))
00641         return;
00642 
00643     dealloc_float(&x, &y, 1);
00644 
00645     if (cull_polygon_plane(&n, x, y, &pl_bot))
00646         return;
00647 
00648     dealloc_float(&x, &y, 1);
00649 
00650     if (cull_polygon_plane(&n, x, y, &pl_top))
00651         return;
00652 
00653     dealloc_float(&x, &y, 1);
00654 
00655     do_round(x, y, n);
00656     n = do_filter(xi, yi, n);
00657 
00658     R_polygon_abs(xi, yi, n);
00659 }
00660 
00661 void D_polygon_cull(const double *x, const double *y, int n)
00662 {
00663     if (!window_set)
00664         D_clip_to_map();
00665 
00666     if (D_is_lat_lon())
00667         do_ll_wrap(x, y, n, polygon_cull);
00668     else
00669         polygon_cull(x, y, n);
00670 }
00671 
00672 static int clip_polygon_plane(int *pn, const double *x, const double *y,
00673                               const struct plane *p)
00674 {
00675     int n = *pn;
00676     double x0 = x[n - 1];
00677     double y0 = y[n - 1];
00678     double d0 = dist_plane(x0, y0, p);
00679     int i, j;
00680 
00681     for (i = j = 0; i < n; i++) {
00682         double x1 = x[i];
00683         double y1 = y[i];
00684         double d1 = dist_plane(x1, y1, p);
00685         int in0 = d0 <= 0;
00686         int in1 = d1 <= 0;
00687 
00688         if (in0 != in1) {       /* edge crossing */
00689             alloc_float(j + 1);
00690             xf[j] = interpolate(x0, x1, d0, d1);
00691             yf[j] = interpolate(y0, y1, d0, d1);
00692             j++;
00693         }
00694 
00695         if (in1) {              /* point inside */
00696             alloc_float(j + 1);
00697             xf[j] = x[i];
00698             yf[j] = y[i];
00699             j++;
00700         }
00701 
00702         x0 = x1;
00703         y0 = y1;
00704         d0 = d1;
00705     }
00706 
00707     *pn = j;
00708 
00709     return (j == 0);
00710 }
00711 
00712 static void polygon_clip(const double *x, const double *y, int n)
00713 {
00714     alloc_float(n + 10);
00715 
00716     if (clip_polygon_plane(&n, x, y, &pl_left))
00717         return;
00718 
00719     dealloc_float(&x, &y, 0);
00720 
00721     if (clip_polygon_plane(&n, x, y, &pl_rite))
00722         return;
00723 
00724     dealloc_float(&x, &y, 1);
00725 
00726     if (clip_polygon_plane(&n, x, y, &pl_bot))
00727         return;
00728 
00729     dealloc_float(&x, &y, 1);
00730 
00731     if (clip_polygon_plane(&n, x, y, &pl_top))
00732         return;
00733 
00734     dealloc_float(&x, &y, 1);
00735 
00736     do_round(x, y, n);
00737     n = do_filter(xi, yi, n);
00738 
00739     R_polygon_abs(xi, yi, n);
00740 }
00741 
00742 void D_polygon_clip(const double *x, const double *y, int n)
00743 {
00744     if (!window_set)
00745         D_clip_to_map();
00746 
00747     if (D_is_lat_lon())
00748         do_ll_wrap(x, y, n, polygon_clip);
00749     else
00750         polygon_clip(x, y, n);
00751 }
00752 
00753 static void box_clip(double x1, double y1, double x2, double y2)
00754 {
00755     double t, b, l, r;
00756     int ti, bi, li, ri;
00757 
00758     l = max(clip.left, min(x1, x2));
00759     r = min(clip.rite, max(x1, x2));
00760     b = max(clip.bot, min(y1, y2));
00761     t = min(clip.top, max(y1, y2));
00762 
00763     li = round(D_u_to_d_col(l));
00764     ri = round(D_u_to_d_col(r));
00765     bi = round(D_u_to_d_row(b));
00766     ti = round(D_u_to_d_row(t));
00767 
00768     R_box_abs(li, ti, ri, bi);
00769 }
00770 
00771 static void box_clip_ll(double x1, double y1, double x2, double y2)
00772 {
00773     double ux0 = clip.left;
00774     double ux1 = clip.rite;
00775     int lo, hi, i;
00776 
00777     x2 = x1 + coerce(x2 - x1);
00778 
00779     lo = -shift_count(ux1 - x1);
00780     hi = shift_count(x2 - ux0);
00781 
00782     for (i = lo; i <= hi; i++)
00783         box_clip(x1 + i * 360, y1, x2 + i * 360, y2);
00784 }
00785 
00786 void D_box_clip(double x1, double y1, double x2, double y2)
00787 {
00788     if (!window_set)
00789         D_clip_to_map();
00790 
00791     if (D_is_lat_lon())
00792         box_clip_ll(x1, y1, x2, y2);
00793     else
00794         box_clip(x1, y1, x2, y2);
00795 }
00796 
00797 void D_move(double x, double y)
00798 {
00799     int xi = round(D_u_to_d_col(x));
00800     int yi = round(D_u_to_d_row(y));
00801 
00802     R_move_abs(xi, yi);
00803 }
00804 
00805 void D_cont(double x, double y)
00806 {
00807     int xi = round(D_u_to_d_col(x));
00808     int yi = round(D_u_to_d_row(y));
00809 
00810     R_cont_abs(xi, yi);
00811 }
00812 
00813 void D_polydots(const double *x, const double *y, int n)
00814 {
00815     do_floor(x, y, n);
00816     n = do_filter(xi, yi, n);
00817     R_polydots_abs(xi, yi, n);
00818 }
00819 
00820 void D_polyline(const double *x, const double *y, int n)
00821 {
00822     do_floor(x, y, n);
00823     n = do_filter(xi, yi, n);
00824     R_polyline_abs(xi, yi, n);
00825 }
00826 
00827 void D_polygon(const double *x, const double *y, int n)
00828 {
00829     do_round(x, y, n);
00830     n = do_filter(xi, yi, n);
00831     R_polygon_abs(xi, yi, n);
00832 }
00833 
00834 void D_box(double x1, double y1, double x2, double y2)
00835 {
00836     double l = min(x1, x2);
00837     double r = max(x1, x2);
00838     double b = min(y1, y2);
00839     double t = max(y1, y2);
00840     int li = round(D_u_to_d_col(l));
00841     int ri = round(D_u_to_d_col(r));
00842     int bi = round(D_u_to_d_row(b));
00843     int ti = round(D_u_to_d_row(t));
00844 
00845     R_box_abs(li, ti, ri, bi);
00846 }
00847 
00848 void D_line_width(double d)
00849 {
00850     int w = round(d);
00851 
00852     if (w < 0)
00853         w = 0;
00854 
00855     R_line_width(w);
00856 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines