GRASS GIS 8 Programmer's Manual  8.5.0dev(2024)-f0497c5f6c
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
5  surface polygons
6
7  GRASS OpenGL gsurf OGSF Library
8
9  For efficiency, intersections are found without respect to which
10  specific triangle edge is intersected, but on a broader sense with
11  the horizontal, vertical, and diagonal seams in the grid, then
12  the intersections are ordered. If quadstrips are used for drawing
13  rather than tmesh, triangulation is not consistent; for diagonal
14  intersections, the proper diagonal to intersect would need to be
15  determined according to the algorithm used by qstrip (look at nearby
16  normals). It may be faster to go ahead and find the intersections
17  with the other diagonals using the same methods, then at sorting
18  time determine which diagonal array to look at for each quad.
19  It would also require a mechanism for throwing out unused intersections
20  with the diagonals during the ordering phase.
21  Do intersections in 2D, fill line structure with 3D pts (maybe calling
22  routine will cache for redrawing). Get Z value by using linear interp
23  between corners.
24
25  - check for easy cases:
26  - single point
27  - colinear with horizontal or vertical edges
28  - colinear with diagonal edges of triangles
29  - calculate three arrays of ordered intersections:
30  - with vertical edges
31  - with horizontal edges
32  - with diagonal edges and interpolate Z, using simple linear interpolation.
33  - eliminate duplicate intersections (need only compare one coord for each)
34  - build ordered set of points.
35
36  Return static pointer to 3D set of points. Max number of intersections
37  will be rows + cols + diags, so should allocate this number to initialize.
38  Let calling routine worry about copying points for caching.
39
40  (C) 1999-2008 by the GRASS Development Team
41
42  This program is free software under the
43  GNU General Public License (>=v2).
44  Read the file COPYING that comes with GRASS
45  for details.
46
47  \author Bill Brown UI GMS Lab
48  \author Doxygenized by Martin Landa <landa.martin gmail.com> (May 2008)
49  */
50
51 #include <stdlib.h>
52
53 #include <grass/ogsf.h>
54 #include <grass/glocale.h>
55
56 #include "gsget.h"
57 #include "rowcol.h"
58 #include "math.h"
59
60 #define DONT_INTERSECT 0
61 #define DO_INTERSECT 1
62 #define COLLINEAR 2
63
64 #define LERP(a, l, h) ((l) + (((h) - (l)) * (a)))
65 #define EQUAL(a, b) (fabs((a) - (b)) < EPSILON)
66 #define ISNODE(p, res) (fmod((double)p, (double)res) < EPSILON)
67
68 #define SAME_SIGNS(a, b) ((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  \brief Initizalize
88
89  \param rows number of rows
90  \param cols number of columns
91
92  \return -1 on failure
93  \return 1 on success
94  */
95 static int drape_line_init(int rows, int cols)
96 {
97  /* use G_calloc() [-> G_fatal_error] instead of calloc ? */
98  if (NULL == (I3d = (Point3 *)calloc(2 * (rows + cols), sizeof(Point3)))) {
99  return (-1);
100  }
101
102  if (NULL == (Vi = (Point3 *)calloc(cols, sizeof(Point3)))) {
103  G_free(I3d);
104
105  return (-1);
106  }
107
108  if (NULL == (Hi = (Point3 *)calloc(rows, sizeof(Point3)))) {
109  G_free(I3d);
110  G_free(Vi);
111
112  return (-1);
113  }
114
115  if (NULL == (Di = (Point3 *)calloc(rows + cols, sizeof(Point3)))) {
116  G_free(I3d);
117  G_free(Vi);
118  G_free(Hi);
119
120  return (-1);
121  }
122
123  return (1);
124 }
125
126 /*!
127  \brief Get segments
128
129  \param gs surface (geosurf)
130  \param bgn begin point
131  \param end end point
132  \param num
133
134  \return pointer to Point3 struct
135  */
136 static Point3 *_gsdrape_get_segments(geosurf *gs, float *bgn, float *end,
137  int *num)
138 {
139  float f[3], l[3];
140  int vi, hi, di;
141  float dir[2], yres, xres;
142
143  xres = VXRES(gs);
144  yres = VYRES(gs);
145
146  vi = hi = di = 0;
147  GS_v2dir(bgn, end, dir);
148
149  if (dir[X]) {
150  vi = get_vert_intersects(gs, bgn, end, dir);
151  }
152
153  if (dir[Y]) {
154  hi = get_horz_intersects(gs, bgn, end, dir);
155  }
156
157  if (!((end[Y] - bgn[Y]) / (end[X] - bgn[X]) == yres / xres)) {
158  di = get_diag_intersects(gs, bgn, end, dir);
159  }
160
161  interp_first_last(gs, bgn, end, f, l);
162
163  *num = order_intersects(gs, f, l, vi, hi, di);
164  /* fills in return values, eliminates dupes (corners) */
165
166  G_debug(5, "_gsdrape_get_segments vi=%d, hi=%d, di=%d, num=%d", vi, hi, di,
167  *num);
168
169  return (I3d);
170 }
171
172 /*!
173  \brief Calculate 2D distance
174
175  \param p1 first point
176  \param p2 second point
177
178  \return distance
179  */
180 static float dist_squared_2d(float *p1, float *p2)
181 {
182  float dx, dy;
183
184  dx = p2[X] - p1[X];
185  dy = p2[Y] - p1[Y];
186
187  return (dx * dx + dy * dy);
188 }
189
190 /*!
192
193  \param gs surface (geosurf)
194
195  \return -1 on failure
196  \return 1 on success
197  */
199 {
200  static int first = 1;
201
202  if (first) {
203  first = 0;
204
205  if (0 > drape_line_init(gs->rows, gs->cols)) {
206  G_warning(_("Unable to process vector map - out of memory"));
207  Ebuf = NULL;
208
209  return (-1);
210  }
211  }
212
213  Ebuf = gs_get_att_typbuff(gs, ATT_TOPO, 0);
214
215  return (1);
216 }
217
218 /*!
219  \brief Check if segment intersect vector region
220
221  Clipping performed:
222  - bgn and end are replaced so that both points are within viewregion
223  - if seg intersects
224
225  \param gs surface (geosurf)
226  \param bgn begin point
227  \param end end point
228
229  \return 0 if segment doesn't intersect the viewregion, or intersects only at
230  corner \return otherwise returns 1
231  */
232 int seg_intersect_vregion(geosurf *gs, float *bgn, float *end)
233 {
234  float *replace, xl, yb, xr, yt, xi, yi;
235  int inside = 0;
236
237  xl = 0.0;
238  xr = VCOL2X(gs, VCOLS(gs));
239  yt = VROW2Y(gs, 0);
240  yb = VROW2Y(gs, VROWS(gs));
241
242  if (in_vregion(gs, bgn)) {
243  replace = end;
244  inside++;
245  }
246
247  if (in_vregion(gs, end)) {
248  replace = bgn;
249  inside++;
250  }
251
252  if (inside == 2) {
253  return (1);
254  }
255  else if (inside) {
256  /* one in & one out - replace gets first intersection */
257  if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yb, xl, yt, &xi,
258  &yi)) {
259  /* left */
260  }
261  else if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xr, yb, xr, yt,
262  &xi, &yi)) {
263  /* right */
264  }
265  else if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yb, xr, yb,
266  &xi, &yi)) {
267  /* bottom */
268  }
269  else if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yt, xr, yt,
270  &xi, &yi)) {
271  /* top */
272  }
273
274  replace[X] = xi;
275  replace[Y] = yi;
276  }
277  else {
278  /* both out - find 2 intersects & replace both */
279  float pt1[2], pt2[2];
280
281  replace = pt1;
282  if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yb, xl, yt, &xi,
283  &yi)) {
284  replace[X] = xi;
285  replace[Y] = yi;
286  replace = pt2;
287  inside++;
288  }
289
290  if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xr, yb, xr, yt, &xi,
291  &yi)) {
292  replace[X] = xi;
293  replace[Y] = yi;
294  replace = pt2;
295  inside++;
296  }
297
298  if (inside < 2) {
299  if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yb, xr, yb,
300  &xi, &yi)) {
301  replace[X] = xi;
302  replace[Y] = yi;
303  replace = pt2;
304  inside++;
305  }
306  }
307
308  if (inside < 2) {
309  if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yt, xr, yt,
310  &xi, &yi)) {
311  replace[X] = xi;
312  replace[Y] = yi;
313  inside++;
314  }
315  }
316
317  if (inside < 2) {
318  return (0); /* no intersect or only 1 point on corner */
319  }
320
321  /* compare dist of intersects to bgn - closest replaces bgn */
322  if (GS_P2distance(bgn, pt1) < GS_P2distance(bgn, pt2)) {
323  bgn[X] = pt1[X];
324  bgn[Y] = pt1[Y];
325  end[X] = pt2[X];
326  end[Y] = pt2[Y];
327  }
328  else {
329  bgn[X] = pt2[X];
330  bgn[Y] = pt2[Y];
331  end[X] = pt1[X];
332  end[Y] = pt1[Y];
333  }
334  }
335
336  return (1);
337 }
338
339 /*!
341
342  \param gs surface (geosurf)
343  \param bgn begin point (x,y)
344  \param end end point (x,y)
345  \param num
346
347  \return pointer to Point3 struct
348  */
349 Point3 *gsdrape_get_segments(geosurf *gs, float *bgn, float *end, int *num)
350 {
352
353  if (!seg_intersect_vregion(gs, bgn, end)) {
354  *num = 0;
355
356  return (NULL);
357  }
358
359  if (CONST_ATT == gs_get_att_src(gs, ATT_TOPO)) {
360  /* will probably want a force_drape option to get all intersects */
361  I3d[0][X] = bgn[X];
362  I3d[0][Y] = bgn[Y];
363  I3d[0][Z] = gs->att[ATT_TOPO].constant;
364  I3d[1][X] = end[X];
365  I3d[1][Y] = end[Y];
366  I3d[1][Z] = gs->att[ATT_TOPO].constant;
367  *num = 2;
368
369  return (I3d);
370  }
371
372  if (bgn[X] == end[X] && bgn[Y] == end[Y]) {
373  float f[3], l[3];
374
375  interp_first_last(gs, bgn, end, f, l);
376  GS_v3eq(I3d[0], f);
377  GS_v3eq(I3d[1], l);
378
379  /* CHANGE (*num = 1) to reflect degenerate line ? */
380  *num = 2;
381
382  return (I3d);
383  }
384
385  Flat = 0;
386  return (_gsdrape_get_segments(gs, bgn, end, num));
387 }
388
389 /*!
390  \brief Get all segments
391
392  \param gs surface (geosurf)
393  \param bgn begin point
394  \param end end point
395  \param num
396
397  \return pointer to Point3 struct
398  */
399 Point3 *gsdrape_get_allsegments(geosurf *gs, float *bgn, float *end, int *num)
400 {
402
403  if (!seg_intersect_vregion(gs, bgn, end)) {
404  *num = 0;
405  return (NULL);
406  }
407
408  if (bgn[X] == end[X] && bgn[Y] == end[Y]) {
409  float f[3], l[3];
410
411  interp_first_last(gs, bgn, end, f, l);
412  GS_v3eq(I3d[0], f);
413  GS_v3eq(I3d[1], l);
414  *num = 2;
415
416  return (I3d);
417  }
418
419  if (CONST_ATT == gs_get_att_src(gs, ATT_TOPO)) {
420  Flat = 1;
421  }
422  else {
423  Flat = 0;
424  }
425
426  return (_gsdrape_get_segments(gs, bgn, end, num));
427 }
428
429 /*!
431
432  \param gs surface (geosurf)
433  \param bgn begin point
434  \param end end point
435  \param f first
436  \param l last
437  */
438 void interp_first_last(geosurf *gs, float *bgn, float *end, Point3 f, Point3 l)
439 {
440  f[X] = bgn[X];
441  f[Y] = bgn[Y];
442
443  l[X] = end[X];
444  l[Y] = end[Y];
445
446  if (Flat) {
447  f[Z] = l[Z] = gs->att[ATT_TOPO].constant;
448  }
449  else {
450  viewcell_tri_interp(gs, Ebuf, f, 0);
451  viewcell_tri_interp(gs, Ebuf, l, 0);
452  }
453
454  return;
455 }
456
457 /*!
459
460  \param gs surface (geosurf)
461  \param pt
462  */
464 {
465  typbuff *buf;
466
467  buf = gs_get_att_typbuff(gs, ATT_TOPO, 0);
468
469  return (viewcell_tri_interp(gs, buf, pt, 0));
470 }
471
472 /*!
474
475  In gsd_surf, tmesh draws polys like so:
476  <pre>
477  --------------
478  | /|
479  | / |
480  | / |
481  | / |
482  | / |
483  | / |
484  | / |
485  | / |
486  | / |
487  | / |
488  | / |
489  |/ |
490  --------------
491  </pre>
492
493  UNLESS the top right or bottom left point is masked, in which case a
494  single triangle with the opposite diagonal is drawn. This case is
495  not yet handled here & should only occur on edges.
496  pt has X & Y coordinates in it, we interpolate Z here
497
498  This could probably be much shorter, but not much faster.
499
500  \return 1 if point is in view region
501  \return otherwise 0 (if masked)
502  */
503 int viewcell_tri_interp(geosurf *gs, typbuff *buf, Point3 pt, int check_mask)
504 {
505  Point3 p1, p2, p3;
506  int offset, drow, dcol, vrow, vcol;
507  float xmax, ymin, ymax, alpha;
508
509  xmax = VCOL2X(gs, VCOLS(gs));
510  ymax = VROW2Y(gs, 0);
511  ymin = VROW2Y(gs, VROWS(gs));
512
515  return (0);
516  }
517  }
518
519  if (pt[X] < 0.0 || pt[Y] > ymax) {
520  /* outside on left or top */
521  return (0);
522  }
523
524  if (pt[Y] < ymin || pt[X] > xmax) {
525  /* outside on bottom or right */
526  return (0);
527  }
528
529  if (CONST_ATT == gs_get_att_src(gs, ATT_TOPO)) {
530  pt[Z] = gs->att[ATT_TOPO].constant;
531
532  return (1);
533  }
534  else if (MAP_ATT != gs_get_att_src(gs, ATT_TOPO)) {
535  return (0);
536  }
537
538  vrow = Y2VROW(gs, pt[Y]);
539  vcol = X2VCOL(gs, pt[X]);
540
541  if (vrow < VROWS(gs) && vcol < VCOLS(gs)) {
542  /*not on bottom or right edge */
543  if (pt[X] > 0.0 && pt[Y] < ymax) {
544  /* not on left or top edge */
545  p1[X] = VCOL2X(gs, vcol + 1);
546  p1[Y] = VROW2Y(gs, vrow);
547  drow = VROW2DROW(gs, vrow);
548  dcol = VCOL2DCOL(gs, vcol + 1);
549  offset = DRC2OFF(gs, drow, dcol);
550  GET_MAPATT(buf, offset, p1[Z]); /* top right */
551
552  p2[X] = VCOL2X(gs, vcol);
553  p2[Y] = VROW2Y(gs, vrow + 1);
554  drow = VROW2DROW(gs, vrow + 1);
555  dcol = VCOL2DCOL(gs, vcol);
556  offset = DRC2OFF(gs, drow, dcol);
557  GET_MAPATT(buf, offset, p2[Z]); /* bottom left */
558
559  if ((pt[X] - p2[X]) / VXRES(gs) > (pt[Y] - p2[Y]) / VYRES(gs)) {
560  /* lower triangle */
561  p3[X] = VCOL2X(gs, vcol + 1);
562  p3[Y] = VROW2Y(gs, vrow + 1);
563  drow = VROW2DROW(gs, vrow + 1);
564  dcol = VCOL2DCOL(gs, vcol + 1);
565  offset = DRC2OFF(gs, drow, dcol);
566  GET_MAPATT(buf, offset, p3[Z]); /* bottom right */
567  }
568  else {
569  /* upper triangle */
570  p3[X] = VCOL2X(gs, vcol);
571  p3[Y] = VROW2Y(gs, vrow);
572  drow = VROW2DROW(gs, vrow);
573  dcol = VCOL2DCOL(gs, vcol);
574  offset = DRC2OFF(gs, drow, dcol);
575  GET_MAPATT(buf, offset, p3[Z]); /* top left */
576  }
577
578  return (Point_on_plane(p1, p2, p3, pt));
579  }
580  else if (pt[X] == 0.0) {
581  /* on left edge */
582  if (pt[Y] < ymax) {
583  vrow = Y2VROW(gs, pt[Y]);
584  drow = VROW2DROW(gs, vrow);
585  offset = DRC2OFF(gs, drow, 0);
586  GET_MAPATT(buf, offset, p1[Z]);
587
588  drow = VROW2DROW(gs, vrow + 1);
589  offset = DRC2OFF(gs, drow, 0);
590  GET_MAPATT(buf, offset, p2[Z]);
591
592  alpha = (VROW2Y(gs, vrow) - pt[Y]) / VYRES(gs);
593  pt[Z] = LERP(alpha, p1[Z], p2[Z]);
594  }
595  else {
596  /* top left corner */
597  GET_MAPATT(buf, 0, pt[Z]);
598  }
599
600  return (1);
601  }
602  else if (pt[Y] == gs->yrange) {
603  /* on top edge, not a corner */
604  vcol = X2VCOL(gs, pt[X]);
605  dcol = VCOL2DCOL(gs, vcol);
606  GET_MAPATT(buf, dcol, p1[Z]);
607
608  dcol = VCOL2DCOL(gs, vcol + 1);
609  GET_MAPATT(buf, dcol, p2[Z]);
610
611  alpha = (pt[X] - VCOL2X(gs, vcol)) / VXRES(gs);
612  pt[Z] = LERP(alpha, p1[Z], p2[Z]);
613
614  return (1);
615  }
616  }
617  else if (vrow == VROWS(gs)) {
618  /* on bottom edge */
619  drow = VROW2DROW(gs, VROWS(gs));
620
621  if (pt[X] > 0.0 && pt[X] < xmax) {
622  /* not a corner */
623  vcol = X2VCOL(gs, pt[X]);
624  dcol = VCOL2DCOL(gs, vcol);
625  offset = DRC2OFF(gs, drow, dcol);
626  GET_MAPATT(buf, offset, p1[Z]);
627
628  dcol = VCOL2DCOL(gs, vcol + 1);
629  offset = DRC2OFF(gs, drow, dcol);
630  GET_MAPATT(buf, offset, p2[Z]);
631
632  alpha = (pt[X] - VCOL2X(gs, vcol)) / VXRES(gs);
633  pt[Z] = LERP(alpha, p1[Z], p2[Z]);
634
635  return (1);
636  }
637  else if (pt[X] == 0.0) {
638  /* bottom left corner */
639  offset = DRC2OFF(gs, drow, 0);
640  GET_MAPATT(buf, offset, pt[Z]);
641
642  return (1);
643  }
644  else {
645  /* bottom right corner */
646  dcol = VCOL2DCOL(gs, VCOLS(gs));
647  offset = DRC2OFF(gs, drow, dcol);
648  GET_MAPATT(buf, offset, pt[Z]);
649
650  return (1);
651  }
652  }
653  else {
654  /* on right edge, not bottom corner */
655  dcol = VCOL2DCOL(gs, VCOLS(gs));
656
657  if (pt[Y] < ymax) {
658  vrow = Y2VROW(gs, pt[Y]);
659  drow = VROW2DROW(gs, vrow);
660  offset = DRC2OFF(gs, drow, dcol);
661  GET_MAPATT(buf, offset, p1[Z]);
662
663  drow = VROW2DROW(gs, vrow + 1);
664  offset = DRC2OFF(gs, drow, dcol);
665  GET_MAPATT(buf, offset, p2[Z]);
666
667  alpha = (VROW2Y(gs, vrow) - pt[Y]) / VYRES(gs);
668  pt[Z] = LERP(alpha, p1[Z], p2[Z]);
669
670  return (1);
671  }
672  else {
673  /* top right corner */
674  GET_MAPATT(buf, dcol, pt[Z]);
675
676  return (1);
677  }
678  }
679
680  return (0);
681 }
682
683 /*!
685
686  \param gs surface (geosurf)
687
688  \return 1
689  \return 0
690  */
691 int in_vregion(geosurf *gs, float *pt)
692 {
693  if (pt[X] >= 0.0 && pt[Y] <= gs->yrange) {
694  if (pt[X] <= VCOL2X(gs, VCOLS(gs))) {
695  return (pt[Y] >= VROW2Y(gs, VROWS(gs)));
696  }
697  }
698
699  return (0);
700 }
701
702 /*!
704
705  After all the intersections between the segment and triangle
706  edges have been found, they are in three lists. (intersections
707  with vertical, horizontal, and diagonal triangle edges)
708
709  Each list is ordered in space from first to last segment points,
710  but now the lists need to be woven together. This routine
711  starts with the first point of the segment and then checks the
712  next point in each list to find the closest, eliminating duplicates
713  along the way and storing the result in I3d.
714
715  \param gs surface (geosurf)
716  \param first first point
717  \param last last point
718  \param vi
719  \param hi
720  \param di
721
722  \return
723  */
724 int order_intersects(geosurf *gs, Point3 first, Point3 last, int vi, int hi,
725  int di)
726 {
727  int num, i, found, cv, ch, cd, cnum;
728  float dv, dh, dd, big, cpoint[2];
729
730  cv = ch = cd = cnum = 0;
731  num = vi + hi + di;
732
733  cpoint[X] = first[X];
734  cpoint[Y] = first[Y];
735
736  if (in_vregion(gs, first)) {
737  I3d[cnum][X] = first[X];
738  I3d[cnum][Y] = first[Y];
739  I3d[cnum][Z] = first[Z];
740  cnum++;
741  }
742
743  /* TODO: big could still be less than first dist */
744  big = gs->yrange * gs->yrange + gs->xrange * gs->xrange; /*BIG distance */
745  dv = dh = dd = big;
746
747  for (i = 0; i < num; i = cv + ch + cd) {
748  if (cv < vi) {
749  dv = dist_squared_2d(Vi[cv], cpoint);
750
751  if (dv < EPSILON) {
752  cv++;
753  continue;
754  }
755  }
756  else {
757  dv = big;
758  }
759
760  if (ch < hi) {
761  dh = dist_squared_2d(Hi[ch], cpoint);
762
763  if (dh < EPSILON) {
764  ch++;
765  continue;
766  }
767  }
768  else {
769  dh = big;
770  }
771
772  if (cd < di) {
773  dd = dist_squared_2d(Di[cd], cpoint);
774
775  if (dd < EPSILON) {
776  cd++;
777  continue;
778  }
779  }
780  else {
781  dd = big;
782  }
783
784  found = 0;
785
786  if (cd < di) {
787  if (dd <= dv && dd <= dh) {
788  found = 1;
789  cpoint[X] = I3d[cnum][X] = Di[cd][X];
790  cpoint[Y] = I3d[cnum][Y] = Di[cd][Y];
791  I3d[cnum][Z] = Di[cd][Z];
792  cnum++;
793
794  if (EQUAL(dd, dv)) {
795  cv++;
796  }
797
798  if (EQUAL(dd, dh)) {
799  ch++;
800  }
801
802  cd++;
803  }
804  }
805
806  if (!found) {
807  if (cv < vi) {
808  if (dv <= dh) {
809  found = 1;
810  cpoint[X] = I3d[cnum][X] = Vi[cv][X];
811  cpoint[Y] = I3d[cnum][Y] = Vi[cv][Y];
812  I3d[cnum][Z] = Vi[cv][Z];
813  cnum++;
814
815  if (EQUAL(dv, dh)) {
816  ch++;
817  }
818
819  cv++;
820  }
821  }
822  }
823
824  if (!found) {
825  if (ch < hi) {
826  cpoint[X] = I3d[cnum][X] = Hi[ch][X];
827  cpoint[Y] = I3d[cnum][Y] = Hi[ch][Y];
828  I3d[cnum][Z] = Hi[ch][Z];
829  cnum++;
830  ch++;
831  }
832  }
833
834  if (i == cv + ch + cd) {
835  G_debug(5, "order_intersects(): stuck on %d", cnum);
836  G_debug(5, "order_intersects(): cv = %d, ch = %d, cd = %d", cv, ch,
837  cd);
838  G_debug(5, "order_intersects(): dv = %f, dh = %f, dd = %f", dv, dh,
839  dd);
840
841  break;
842  }
843  }
844
845  if (EQUAL(last[X], cpoint[X]) && EQUAL(last[Y], cpoint[Y])) {
846  return (cnum);
847  }
848
849  if (in_vregion(gs, last)) {
850  /* TODO: check for last point on corner ? */
851  I3d[cnum][X] = last[X];
852  I3d[cnum][Y] = last[Y];
853  I3d[cnum][Z] = last[Z];
854  ++cnum;
855  }
856
857  return (cnum);
858 }
859
860 /*!
862
863  \todo For consistency, need to decide how last row & last column are
864  displayed - would it look funny to always draw last row/col with
865  finer resolution if necessary, or would it be better to only show
866  full rows/cols?
867
869
870  \param gs surface (geosurf)
871  \param bgn begin point
872  \param end end point
873  \param dir direction
874
875  \return
876  */
877 int get_vert_intersects(geosurf *gs, float *bgn, float *end, float *dir)
878 {
879  int fcol, lcol, incr, hits, num, offset, drow1, drow2;
880  float xl, yb, xr, yt, z1, z2, alpha;
881  float yres, xi, yi;
882  int bgncol, endcol, cols, rows;
883
884  yres = VYRES(gs);
885  cols = VCOLS(gs);
886  rows = VROWS(gs);
887
888  bgncol = X2VCOL(gs, bgn[X]);
889  endcol = X2VCOL(gs, end[X]);
890
891  if (bgncol > cols && endcol > cols) {
892  return 0;
893  }
894
895  if (bgncol == endcol) {
896  return 0;
897  }
898
899  fcol = dir[X] > 0 ? bgncol + 1 : bgncol;
900  lcol = dir[X] > 0 ? endcol : endcol + 1;
901
902  /* assuming only showing FULL cols */
903  incr = lcol - fcol > 0 ? 1 : -1;
904
905  while (fcol > cols || fcol < 0) {
906  fcol += incr;
907  }
908
909  while (lcol > cols || lcol < 0) {
910  lcol -= incr;
911  }
912
913  num = abs(lcol - fcol) + 1;
914
915  yb = gs->yrange - (yres * rows) - EPSILON;
916  yt = gs->yrange + EPSILON;
917
918  for (hits = 0; hits < num; hits++) {
919  xl = xr = VCOL2X(gs, fcol);
920
921  if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yt, xr, yb, &xi,
922  &yi)) {
923  Vi[hits][X] = xi;
924  Vi[hits][Y] = yi;
925
926  /* find data rows */
927  if (Flat) {
928  Vi[hits][Z] = gs->att[ATT_TOPO].constant;
929  }
930  else {
931  drow1 = Y2VROW(gs, Vi[hits][Y]) * gs->y_mod;
932  drow2 = (1 + Y2VROW(gs, Vi[hits][Y])) * gs->y_mod;
933
934  if (drow2 >= gs->rows) {
935  drow2 = gs->rows - 1; /*bottom edge */
936  }
937
938  alpha = ((gs->yrange - drow1 * gs->yres) - Vi[hits][Y]) / yres;
939
940  offset = DRC2OFF(gs, drow1, fcol * gs->x_mod);
941  GET_MAPATT(Ebuf, offset, z1);
942  offset = DRC2OFF(gs, drow2, fcol * gs->x_mod);
943  GET_MAPATT(Ebuf, offset, z2);
944  Vi[hits][Z] = LERP(alpha, z1, z2);
945  }
946  }
947
948  /* if they don't intersect, something's wrong! */
949  /* should only happen on endpoint, so it will be added later */
950  else {
951  hits--;
952  num--;
953  }
954
955  fcol += incr;
956  }
957
958  return (hits);
959 }
960
961 /*!
962  \brief Get horizontal intersects
963
964  \param gs surface (geosurf)
965  \param bgn begin point
966  \param end end point
967  \param dir
968
969  \return number of intersects
970  */
971 int get_horz_intersects(geosurf *gs, float *bgn, float *end, float *dir)
972 {
973  int frow, lrow, incr, hits, num, offset, dcol1, dcol2;
974  float xl, yb, xr, yt, z1, z2, alpha;
975  float xres, xi, yi;
976  int bgnrow, endrow, rows, cols;
977
978  xres = VXRES(gs);
979  cols = VCOLS(gs);
980  rows = VROWS(gs);
981
982  bgnrow = Y2VROW(gs, bgn[Y]);
983  endrow = Y2VROW(gs, end[Y]);
984  if (bgnrow == endrow) {
985  return 0;
986  }
987
988  if (bgnrow > rows && endrow > rows) {
989  return 0;
990  }
991
992  frow = dir[Y] > 0 ? bgnrow : bgnrow + 1;
993  lrow = dir[Y] > 0 ? endrow + 1 : endrow;
994
995  /* assuming only showing FULL rows */
996  incr = lrow - frow > 0 ? 1 : -1;
997
998  while (frow > rows || frow < 0) {
999  frow += incr;
1000  }
1001
1002  while (lrow > rows || lrow < 0) {
1003  lrow -= incr;
1004  }
1005
1006  num = abs(lrow - frow) + 1;
1007
1008  xl = 0.0 - EPSILON;
1009  xr = xres * cols + EPSILON;
1010
1011  for (hits = 0; hits < num; hits++) {
1012  yb = yt = VROW2Y(gs, frow);
1013
1014  if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yt, xr, yb, &xi,
1015  &yi)) {
1016  Hi[hits][X] = xi;
1017  Hi[hits][Y] = yi;
1018
1019  /* find data cols */
1020  if (Flat) {
1021  Hi[hits][Z] = gs->att[ATT_TOPO].constant;
1022  }
1023  else {
1024  dcol1 = X2VCOL(gs, Hi[hits][X]) * gs->x_mod;
1025  dcol2 = (1 + X2VCOL(gs, Hi[hits][X])) * gs->x_mod;
1026
1027  if (dcol2 >= gs->cols) {
1028  dcol2 = gs->cols - 1; /* right edge */
1029  }
1030
1031  alpha = (Hi[hits][X] - (dcol1 * gs->xres)) / xres;
1032
1033  offset = DRC2OFF(gs, frow * gs->y_mod, dcol1);
1034  GET_MAPATT(Ebuf, offset, z1);
1035  offset = DRC2OFF(gs, frow * gs->y_mod, dcol2);
1036  GET_MAPATT(Ebuf, offset, z2);
1037  Hi[hits][Z] = LERP(alpha, z1, z2);
1038  }
1039  }
1040
1041  /* if they don't intersect, something's wrong! */
1042  /* should only happen on endpoint, so it will be added later */
1043  else {
1044  hits--;
1045  num--;
1046  }
1047
1048  frow += incr;
1049  }
1050
1051  return (hits);
1052 }
1053
1054 /*!
1055  \brief Get diagonal intersects
1056
1058
1059  \param gs surface (geosurf)
1060  \param bgn begin point
1061  \param end end point
1062  \param dir ? (unused)
1063
1064  \return number of intersects
1065  */
1066 int get_diag_intersects(geosurf *gs, float *bgn, float *end, float *dir UNUSED)
1067 {
1068  int fdig, ldig, incr, hits, num, offset;
1069  int vrow, vcol, drow1, drow2, dcol1, dcol2;
1070  float xl, yb, xr, yt, z1, z2, alpha;
1071  float xres, yres, xi, yi, dx, dy;
1072  int diags, cols, rows, lower;
1073  Point3 pt;
1074
1075  xres = VXRES(gs);
1076  yres = VYRES(gs);
1077  cols = VCOLS(gs);
1078  rows = VROWS(gs);
1079  diags = rows + cols; /* -1 ? */
1080
1081  /* determine upper/lower triangle for last */
1082  vrow = Y2VROW(gs, end[Y]);
1083  vcol = X2VCOL(gs, end[X]);
1084  pt[X] = VCOL2X(gs, vcol);
1085  pt[Y] = VROW2Y(gs, vrow + 1);
1086  lower = ((end[X] - pt[X]) / xres > (end[Y] - pt[Y]) / yres);
1087  ldig = lower ? vrow + vcol + 1 : vrow + vcol;
1088
1089  /* determine upper/lower triangle for first */
1090  vrow = Y2VROW(gs, bgn[Y]);
1091  vcol = X2VCOL(gs, bgn[X]);
1092  pt[X] = VCOL2X(gs, vcol);
1093  pt[Y] = VROW2Y(gs, vrow + 1);
1094  lower = ((bgn[X] - pt[X]) / xres > (bgn[Y] - pt[Y]) / yres);
1095  fdig = lower ? vrow + vcol + 1 : vrow + vcol;
1096
1097  /* adjust according to direction */
1098  if (ldig > fdig) {
1099  fdig++;
1100  }
1101
1102  if (fdig > ldig) {
1103  ldig++;
1104  }
1105
1106  incr = ldig - fdig > 0 ? 1 : -1;
1107
1108  while (fdig > diags || fdig < 0) {
1109  fdig += incr;
1110  }
1111
1112  while (ldig > diags || ldig < 0) {
1113  ldig -= incr;
1114  }
1115
1116  num = abs(ldig - fdig) + 1;
1117
1118  for (hits = 0; hits < num; hits++) {
1119  yb = gs->yrange - (yres * (fdig < rows ? fdig : rows)) - EPSILON;
1120  xl = VCOL2X(gs, (fdig < rows ? 0 : fdig - rows)) - EPSILON;
1121  yt = gs->yrange - (yres * (fdig < cols ? 0 : fdig - cols)) + EPSILON;
1122  xr = VCOL2X(gs, (fdig < cols ? fdig : cols)) + EPSILON;
1123
1124  if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yb, xr, yt, &xi,
1125  &yi)) {
1126  Di[hits][X] = xi;
1127  Di[hits][Y] = yi;
1128
1129  if (ISNODE(xi, xres)) {
1130  /* then it's also a ynode */
1131  num--;
1132  hits--;
1133  continue;
1134  }
1135
1136  /* find data rows */
1137  drow1 = Y2VROW(gs, Di[hits][Y]) * gs->y_mod;
1138  drow2 = (1 + Y2VROW(gs, Di[hits][Y])) * gs->y_mod;
1139
1140  if (drow2 >= gs->rows) {
1141  drow2 = gs->rows - 1; /* bottom edge */
1142  }
1143
1144  /* find data cols */
1145  if (Flat) {
1146  Di[hits][Z] = gs->att[ATT_TOPO].constant;
1147  }
1148  else {
1149  dcol1 = X2VCOL(gs, Di[hits][X]) * gs->x_mod;
1150  dcol2 = (1 + X2VCOL(gs, Di[hits][X])) * gs->x_mod;
1151
1152  if (dcol2 >= gs->cols) {
1153  dcol2 = gs->cols - 1; /* right edge */
1154  }
1155
1156  dx = DCOL2X(gs, dcol2) - Di[hits][X];
1157  dy = DROW2Y(gs, drow1) - Di[hits][Y];
1158  alpha =
1159  sqrt(dx * dx + dy * dy) / sqrt(xres * xres + yres * yres);
1160
1161  offset = DRC2OFF(gs, drow1, dcol2);
1162  GET_MAPATT(Ebuf, offset, z1);
1163  offset = DRC2OFF(gs, drow2, dcol1);
1164  GET_MAPATT(Ebuf, offset, z2);
1165  Di[hits][Z] = LERP(alpha, z1, z2);
1166  }
1167  }
1168
1169  /* if they don't intersect, something's wrong! */
1170  /* should only happen on endpoint, so it will be added later */
1171  else {
1172  hits--;
1173  num--;
1174  }
1175
1176  fdig += incr;
1177  }
1178
1179  return (hits);
1180 }
1181
1182 /*!
1183  \brief Line intersect
1184
1186  Modified for floating point: Bill Brown
1187
1188  This function computes whether two line segments,
1189  respectively joining the input points (x1,y1) -- (x2,y2)
1190  and the input points (x3,y3) -- (x4,y4) intersect.
1191  If the lines intersect, the output variables x, y are
1192  set to coordinates of the point of intersection.
1193
1194  \param x1,y1,x2,y2 coordinates of endpoints of one segment
1195  \param x3,y3,x4,y4 coordinates of endpoints of other segment
1196  \param[out] x,y coordinates of intersection point
1197
1198  \return 0 no intersection
1199  \return 1 intersect
1200  \return 2 collinear
1201  */
1202 int segs_intersect(float x1, float y1, float x2, float y2, float x3, float y3,
1203  float x4, float y4, float *x, float *y)
1204 {
1205  float a1, a2, b1, b2, c1, c2; /* Coefficients of line eqns. */
1206  float r1, r2, r3, r4; /* 'Sign' values */
1207  float denom, /* offset, */ num; /* Intermediate values */
1208
1209  /* Compute a1, b1, c1, where line joining points 1 and 2
1210  * is "a1 x + b1 y + c1 = 0".
1211  */
1212  a1 = y2 - y1;
1213  b1 = x1 - x2;
1214  c1 = x2 * y1 - x1 * y2;
1215
1216  /* Compute r3 and r4.
1217  */
1218  r3 = a1 * x3 + b1 * y3 + c1;
1219  r4 = a1 * x4 + b1 * y4 + c1;
1220
1221  /* Check signs of r3 and r4. If both point 3 and point 4 lie on
1222  * same side of line 1, the line segments do not intersect.
1223  */
1224
1225  if (!EQUAL(r3, 0.0) && !EQUAL(r4, 0.0) && SAME_SIGNS(r3, r4)) {
1226  return (DONT_INTERSECT);
1227  }
1228
1229  /* Compute a2, b2, c2 */
1230  a2 = y4 - y3;
1231  b2 = x3 - x4;
1232  c2 = x4 * y3 - x3 * y4;
1233
1234  /* Compute r1 and r2 */
1235  r1 = a2 * x1 + b2 * y1 + c2;
1236  r2 = a2 * x2 + b2 * y2 + c2;
1237
1238  /* Check signs of r1 and r2. If both point 1 and point 2 lie
1239  * on same side of second line segment, the line segments do
1240  * not intersect.
1241  */
1242
1243  if (!EQUAL(r1, 0.0) && !EQUAL(r2, 0.0) && SAME_SIGNS(r1, r2)) {
1244  return (DONT_INTERSECT);
1245  }
1246
1247  /* Line segments intersect: compute intersection point.
1248  */
1249  denom = a1 * b2 - a2 * b1;
1250
1251  if (denom == 0) {
1252  return (COLLINEAR);
1253  }
1254
1255  /* offset = denom < 0 ? -denom / 2 : denom / 2; */
1256
1257  /* The denom/2 is to get rounding instead of truncating. It
1258  * is added or subtracted to the numerator, depending upon the
1259  * sign of the numerator.
1260  */
1261  num = b1 * c2 - b2 * c1;
1262
1263  *x = num / denom;
1264
1265  num = a2 * c1 - a1 * c2;
1266  *y = num / denom;
1267
1268  return (DO_INTERSECT);
1269 }
1270
1271 /*!
1272  \brief Check if point is on plane
1273
1274  Plane defined by three points here; user fills in unk[X] & unk[Y]
1275
1276  \param p1,p2,p3 points defining plane
1277  \param unk point
1278
1279  \return 1 point on plane
1280  \return 0 point not on plane
1281  */
1283 {
1284  float plane[4];
1285
1286  P3toPlane(p1, p2, p3, plane);
1287
1288  return (XY_intersect_plane(unk, plane));
1289 }
1290
1291 /*!
1292  \brief Check for intersection (point and plane)
1293
1294  Ax + By + Cz + D = 0, so z = (Ax + By + D) / -C
1295
1296  User fills in intersect[X] & intersect[Y]
1297
1298  \param[out] intersect intersect coordinates
1299  \param plane plane definition
1300
1301  \return 0 doesn't intersect
1302  \return 1 intesects
1303  */
1304 int XY_intersect_plane(float *intersect, float *plane)
1305 {
1306  float x, y;
1307
1308  if (!plane[Z]) {
1309  return (0); /* doesn't intersect */
1310  }
1311
1312  x = intersect[X];
1313  y = intersect[Y];
1314  intersect[Z] = (plane[X] * x + plane[Y] * y + plane[W]) / -plane[Z];
1315
1316  return (1);
1317 }
1318
1319 /*!
1320  \brief Define plane
1321
1322  \param p1,p2,p3 three point on plane
1323  \param[out] plane plane definition
1324
1325  \return 1
1326  */
1327 int P3toPlane(Point3 p1, Point3 p2, Point3 p3, float *plane)
1328 {
1329  Point3 v1, v2, norm;
1330
1331  v1[X] = p1[X] - p3[X];
1332  v1[Y] = p1[Y] - p3[Y];
1333  v1[Z] = p1[Z] - p3[Z];
1334
1335  v2[X] = p2[X] - p3[X];
1336  v2[Y] = p2[Y] - p3[Y];
1337  v2[Z] = p2[Z] - p3[Z];
1338
1339  V3Cross(v1, v2, norm);
1340
1341  plane[X] = norm[X];
1342  plane[Y] = norm[Y];
1343  plane[Z] = norm[Z];
1344  plane[W] = -p3[X] * norm[X] - p3[Y] * norm[Y] - p3[Z] * norm[Z];
1345
1346  return (1);
1347 }
1348
1349 /*!
1350  \brief Get cross product
1351
1352  \param a,b,c
1353
1354  \return cross product c = a cross b
1355  */
1357 {
1358  c[X] = (a[Y] * b[Z]) - (a[Z] * b[Y]);
1359  c[Y] = (a[Z] * b[X]) - (a[X] * b[Z]);
1360  c[Z] = (a[X] * b[Y]) - (a[Y] * b[X]);
1361
1362  return (1);
1363 }
#define EPSILON
Definition: blas_level_2.c:26
#define NULL
Definition: ccmath.h:32
void G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:150
void G_warning(const char *,...) __attribute__((format(printf
int G_debug(int, const char *,...) __attribute__((format(printf
Definition: gs.c:1314
float GS_P2distance(float *, float *)
Calculate distance in plane.
Definition: gs_util.c:160
int gs_get_att_src(geosurf *, int)
Get attribute source.
Definition: gs.c:656
void GS_v3eq(float *, float *)
Copy vector values.
Definition: gs_util.c:178
typbuff * gs_get_att_typbuff(geosurf *, int, int)
Get attribute data buffer.
Definition: gs.c:681
void GS_v2dir(float *, float *, float *)
Get a normalized direction from v1 to v2, store in v3 (2D)
Definition: gs_util.c:382
#define UNUSED
A macro for an attribute, if attached to a variable, indicating that the variable is not used.
Definition: gis.h:47
#define _(str)
Definition: glocale.h:10
#define COLLINEAR
Definition: gsdrape.c:62
int P3toPlane(Point3 p1, Point3 p2, Point3 p3, float *plane)
Define plane.
Definition: gsdrape.c:1327
#define ISNODE(p, res)
Definition: gsdrape.c:66
int in_vregion(geosurf *gs, float *pt)
Definition: gsdrape.c:691
#define SAME_SIGNS(a, b)
Definition: gsdrape.c:68
Point3 * gsdrape_get_segments(geosurf *gs, float *bgn, float *end, int *num)
Definition: gsdrape.c:349
int order_intersects(geosurf *gs, Point3 first, Point3 last, int vi, int hi, int di)
Definition: gsdrape.c:724
int _viewcell_tri_interp(geosurf *gs, Point3 pt)
Definition: gsdrape.c:463
Point3 * gsdrape_get_allsegments(geosurf *gs, float *bgn, float *end, int *num)
Get all segments.
Definition: gsdrape.c:399
int gsdrape_set_surface(geosurf *gs)
Definition: gsdrape.c:198
#define EQUAL(a, b)
Definition: gsdrape.c:65
void interp_first_last(geosurf *gs, float *bgn, float *end, Point3 f, Point3 l)
Definition: gsdrape.c:438
#define DONT_INTERSECT
Definition: gsdrape.c:60
int V3Cross(Point3 a, Point3 b, Point3 c)
Get cross product.
Definition: gsdrape.c:1356
int viewcell_tri_interp(geosurf *gs, typbuff *buf, Point3 pt, int check_mask)
Definition: gsdrape.c:503
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:1202
int get_horz_intersects(geosurf *gs, float *bgn, float *end, float *dir)
Get horizontal intersects.
Definition: gsdrape.c:971
#define LERP(a, l, h)
Definition: gsdrape.c:64
int get_vert_intersects(geosurf *gs, float *bgn, float *end, float *dir)
Definition: gsdrape.c:877
int seg_intersect_vregion(geosurf *gs, float *bgn, float *end)
Check if segment intersect vector region.
Definition: gsdrape.c:232
#define DO_INTERSECT
Definition: gsdrape.c:61
int get_diag_intersects(geosurf *gs, float *bgn, float *end, float *dir UNUSED)
Get diagonal intersects.
Definition: gsdrape.c:1066
int XY_intersect_plane(float *intersect, float *plane)
Check for intersection (point and plane)
Definition: gsdrape.c:1304
int Point_on_plane(Point3 p1, Point3 p2, Point3 p3, Point3 unk)
Check if point is on plane.
Definition: gsdrape.c:1282
#define GET_MAPATT(buff, offset, att)
Definition: gsget.h:29
#define X
Definition: ogsf.h:140
#define ATT_TOPO
Definition: ogsf.h:75
float Point3[3]
Definition: ogsf.h:205
#define Z
Definition: ogsf.h:142
#define W
Definition: ogsf.h:143
#define Y
Definition: ogsf.h:141
#define MAP_ATT
Definition: ogsf.h:85
#define CONST_ATT
Definition: ogsf.h:86
double b
Definition: r_raster.c:39
double l
Definition: r_raster.c:39
#define VYRES(gs)
Definition: rowcol.h:10
#define Y2VROW(gs, py)
Definition: rowcol.h:27
#define VXRES(gs)
Definition: rowcol.h:9
#define VCOL2X(gs, vcol)
Definition: rowcol.h:40
#define VCOLS(gs)
Definition: rowcol.h:14
#define VROWS(gs)
Definition: rowcol.h:13
#define DRC2OFF(gs, drow, dcol)
Definition: rowcol.h:17
#define DCOL2X(gs, dcol)
Definition: rowcol.h:36
#define VROW2Y(gs, vrow)
Definition: rowcol.h:39
#define VROW2DROW(gs, vrow)
Definition: rowcol.h:31
#define X2VCOL(gs, px)
Definition: rowcol.h:28
#define VCOL2DCOL(gs, vcol)
Definition: rowcol.h:32
#define DROW2Y(gs, drow)
Definition: rowcol.h:35
if(!(yy_init))
Definition: sqlp.yy.c:775
Definition: ogsf.h:256
int rows
Definition: ogsf.h:258
double yres
Definition: ogsf.h:264
float yrange
Definition: ogsf.h:268
int cols
Definition: ogsf.h:258
float xrange
Definition: ogsf.h:268
gsurf_att att[MAX_ATTS]
Definition: ogsf.h:259
int y_mod
Definition: ogsf.h:270
double xres
Definition: ogsf.h:264
int x_mod
Definition: ogsf.h:270
float constant
Definition: ogsf.h:250
Definition: clip.h:6
Definition: ogsf.h:208
#define x