GRASS GIS 8 Programmer's Manual  8.5.0dev(2024)-826ab62d78
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);
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)
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