GRASS GIS 8 Programmer's Manual 8.4.1(2025)-45ca3179ab
Loading...
Searching...
No Matches
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 */
52int 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 */
192int 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 */
384int 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 */
469void gs_get_databounds_planes(Point4 *planes)
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 */
529int gs_setlos_enterdata(Point3 *los)
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
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 level, const char *msg,...)
Print debugging message.
Definition debug.c:66
double b
double t
int GS_get_zrange(float *min, float *max, int doexag)
Get z-extent for all loaded surfaces.
Definition gs2.c:2685
int gs_get_xrange(float *min, float *max)
Get x-range.
Definition gs.c:1124
geosurf * gs_get_surf(int id)
Get geosurf struct.
Definition gs.c:63
int gs_get_yrange(float *min, float *max)
Get y-range.
Definition gs.c:1162
typbuff * gs_get_att_typbuff(geosurf *gs, int desc, int to_write)
Get attribute data buffer.
Definition gs.c:681
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
void GS_v3mult(float *v1, float k)
Multiple vectors.
Definition gs_util.c:229
float GS_distance(float *from, float *to)
Calculate distance.
Definition gs_util.c:141
void GS_v3add(float *v1, float *v2)
Sum vectors.
Definition gs_util.c:195
void GS_v3eq(float *v1, float *v2)
Copy vector values.
Definition gs_util.c:178
int GS_v3dir(float *v1, float *v2, float *v3)
Get a normalized direction from v1 to v2, store in v3.
Definition gs_util.c:351
int gsd_get_cplanes(Point4 *planes)
Get cplaces.
Definition gsd_cplane.c:162
int viewcell_tri_interp(geosurf *gs, typbuff *buf, Point3 pt, int check_mask)
ADD.
Definition gsdrape.c:503
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:1202
Point3 * gsdrape_get_allsegments(geosurf *gs, float *bgn, float *end, int *num)
Get all segments.
Definition gsdrape.c:399
#define X(j)
#define Y(j)