GRASS GIS 7 Programmer's Manual  7.5.svn(2018)-r72729
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 /*!
192  \brief ADD
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 /*!
341  \brief ADD
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 /*!
433  \brief ADD
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 /*!
462  \brief ADD
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 /*!
477  \brief ADD
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  */
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 
688 /*!
689  \brief ADD
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 /*!
708  \brief ADD
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 /*!
866  \brief ADD
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 
873  Colinear already eliminated
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 
1065  Colinear already eliminated
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 
1193  Author: Mukesh Prasad
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)
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: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)
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
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)
ADD.
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)
ADD.
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)
ADD.
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)
ADD.
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