GRASS GIS 7 Programmer's Manual  7.7.svn(2018)-r73679
gsdrape.c
Go to the documentation of this file.
1 /*!
2  \file lib/ogsf/gsdrape.c
3
4  \brief OGSF library - functions to intersect line segments with edges of surface polygons
5
6  GRASS OpenGL gsurf OGSF Library
7
8  For efficiency, intersections are found without respect to which
9  specific triangle edge is intersected, but on a broader sense with
10  the horizontal, vertical, and diagonal seams in the grid, then
11  the intersections are ordered. If quadstrips are used for drawing
12  rather than tmesh, triangulation is not consistant; for diagonal
13  intersections, the proper diagonal to intersect would need to be
14  determined according to the algorithm used by qstrip (look at nearby
15  normals). It may be faster to go ahead and find the intersections
16  with the other diagonals using the same methods, then at sorting
17  time determine which diagonal array to look at for each quad.
18  It would also require a mechanism for throwing out unused intersections
19  with the diagonals during the ordering phase.
20  Do intersections in 2D, fill line structure with 3D pts (maybe calling
21  routine will cache for redrawing). Get Z value by using linear interp
22  between corners.
23
24  - check for easy cases:
25  - single point
26  - colinear with horizontal or vertical edges
27  - colinear with diagonal edges of triangles
28  - calculate three arrays of ordered intersections:
29  - with vertical edges
30  - with horizontal edges
31  - with diagonal edges and interpolate Z, using simple linear interpolation.
32  - eliminate duplicate intersections (need only compare one coord for each)
33  - build ordered set of points.
34
35  Return static pointer to 3D set of points. Max number of intersections
36  will be rows + cols + diags, so should allocate this number to initialize.
37  Let calling routine worry about copying points for caching.
38
39  (C) 1999-2008 by the GRASS Development Team
40
41  This program is free software under the
42  GNU General Public License (>=v2).
43  Read the file COPYING that comes with GRASS
44  for details.
45
46  \author Bill Brown UI GMS Lab
47  \author Doxygenized by Martin Landa <landa.martin gmail.com> (May 2008)
48  */
49
50 #include <stdlib.h>
51
52 #include <grass/ogsf.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
87 /*!
88  \brief Initizalize
89
90  \param rows number of rows
91  \param cols number of columns
92
93  \return -1 on failure
94  \return 1 on success
95  */
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
127 /*!
128  \brief Get segments
129
130  \param gs surface (geosurf)
131  \param bgn begin point
132  \param end end point
133  \param num
134
135  \return pointer to Point3 struct
136  */
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
173 /*!
174  \brief Calculate 2D distance
175
176  \param p1 first point
177  \param p2 second point
178
179  \return distance
180  */
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
191 /*!
193
194  \param gs surface (geosurf)
195
196  \return -1 on failure
197  \return 1 on success
198  */
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
219 /*!
220  \brief Check if segment intersect vector region
221
222  Clipping performed:
223  - bgn and end are replaced so that both points are within viewregion
224  - if seg intersects
225
226  \param gs surface (geosurf)
227  \param bgn begin point
228  \param end end point
229
230  \return 0 if segment doesn't intersect the viewregion, or intersects only at corner
231  \return otherwise returns 1
232  */
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
340 /*!
342
343  \param gs surface (geosurf)
344  \param bgn begin point (x,y)
345  \param end end point (x,y)
346  \param num
347
348  \return pointer to Point3 struct
349  */
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
391 /*!
392  \brief Get all segments
393
394  \param gs surface (geosurf)
395  \param bgn begin point
396  \param end end point
397  \param num
398
399  \return pointer to Point3 struct
400  */
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
432 /*!
434
435  \param gs surface (geosurf)
436  \param bgn begin point
437  \param end end point
438  \param f first
439  \param l last
440  */
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
461 /*!
463
464  \param gs surface (geosurf)
465  \param pt
466  */
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
476 /*!
478
479  In gsd_surf, tmesh draws polys like so:
480  <pre>
481  --------------
482  | /|
483  | / |
484  | / |
485  | / |
486  | / |
487  | / |
488  | / |
489  | / |
490  | / |
491  | / |
492  | / |
493  |/ |
494  --------------
495  </pre>
496
497  UNLESS the top right or bottom left point is masked, in which case a
498  single triangle with the opposite diagonal is drawn. This case is
499  not yet handled here & should only occur on edges.
500  pt has X & Y coordinates in it, we interpolate Z here
501
502  This could probably be much shorter, but not much faster.
503
504  \return 1 if point is in view region
505  \return otherwise 0 (if masked)
506  */
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
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
688 /*!
690
691  \param gs surface (geosurf)
692
693  \return 1
694  \return 0
695  */
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
707 /*!
709
710  After all the intersections between the segment and triangle
711  edges have been found, they are in three lists. (intersections
712  with vertical, horizontal, and diagonal triangle edges)
713
714  Each list is ordered in space from first to last segment points,
715  but now the lists need to be woven together. This routine
716  starts with the first point of the segment and then checks the
717  next point in each list to find the closest, eliminating duplicates
718  along the way and storing the result in I3d.
719
720  \param gs surface (geosurf)
721  \param first first point
722  \param last last point
723  \param vi
724  \param hi
725  \param di
726
727  \return
728  */
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
865 /*!
867
868  \todo For consistancy, need to decide how last row & last column are
869  displayed - would it look funny to always draw last row/col with
870  finer resolution if necessary, or would it be better to only show
871  full rows/cols?
872
874
875  \param gs surface (geosurf)
876  \param bgn begin point
877  \param end end point
878  \param dir direction
879
880  \return
881  */
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
968 /*!
969  \brief Get horizontal intersects
970
971  \param gs surface (geosurf)
972  \param bgn begin point
973  \param end end point
974  \param dir
975
976  \return number of intersects
977  */
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
1062 /*!
1063  \brief Get diagonal intersects
1064
1066
1067  \param gs surface (geosurf)
1068  \param bgn begin point
1069  \param end end point
1070  \param dir ? (unused)
1071
1072  \return number of intersects
1073  */
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
1190 /*!
1191  \brief Line intersect
1192
1194  Modified for floating point: Bill Brown
1195
1196  This function computes whether two line segments,
1197  respectively joining the input points (x1,y1) -- (x2,y2)
1198  and the input points (x3,y3) -- (x4,y4) intersect.
1199  If the lines intersect, the output variables x, y are
1200  set to coordinates of the point of intersection.
1201
1202  \param x1,y1,x2,y2 coordinates of endpoints of one segment
1203  \param x3,y3,x4,y4 coordinates of endpoints of other segment
1204  \param[out] x,y coordinates of intersection point
1205
1206  \return 0 no intersection
1207  \return 1 intersect
1208  \return 2 collinear
1209  */
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
1279 /*!
1280  \brief Check if point is on plane
1281
1282  Plane defined by three points here; user fills in unk[X] & unk[Y]
1283
1284  \param p1,p2,p3 points defining plane
1285  \param unk point
1286
1287  \return 1 point on plane
1288  \return 0 point not on plane
1289  */
1291 {
1292  float plane[4];
1293
1294  P3toPlane(p1, p2, p3, plane);
1295
1296  return (XY_intersect_plane(unk, plane));
1297 }
1298
1299 /*!
1300  \brief Check for intersection (point and plane)
1301
1302  Ax + By + Cz + D = 0, so z = (Ax + By + D) / -C
1303
1304  User fills in intersect[X] & intersect[Y]
1305
1306  \param[out] intersect intersect coordinates
1307  \param plane plane definition
1308
1309  \return 0 doesn't intersect
1310  \return 1 intesects
1311  */
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
1327 /*!
1328  \brief Define plane
1329
1330  \param p1,p2,p3 three point on plane
1331  \param[out] plane plane defintion
1332
1333  \return 1
1334  */
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
1358 /*!
1359  \brief Get cross product
1360
1361  \param a,b,c
1362
1363  \return cross product c = a cross b
1364  */
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:9
float xrange
Definition: ogsf.h:269
void G_free(void *buf)
Free allocated memory.
Definition: gis/alloc.c:149
#define DRC2OFF(gs, drow, dcol)
Definition: rowcol.h:17
#define VCOL2X(gs, vcol)
Definition: rowcol.h:40
if(!YY_CURRENT_BUFFER)
Definition: sqlp.yy.c:803
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)
Definition: gsdrape.c:696
int gsdrape_set_surface(geosurf *gs)
Definition: gsdrape.c:199
void interp_first_last(geosurf *gs, float *bgn, float *end, Point3 f, Point3 l)
Definition: gsdrape.c:441
Definition: gs.c:1317
#define VROWS(gs)
Definition: rowcol.h:13
int rows
Definition: ogsf.h:260
#define VCOLS(gs)
Definition: rowcol.h:14
#define ATT_TOPO
Definition: ogsf.h:73
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:27
int viewcell_tri_interp(geosurf *gs, typbuff *buf, Point3 pt, int check_mask)
Definition: gsdrape.c:507
int seg_intersect_vregion(geosurf *gs, float *bgn, float *end)
Check if segment intersect vector region.
Definition: gsdrape.c:233
float Point3[3]
Definition: ogsf.h:201
#define NULL
Definition: ccmath.h:32
#define x
int get_vert_intersects(geosurf *gs, float *bgn, float *end, float *dir)
Definition: gsdrape.c:882
#define VYRES(gs)
Definition: rowcol.h:10
#define W
Definition: ogsf.h:140
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 VROW2Y(gs, vrow)
Definition: rowcol.h:39
int cols
Definition: ogsf.h:260
double l
Definition: r_raster.c:39
#define DONT_INTERSECT
Definition: gsdrape.c:59
float yrange
Definition: ogsf.h:269
double yres
Definition: ogsf.h:265
typbuff * gs_get_att_typbuff(geosurf *gs, int desc, int to_write)
Get attribute data buffer.
Definition: gs.c:681
int x_mod
Definition: ogsf.h:271
double b
Definition: r_raster.c:39
#define Y2VROW(gs, py)
Definition: rowcol.h:27
#define X2VCOL(gs, px)
Definition: rowcol.h:28
#define DCOL2X(gs, dcol)
Definition: rowcol.h:36
#define ISNODE(p, res)
Definition: gsdrape.c:65
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: debug.c:65
Point3 * gsdrape_get_allsegments(geosurf *gs, float *bgn, float *end, int *num)
Get all segments.
Definition: gsdrape.c:401
#define MAP_ATT
Definition: ogsf.h:83
int y_mod
Definition: ogsf.h:271
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
#define Z
Definition: ogsf.h:139
int Point_on_plane(Point3 p1, Point3 p2, Point3 p3, Point3 unk)
Check if point is on plane.
Definition: gsdrape.c:1290
gsurf_att att[MAX_ATTS]
Definition: ogsf.h:261
int _viewcell_tri_interp(geosurf *gs, Point3 pt)
Definition: gsdrape.c:467
#define Y
Definition: ogsf.h:138
int P3toPlane(Point3 p1, Point3 p2, Point3 p3, float *plane)
Define plane.
Definition: gsdrape.c:1335
Point3 * gsdrape_get_segments(geosurf *gs, float *bgn, float *end, int *num)
Definition: gsdrape.c:350
Definition: ogsf.h:204
void GS_v3eq(float *v1, float *v2)
Copy vector values.
Definition: gs_util.c:178
int gs_get_att_src(geosurf *gs, int desc)
Get attribute source.
Definition: gs.c:656
#define DROW2Y(gs, drow)
Definition: rowcol.h:35
#define VCOL2DCOL(gs, vcol)
Definition: rowcol.h:32
#define _(str)
Definition: glocale.h:13
#define X
Definition: ogsf.h:137
long num
Definition: raster3d/cats.c:85
#define COLLINEAR
Definition: gsdrape.c:61
int order_intersects(geosurf *gs, Point3 first, Point3 last, int vi, int hi, int di)
Definition: gsdrape.c:729
Definition: clip.h:7
#define VROW2DROW(gs, vrow)
Definition: rowcol.h:31
double xres
Definition: ogsf.h:265
int V3Cross(Point3 a, Point3 b, Point3 c)
Get cross product.
Definition: gsdrape.c:1365
#define CONST_ATT
Definition: ogsf.h:84
#define DO_INTERSECT
Definition: gsdrape.c:60
float GS_P2distance(float *from, float *to)
Calculate distance in plane.
Definition: gs_util.c:160
Definition: ogsf.h:257
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition: gis/error.c:204
float constant
Definition: ogsf.h:251
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 EQUAL
Definition: sqlp.tab.c:159
#define GET_MAPATT(buff, offset, att)
Definition: gsget.h:27