Loading...
Searching...
No Matches
gs_query.c
Go to the documentation of this file.
/*!
\file lib/ogsf/gs_query.c
\brief OGSF library - query (lower level functions)
GRASS OpenGL gsurf OGSF Library
(C) 1999-2008 by the GRASS Development Team
This program is free software under the
GNU General Public License (>=v2).
Read the file COPYING that comes with GRASS
for details.
\author Bill Brown USACERL (January 1994)
\author Doxygenized by Martin Landa <landa.martin gmail.com> (May 2008)
*/
#include <math.h>
#include <grass/gis.h>
#include <grass/ogsf.h>
/*!
\brief Values needed for Ray-Convex Polyhedron Intersection Test below
originally by Eric Haines, erich@eye.com
*/
#ifndef HUGE_VAL
#define HUGE_VAL 1.7976931348623157e+308
#endif
/* return codes */
#define MISSED 0
#define FRONTFACE 1
#define BACKFACE -1
/* end Ray-Convex Polyhedron Intersection Test values */
/*!
\brief Crude method of intersecting line of sight with closest part of
surface.
Uses los vector to determine the point of first intersection
which is returned in point. Returns 0 if los doesn't intersect.
\param surfid surface id
\param los should be in surf-world coordinates
\param[out] point intersect point (real)
\return 0 on failure
\return 1 on success
*/
int gs_los_intersect1(int surfid, float (*los)[3], float *point)
{
float dx, dy, dz, u_d[3];
float a[3], incr, min_incr, tlen, len;
int outside, above, below, edge, istep;
float b[3];
geosurf *gs;
typbuff *buf;
G_debug(3, "gs_los_intersect1():");
if (NULL == (gs = gs_get_surf(surfid))) {
return (0);
}
if (0 == GS_v3dir(los[FROM], los[TO], u_d)) {
return (0);
}
buf = gs_get_att_typbuff(gs, ATT_TOPO, 0);
istep = edge = below = 0;
len = 0.0;
tlen = GS_distance(los[FROM], los[TO]);
incr = tlen / 1000.0;
min_incr = incr / 1000.0;
dx = incr * u_d[X];
dy = incr * u_d[Y];
dz = incr * u_d[Z];
a[X] = los[FROM][X];
a[Y] = los[FROM][Y];
a[Z] = los[FROM][Z];
b[X] = a[X] - gs->x_trans;
b[Y] = a[Y] - gs->y_trans;
if (viewcell_tri_interp(gs, buf, b, 0)) {
/* expects surface coords */
b[Z] += gs->z_trans;
if (a[Z] < b[Z]) {
/* viewing from below surface */
/* don't use this method
fprintf(stderr,"view from below\n");
below = 1;
*/
return (0);
}
}
while (incr > min_incr) {
outside = 0;
above = 0;
b[X] = a[X] - gs->x_trans;
b[Y] = a[Y] - gs->y_trans;
if (viewcell_tri_interp(gs, buf, b, 0)) {
/* ignores masks */
b[Z] += gs->z_trans;
above = (a[Z] > b[Z]);
}
else {
outside = 1;
if (istep > 10) {
edge = 1;
below = 1;
}
}
while (outside || above) {
a[X] += dx;
a[Y] += dy;
a[Z] += dz;
len += incr;
outside = 0;
above = 0;
b[X] = a[X] - gs->x_trans;
b[Y] = a[Y] - gs->y_trans;
if (viewcell_tri_interp(gs, buf, b, 0)) {
b[Z] += gs->z_trans;
above = (a[Z] > b[Z]);
}
else {
outside = 1;
}
if (len > tlen) {
return 0; /* over surface */ /* under surface */
}
}
/* could look for spikes here - see if any data points along
shadow of line on surf go above los */
/* back up several spots? */
a[X] -= (1.0 * dx);
a[Y] -= (1.0 * dy);
a[Z] -= (1.0 * dz);
incr /= 2.0;
++istep;
dx = incr * u_d[X];
dy = incr * u_d[Y];
dz = incr * u_d[Z];
}
if ((edge) && (b[Z] - (a[Z] + dz * 2.0) > incr * u_d[Z])) {
G_debug(3, " looking under surface");
return 0;
}
point[X] = b[X];
point[Y] = b[Y];
point[Z] = b[Z] - gs->z_trans;
return (1);
}
/*!
\brief Crude method of intersecting line of sight with closest part of
surface.
This version uses the shadow of the los projected down to
the surface to generate a line_on_surf, then follows each
point in that line until the los intersects it.
\param surfid surface id
\param los should be in surf-world coordinates
\param[out] point intersect point (real)
\return 0 on failure
\return 1 on success
*/
int gs_los_intersect(int surfid, float **los, float *point)
{
double incr;
float p1, p2, u_d[3];
int above, ret, num, i, usedx;
float a[3], b[3];
float bgn[3], end[3], a1[3];
geosurf *gs;
typbuff *buf;
Point3 *points;
G_debug(3, "gs_los_intersect");
if (NULL == (gs = gs_get_surf(surfid))) {
return (0);
}
if (0 == GS_v3dir(los[FROM], los[TO], u_d)) {
return (0);
}
buf = gs_get_att_typbuff(gs, ATT_TOPO, 0);
GS_v3eq(bgn, los[FROM]);
GS_v3eq(end, los[TO]);
bgn[X] -= gs->x_trans;
bgn[Y] -= gs->y_trans;
end[X] -= gs->x_trans;
end[Y] -= gs->y_trans;
/* trans? */
points = gsdrape_get_allsegments(gs, bgn, end, &num);
/* DEBUG
{
float t1[3], t2[3];
t1[X] = los[FROM][X] ;
t1[Y] = los[FROM][Y] ;
t2[X] = los[TO][X] ;
t2[Y] = los[TO][Y] ;
GS_set_draw(GSD_FRONT);
gsd_pushmatrix();
gsd_do_scale(1);
gsd_translate(gs->x_trans, gs->y_trans, gs->z_trans);
gsd_linewidth(1);
gsd_color_func(GS_default_draw_color());
gsd_line_onsurf(gs, t1, t2);
gsd_popmatrix();
GS_set_draw(GSD_BACK);
gsd_flush();
}
fprintf(stderr,"%d points to check\n", num);
fprintf(stderr,"point0 = %.6lf %.6lf %.6lf FT =%.6lf %.6lf %.6lf\n",
points[0][X],points[0][Y],points[0][Z],
los[FROM][X],los[FROM][Y],los[FROM][Z]);
fprintf(stderr,"incr1 = %.6lf: %.6lf %.6lf
%.6lf\n",incr,u_d[X],u_d[Y],u_d[Z]); fprintf(stderr,"first point below
surf\n"); fprintf(stderr,"incr2 = %f\n", (float)incr);
fprintf(stderr,"(%d/%d) %f > %f\n", i,num, a[Z], points[i][Z]);
fprintf(stderr,"incr3 = %f\n", (float)incr);
fprintf(stderr,"all points above surf\n");
*/
if (num < 2) {
G_debug(3, " %d points to check", num);
return (0);
}
/* use larger of deltas for better precision */
usedx = (fabs(u_d[X]) > fabs(u_d[Y]));
if (usedx) {
incr = ((points[0][X] - (los[FROM][X] - gs->x_trans)) / u_d[X]);
}
else if (u_d[Y]) {
incr = ((points[0][Y] - (los[FROM][Y] - gs->y_trans)) / u_d[Y]);
}
else {
point[X] = los[FROM][X] - gs->x_trans;
point[Y] = los[FROM][Y] - gs->y_trans;
return (viewcell_tri_interp(gs, buf, point, 1));
}
/* DEBUG
fprintf(stderr,"-----------------------------\n");
fprintf(stderr,"%d points to check\n", num);
fprintf(stderr,"incr1 = %.6lf: %.9f %.9f
%.9f\n",incr,u_d[X],u_d[Y],u_d[Z]); fprintf(stderr,
"\tpoint0 = %.6f %.6f %.6f\n\tFT = %.6f %.6f %.6f\n\tpoint%d = %.6f
%.6f\n", points[0][X],points[0][Y],points[0][Z],
los[FROM][X],los[FROM][Y],los[FROM][Z],
num-1, points[num-1][X],points[num-1][Y]);
*/
/* This should bring us right above (or below) the first point */
a[X] = los[FROM][X] + incr * u_d[X] - gs->x_trans;
a[Y] = los[FROM][Y] + incr * u_d[Y] - gs->y_trans;
a[Z] = los[FROM][Z] + incr * u_d[Z] - gs->z_trans;
if (a[Z] < points[0][Z]) {
/* viewing from below surface */
/* don't use this method */
/* DEBUG
fprintf(stderr,"first point below surf\n");
fprintf(stderr,"aZ= %.6f point0 = %.6f %.6f %.6f FT =%.6f %.6f
%.6f\n", a[Z], points[0][X],points[0][Y],points[0][Z],
los[FROM][X],los[FROM][Y],los[FROM][Z]);
*/
return (0);
}
GS_v3eq(a1, a);
GS_v3eq(b, a);
for (i = 1; i < num; i++) {
if (usedx) {
incr = ((points[i][X] - a1[X]) / u_d[X]);
}
else {
incr = ((points[i][Y] - a1[Y]) / u_d[Y]);
}
a[X] = a1[X] + (incr * u_d[X]);
a[Y] = a1[Y] + (incr * u_d[Y]);
a[Z] = a1[Z] + (incr * u_d[Z]);
above = (a[Z] >= points[i][Z]);
if (above) {
GS_v3eq(b, a);
continue;
}
/*
* Now we know b[Z] is above points[i-1]
* and a[Z] is below points[i]
* Since there should only be one polygon along this seg,
* just interpolate to intersect
*/
if (usedx) {
incr = ((a[X] - b[X]) / u_d[X]);
}
else {
incr = ((a[Y] - b[Y]) / u_d[Y]);
}
if (1 == (ret = segs_intersect(1.0, points[i][Z], 0.0, points[i - 1][Z],
1.0, a[Z], 0.0, b[Z], &p1, &p2))) {
point[X] = points[i - 1][X] + (u_d[X] * incr * p1);
point[Y] = points[i - 1][Y] + (u_d[Y] * incr * p1);
point[Z] = p2;
return (1);
}
G_debug(3, " line of sight error %d", ret);
return 0;
}
/* over surface */
return 0;
}
/*!
\brief Ray-Convex Polyhedron Intersection Test
Originally by Eric Haines, erich@eye.com
This test checks the ray against each face of a polyhedron, checking whether
the set of intersection points found for each ray-plane intersection
overlaps the previous intersection results. If there is no overlap (i.e.
no line segment along the ray that is inside the polyhedron), then the
ray misses and returns 0; else 1 is returned if the ray is entering the
polyhedron, -1 if the ray originates inside the polyhedron. If there is
an intersection, the distance and the number of the face hit is returned.
\param org,dir origin and direction of ray
\param tmax maximum useful distance along ray
\param phdrn list of planes in convex polyhedron
\param ph_num number of planes in convex polyhedron
\param[out] tresult distance of intersection along ray
\param[out] pn number of face hit (0 to ph_num-1)
\return FACE code
*/
int RayCvxPolyhedronInt(Point3 org, Point3 dir, double tmax, Point4 *phdrn,
int ph_num, double *tresult, int *pn)
{
double tnear, tfar, t, vn, vd;
int fnorm_num, bnorm_num; /* front/back face # hit */
tnear = -HUGE_VAL;
tfar = tmax;
/* Test each plane in polyhedron */
for (; ph_num--;) {
/* Compute intersection point T and sidedness */
vd = DOT3(dir, phdrn[ph_num]);
vn = DOT3(org, phdrn[ph_num]) + phdrn[ph_num][W];
if (vd == 0.0) {
/* ray is parallel to plane - check if ray origin is inside plane's
half-space */
if (vn > 0.0) {
/* ray origin is outside half-space */
return (MISSED);
}
}
else {
/* ray not parallel - get distance to plane */
t = -vn / vd;
if (vd < 0.0) {
/* front face - T is a near point */
if (t > tfar) {
return (MISSED);
}
if (t > tnear) {
/* hit near face, update normal */
fnorm_num = ph_num;
tnear = t;
}
}
else {
/* back face - T is a far point */
if (t < tnear) {
return (MISSED);
}
if (t < tfar) {
/* hit far face, update normal */
bnorm_num = ph_num;
tfar = t;
}
}
}
}
/* survived all tests */
/* Note: if ray originates on polyhedron, may want to change 0.0 to some
* epsilon to avoid intersecting the originating face.
*/
if (tnear >= 0.0) {
/* outside, hitting front face */
*tresult = tnear;
*pn = fnorm_num;
return (FRONTFACE);
}
else {
if (tfar < tmax) {
/* inside, hitting back face */
*tresult = tfar;
*pn = bnorm_num;
return (BACKFACE);
}
else {
/* inside, but back face beyond tmax */
return (MISSED);
}
}
}
/*!
\brief Get data bounds for plane
\param[out] planes
*/
void gs_get_databounds_planes(Point4 *planes)
{
float n, s, w, e, b, t;
Point3 tlfront, brback;
GS_get_zrange(&b, &t, 0);
gs_get_xrange(&w, &e);
gs_get_yrange(&s, &n);
tlfront[X] = tlfront[Y] = 0.0;
tlfront[Z] = t;
brback[X] = e - w;
brback[Y] = n - s;
brback[Z] = b;
/* top */
planes[0][X] = planes[0][Y] = 0.0;
planes[0][Z] = 1.0;
planes[0][W] = -(DOT3(planes[0], tlfront));
/* bottom */
planes[1][X] = planes[1][Y] = 0.0;
planes[1][Z] = -1.0;
planes[1][W] = -(DOT3(planes[1], brback));
/* left */
planes[2][Y] = planes[2][Z] = 0.0;
planes[2][X] = -1.0;
planes[2][W] = -(DOT3(planes[2], tlfront));
/* right */
planes[3][Y] = planes[3][Z] = 0.0;
planes[3][X] = 1.0;
planes[3][W] = -(DOT3(planes[3], brback));
/* front */
planes[4][X] = planes[4][Z] = 0.0;
planes[4][Y] = -1.0;
planes[4][W] = -(DOT3(planes[4], tlfront));
/* back */
planes[5][X] = planes[5][Z] = 0.0;
planes[5][Y] = 1.0;
planes[5][W] = -(DOT3(planes[5], brback));
return;
}
/*!
Gets all current cutting planes & data bounding planes
Intersects los with resulting convex polyhedron, then replaces los[FROM] with
first point on ray inside data.
\param[out] los
\return 0 on failure
\return 1 on success
*/
int gs_setlos_enterdata(Point3 *los)
{
Point4 planes[12]; /* MAX_CPLANES + 6 - should define this */
Point3 dir;
double dist, maxdist;
int num, ret, retp; /* might want to tell if retp is a clipping plane */
gs_get_databounds_planes(planes);
num = gsd_get_cplanes(planes + 6);
GS_v3dir(los[FROM], los[TO], dir);
maxdist = GS_distance(los[FROM], los[TO]);
ret = RayCvxPolyhedronInt(los[0], dir, maxdist, planes, num + 6, &dist,
&retp);
if (ret == MISSED) {
return (0);
}
if (ret == FRONTFACE) {
GS_v3mult(dir, (float)dist);
GS_v3add(los[FROM], dir);
}
return (1);
}
/***********************************************************************/
/* DEBUG ****
void pr_plane(int pnum)
{
switch(pnum)
{
case 0:
fprintf(stderr,"top plane");
break;
case 1:
fprintf(stderr,"bottom plane");
break;
case 2:
fprintf(stderr,"left plane");
break;
case 3:
fprintf(stderr,"right plane");
break;
case 4:
fprintf(stderr,"front plane");
break;
case 5:
fprintf(stderr,"back plane");
break;
default:
fprintf(stderr,"clipping plane %d", 6 - pnum);
break;
}
return;
}
******* */
NULL
#define NULL
Definition ccmath.h:32
G_debug
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition debug.c:66
b
double b
Definition driver/set_window.c:5
t
double t
Definition driver/set_window.c:5
GS_get_zrange
int GS_get_zrange(float *min, float *max, int doexag)
Get z-extent for all loaded surfaces.
Definition gs2.c:2688
gs_get_xrange
int gs_get_xrange(float *min, float *max)
Get x-range.
Definition gs.c:1124
gs_get_surf
geosurf * gs_get_surf(int id)
Get geosurf struct.
Definition gs.c:63
gs_get_yrange
int gs_get_yrange(float *min, float *max)
Get y-range.
Definition gs.c:1162
gs_get_att_typbuff
typbuff * gs_get_att_typbuff(geosurf *gs, int desc, int to_write)
Get attribute data buffer.
Definition gs.c:681
gs_get_databounds_planes
void gs_get_databounds_planes(Point4 *planes)
Get data bounds for plane.
Definition gs_query.c:469
BACKFACE
#define BACKFACE
Definition gs_query.c:35
MISSED
#define MISSED
Definition gs_query.c:33
gs_los_intersect1
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
FRONTFACE
#define FRONTFACE
Definition gs_query.c:34
gs_setlos_enterdata
int gs_setlos_enterdata(Point3 *los)
Definition gs_query.c:529
gs_los_intersect
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
HUGE_VAL
#define HUGE_VAL
Values needed for Ray-Convex Polyhedron Intersection Test below originally by Eric Haines,...
Definition gs_query.c:29
RayCvxPolyhedronInt
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
GS_v3mult
void GS_v3mult(float *v1, float k)
Multiple vectors.
Definition gs_util.c:229
GS_distance
float GS_distance(float *from, float *to)
Calculate distance.
Definition gs_util.c:141
GS_v3add
void GS_v3add(float *v1, float *v2)
Sum vectors.
Definition gs_util.c:195
GS_v3eq
void GS_v3eq(float *v1, float *v2)
Copy vector values.
Definition gs_util.c:178
GS_v3dir
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
gsd_get_cplanes
int gsd_get_cplanes(Point4 *planes)
Get cplaces.
Definition gsd_cplane.c:162
viewcell_tri_interp
int viewcell_tri_interp(geosurf *gs, typbuff *buf, Point3 pt, int check_mask)
ADD.
Definition gsdrape.c:509
segs_intersect
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
gsdrape_get_allsegments
Point3 * gsdrape_get_allsegments(geosurf *gs, float *bgn, float *end, int *num)
Get all segments.
Definition gsdrape.c:400
X
#define X(j)
Y
#define Y(j)