GRASS Programmer's Manual  6.4.3(2013)-r
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Macros Pages
plot.c
Go to the documentation of this file.
1 
2 /*****************************************************************
3  * Plot lines and filled polygons. Input space is database window.
4  * Output space and output functions are user defined.
5  * Converts input east,north lines and polygons to output x,y
6  * and calls user supplied line drawing routines to do the plotting.
7  *
8  * Handles global wrap-around for lat-lon databases.
9  *
10  * Does not perform window clipping.
11  * Clipping must be done by the line draw routines supplied by the user.
12  *
13  * Note:
14  * Hopefully, cartographic style projection plotting will be added later.
15  *******************************************************************/
16 #include <stdlib.h>
17 #include <math.h>
18 #include <grass/gis.h>
19 
20 static double xconv, yconv;
21 static double left, right, top, bottom;
22 static int ymin, ymax;
23 static struct Cell_head window;
24 static int fastline(double, double, double, double);
25 static int slowline(double, double, double, double);
26 static int plot_line(double, double, double, double, int (*)());
27 static double nearest(double, double);
28 static int edge(double, double, double, double);
29 static int edge_point(double, int);
30 
31 #define POINT struct point
33  double x;
34  int y;
35 };
36 static int edge_order(const void *, const void *);
37 static int row_solid_fill(int, double, double);
38 static int row_dotted_fill(int, double, double);
39 static int dotted_fill_gap = 2;
40 static int ifloor(double);
41 static int iceil(double);
42 static int (*row_fill) () = row_solid_fill;
43 static int (*move) (int, int);
44 static int (*cont) (int, int);
45 
60 /*
61  * G_setup_plot (t, b, l, r, Move, Cont)
62  * double t, b, l, r;
63  * int (*Move)(), (*Cont)();
64  *
65  * initialize the plotting capability.
66  * t,b,l,r: top, bottom, left, right of the output x,y coordinate space.
67  * Move,Cont: subroutines that will draw lines in x,y space.
68  * Move(x,y) move to x,y (no draw)
69  * Cont(x,y) draw from previous position to x,y
70  * Notes:
71  * Cont() is responsible for clipping.
72  * The t,b,l,r are only used to compute coordinate transformations.
73  * The input space is assumed to be the current GRASS window.
74  */
75 
97 int G_setup_plot(double t, double b, double l, double r,
98  int (*Move) (int, int), int (*Cont) (int, int))
99 {
101 
102  left = l;
103  right = r;
104  top = t;
105  bottom = b;
106 
107  xconv = (right - left) / (window.east - window.west);
108  yconv = (bottom - top) / (window.north - window.south);
109 
110  if (top < bottom) {
111  ymin = iceil(top);
112  ymax = ifloor(bottom);
113  }
114  else {
115  ymin = iceil(bottom);
116  ymax = ifloor(top);
117  }
118 
119  move = Move;
120  cont = Cont;
121 
122  return 0;
123 }
124 
136 int G_setup_fill(int gap)
137 {
138  if (gap > 0) {
139  row_fill = row_dotted_fill;
140  dotted_fill_gap = gap + 1;
141  }
142  else
143  row_fill = row_solid_fill;
144 
145  return 0;
146 }
147 
148 #define X(e) (left + xconv * ((e) - window.west))
149 #define Y(n) (top + yconv * (window.north - (n)))
150 
151 #define EAST(x) (window.west + ((x)-left)/xconv)
152 #define NORTH(y) (window.north - ((y)-top)/yconv)
153 
154 
168 int G_plot_where_xy(double east, double north, int *x, int *y)
169 {
170  *x = ifloor(X(G_adjust_easting(east, &window)) + 0.5);
171  *y = ifloor(Y(north) + 0.5);
172 
173  return 0;
174 }
175 
176 
190 int G_plot_where_en(int x, int y, double *east, double *north)
191 {
192  *east = G_adjust_easting(EAST(x), &window);
193  *north = NORTH(y);
194 
195  return 0;
196 }
197 
198 int G_plot_point(double east, double north)
199 {
200  int x, y;
201 
202  G_plot_where_xy(east, north, &x, &y);
203  move(x, y);
204  cont(x, y);
205 
206  return 0;
207 }
208 
209 /*
210  * Line in map coordinates is plotted in output x,y coordinates
211  * This routine handles global wrap-around for lat-long databses.
212  *
213  */
214 
230 int G_plot_line(double east1, double north1, double east2, double north2)
231 {
232  return plot_line(east1, north1, east2, north2, fastline);
233 }
234 
235 int G_plot_line2(double east1, double north1, double east2, double north2)
236 {
237  return plot_line(east1, north1, east2, north2, slowline);
238 }
239 
240 /* fastline converts double rows/cols to ints then plots
241  * this is ok for graphics, but not the best for vector to raster
242  */
243 
244 static int fastline(double x1, double y1, double x2, double y2)
245 {
246  move(ifloor(x1 + 0.5), ifloor(y1 + 0.5));
247  cont(ifloor(x2 + 0.5), ifloor(y2 + 0.5));
248 
249  return 0;
250 }
251 
252 /* NOTE (shapiro):
253  * I think the adding of 0.5 in slowline is not correct
254  * the output window (left, right, top, bottom) should already
255  * be adjusted for this: left=-0.5; right = window.cols-0.5;
256  */
257 
258 static int slowline(double x1, double y1, double x2, double y2)
259 {
260  double dx, dy;
261  double m, b;
262  int xstart, xstop, ystart, ystop;
263 
264  dx = x2 - x1;
265  dy = y2 - y1;
266 
267  if (fabs(dx) > fabs(dy)) {
268  m = dy / dx;
269  b = y1 - m * x1;
270 
271  if (x1 > x2) {
272  xstart = iceil(x2 - 0.5);
273  xstop = ifloor(x1 + 0.5);
274  }
275  else {
276  xstart = iceil(x1 - 0.5);
277  xstop = ifloor(x2 + 0.5);
278  }
279  if (xstart <= xstop) {
280  ystart = ifloor(m * xstart + b + 0.5);
281  move(xstart, ystart);
282  while (xstart <= xstop) {
283  cont(xstart++, ystart);
284  ystart = ifloor(m * xstart + b + 0.5);
285  }
286  }
287  }
288  else {
289  if (dx == dy) /* they both might be 0 */
290  m = 1.;
291  else
292  m = dx / dy;
293  b = x1 - m * y1;
294 
295  if (y1 > y2) {
296  ystart = iceil(y2 - 0.5);
297  ystop = ifloor(y1 + 0.5);
298  }
299  else {
300  ystart = iceil(y1 - 0.5);
301  ystop = ifloor(y2 + 0.5);
302  }
303  if (ystart <= ystop) {
304  xstart = ifloor(m * ystart + b + 0.5);
305  move(xstart, ystart);
306  while (ystart <= ystop) {
307  cont(xstart, ystart++);
308  xstart = ifloor(m * ystart + b + 0.5);
309  }
310  }
311  }
312 
313  return 0;
314 }
315 
316 static int plot_line(double east1, double north1, double east2, double north2,
317  int (*line) (double, double, double, double))
318 {
319  double x1, x2, y1, y2;
320 
321  y1 = Y(north1);
322  y2 = Y(north2);
323 
324  if (window.proj == PROJECTION_LL) {
325  if (east1 > east2)
326  while ((east1 - east2) > 180)
327  east2 += 360;
328  else if (east2 > east1)
329  while ((east2 - east1) > 180)
330  east1 += 360;
331  while (east1 > window.east) {
332  east1 -= 360.0;
333  east2 -= 360.0;
334  }
335  while (east1 < window.west) {
336  east1 += 360.0;
337  east2 += 360.0;
338  }
339  x1 = X(east1);
340  x2 = X(east2);
341 
342  line(x1, y1, x2, y2);
343 
344  if (east2 > window.east || east2 < window.west) {
345  while (east2 > window.east) {
346  east1 -= 360.0;
347  east2 -= 360.0;
348  }
349  while (east2 < window.west) {
350  east1 += 360.0;
351  east2 += 360.0;
352  }
353  x1 = X(east1);
354  x2 = X(east2);
355  line(x1, y1, x2, y2);
356  }
357  }
358  else {
359  x1 = X(east1);
360  x2 = X(east2);
361  line(x1, y1, x2, y2);
362  }
363 
364  return 0;
365 }
366 
367 /*
368  * G_plot_polygon (x, y, n)
369  *
370  * double *x x coordinates of vertices
371  * double *y y coordinates of vertices
372  * int n number of verticies
373  *
374  * polygon fill from map coordinate space to plot x,y space.
375  * for lat-lon, handles global wrap-around as well as polar polygons.
376  *
377  * returns 0 ok, 2 n<3, -1 weird internal error, 1 no memory
378  */
379 
380 static POINT *P;
381 static int np;
382 static int npalloc = 0;
383 
384 #define OK 0
385 #define TOO_FEW_EDGES 2
386 #define NO_MEMORY 1
387 #define OUT_OF_SYNC -1
388 
389 static double nearest(double e0, double e1)
390 {
391  while (e0 - e1 > 180)
392  e1 += 360.0;
393  while (e1 - e0 > 180)
394  e1 -= 360.0;
395 
396  return e1;
397 }
398 
399 
412 int G_plot_polygon(const double *x, const double *y, int n)
413 {
414  int i;
415  int pole;
416  double x0, x1;
417  double y0, y1;
418  double shift, E, W = 0L;
419  double e0, e1;
420  int shift1, shift2;
421 
422  if (n < 3)
423  return TOO_FEW_EDGES;
424 
425  /* traverse the perimeter */
426 
427  np = 0;
428  shift1 = 0;
429 
430  /* global wrap-around for lat-lon, part1 */
431  if (window.proj == PROJECTION_LL) {
432  /*
433  pole = G_pole_in_polygon(x,y,n);
434  */
435  pole = 0;
436 
437  e0 = x[n - 1];
438  E = W = e0;
439 
440  x0 = X(e0);
441  y0 = Y(y[n - 1]);
442 
443  if (pole && !edge(x0, y0, x0, Y(90.0 * pole)))
444  return NO_MEMORY;
445 
446  for (i = 0; i < n; i++) {
447  e1 = nearest(e0, x[i]);
448  if (e1 > E)
449  E = e1;
450  if (e1 < W)
451  W = e1;
452 
453  x1 = X(e1);
454  y1 = Y(y[i]);
455 
456  if (!edge(x0, y0, x1, y1))
457  return NO_MEMORY;
458 
459  x0 = x1;
460  y0 = y1;
461  e0 = e1;
462  }
463  if (pole && !edge(x0, y0, x0, Y(90.0 * pole)))
464  return NO_MEMORY;
465 
466  shift = 0; /* shift into window */
467  while (E + shift > window.east)
468  shift -= 360.0;
469  while (E + shift < window.west)
470  shift += 360.0;
471  shift1 = X(x[n - 1] + shift) - X(x[n - 1]);
472  }
473  else {
474  x0 = X(x[n - 1]);
475  y0 = Y(y[n - 1]);
476 
477  for (i = 0; i < n; i++) {
478  x1 = X(x[i]);
479  y1 = Y(y[i]);
480  if (!edge(x0, y0, x1, y1))
481  return NO_MEMORY;
482  x0 = x1;
483  y0 = y1;
484  }
485  }
486 
487  /* check if perimeter has odd number of points */
488  if (np % 2)
489  return OUT_OF_SYNC;
490 
491  /* sort the edge points by col(x) and then by row(y) */
492  qsort(P, np, sizeof(POINT), &edge_order);
493 
494  /* plot */
495  for (i = 1; i < np; i += 2) {
496  if (P[i].y != P[i - 1].y)
497  return OUT_OF_SYNC;
498  row_fill(P[i].y, P[i - 1].x + shift1, P[i].x + shift1);
499  }
500  if (window.proj == PROJECTION_LL) { /* now do wrap-around, part 2 */
501  shift = 0;
502  while (W + shift < window.west)
503  shift += 360.0;
504  while (W + shift > window.east)
505  shift -= 360.0;
506  shift2 = X(x[n - 1] + shift) - X(x[n - 1]);
507  if (shift2 != shift1) {
508  for (i = 1; i < np; i += 2) {
509  row_fill(P[i].y, P[i - 1].x + shift2, P[i].x + shift2);
510  }
511  }
512  }
513  return OK;
514 }
515 
516 /*
517  * G_plot_area (xs, ys, rpnts, rings)
518  * double **xs; -- pointer to pointer for X's
519  * double **ys; -- pointer to pointer for Y's
520  * int *rpnts; -- array of ints w/ num points per ring
521  * int rings; -- number of rings
522  *
523  * Essentially a copy of G_plot_polygon, with minor mods to
524  * handle a set of polygons. return values are the same.
525  */
526 
542 int G_plot_area(double *const *xs, double *const *ys, int *rpnts, int rings)
543 {
544  int i, j, n;
545  int pole;
546  double x0, x1, *x;
547  double y0, y1, *y;
548  double shift, E, W = 0L;
549  double e0, e1;
550  int *shift1 = NULL, shift2;
551 
552  /* traverse the perimeter */
553 
554  np = 0;
555  shift1 = (int *)G_calloc(sizeof(int), rings);
556 
557  for (j = 0; j < rings; j++) {
558  n = rpnts[j];
559 
560  if (n < 3)
561  return TOO_FEW_EDGES;
562 
563  x = xs[j];
564  y = ys[j];
565 
566  /* global wrap-around for lat-lon, part1 */
567  if (window.proj == PROJECTION_LL) {
568  /*
569  pole = G_pole_in_polygon(x,y,n);
570  */
571  pole = 0;
572 
573  e0 = x[n - 1];
574  E = W = e0;
575 
576  x0 = X(e0);
577  y0 = Y(y[n - 1]);
578 
579  if (pole && !edge(x0, y0, x0, Y(90.0 * pole)))
580  return NO_MEMORY;
581 
582  for (i = 0; i < n; i++) {
583  e1 = nearest(e0, x[i]);
584  if (e1 > E)
585  E = e1;
586  if (e1 < W)
587  W = e1;
588 
589  x1 = X(e1);
590  y1 = Y(y[i]);
591 
592  if (!edge(x0, y0, x1, y1))
593  return NO_MEMORY;
594 
595  x0 = x1;
596  y0 = y1;
597  e0 = e1;
598  }
599  if (pole && !edge(x0, y0, x0, Y(90.0 * pole)))
600  return NO_MEMORY;
601 
602  shift = 0; /* shift into window */
603  while (E + shift > window.east)
604  shift -= 360.0;
605  while (E + shift < window.west)
606  shift += 360.0;
607  shift1[j] = X(x[n - 1] + shift) - X(x[n - 1]);
608  }
609  else {
610  x0 = X(x[n - 1]);
611  y0 = Y(y[n - 1]);
612 
613  for (i = 0; i < n; i++) {
614  x1 = X(x[i]);
615  y1 = Y(y[i]);
616  if (!edge(x0, y0, x1, y1))
617  return NO_MEMORY;
618  x0 = x1;
619  y0 = y1;
620  }
621  }
622  } /* for() */
623 
624  /* check if perimeter has odd number of points */
625  if (np % 2)
626  return OUT_OF_SYNC;
627 
628  /* sort the edge points by col(x) and then by row(y) */
629  qsort(P, np, sizeof(POINT), &edge_order);
630 
631  /* plot */
632  for (j = 0; j < rings; j++) {
633  for (i = 1; i < np; i += 2) {
634  if (P[i].y != P[i - 1].y)
635  return OUT_OF_SYNC;
636  row_fill(P[i].y, P[i - 1].x + shift1[j], P[i].x + shift1[j]);
637  }
638  if (window.proj == PROJECTION_LL) { /* now do wrap-around, part 2 */
639  n = rpnts[j];
640  x = xs[j];
641  y = ys[j];
642 
643  shift = 0;
644  while (W + shift < window.west)
645  shift += 360.0;
646  while (W + shift > window.east)
647  shift -= 360.0;
648  shift2 = X(x[n - 1] + shift) - X(x[n - 1]);
649  if (shift2 != shift1[j]) {
650  for (i = 1; i < np; i += 2) {
651  row_fill(P[i].y, P[i - 1].x + shift2, P[i].x + shift2);
652  }
653  }
654  }
655  }
656  G_free(shift1);
657  return OK;
658 
659 }
660 
661 static int edge(double x0, double y0, double x1, double y1)
662 {
663  register double m;
664  double dy, x;
665  int ystart, ystop;
666 
667 
668  /* tolerance to avoid FPE */
669  dy = y0 - y1;
670  if (fabs(dy) < 1e-10)
671  return 1;
672 
673  m = (x0 - x1) / dy;
674 
675  if (y0 < y1) {
676  ystart = iceil(y0);
677  ystop = ifloor(y1);
678  if (ystop == y1)
679  ystop--; /* if line stops at row center, don't include point */
680  }
681  else {
682  ystart = iceil(y1);
683  ystop = ifloor(y0);
684  if (ystop == y0)
685  ystop--; /* if line stops at row center, don't include point */
686  }
687  if (ystart > ystop)
688  return 1; /* does not cross center line of row */
689 
690  x = m * (ystart - y0) + x0;
691  while (ystart <= ystop) {
692  if (!edge_point(x, ystart++))
693  return 0;
694  x += m;
695  }
696  return 1;
697 }
698 
699 static int edge_point(double x, int y)
700 {
701 
702  if (y < ymin || y > ymax)
703  return 1;
704  if (np >= npalloc) {
705  if (npalloc > 0) {
706  npalloc *= 2;
707  P = (POINT *) G_realloc(P, npalloc * sizeof(POINT));
708  }
709  else {
710  npalloc = 32;
711  P = (POINT *) G_malloc(npalloc * sizeof(POINT));
712  }
713  if (P == NULL) {
714  npalloc = 0;
715  return 0;
716  }
717  }
718  P[np].x = x;
719  P[np++].y = y;
720  return 1;
721 }
722 
723 static int edge_order(const void *aa, const void *bb)
724 {
725  const struct point *a = aa, *b = bb;
726 
727  if (a->y < b->y)
728  return (-1);
729  if (a->y > b->y)
730  return (1);
731 
732  if (a->x < b->x)
733  return (-1);
734  if (a->x > b->x)
735  return (1);
736 
737  return (0);
738 }
739 
740 static int row_solid_fill(int y, double x1, double x2)
741 {
742  int i1, i2;
743 
744  i1 = iceil(x1);
745  i2 = ifloor(x2);
746  if (i1 <= i2) {
747  move(i1, y);
748  cont(i2, y);
749  }
750 
751  return 0;
752 }
753 
754 static int row_dotted_fill(int y, double x1, double x2)
755 {
756  int i1, i2, i;
757 
758  if (y != iceil(y / dotted_fill_gap) * dotted_fill_gap)
759  return 0;
760 
761  i1 = iceil(x1 / dotted_fill_gap) * dotted_fill_gap;
762  i2 = ifloor(x2);
763  if (i1 <= i2) {
764  for (i = i1; i <= i2; i += dotted_fill_gap) {
765  move(i, y);
766  cont(i, y);
767  }
768  }
769 
770  return 0;
771 }
772 
773 static int ifloor(double x)
774 {
775  int i;
776 
777  i = (int)x;
778  if (i > x)
779  i--;
780  return i;
781 }
782 
783 static int iceil(double x)
784 {
785  int i;
786 
787  i = (int)x;
788  if (i < x)
789  i++;
790  return i;
791 }
792 
793 /*
794  * G_plot_fx(e1,e2)
795  *
796  * plot f(x) from x=e1 to x=e2
797  */
798 
799 
811 int G_plot_fx(double (*f) (double), double east1, double east2)
812 {
813  double east, north, north1;
814  double incr;
815 
816 
817  incr = fabs(1.0 / xconv);
818 
819  east = east1;
820  north = f(east1);
821 
822  if (east1 > east2) {
823  while ((east1 -= incr) > east2) {
824  north1 = f(east1);
825  G_plot_line(east, north, east1, north1);
826  north = north1;
827  east = east1;
828  }
829  }
830  else {
831  while ((east1 += incr) < east2) {
832  north1 = f(east1);
833  G_plot_line(east, north, east1, north1);
834  north = north1;
835  east = east1;
836  }
837  }
838  G_plot_line(east, north, east2, f(east2));
839 
840  return 0;
841 }