GRASS GIS 7 Programmer's Manual  7.7.svn(2019)-r74095
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)) {
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);
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
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
int GS_get_zrange(float *min, float *max, int doexag)
Get z-extent for all loaded surfaces.
Definition: gs2.c:2689
#define ATT_TOPO
Definition: ogsf.h:73
int GS_v3dir(float *v1, float *v2, float *v3)
Get a normalized direction from v1 to v2, store in v3.
Definition: gs_util.c:353
#define FRONTFACE
Definition: gs_query.c:34
int viewcell_tri_interp(geosurf *gs, typbuff *buf, Point3 pt, int check_mask)
Definition: gsdrape.c:507
Sum vectors.
Definition: gs_util.c:195
float Point3[3]
Definition: ogsf.h:201
int gs_get_yrange(float *min, float *max)
Get y-range.
Definition: gs.c:1164
#define NULL
Definition: ccmath.h:32
#define W
Definition: ogsf.h:140
double t
Definition: r_raster.c:39
typbuff * gs_get_att_typbuff(geosurf *gs, int desc, int to_write)
Get attribute data buffer.
Definition: gs.c:681
double b
Definition: r_raster.c:39
float x_trans
Definition: ogsf.h:267
void GS_v3mult(float *v1, float k)
Multiple vectors.
Definition: gs_util.c:229
int gs_setlos_enterdata(Point3 *los)
Definition: gs_query.c:529
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: debug.c:65
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
Point3 * gsdrape_get_allsegments(geosurf *gs, float *bgn, float *end, int *num)
Get all segments.
Definition: gsdrape.c:401
#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 gs_get_xrange(float *min, float *max)
Get x-range.
Definition: gs.c:1126
float y_trans
Definition: ogsf.h:267
#define Y
Definition: ogsf.h:138
geosurf * gs_get_surf(int id)
Get geosurf struct.
Definition: gs.c:62
#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
void GS_v3eq(float *v1, float *v2)
Copy vector values.
Definition: gs_util.c:178
int gsd_get_cplanes(Point4 *planes)
Get cplaces.
Definition: gsd_cplane.c:162
#define TO
Definition: ogsf.h:142
#define MISSED
Definition: gs_query.c:33
#define X
Definition: ogsf.h:137
long num
Definition: raster3d/cats.c:85
float GS_distance(float *from, float *to)
Calculate distance.
Definition: gs_util.c:141
float z_trans
Definition: ogsf.h:267
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
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