GRASS Programmer's Manual  6.5.svn(2014)-r66266
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
gsdrape.c
Go to the documentation of this file.
1 
50 #include <stdlib.h>
51 
52 #include <grass/gstypes.h>
53 #include <grass/glocale.h>
54 
55 #include "gsget.h"
56 #include "rowcol.h"
57 #include "math.h"
58 
59 #define DONT_INTERSECT 0
60 #define DO_INTERSECT 1
61 #define COLLINEAR 2
62 
63 #define LERP(a,l,h) ((l)+(((h)-(l))*(a)))
64 #define EQUAL(a,b) (fabs((a)-(b))<EPSILON)
65 #define ISNODE(p,res) (fmod((double)p,(double)res)<EPSILON)
66 
67 #define SAME_SIGNS( a, b ) \
68  ((a >= 0 && b >= 0) || (a < 0 && b < 0))
69 
70 static int drape_line_init(int, int);
71 static Point3 *_gsdrape_get_segments(geosurf *, float *, float *, int *);
72 static float dist_squared_2d(float *, float *);
73 
74 /* array of points to be returned */
75 static Point3 *I3d;
76 
77 /* make dependent on resolution? */
78 static float EPSILON = 0.000001;
79 
80 /*vertical, horizontal, & diagonal intersections */
81 static Point3 *Vi, *Hi, *Di;
82 
83 static typbuff *Ebuf; /* elevation buffer */
84 static int Flat;
85 
86 
96 static int drape_line_init(int rows, int cols)
97 {
98  /* use G_calloc() [-> G_fatal_error] instead of calloc ? */
99  if (NULL == (I3d = (Point3 *) calloc(2 * (rows + cols), sizeof(Point3)))) {
100  return (-1);
101  }
102 
103  if (NULL == (Vi = (Point3 *) calloc(cols, sizeof(Point3)))) {
104  G_free(I3d);
105 
106  return (-1);
107  }
108 
109  if (NULL == (Hi = (Point3 *) calloc(rows, sizeof(Point3)))) {
110  G_free(I3d);
111  G_free(Vi);
112 
113  return (-1);
114  }
115 
116  if (NULL == (Di = (Point3 *) calloc(rows + cols, sizeof(Point3)))) {
117  G_free(I3d);
118  G_free(Vi);
119  G_free(Hi);
120 
121  return (-1);
122  }
123 
124  return (1);
125 }
126 
137 static Point3 *_gsdrape_get_segments(geosurf * gs, float *bgn, float *end,
138  int *num)
139 {
140  float f[3], l[3];
141  int vi, hi, di;
142  float dir[2], yres, xres;
143 
144  xres = VXRES(gs);
145  yres = VYRES(gs);
146 
147  vi = hi = di = 0;
148  GS_v2dir(bgn, end, dir);
149 
150  if (dir[X]) {
151  vi = get_vert_intersects(gs, bgn, end, dir);
152  }
153 
154  if (dir[Y]) {
155  hi = get_horz_intersects(gs, bgn, end, dir);
156  }
157 
158  if (!((end[Y] - bgn[Y]) / (end[X] - bgn[X]) == yres / xres)) {
159  di = get_diag_intersects(gs, bgn, end, dir);
160  }
161 
162  interp_first_last(gs, bgn, end, f, l);
163 
164  *num = order_intersects(gs, f, l, vi, hi, di);
165  /* fills in return values, eliminates dupes (corners) */
166 
167  G_debug(5, "_gsdrape_get_segments vi=%d, hi=%d, di=%d, num=%d",
168  vi, hi, di, *num);
169 
170  return (I3d);
171 }
172 
181 static float dist_squared_2d(float *p1, float *p2)
182 {
183  float dx, dy;
184 
185  dx = p2[X] - p1[X];
186  dy = p2[Y] - p1[Y];
187 
188  return (dx * dx + dy * dy);
189 }
190 
199 int gsdrape_set_surface(geosurf * gs)
200 {
201  static int first = 1;
202 
203  if (first) {
204  first = 0;
205 
206  if (0 > drape_line_init(gs->rows, gs->cols)) {
207  G_warning(_("Unable to process vector map - out of memory"));
208  Ebuf = NULL;
209 
210  return (-1);
211  }
212  }
213 
214  Ebuf = gs_get_att_typbuff(gs, ATT_TOPO, 0);
215 
216  return (1);
217 }
218 
233 int seg_intersect_vregion(geosurf * gs, float *bgn, float *end)
234 {
235  float *replace, xl, yb, xr, yt, xi, yi;
236  int inside = 0;
237 
238  xl = 0.0;
239  xr = VCOL2X(gs, VCOLS(gs));
240  yt = VROW2Y(gs, 0);
241  yb = VROW2Y(gs, VROWS(gs));
242 
243  if (in_vregion(gs, bgn)) {
244  replace = end;
245  inside++;
246  }
247 
248  if (in_vregion(gs, end)) {
249  replace = bgn;
250  inside++;
251  }
252 
253  if (inside == 2) {
254  return (1);
255  }
256  else if (inside) {
257  /* one in & one out - replace gets first intersection */
258  if (segs_intersect
259  (bgn[X], bgn[Y], end[X], end[Y], xl, yb, xl, yt, &xi, &yi)) {
260  /* left */
261  }
262  else if (segs_intersect
263  (bgn[X], bgn[Y], end[X], end[Y], xr, yb, xr, yt, &xi, &yi)) {
264  /* right */
265  }
266  else if (segs_intersect
267  (bgn[X], bgn[Y], end[X], end[Y], xl, yb, xr, yb, &xi, &yi)) {
268  /* bottom */
269  }
270  else if (segs_intersect
271  (bgn[X], bgn[Y], end[X], end[Y], xl, yt, xr, yt, &xi, &yi)) {
272  /* top */
273  }
274 
275  replace[X] = xi;
276  replace[Y] = yi;
277  }
278  else {
279  /* both out - find 2 intersects & replace both */
280  float pt1[2], pt2[2];
281 
282  replace = pt1;
283  if (segs_intersect
284  (bgn[X], bgn[Y], end[X], end[Y], xl, yb, xl, yt, &xi, &yi)) {
285  replace[X] = xi;
286  replace[Y] = yi;
287  replace = pt2;
288  inside++;
289  }
290 
291  if (segs_intersect
292  (bgn[X], bgn[Y], end[X], end[Y], xr, yb, xr, yt, &xi, &yi)) {
293  replace[X] = xi;
294  replace[Y] = yi;
295  replace = pt2;
296  inside++;
297  }
298 
299  if (inside < 2) {
300  if (segs_intersect
301  (bgn[X], bgn[Y], end[X], end[Y], xl, yb, xr, yb, &xi, &yi)) {
302  replace[X] = xi;
303  replace[Y] = yi;
304  replace = pt2;
305  inside++;
306  }
307  }
308 
309  if (inside < 2) {
310  if (segs_intersect
311  (bgn[X], bgn[Y], end[X], end[Y], xl, yt, xr, yt, &xi, &yi)) {
312  replace[X] = xi;
313  replace[Y] = yi;
314  inside++;
315  }
316  }
317 
318  if (inside < 2) {
319  return (0); /* no intersect or only 1 point on corner */
320  }
321 
322  /* compare dist of intersects to bgn - closest replaces bgn */
323  if (GS_P2distance(bgn, pt1) < GS_P2distance(bgn, pt2)) {
324  bgn[X] = pt1[X];
325  bgn[Y] = pt1[Y];
326  end[X] = pt2[X];
327  end[Y] = pt2[Y];
328  }
329  else {
330  bgn[X] = pt2[X];
331  bgn[Y] = pt2[Y];
332  end[X] = pt1[X];
333  end[Y] = pt1[Y];
334  }
335  }
336 
337  return (1);
338 }
339 
350 Point3 *gsdrape_get_segments(geosurf * gs, float *bgn, float *end, int *num)
351 {
353 
354  if (!seg_intersect_vregion(gs, bgn, end)) {
355  *num = 0;
356 
357  return (NULL);
358  }
359 
360  if (CONST_ATT == gs_get_att_src(gs, ATT_TOPO)) {
361  /* will probably want a force_drape option to get all intersects */
362  I3d[0][X] = bgn[X];
363  I3d[0][Y] = bgn[Y];
364  I3d[0][Z] = gs->att[ATT_TOPO].constant;
365  I3d[1][X] = end[X];
366  I3d[1][Y] = end[Y];
367  I3d[1][Z] = gs->att[ATT_TOPO].constant;
368  *num = 2;
369 
370  return (I3d);
371  }
372 
373  if (bgn[X] == end[X] && bgn[Y] == end[Y]) {
374  float f[3], l[3];
375 
376  interp_first_last(gs, bgn, end, f, l);
377  GS_v3eq(I3d[0], f);
378  GS_v3eq(I3d[1], l);
379 
380  /* CHANGE (*num = 1) to reflect degenerate line ? */
381  *num = 2;
382 
383  return (I3d);
384  }
385 
386  Flat = 0;
387  return (_gsdrape_get_segments(gs, bgn, end, num));
388 }
389 
390 
401 Point3 *gsdrape_get_allsegments(geosurf * gs, float *bgn, float *end,
402  int *num)
403 {
405 
406  if (!seg_intersect_vregion(gs, bgn, end)) {
407  *num = 0;
408  return (NULL);
409  }
410 
411  if (bgn[X] == end[X] && bgn[Y] == end[Y]) {
412  float f[3], l[3];
413 
414  interp_first_last(gs, bgn, end, f, l);
415  GS_v3eq(I3d[0], f);
416  GS_v3eq(I3d[1], l);
417  *num = 2;
418 
419  return (I3d);
420  }
421 
422  if (CONST_ATT == gs_get_att_src(gs, ATT_TOPO)) {
423  Flat = 1;
424  }
425  else {
426  Flat = 0;
427  }
428 
429  return (_gsdrape_get_segments(gs, bgn, end, num));
430 }
431 
441 void interp_first_last(geosurf * gs, float *bgn, float *end, Point3 f,
442  Point3 l)
443 {
444  f[X] = bgn[X];
445  f[Y] = bgn[Y];
446 
447  l[X] = end[X];
448  l[Y] = end[Y];
449 
450  if (Flat) {
451  f[Z] = l[Z] = gs->att[ATT_TOPO].constant;
452  }
453  else {
454  viewcell_tri_interp(gs, Ebuf, f, 0);
455  viewcell_tri_interp(gs, Ebuf, l, 0);
456  }
457 
458  return;
459 }
460 
467 int _viewcell_tri_interp(geosurf * gs, Point3 pt)
468 {
469  typbuff *buf;
470 
471  buf = gs_get_att_typbuff(gs, ATT_TOPO, 0);
472 
473  return (viewcell_tri_interp(gs, buf, pt, 0));
474 }
475 
507 int viewcell_tri_interp(geosurf * gs, typbuff * buf, Point3 pt,
508  int check_mask)
509 {
510  Point3 p1, p2, p3;
511  int offset, drow, dcol, vrow, vcol;
512  float xmax, ymin, ymax, alpha;
513 
514  xmax = VCOL2X(gs, VCOLS(gs));
515  ymax = VROW2Y(gs, 0);
516  ymin = VROW2Y(gs, VROWS(gs));
517 
518  if (check_mask) {
519  if (gs_point_is_masked(gs, pt)) {
520  return (0);
521  }
522  }
523 
524  if (pt[X] < 0.0 || pt[Y] > ymax) {
525  /* outside on left or top */
526  return (0);
527  }
528 
529  if (pt[Y] < ymin || pt[X] > xmax) {
530  /* outside on bottom or right */
531  return (0);
532  }
533 
534  if (CONST_ATT == gs_get_att_src(gs, ATT_TOPO)) {
535  pt[Z] = gs->att[ATT_TOPO].constant;
536 
537  return (1);
538  }
539  else if (MAP_ATT != gs_get_att_src(gs, ATT_TOPO)) {
540  return (0);
541  }
542 
543  vrow = Y2VROW(gs, pt[Y]);
544  vcol = X2VCOL(gs, pt[X]);
545 
546  if (vrow < VROWS(gs) && vcol < VCOLS(gs)) {
547  /*not on bottom or right edge */
548  if (pt[X] > 0.0 && pt[Y] < ymax) {
549  /* not on left or top edge */
550  p1[X] = VCOL2X(gs, vcol + 1);
551  p1[Y] = VROW2Y(gs, vrow);
552  drow = VROW2DROW(gs, vrow);
553  dcol = VCOL2DCOL(gs, vcol + 1);
554  offset = DRC2OFF(gs, drow, dcol);
555  GET_MAPATT(buf, offset, p1[Z]); /* top right */
556 
557  p2[X] = VCOL2X(gs, vcol);
558  p2[Y] = VROW2Y(gs, vrow + 1);
559  drow = VROW2DROW(gs, vrow + 1);
560  dcol = VCOL2DCOL(gs, vcol);
561  offset = DRC2OFF(gs, drow, dcol);
562  GET_MAPATT(buf, offset, p2[Z]); /* bottom left */
563 
564  if ((pt[X] - p2[X]) / VXRES(gs) > (pt[Y] - p2[Y]) / VYRES(gs)) {
565  /* lower triangle */
566  p3[X] = VCOL2X(gs, vcol + 1);
567  p3[Y] = VROW2Y(gs, vrow + 1);
568  drow = VROW2DROW(gs, vrow + 1);
569  dcol = VCOL2DCOL(gs, vcol + 1);
570  offset = DRC2OFF(gs, drow, dcol);
571  GET_MAPATT(buf, offset, p3[Z]); /* bottom right */
572  }
573  else {
574  /* upper triangle */
575  p3[X] = VCOL2X(gs, vcol);
576  p3[Y] = VROW2Y(gs, vrow);
577  drow = VROW2DROW(gs, vrow);
578  dcol = VCOL2DCOL(gs, vcol);
579  offset = DRC2OFF(gs, drow, dcol);
580  GET_MAPATT(buf, offset, p3[Z]); /* top left */
581  }
582 
583  return (Point_on_plane(p1, p2, p3, pt));
584  }
585  else if (pt[X] == 0.0) {
586  /* on left edge */
587  if (pt[Y] < ymax) {
588  vrow = Y2VROW(gs, pt[Y]);
589  drow = VROW2DROW(gs, vrow);
590  offset = DRC2OFF(gs, drow, 0);
591  GET_MAPATT(buf, offset, p1[Z]);
592 
593  drow = VROW2DROW(gs, vrow + 1);
594  offset = DRC2OFF(gs, drow, 0);
595  GET_MAPATT(buf, offset, p2[Z]);
596 
597  alpha = (VROW2Y(gs, vrow) - pt[Y]) / VYRES(gs);
598  pt[Z] = LERP(alpha, p1[Z], p2[Z]);
599  }
600  else {
601  /* top left corner */
602  GET_MAPATT(buf, 0, pt[Z]);
603  }
604 
605  return (1);
606  }
607  else if (pt[Y] == gs->yrange) {
608  /* on top edge, not a corner */
609  vcol = X2VCOL(gs, pt[X]);
610  dcol = VCOL2DCOL(gs, vcol);
611  GET_MAPATT(buf, dcol, p1[Z]);
612 
613  dcol = VCOL2DCOL(gs, vcol + 1);
614  GET_MAPATT(buf, dcol, p2[Z]);
615 
616  alpha = (pt[X] - VCOL2X(gs, vcol)) / VXRES(gs);
617  pt[Z] = LERP(alpha, p1[Z], p2[Z]);
618 
619  return (1);
620  }
621  }
622  else if (vrow == VROWS(gs)) {
623  /* on bottom edge */
624  drow = VROW2DROW(gs, VROWS(gs));
625 
626  if (pt[X] > 0.0 && pt[X] < xmax) {
627  /* not a corner */
628  vcol = X2VCOL(gs, pt[X]);
629  dcol = VCOL2DCOL(gs, vcol);
630  offset = DRC2OFF(gs, drow, dcol);
631  GET_MAPATT(buf, offset, p1[Z]);
632 
633  dcol = VCOL2DCOL(gs, vcol + 1);
634  offset = DRC2OFF(gs, drow, dcol);
635  GET_MAPATT(buf, offset, p2[Z]);
636 
637  alpha = (pt[X] - VCOL2X(gs, vcol)) / VXRES(gs);
638  pt[Z] = LERP(alpha, p1[Z], p2[Z]);
639 
640  return (1);
641  }
642  else if (pt[X] == 0.0) {
643  /* bottom left corner */
644  offset = DRC2OFF(gs, drow, 0);
645  GET_MAPATT(buf, offset, pt[Z]);
646 
647  return (1);
648  }
649  else {
650  /* bottom right corner */
651  dcol = VCOL2DCOL(gs, VCOLS(gs));
652  offset = DRC2OFF(gs, drow, dcol);
653  GET_MAPATT(buf, offset, pt[Z]);
654 
655  return (1);
656  }
657  }
658  else {
659  /* on right edge, not bottom corner */
660  dcol = VCOL2DCOL(gs, VCOLS(gs));
661 
662  if (pt[Y] < ymax) {
663  vrow = Y2VROW(gs, pt[Y]);
664  drow = VROW2DROW(gs, vrow);
665  offset = DRC2OFF(gs, drow, dcol);
666  GET_MAPATT(buf, offset, p1[Z]);
667 
668  drow = VROW2DROW(gs, vrow + 1);
669  offset = DRC2OFF(gs, drow, dcol);
670  GET_MAPATT(buf, offset, p2[Z]);
671 
672  alpha = (VROW2Y(gs, vrow) - pt[Y]) / VYRES(gs);
673  pt[Z] = LERP(alpha, p1[Z], p2[Z]);
674 
675  return (1);
676  }
677  else {
678  /* top right corner */
679  GET_MAPATT(buf, dcol, pt[Z]);
680 
681  return (1);
682  }
683  }
684 
685  return (0);
686 }
687 
696 int in_vregion(geosurf * gs, float *pt)
697 {
698  if (pt[X] >= 0.0 && pt[Y] <= gs->yrange) {
699  if (pt[X] <= VCOL2X(gs, VCOLS(gs))) {
700  return (pt[Y] >= VROW2Y(gs, VROWS(gs)));
701  }
702  }
703 
704  return (0);
705 }
706 
729 int order_intersects(geosurf * gs, Point3 first, Point3 last, int vi, int hi,
730  int di)
731 {
732  int num, i, found, cv, ch, cd, cnum;
733  float dv, dh, dd, big, cpoint[2];
734 
735  cv = ch = cd = cnum = 0;
736  num = vi + hi + di;
737 
738  cpoint[X] = first[X];
739  cpoint[Y] = first[Y];
740 
741  if (in_vregion(gs, first)) {
742  I3d[cnum][X] = first[X];
743  I3d[cnum][Y] = first[Y];
744  I3d[cnum][Z] = first[Z];
745  cnum++;
746  }
747 
748  /* TODO: big could still be less than first dist */
749  big = gs->yrange * gs->yrange + gs->xrange * gs->xrange; /*BIG distance */
750  dv = dh = dd = big;
751 
752  for (i = 0; i < num; i = cv + ch + cd) {
753  if (cv < vi) {
754  dv = dist_squared_2d(Vi[cv], cpoint);
755 
756  if (dv < EPSILON) {
757  cv++;
758  continue;
759  }
760  }
761  else {
762  dv = big;
763  }
764 
765  if (ch < hi) {
766  dh = dist_squared_2d(Hi[ch], cpoint);
767 
768  if (dh < EPSILON) {
769  ch++;
770  continue;
771  }
772  }
773  else {
774  dh = big;
775  }
776 
777  if (cd < di) {
778  dd = dist_squared_2d(Di[cd], cpoint);
779 
780  if (dd < EPSILON) {
781  cd++;
782  continue;
783  }
784  }
785  else {
786  dd = big;
787  }
788 
789  found = 0;
790 
791  if (cd < di) {
792  if (dd <= dv && dd <= dh) {
793  found = 1;
794  cpoint[X] = I3d[cnum][X] = Di[cd][X];
795  cpoint[Y] = I3d[cnum][Y] = Di[cd][Y];
796  I3d[cnum][Z] = Di[cd][Z];
797  cnum++;
798 
799  if (EQUAL(dd, dv)) {
800  cv++;
801  }
802 
803  if (EQUAL(dd, dh)) {
804  ch++;
805  }
806 
807  cd++;
808  }
809  }
810 
811  if (!found) {
812  if (cv < vi) {
813  if (dv <= dh) {
814  found = 1;
815  cpoint[X] = I3d[cnum][X] = Vi[cv][X];
816  cpoint[Y] = I3d[cnum][Y] = Vi[cv][Y];
817  I3d[cnum][Z] = Vi[cv][Z];
818  cnum++;
819 
820  if (EQUAL(dv, dh)) {
821  ch++;
822  }
823 
824  cv++;
825  }
826  }
827  }
828 
829  if (!found) {
830  if (ch < hi) {
831  cpoint[X] = I3d[cnum][X] = Hi[ch][X];
832  cpoint[Y] = I3d[cnum][Y] = Hi[ch][Y];
833  I3d[cnum][Z] = Hi[ch][Z];
834  cnum++;
835  ch++;
836  }
837  }
838 
839  if (i == cv + ch + cd) {
840  G_debug(5, "order_intersects(): stuck on %d", cnum);
841  G_debug(5, "order_intersects(): cv = %d, ch = %d, cd = %d", cv,
842  ch, cd);
843  G_debug(5, "order_intersects(): dv = %f, dh = %f, dd = %f", dv,
844  dh, dd);
845 
846  break;
847  }
848  }
849 
850  if (EQUAL(last[X], cpoint[X]) && EQUAL(last[Y], cpoint[Y])) {
851  return (cnum);
852  }
853 
854  if (in_vregion(gs, last)) {
855  /* TODO: check for last point on corner ? */
856  I3d[cnum][X] = last[X];
857  I3d[cnum][Y] = last[Y];
858  I3d[cnum][Z] = last[Z];
859  ++cnum;
860  }
861 
862  return (cnum);
863 }
864 
882 int get_vert_intersects(geosurf * gs, float *bgn, float *end, float *dir)
883 {
884  int fcol, lcol, incr, hits, num, offset, drow1, drow2;
885  float xl, yb, xr, yt, z1, z2, alpha;
886  float xres, yres, xi, yi;
887  int bgncol, endcol, cols, rows;
888 
889  xres = VXRES(gs);
890  yres = VYRES(gs);
891  cols = VCOLS(gs);
892  rows = VROWS(gs);
893 
894  bgncol = X2VCOL(gs, bgn[X]);
895  endcol = X2VCOL(gs, end[X]);
896 
897  if (bgncol > cols && endcol > cols) {
898  return 0;
899  }
900 
901  if (bgncol == endcol) {
902  return 0;
903  }
904 
905  fcol = dir[X] > 0 ? bgncol + 1 : bgncol;
906  lcol = dir[X] > 0 ? endcol : endcol + 1;
907 
908  /* assuming only showing FULL cols */
909  incr = lcol - fcol > 0 ? 1 : -1;
910 
911  while (fcol > cols || fcol < 0) {
912  fcol += incr;
913  }
914 
915  while (lcol > cols || lcol < 0) {
916  lcol -= incr;
917  }
918 
919  num = abs(lcol - fcol) + 1;
920 
921  yb = gs->yrange - (yres * rows) - EPSILON;
922  yt = gs->yrange + EPSILON;
923 
924  for (hits = 0; hits < num; hits++) {
925  xl = xr = VCOL2X(gs, fcol);
926 
927  if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yt, xr, yb,
928  &xi, &yi)) {
929  Vi[hits][X] = xi;
930  Vi[hits][Y] = yi;
931 
932  /* find data rows */
933  if (Flat) {
934  Vi[hits][Z] = gs->att[ATT_TOPO].constant;
935  }
936  else {
937  drow1 = Y2VROW(gs, Vi[hits][Y]) * gs->y_mod;
938  drow2 = (1 + Y2VROW(gs, Vi[hits][Y])) * gs->y_mod;
939 
940  if (drow2 >= gs->rows) {
941  drow2 = gs->rows - 1; /*bottom edge */
942  }
943 
944  alpha =
945  ((gs->yrange - drow1 * gs->yres) - Vi[hits][Y]) / yres;
946 
947  offset = DRC2OFF(gs, drow1, fcol * gs->x_mod);
948  GET_MAPATT(Ebuf, offset, z1);
949  offset = DRC2OFF(gs, drow2, fcol * gs->x_mod);
950  GET_MAPATT(Ebuf, offset, z2);
951  Vi[hits][Z] = LERP(alpha, z1, z2);
952  }
953  }
954 
955  /* if they don't intersect, something's wrong! */
956  /* should only happen on endpoint, so it will be added later */
957  else {
958  hits--;
959  num--;
960  }
961 
962  fcol += incr;
963  }
964 
965  return (hits);
966 }
967 
978 int get_horz_intersects(geosurf * gs, float *bgn, float *end, float *dir)
979 {
980  int frow, lrow, incr, hits, num, offset, dcol1, dcol2;
981  float xl, yb, xr, yt, z1, z2, alpha;
982  float xres, yres, xi, yi;
983  int bgnrow, endrow, rows, cols;
984 
985  xres = VXRES(gs);
986  yres = VYRES(gs);
987  cols = VCOLS(gs);
988  rows = VROWS(gs);
989 
990  bgnrow = Y2VROW(gs, bgn[Y]);
991  endrow = Y2VROW(gs, end[Y]);
992  if (bgnrow == endrow) {
993  return 0;
994  }
995 
996  if (bgnrow > rows && endrow > rows) {
997  return 0;
998  }
999 
1000  frow = dir[Y] > 0 ? bgnrow : bgnrow + 1;
1001  lrow = dir[Y] > 0 ? endrow + 1 : endrow;
1002 
1003  /* assuming only showing FULL rows */
1004  incr = lrow - frow > 0 ? 1 : -1;
1005 
1006  while (frow > rows || frow < 0) {
1007  frow += incr;
1008  }
1009 
1010  while (lrow > rows || lrow < 0) {
1011  lrow -= incr;
1012  }
1013 
1014  num = abs(lrow - frow) + 1;
1015 
1016  xl = 0.0 - EPSILON;
1017  xr = xres * cols + EPSILON;
1018 
1019  for (hits = 0; hits < num; hits++) {
1020  yb = yt = VROW2Y(gs, frow);
1021 
1022  if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yt, xr, yb,
1023  &xi, &yi)) {
1024  Hi[hits][X] = xi;
1025  Hi[hits][Y] = yi;
1026 
1027  /* find data cols */
1028  if (Flat) {
1029  Hi[hits][Z] = gs->att[ATT_TOPO].constant;
1030  }
1031  else {
1032  dcol1 = X2VCOL(gs, Hi[hits][X]) * gs->x_mod;
1033  dcol2 = (1 + X2VCOL(gs, Hi[hits][X])) * gs->x_mod;
1034 
1035  if (dcol2 >= gs->cols) {
1036  dcol2 = gs->cols - 1; /* right edge */
1037  }
1038 
1039  alpha = (Hi[hits][X] - (dcol1 * gs->xres)) / xres;
1040 
1041  offset = DRC2OFF(gs, frow * gs->y_mod, dcol1);
1042  GET_MAPATT(Ebuf, offset, z1);
1043  offset = DRC2OFF(gs, frow * gs->y_mod, dcol2);
1044  GET_MAPATT(Ebuf, offset, z2);
1045  Hi[hits][Z] = LERP(alpha, z1, z2);
1046  }
1047  }
1048 
1049  /* if they don't intersect, something's wrong! */
1050  /* should only happen on endpoint, so it will be added later */
1051  else {
1052  hits--;
1053  num--;
1054  }
1055 
1056  frow += incr;
1057  }
1058 
1059  return (hits);
1060 }
1061 
1074 int get_diag_intersects(geosurf * gs, float *bgn, float *end, float *dir)
1075 {
1076  int fdig, ldig, incr, hits, num, offset;
1077  int vrow, vcol, drow1, drow2, dcol1, dcol2;
1078  float xl, yb, xr, yt, z1, z2, alpha;
1079  float xres, yres, xi, yi, dx, dy;
1080  int diags, cols, rows, lower;
1081  Point3 pt;
1082 
1083  xres = VXRES(gs);
1084  yres = VYRES(gs);
1085  cols = VCOLS(gs);
1086  rows = VROWS(gs);
1087  diags = rows + cols; /* -1 ? */
1088 
1089  /* determine upper/lower triangle for last */
1090  vrow = Y2VROW(gs, end[Y]);
1091  vcol = X2VCOL(gs, end[X]);
1092  pt[X] = VCOL2X(gs, vcol);
1093  pt[Y] = VROW2Y(gs, vrow + 1);
1094  lower = ((end[X] - pt[X]) / xres > (end[Y] - pt[Y]) / yres);
1095  ldig = lower ? vrow + vcol + 1 : vrow + vcol;
1096 
1097  /* determine upper/lower triangle for first */
1098  vrow = Y2VROW(gs, bgn[Y]);
1099  vcol = X2VCOL(gs, bgn[X]);
1100  pt[X] = VCOL2X(gs, vcol);
1101  pt[Y] = VROW2Y(gs, vrow + 1);
1102  lower = ((bgn[X] - pt[X]) / xres > (bgn[Y] - pt[Y]) / yres);
1103  fdig = lower ? vrow + vcol + 1 : vrow + vcol;
1104 
1105  /* adjust according to direction */
1106  if (ldig > fdig) {
1107  fdig++;
1108  }
1109 
1110  if (fdig > ldig) {
1111  ldig++;
1112  }
1113 
1114  incr = ldig - fdig > 0 ? 1 : -1;
1115 
1116  while (fdig > diags || fdig < 0) {
1117  fdig += incr;
1118  }
1119 
1120  while (ldig > diags || ldig < 0) {
1121  ldig -= incr;
1122  }
1123 
1124  num = abs(ldig - fdig) + 1;
1125 
1126  for (hits = 0; hits < num; hits++) {
1127  yb = gs->yrange - (yres * (fdig < rows ? fdig : rows)) - EPSILON;
1128  xl = VCOL2X(gs, (fdig < rows ? 0 : fdig - rows)) - EPSILON;
1129  yt = gs->yrange - (yres * (fdig < cols ? 0 : fdig - cols)) + EPSILON;
1130  xr = VCOL2X(gs, (fdig < cols ? fdig : cols)) + EPSILON;
1131 
1132  if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yb, xr, yt,
1133  &xi, &yi)) {
1134  Di[hits][X] = xi;
1135  Di[hits][Y] = yi;
1136 
1137  if (ISNODE(xi, xres)) {
1138  /* then it's also a ynode */
1139  num--;
1140  hits--;
1141  continue;
1142  }
1143 
1144  /* find data rows */
1145  drow1 = Y2VROW(gs, Di[hits][Y]) * gs->y_mod;
1146  drow2 = (1 + Y2VROW(gs, Di[hits][Y])) * gs->y_mod;
1147 
1148  if (drow2 >= gs->rows) {
1149  drow2 = gs->rows - 1; /* bottom edge */
1150  }
1151 
1152  /* find data cols */
1153  if (Flat) {
1154  Di[hits][Z] = gs->att[ATT_TOPO].constant;
1155  }
1156  else {
1157  dcol1 = X2VCOL(gs, Di[hits][X]) * gs->x_mod;
1158  dcol2 = (1 + X2VCOL(gs, Di[hits][X])) * gs->x_mod;
1159 
1160  if (dcol2 >= gs->cols) {
1161  dcol2 = gs->cols - 1; /* right edge */
1162  }
1163 
1164  dx = DCOL2X(gs, dcol2) - Di[hits][X];
1165  dy = DROW2Y(gs, drow1) - Di[hits][Y];
1166  alpha =
1167  sqrt(dx * dx + dy * dy) / sqrt(xres * xres + yres * yres);
1168 
1169  offset = DRC2OFF(gs, drow1, dcol2);
1170  GET_MAPATT(Ebuf, offset, z1);
1171  offset = DRC2OFF(gs, drow2, dcol1);
1172  GET_MAPATT(Ebuf, offset, z2);
1173  Di[hits][Z] = LERP(alpha, z1, z2);
1174  }
1175  }
1176 
1177  /* if they don't intersect, something's wrong! */
1178  /* should only happen on endpoint, so it will be added later */
1179  else {
1180  hits--;
1181  num--;
1182  }
1183 
1184  fdig += incr;
1185  }
1186 
1187  return (hits);
1188 }
1189 
1210 int segs_intersect(float x1, float y1, float x2, float y2, float x3, float y3,
1211  float x4, float y4, float *x, float *y)
1212 {
1213  float a1, a2, b1, b2, c1, c2; /* Coefficients of line eqns. */
1214  float r1, r2, r3, r4; /* 'Sign' values */
1215  float denom, offset, num; /* Intermediate values */
1216 
1217  /* Compute a1, b1, c1, where line joining points 1 and 2
1218  * is "a1 x + b1 y + c1 = 0".
1219  */
1220  a1 = y2 - y1;
1221  b1 = x1 - x2;
1222  c1 = x2 * y1 - x1 * y2;
1223 
1224  /* Compute r3 and r4.
1225  */
1226  r3 = a1 * x3 + b1 * y3 + c1;
1227  r4 = a1 * x4 + b1 * y4 + c1;
1228 
1229  /* Check signs of r3 and r4. If both point 3 and point 4 lie on
1230  * same side of line 1, the line segments do not intersect.
1231  */
1232 
1233  if (!EQUAL(r3, 0.0) && !EQUAL(r4, 0.0) && SAME_SIGNS(r3, r4)) {
1234  return (DONT_INTERSECT);
1235  }
1236 
1237  /* Compute a2, b2, c2 */
1238  a2 = y4 - y3;
1239  b2 = x3 - x4;
1240  c2 = x4 * y3 - x3 * y4;
1241 
1242  /* Compute r1 and r2 */
1243  r1 = a2 * x1 + b2 * y1 + c2;
1244  r2 = a2 * x2 + b2 * y2 + c2;
1245 
1246  /* Check signs of r1 and r2. If both point 1 and point 2 lie
1247  * on same side of second line segment, the line segments do
1248  * not intersect.
1249  */
1250 
1251  if (!EQUAL(r1, 0.0) && !EQUAL(r2, 0.0) && SAME_SIGNS(r1, r2)) {
1252  return (DONT_INTERSECT);
1253  }
1254 
1255  /* Line segments intersect: compute intersection point.
1256  */
1257  denom = a1 * b2 - a2 * b1;
1258 
1259  if (denom == 0) {
1260  return (COLLINEAR);
1261  }
1262 
1263  offset = denom < 0 ? -denom / 2 : denom / 2;
1264 
1265  /* The denom/2 is to get rounding instead of truncating. It
1266  * is added or subtracted to the numerator, depending upon the
1267  * sign of the numerator.
1268  */
1269  num = b1 * c2 - b2 * c1;
1270 
1271  *x = num / denom;
1272 
1273  num = a2 * c1 - a1 * c2;
1274  *y = num / denom;
1275 
1276  return (DO_INTERSECT);
1277 }
1278 
1290 int Point_on_plane(Point3 p1, Point3 p2, Point3 p3, Point3 unk)
1291 {
1292  float plane[4];
1293 
1294  P3toPlane(p1, p2, p3, plane);
1295 
1296  return (XY_intersect_plane(unk, plane));
1297 }
1298 
1312 int XY_intersect_plane(float *intersect, float *plane)
1313 {
1314  float x, y;
1315 
1316  if (!plane[Z]) {
1317  return (0); /* doesn't intersect */
1318  }
1319 
1320  x = intersect[X];
1321  y = intersect[Y];
1322  intersect[Z] = (plane[X] * x + plane[Y] * y + plane[W]) / -plane[Z];
1323 
1324  return (1);
1325 }
1326 
1335 int P3toPlane(Point3 p1, Point3 p2, Point3 p3, float *plane)
1336 {
1337  Point3 v1, v2, norm;
1338 
1339  v1[X] = p1[X] - p3[X];
1340  v1[Y] = p1[Y] - p3[Y];
1341  v1[Z] = p1[Z] - p3[Z];
1342 
1343  v2[X] = p2[X] - p3[X];
1344  v2[Y] = p2[Y] - p3[Y];
1345  v2[Z] = p2[Z] - p3[Z];
1346 
1347  V3Cross(v1, v2, norm);
1348 
1349  plane[X] = norm[X];
1350  plane[Y] = norm[Y];
1351  plane[Z] = norm[Z];
1352  plane[W] = -p3[X] * norm[X] - p3[Y] * norm[Y] - p3[Z] * norm[Z];
1353 
1354  return (1);
1355 }
1356 
1357 
1365 int V3Cross(Point3 a, Point3 b, Point3 c)
1366 {
1367  c[X] = (a[Y] * b[Z]) - (a[Z] * b[Y]);
1368  c[Y] = (a[Z] * b[X]) - (a[X] * b[Z]);
1369  c[Z] = (a[X] * b[Y]) - (a[Y] * b[X]);
1370 
1371  return (1);
1372 }
#define VXRES(gs)
Definition: rowcol.h:7
void G_free(void *buf)
Free allocated memory.
Definition: gis/alloc.c:142
int l
Definition: dataquad.c:292
float b
Definition: named_colr.c:8
double xmax
Definition: dataquad.c:293
#define DRC2OFF(gs, drow, dcol)
Definition: rowcol.h:15
#define VCOL2X(gs, vcol)
Definition: rowcol.h:38
int XY_intersect_plane(float *intersect, float *plane)
Check for intersection (point and plane)
Definition: gsdrape.c:1312
int in_vregion(geosurf *gs, float *pt)
ADD.
Definition: gsdrape.c:696
int gsdrape_set_surface(geosurf *gs)
ADD.
Definition: gsdrape.c:199
void interp_first_last(geosurf *gs, float *bgn, float *end, Point3 f, Point3 l)
ADD.
Definition: gsdrape.c:441
int gs_point_is_masked(geosurf *gs, float *pt)
Check if point is masked.
Definition: gs.c:1317
#define VROWS(gs)
Definition: rowcol.h:11
#define VCOLS(gs)
Definition: rowcol.h:12
int get_diag_intersects(geosurf *gs, float *bgn, float *end, float *dir)
Get diagonal intersects.
Definition: gsdrape.c:1074
#define EPSILON
Definition: blas_level_2.c:28
long num
Definition: g3dcats.c:93
#define Y(x)
Definition: display/draw.c:246
int viewcell_tri_interp(geosurf *gs, typbuff *buf, Point3 pt, int check_mask)
ADD.
Definition: gsdrape.c:507
int seg_intersect_vregion(geosurf *gs, float *bgn, float *end)
Check if segment intersect vector region.
Definition: gsdrape.c:233
#define EQUAL
Definition: y.tab.c:159
#define X(y)
Definition: display/draw.c:248
double ymin
Definition: dataquad.c:293
int y
Definition: plot.c:34
int get_vert_intersects(geosurf *gs, float *bgn, float *end, float *dir)
ADD.
Definition: gsdrape.c:882
void GS_v2dir(float *v1, float *v2, float *v3)
Get a normalized direction from v1 to v2, store in v3 (2D)
Definition: GS_util.c:385
#define VYRES(gs)
Definition: rowcol.h:8
#define VROW2Y(gs, vrow)
Definition: rowcol.h:37
#define DONT_INTERSECT
Definition: gsdrape.c:59
typbuff * gs_get_att_typbuff(geosurf *gs, int desc, int to_write)
Get attribute data buffer.
Definition: gs.c:681
#define Y2VROW(gs, py)
Definition: rowcol.h:25
#define X2VCOL(gs, px)
Definition: rowcol.h:26
#define DCOL2X(gs, dcol)
Definition: rowcol.h:34
#define ISNODE(p, res)
Definition: gsdrape.c:65
Point3 * gsdrape_get_allsegments(geosurf *gs, float *bgn, float *end, int *num)
Get all segments.
Definition: gsdrape.c:401
int get_horz_intersects(geosurf *gs, float *bgn, float *end, float *dir)
Get horizontal intersects.
Definition: gsdrape.c:978
#define LERP(a, l, h)
Definition: gsdrape.c:63
#define SAME_SIGNS(a, b)
Definition: gsdrape.c:67
double ymax
Definition: dataquad.c:293
int first
Definition: form/open.c:25
int Point_on_plane(Point3 p1, Point3 p2, Point3 p3, Point3 unk)
Check if point is on plane.
Definition: gsdrape.c:1290
if(!YY_CURRENT_BUFFER)
Definition: lex.yy.c:799
int _viewcell_tri_interp(geosurf *gs, Point3 pt)
ADD.
Definition: gsdrape.c:467
char buf[GNAME_MAX+sizeof(G3D_DIRECTORY)+2]
Definition: g3drange.c:62
int P3toPlane(Point3 p1, Point3 p2, Point3 p3, float *plane)
Define plane.
Definition: gsdrape.c:1335
return NULL
Definition: dbfopen.c:1394
Point3 * gsdrape_get_segments(geosurf *gs, float *bgn, float *end, int *num)
ADD.
Definition: gsdrape.c:350
void GS_v3eq(float *v1, float *v2)
Copy vector values.
Definition: GS_util.c:178
G_warning("category support for [%s] in mapset [%s] %s", name, mapset, type)
tuple cols
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: gis/debug.c:51
int gs_get_att_src(geosurf *gs, int desc)
Get attribute source.
Definition: gs.c:656
#define DROW2Y(gs, drow)
Definition: rowcol.h:33
#define VCOL2DCOL(gs, vcol)
Definition: rowcol.h:30
#define COLLINEAR
Definition: gsdrape.c:61
int order_intersects(geosurf *gs, Point3 first, Point3 last, int vi, int hi, int di)
ADD.
Definition: gsdrape.c:729
Definition: draw2.c:52
float GS_P2distance(float *from, float *to)
Calculate distance in plane.
Definition: GS_util.c:160
#define VROW2DROW(gs, vrow)
Definition: rowcol.h:29
int V3Cross(Point3 a, Point3 b, Point3 c)
Get cross product.
Definition: gsdrape.c:1365
#define DO_INTERSECT
Definition: gsdrape.c:60
int segs_intersect(float x1, float y1, float x2, float y2, float x3, float y3, float x4, float y4, float *x, float *y)
Line intersect.
Definition: gsdrape.c:1210
#define GET_MAPATT(buff, offset, att)
Definition: gsget.h:27