GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
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 /*!
40  \brief Crude method of intersecting line of sight with closest part of 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 surface.
179 
180  This version uses the shadow of the los projected down to
181  the surface to generate a line_on_surf, then follows each
182  point in that line until the los intersects it.
183 
184  \param surfid surface id
185  \param los should be in surf-world coordinates
186  \param[out] point intersect point (real)
187 
188  \return 0 on failure
189  \return 1 on success
190  */
191 int gs_los_intersect(int surfid, float **los, float *point)
192 {
193  double incr;
194  float p1, p2, u_d[3];
195  int above, ret, num, i, usedx;
196  float a[3], b[3];
197  float bgn[3], end[3], a1[3];
198  geosurf *gs;
199  typbuff *buf;
200  Point3 *points;
201 
202  G_debug(3, "gs_los_intersect");
203 
204  if (NULL == (gs = gs_get_surf(surfid))) {
205  return (0);
206  }
207 
208  if (0 == GS_v3dir(los[FROM], los[TO], u_d)) {
209  return (0);
210  }
211 
212  buf = gs_get_att_typbuff(gs, ATT_TOPO, 0);
213 
214  GS_v3eq(bgn, los[FROM]);
215  GS_v3eq(end, los[TO]);
216 
217  bgn[X] -= gs->x_trans;
218  bgn[Y] -= gs->y_trans;
219 
220  end[X] -= gs->x_trans;
221  end[Y] -= gs->y_trans;
222 
223  /* trans? */
224  points = gsdrape_get_allsegments(gs, bgn, end, &num);
225 
226  /* DEBUG
227  {
228  float t1[3], t2[3];
229 
230  t1[X] = los[FROM][X] ;
231  t1[Y] = los[FROM][Y] ;
232 
233  t2[X] = los[TO][X] ;
234  t2[Y] = los[TO][Y] ;
235 
236  GS_set_draw(GSD_FRONT);
237  gsd_pushmatrix();
238  gsd_do_scale(1);
239  gsd_translate(gs->x_trans, gs->y_trans, gs->z_trans);
240  gsd_linewidth(1);
241  gsd_color_func(GS_default_draw_color());
242  gsd_line_onsurf(gs, t1, t2);
243  gsd_popmatrix();
244  GS_set_draw(GSD_BACK);
245  gsd_flush();
246  }
247  fprintf(stderr,"%d points to check\n", num);
248  fprintf(stderr,"point0 = %.6lf %.6lf %.6lf FT =%.6lf %.6lf %.6lf\n",
249  points[0][X],points[0][Y],points[0][Z],
250  los[FROM][X],los[FROM][Y],los[FROM][Z]);
251  fprintf(stderr,"incr1 = %.6lf: %.6lf %.6lf %.6lf\n",incr,u_d[X],u_d[Y],u_d[Z]);
252  fprintf(stderr,"first point below surf\n");
253  fprintf(stderr,"incr2 = %f\n", (float)incr);
254  fprintf(stderr,"(%d/%d) %f > %f\n", i,num, a[Z], points[i][Z]);
255  fprintf(stderr,"incr3 = %f\n", (float)incr);
256  fprintf(stderr,"all points above surf\n");
257  */
258 
259  if (num < 2) {
260  G_debug(3, " %d points to check", num);
261 
262  return (0);
263  }
264 
265  /* use larger of deltas for better precision */
266  usedx = (fabs(u_d[X]) > fabs(u_d[Y]));
267  if (usedx) {
268  incr = ((points[0][X] - (los[FROM][X] - gs->x_trans)) / u_d[X]);
269  }
270  else if (u_d[Y]) {
271  incr = ((points[0][Y] - (los[FROM][Y] - gs->y_trans)) / u_d[Y]);
272  }
273  else {
274  point[X] = los[FROM][X] - gs->x_trans;
275  point[Y] = los[FROM][Y] - gs->y_trans;
276 
277  return (viewcell_tri_interp(gs, buf, point, 1));
278  }
279 
280  /* DEBUG
281  fprintf(stderr,"-----------------------------\n");
282  fprintf(stderr,"%d points to check\n", num);
283  fprintf(stderr,"incr1 = %.6lf: %.9f %.9f %.9f\n",incr,u_d[X],u_d[Y],u_d[Z]);
284  fprintf(stderr,
285  "\tpoint0 = %.6f %.6f %.6f\n\tFT = %.6f %.6f %.6f\n\tpoint%d = %.6f %.6f\n",
286  points[0][X],points[0][Y],points[0][Z],
287  los[FROM][X],los[FROM][Y],los[FROM][Z],
288  num-1, points[num-1][X],points[num-1][Y]);
289  */
290 
291  /* This should bring us right above (or below) the first point */
292  a[X] = los[FROM][X] + incr * u_d[X] - gs->x_trans;
293  a[Y] = los[FROM][Y] + incr * u_d[Y] - gs->y_trans;
294  a[Z] = los[FROM][Z] + incr * u_d[Z] - gs->z_trans;
295 
296  if (a[Z] < points[0][Z]) {
297  /* viewing from below surface */
298  /* don't use this method */
299  /* DEBUG
300  fprintf(stderr,"first point below surf\n");
301  fprintf(stderr,"aZ= %.6f point0 = %.6f %.6f %.6f FT =%.6f %.6f %.6f\n",
302  a[Z], points[0][X],points[0][Y],points[0][Z],
303  los[FROM][X],los[FROM][Y],los[FROM][Z]);
304  */
305  return (0);
306  }
307 
308  GS_v3eq(a1, a);
309  GS_v3eq(b, a);
310 
311  for (i = 1; i < num; i++) {
312  if (usedx) {
313  incr = ((points[i][X] - a1[X]) / u_d[X]);
314  }
315  else {
316  incr = ((points[i][Y] - a1[Y]) / u_d[Y]);
317  }
318 
319  a[X] = a1[X] + (incr * u_d[X]);
320  a[Y] = a1[Y] + (incr * u_d[Y]);
321  a[Z] = a1[Z] + (incr * u_d[Z]);
322  above = (a[Z] >= points[i][Z]);
323 
324  if (above) {
325  GS_v3eq(b, a);
326  continue;
327  }
328 
329  /*
330  * Now we know b[Z] is above points[i-1]
331  * and a[Z] is below points[i]
332  * Since there should only be one polygon along this seg,
333  * just interpolate to intersect
334  */
335 
336  if (usedx) {
337  incr = ((a[X] - b[X]) / u_d[X]);
338  }
339  else {
340  incr = ((a[Y] - b[Y]) / u_d[Y]);
341  }
342 
343  if (1 == (ret = segs_intersect(1.0, points[i][Z],
344  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 nunber 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,
542  planes, num + 6, &dist, &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  ******* */
float Point4[4]
Definition: ogsf.h:200
geosurf * gs_get_surf(int)
Get geosurf struct.
Definition: gs.c:62
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:191
#define ATT_TOPO
Definition: ogsf.h:73
int segs_intersect(float, float, float, float, float, float, float, float, float *, float *)
Line intersect.
Definition: gsdrape.c:1210
#define FRONTFACE
Definition: gs_query.c:34
void GS_v3add(float *, float *)
Sum vectors.
Definition: gs_util.c:195
Point3 * gsdrape_get_allsegments(geosurf *, float *, float *, int *)
Get all segments.
Definition: gsdrape.c:401
float Point3[3]
Definition: ogsf.h:201
#define NULL
Definition: ccmath.h:32
int gs_get_xrange(float *, float *)
Get x-range.
Definition: gs.c:1126
#define W
Definition: ogsf.h:140
float GS_distance(float *, float *)
Calculate distance.
Definition: gs_util.c:141
double t
Definition: r_raster.c:39
double b
Definition: r_raster.c:39
float x_trans
Definition: ogsf.h:267
int gs_setlos_enterdata(Point3 *los)
Definition: gs_query.c:529
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 DOT3(a, b)
Definition: ogsf.h:176
void gs_get_databounds_planes(Point4 *planes)
Get data bounds for plane.
Definition: gs_query.c:469
#define FROM
Definition: ogsf.h:141
#define Z
Definition: ogsf.h:139
int viewcell_tri_interp(geosurf *, typbuff *, Point3, int)
ADD.
Definition: gsdrape.c:507
float y_trans
Definition: ogsf.h:267
typbuff * gs_get_att_typbuff(geosurf *, int, int)
Get attribute data buffer.
Definition: gs.c:681
#define Y
Definition: ogsf.h:138
int gs_get_yrange(float *, float *)
Get y-range.
Definition: gs.c:1164
int gsd_get_cplanes(Point4 *)
Get cplaces.
Definition: gsd_cplane.c:162
void GS_v3eq(float *, float *)
Copy vector values.
Definition: gs_util.c:178
#define HUGE_VAL
Values needed for Ray-Convex Polyhedron Intersection Test below originally by Eric Haines...
Definition: gs_query.c:29
#define BACKFACE
Definition: gs_query.c:35
Definition: ogsf.h:204
#define TO
Definition: ogsf.h:142
#define MISSED
Definition: gs_query.c:33
int GS_v3dir(float *, float *, float *)
Get a normalized direction from v1 to v2, store in v3.
Definition: gs_util.c:353
#define X
Definition: ogsf.h:137
int GS_get_zrange(float *, float *, int)
Get z-extent for all loaded surfaces.
Definition: gs2.c:2689
float z_trans
Definition: ogsf.h:267
int G_debug(int, const char *,...) __attribute__((format(printf
Definition: ogsf.h:257
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
void GS_v3mult(float *, float)
Multiple vectors.
Definition: gs_util.c:229