GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-9f2952cf23
gs_query.c
Go to the documentation of this file.
1 /*!
2  \file lib/ogsf/gs_query.c
3 
4  \brief OGSF library - query (lower level functions)
5 
6  GRASS OpenGL gsurf OGSF Library
7 
8  (C) 1999-2008 by the GRASS Development Team
9 
10  This program is free software under the
11  GNU General Public License (>=v2).
12  Read the file COPYING that comes with GRASS
13  for details.
14 
15  \author Bill Brown USACERL (January 1994)
16  \author Doxygenized by Martin Landa <landa.martin gmail.com> (May 2008)
17  */
18 
19 #include <math.h>
20 
21 #include <grass/gis.h>
22 #include <grass/ogsf.h>
23 
24 /*!
25  \brief Values needed for Ray-Convex Polyhedron Intersection Test below
26  originally by Eric Haines, erich@eye.com
27  */
28 #ifndef HUGE_VAL
29 #define HUGE_VAL 1.7976931348623157e+308
30 #endif
31 
32 /* return codes */
33 #define MISSED 0
34 #define FRONTFACE 1
35 #define BACKFACE -1
36 /* end Ray-Convex Polyhedron Intersection Test values */
37 
38 /*!
39  \brief Crude method of intersecting line of sight with closest part of
40  surface.
41 
42  Uses los vector to determine the point of first intersection
43  which is returned in point. Returns 0 if los doesn't intersect.
44 
45  \param surfid surface id
46  \param los should be in surf-world coordinates
47  \param[out] point intersect point (real)
48 
49  \return 0 on failure
50  \return 1 on success
51  */
52 int gs_los_intersect1(int surfid, float (*los)[3], float *point)
53 {
54  float dx, dy, dz, u_d[3];
55  float a[3], incr, min_incr, tlen, len;
56  int outside, above, below, edge, istep;
57  float b[3];
58  geosurf *gs;
59  typbuff *buf;
60 
61  G_debug(3, "gs_los_intersect1():");
62 
63  if (NULL == (gs = gs_get_surf(surfid))) {
64  return (0);
65  }
66 
67  if (0 == GS_v3dir(los[FROM], los[TO], u_d)) {
68  return (0);
69  }
70 
71  buf = gs_get_att_typbuff(gs, ATT_TOPO, 0);
72 
73  istep = edge = below = 0;
74 
75  len = 0.0;
76  tlen = GS_distance(los[FROM], los[TO]);
77 
78  incr = tlen / 1000.0;
79  min_incr = incr / 1000.0;
80 
81  dx = incr * u_d[X];
82  dy = incr * u_d[Y];
83  dz = incr * u_d[Z];
84 
85  a[X] = los[FROM][X];
86  a[Y] = los[FROM][Y];
87  a[Z] = los[FROM][Z];
88 
89  b[X] = a[X] - gs->x_trans;
90  b[Y] = a[Y] - gs->y_trans;
91 
92  if (viewcell_tri_interp(gs, buf, b, 0)) {
93  /* expects surface coords */
94  b[Z] += gs->z_trans;
95 
96  if (a[Z] < b[Z]) {
97  /* viewing from below surface */
98  /* don't use this method
99  fprintf(stderr,"view from below\n");
100  below = 1;
101  */
102 
103  return (0);
104  }
105  }
106 
107  while (incr > min_incr) {
108  outside = 0;
109  above = 0;
110  b[X] = a[X] - gs->x_trans;
111  b[Y] = a[Y] - gs->y_trans;
112 
113  if (viewcell_tri_interp(gs, buf, b, 0)) {
114  /* ignores masks */
115  b[Z] += gs->z_trans;
116  above = (a[Z] > b[Z]);
117  }
118  else {
119  outside = 1;
120 
121  if (istep > 10) {
122  edge = 1;
123  below = 1;
124  }
125  }
126 
127  while (outside || above) {
128  a[X] += dx;
129  a[Y] += dy;
130  a[Z] += dz;
131  len += incr;
132  outside = 0;
133  above = 0;
134  b[X] = a[X] - gs->x_trans;
135  b[Y] = a[Y] - gs->y_trans;
136 
137  if (viewcell_tri_interp(gs, buf, b, 0)) {
138  b[Z] += gs->z_trans;
139  above = (a[Z] > b[Z]);
140  }
141  else {
142  outside = 1;
143  }
144 
145  if (len > tlen) {
146  return 0; /* over surface */ /* under surface */
147  }
148  }
149 
150  /* could look for spikes here - see if any data points along
151  shadow of line on surf go above los */
152 
153  /* back up several spots? */
154  a[X] -= (1.0 * dx);
155  a[Y] -= (1.0 * dy);
156  a[Z] -= (1.0 * dz);
157  incr /= 2.0;
158  ++istep;
159  dx = incr * u_d[X];
160  dy = incr * u_d[Y];
161  dz = incr * u_d[Z];
162  }
163 
164  if ((edge) && (b[Z] - (a[Z] + dz * 2.0) > incr * u_d[Z])) {
165  G_debug(3, " looking under surface");
166 
167  return 0;
168  }
169 
170  point[X] = b[X];
171  point[Y] = b[Y];
172  point[Z] = b[Z] - gs->z_trans;
173 
174  return (1);
175 }
176 
177 /*!
178  \brief Crude method of intersecting line of sight with closest part of
179  surface.
180 
181  This version uses the shadow of the los projected down to
182  the surface to generate a line_on_surf, then follows each
183  point in that line until the los intersects it.
184 
185  \param surfid surface id
186  \param los should be in surf-world coordinates
187  \param[out] point intersect point (real)
188 
189  \return 0 on failure
190  \return 1 on success
191  */
192 int gs_los_intersect(int surfid, float **los, float *point)
193 {
194  double incr;
195  float p1, p2, u_d[3];
196  int above, ret, num, i, usedx;
197  float a[3], b[3];
198  float bgn[3], end[3], a1[3];
199  geosurf *gs;
200  typbuff *buf;
201  Point3 *points;
202 
203  G_debug(3, "gs_los_intersect");
204 
205  if (NULL == (gs = gs_get_surf(surfid))) {
206  return (0);
207  }
208 
209  if (0 == GS_v3dir(los[FROM], los[TO], u_d)) {
210  return (0);
211  }
212 
213  buf = gs_get_att_typbuff(gs, ATT_TOPO, 0);
214 
215  GS_v3eq(bgn, los[FROM]);
216  GS_v3eq(end, los[TO]);
217 
218  bgn[X] -= gs->x_trans;
219  bgn[Y] -= gs->y_trans;
220 
221  end[X] -= gs->x_trans;
222  end[Y] -= gs->y_trans;
223 
224  /* trans? */
225  points = gsdrape_get_allsegments(gs, bgn, end, &num);
226 
227  /* DEBUG
228  {
229  float t1[3], t2[3];
230 
231  t1[X] = los[FROM][X] ;
232  t1[Y] = los[FROM][Y] ;
233 
234  t2[X] = los[TO][X] ;
235  t2[Y] = los[TO][Y] ;
236 
237  GS_set_draw(GSD_FRONT);
238  gsd_pushmatrix();
239  gsd_do_scale(1);
240  gsd_translate(gs->x_trans, gs->y_trans, gs->z_trans);
241  gsd_linewidth(1);
242  gsd_color_func(GS_default_draw_color());
243  gsd_line_onsurf(gs, t1, t2);
244  gsd_popmatrix();
245  GS_set_draw(GSD_BACK);
246  gsd_flush();
247  }
248  fprintf(stderr,"%d points to check\n", num);
249  fprintf(stderr,"point0 = %.6lf %.6lf %.6lf FT =%.6lf %.6lf %.6lf\n",
250  points[0][X],points[0][Y],points[0][Z],
251  los[FROM][X],los[FROM][Y],los[FROM][Z]);
252  fprintf(stderr,"incr1 = %.6lf: %.6lf %.6lf
253  %.6lf\n",incr,u_d[X],u_d[Y],u_d[Z]); fprintf(stderr,"first point below
254  surf\n"); fprintf(stderr,"incr2 = %f\n", (float)incr);
255  fprintf(stderr,"(%d/%d) %f > %f\n", i,num, a[Z], points[i][Z]);
256  fprintf(stderr,"incr3 = %f\n", (float)incr);
257  fprintf(stderr,"all points above surf\n");
258  */
259 
260  if (num < 2) {
261  G_debug(3, " %d points to check", num);
262 
263  return (0);
264  }
265 
266  /* use larger of deltas for better precision */
267  usedx = (fabs(u_d[X]) > fabs(u_d[Y]));
268  if (usedx) {
269  incr = ((points[0][X] - (los[FROM][X] - gs->x_trans)) / u_d[X]);
270  }
271  else if (u_d[Y]) {
272  incr = ((points[0][Y] - (los[FROM][Y] - gs->y_trans)) / u_d[Y]);
273  }
274  else {
275  point[X] = los[FROM][X] - gs->x_trans;
276  point[Y] = los[FROM][Y] - gs->y_trans;
277 
278  return (viewcell_tri_interp(gs, buf, point, 1));
279  }
280 
281  /* DEBUG
282  fprintf(stderr,"-----------------------------\n");
283  fprintf(stderr,"%d points to check\n", num);
284  fprintf(stderr,"incr1 = %.6lf: %.9f %.9f
285  %.9f\n",incr,u_d[X],u_d[Y],u_d[Z]); fprintf(stderr,
286  "\tpoint0 = %.6f %.6f %.6f\n\tFT = %.6f %.6f %.6f\n\tpoint%d = %.6f
287  %.6f\n", points[0][X],points[0][Y],points[0][Z],
288  los[FROM][X],los[FROM][Y],los[FROM][Z],
289  num-1, points[num-1][X],points[num-1][Y]);
290  */
291 
292  /* This should bring us right above (or below) the first point */
293  a[X] = los[FROM][X] + incr * u_d[X] - gs->x_trans;
294  a[Y] = los[FROM][Y] + incr * u_d[Y] - gs->y_trans;
295  a[Z] = los[FROM][Z] + incr * u_d[Z] - gs->z_trans;
296 
297  if (a[Z] < points[0][Z]) {
298  /* viewing from below surface */
299  /* don't use this method */
300  /* DEBUG
301  fprintf(stderr,"first point below surf\n");
302  fprintf(stderr,"aZ= %.6f point0 = %.6f %.6f %.6f FT =%.6f %.6f
303  %.6f\n", a[Z], points[0][X],points[0][Y],points[0][Z],
304  los[FROM][X],los[FROM][Y],los[FROM][Z]);
305  */
306  return (0);
307  }
308 
309  GS_v3eq(a1, a);
310  GS_v3eq(b, a);
311 
312  for (i = 1; i < num; i++) {
313  if (usedx) {
314  incr = ((points[i][X] - a1[X]) / u_d[X]);
315  }
316  else {
317  incr = ((points[i][Y] - a1[Y]) / u_d[Y]);
318  }
319 
320  a[X] = a1[X] + (incr * u_d[X]);
321  a[Y] = a1[Y] + (incr * u_d[Y]);
322  a[Z] = a1[Z] + (incr * u_d[Z]);
323  above = (a[Z] >= points[i][Z]);
324 
325  if (above) {
326  GS_v3eq(b, a);
327  continue;
328  }
329 
330  /*
331  * Now we know b[Z] is above points[i-1]
332  * and a[Z] is below points[i]
333  * Since there should only be one polygon along this seg,
334  * just interpolate to intersect
335  */
336 
337  if (usedx) {
338  incr = ((a[X] - b[X]) / u_d[X]);
339  }
340  else {
341  incr = ((a[Y] - b[Y]) / u_d[Y]);
342  }
343 
344  if (1 == (ret = segs_intersect(1.0, points[i][Z], 0.0, points[i - 1][Z],
345  1.0, a[Z], 0.0, b[Z], &p1, &p2))) {
346  point[X] = points[i - 1][X] + (u_d[X] * incr * p1);
347  point[Y] = points[i - 1][Y] + (u_d[Y] * incr * p1);
348  point[Z] = p2;
349 
350  return (1);
351  }
352 
353  G_debug(3, " line of sight error %d", ret);
354 
355  return 0;
356  }
357 
358  /* over surface */
359  return 0;
360 }
361 
362 /*!
363  \brief Ray-Convex Polyhedron Intersection Test
364 
365  Originally by Eric Haines, erich@eye.com
366 
367  This test checks the ray against each face of a polyhedron, checking whether
368  the set of intersection points found for each ray-plane intersection
369  overlaps the previous intersection results. If there is no overlap (i.e.
370  no line segment along the ray that is inside the polyhedron), then the
371  ray misses and returns 0; else 1 is returned if the ray is entering the
372  polyhedron, -1 if the ray originates inside the polyhedron. If there is
373  an intersection, the distance and the number of the face hit is returned.
374 
375  \param org,dir origin and direction of ray
376  \param tmax maximum useful distance along ray
377  \param phdrn list of planes in convex polyhedron
378  \param ph_num number of planes in convex polyhedron
379  \param[out] tresult distance of intersection along ray
380  \param[out] pn number of face hit (0 to ph_num-1)
381 
382  \return FACE code
383  */
384 int RayCvxPolyhedronInt(Point3 org, Point3 dir, double tmax, Point4 *phdrn,
385  int ph_num, double *tresult, int *pn)
386 {
387  double tnear, tfar, t, vn, vd;
388  int fnorm_num, bnorm_num; /* front/back face # hit */
389 
390  tnear = -HUGE_VAL;
391  tfar = tmax;
392 
393  /* Test each plane in polyhedron */
394  for (; ph_num--;) {
395  /* Compute intersection point T and sidedness */
396  vd = DOT3(dir, phdrn[ph_num]);
397  vn = DOT3(org, phdrn[ph_num]) + phdrn[ph_num][W];
398 
399  if (vd == 0.0) {
400  /* ray is parallel to plane - check if ray origin is inside plane's
401  half-space */
402  if (vn > 0.0) {
403  /* ray origin is outside half-space */
404  return (MISSED);
405  }
406  }
407  else {
408  /* ray not parallel - get distance to plane */
409  t = -vn / vd;
410 
411  if (vd < 0.0) {
412  /* front face - T is a near point */
413  if (t > tfar) {
414  return (MISSED);
415  }
416 
417  if (t > tnear) {
418  /* hit near face, update normal */
419  fnorm_num = ph_num;
420  tnear = t;
421  }
422  }
423  else {
424  /* back face - T is a far point */
425  if (t < tnear) {
426  return (MISSED);
427  }
428 
429  if (t < tfar) {
430  /* hit far face, update normal */
431  bnorm_num = ph_num;
432  tfar = t;
433  }
434  }
435  }
436  }
437 
438  /* survived all tests */
439  /* Note: if ray originates on polyhedron, may want to change 0.0 to some
440  * epsilon to avoid intersecting the originating face.
441  */
442  if (tnear >= 0.0) {
443  /* outside, hitting front face */
444  *tresult = tnear;
445  *pn = fnorm_num;
446 
447  return (FRONTFACE);
448  }
449  else {
450  if (tfar < tmax) {
451  /* inside, hitting back face */
452  *tresult = tfar;
453  *pn = bnorm_num;
454 
455  return (BACKFACE);
456  }
457  else {
458  /* inside, but back face beyond tmax */
459  return (MISSED);
460  }
461  }
462 }
463 
464 /*!
465  \brief Get data bounds for plane
466 
467  \param[out] planes
468  */
470 {
471  float n, s, w, e, b, t;
472  Point3 tlfront, brback;
473 
474  GS_get_zrange(&b, &t, 0);
475  gs_get_xrange(&w, &e);
476  gs_get_yrange(&s, &n);
477 
478  tlfront[X] = tlfront[Y] = 0.0;
479  tlfront[Z] = t;
480 
481  brback[X] = e - w;
482  brback[Y] = n - s;
483  brback[Z] = b;
484 
485  /* top */
486  planes[0][X] = planes[0][Y] = 0.0;
487  planes[0][Z] = 1.0;
488  planes[0][W] = -(DOT3(planes[0], tlfront));
489 
490  /* bottom */
491  planes[1][X] = planes[1][Y] = 0.0;
492  planes[1][Z] = -1.0;
493  planes[1][W] = -(DOT3(planes[1], brback));
494 
495  /* left */
496  planes[2][Y] = planes[2][Z] = 0.0;
497  planes[2][X] = -1.0;
498  planes[2][W] = -(DOT3(planes[2], tlfront));
499 
500  /* right */
501  planes[3][Y] = planes[3][Z] = 0.0;
502  planes[3][X] = 1.0;
503  planes[3][W] = -(DOT3(planes[3], brback));
504 
505  /* front */
506  planes[4][X] = planes[4][Z] = 0.0;
507  planes[4][Y] = -1.0;
508  planes[4][W] = -(DOT3(planes[4], tlfront));
509 
510  /* back */
511  planes[5][X] = planes[5][Z] = 0.0;
512  planes[5][Y] = 1.0;
513  planes[5][W] = -(DOT3(planes[5], brback));
514 
515  return;
516 }
517 
518 /*!
519  Gets all current cutting planes & data bounding planes
520 
521  Intersects los with resulting convex polyhedron, then replaces los[FROM] with
522  first point on ray inside data.
523 
524  \param[out] los
525 
526  \return 0 on failure
527  \return 1 on success
528  */
530 {
531  Point4 planes[12]; /* MAX_CPLANES + 6 - should define this */
532  Point3 dir;
533  double dist, maxdist;
534  int num, ret, retp; /* might want to tell if retp is a clipping plane */
535 
536  gs_get_databounds_planes(planes);
537  num = gsd_get_cplanes(planes + 6);
538  GS_v3dir(los[FROM], los[TO], dir);
539  maxdist = GS_distance(los[FROM], los[TO]);
540 
541  ret = RayCvxPolyhedronInt(los[0], dir, maxdist, planes, num + 6, &dist,
542  &retp);
543 
544  if (ret == MISSED) {
545  return (0);
546  }
547 
548  if (ret == FRONTFACE) {
549  GS_v3mult(dir, (float)dist);
550  GS_v3add(los[FROM], dir);
551  }
552 
553  return (1);
554 }
555 
556 /***********************************************************************/
557 /* DEBUG ****
558  void pr_plane(int pnum)
559  {
560  switch(pnum)
561  {
562  case 0:
563  fprintf(stderr,"top plane");
564 
565  break;
566  case 1:
567  fprintf(stderr,"bottom plane");
568 
569  break;
570  case 2:
571  fprintf(stderr,"left plane");
572 
573  break;
574  case 3:
575  fprintf(stderr,"right plane");
576 
577  break;
578  case 4:
579  fprintf(stderr,"front plane");
580 
581  break;
582  case 5:
583  fprintf(stderr,"back plane");
584 
585  break;
586  default:
587  fprintf(stderr,"clipping plane %d", 6 - pnum);
588 
589  break;
590  }
591 
592  return;
593  }
594  ******* */
#define NULL
Definition: ccmath.h:32
int G_debug(int, const char *,...) __attribute__((format(printf
void GS_v3mult(float *, float)
Multiple vectors.
Definition: gs_util.c:229
int gs_get_yrange(float *, float *)
Get y-range.
Definition: gs.c:1162
int GS_get_zrange(float *, float *, int)
Get z-extent for all loaded surfaces.
Definition: gs2.c:2685
int gsd_get_cplanes(Point4 *)
Get cplaces.
Definition: gsd_cplane.c:162
void GS_v3add(float *, float *)
Sum vectors.
Definition: gs_util.c:195
int segs_intersect(float, float, float, float, float, float, float, float, float *, float *)
Line intersect.
Definition: gsdrape.c:1202
void GS_v3eq(float *, float *)
Copy vector values.
Definition: gs_util.c:178
int gs_get_xrange(float *, float *)
Get x-range.
Definition: gs.c:1124
Point3 * gsdrape_get_allsegments(geosurf *, float *, float *, int *)
Get all segments.
Definition: gsdrape.c:399
int viewcell_tri_interp(geosurf *, typbuff *, Point3, int)
ADD.
Definition: gsdrape.c:503
typbuff * gs_get_att_typbuff(geosurf *, int, int)
Get attribute data buffer.
Definition: gs.c:681
int GS_v3dir(float *, float *, float *)
Get a normalized direction from v1 to v2, store in v3.
Definition: gs_util.c:351
float GS_distance(float *, float *)
Calculate distance.
Definition: gs_util.c:141
geosurf * gs_get_surf(int)
Get geosurf struct.
Definition: gs.c:63
void gs_get_databounds_planes(Point4 *planes)
Get data bounds for plane.
Definition: gs_query.c:469
#define BACKFACE
Definition: gs_query.c:35
#define MISSED
Definition: gs_query.c:33
int gs_los_intersect1(int surfid, float(*los)[3], float *point)
Crude method of intersecting line of sight with closest part of surface.
Definition: gs_query.c:52
#define FRONTFACE
Definition: gs_query.c:34
int gs_setlos_enterdata(Point3 *los)
Definition: gs_query.c:529
int gs_los_intersect(int surfid, float **los, float *point)
Crude method of intersecting line of sight with closest part of surface.
Definition: gs_query.c:192
#define HUGE_VAL
Values needed for Ray-Convex Polyhedron Intersection Test below originally by Eric Haines,...
Definition: gs_query.c:29
int RayCvxPolyhedronInt(Point3 org, Point3 dir, double tmax, Point4 *phdrn, int ph_num, double *tresult, int *pn)
Ray-Convex Polyhedron Intersection Test.
Definition: gs_query.c:384
#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
float Point4[4]
Definition: ogsf.h:204
#define DOT3(a, b)
Definition: ogsf.h:180
#define TO
Definition: ogsf.h:145
double b
Definition: r_raster.c:39
double t
Definition: r_raster.c:39
@ FROM
Definition: sqlp.tab.h:67
Definition: ogsf.h:256
float x_trans
Definition: ogsf.h:266
float z_trans
Definition: ogsf.h:266
float y_trans
Definition: ogsf.h:266
Definition: ogsf.h:208