21 #include <grass/gis.h>
22 #include <grass/gstypes.h>
29 #define HUGE_VAL 1.7976931348623157e+308
54 float dx, dy, dz, u_d[3];
55 float a[3], incr, min_incr, tlen, len;
56 int outside, above, below, edge, istep;
61 G_debug(3,
"gs_los_intersect1():");
73 istep = edge = below = 0;
79 min_incr = incr / 1000.0;
89 b[
X] = a[
X] - gs->x_trans;
90 b[
Y] = a[
Y] - gs->y_trans;
107 while (incr > min_incr) {
110 b[
X] = a[
X] - gs->x_trans;
111 b[
Y] = a[
Y] - gs->y_trans;
116 above = (a[Z] > b[Z]);
127 while (outside || above) {
134 b[
X] = a[
X] - gs->x_trans;
135 b[
Y] = a[
Y] - gs->y_trans;
139 above = (a[Z] > b[Z]);
164 if ((edge) && (b[Z] - (a[Z] + dz * 2.0) > incr * u_d[Z])) {
165 G_debug(3,
" looking under surface");
172 point[Z] = b[Z] - gs->z_trans;
194 float p1, p2, u_d[3];
195 int above, ret, num, i, usedx;
197 float bgn[3], end[3], a1[3];
202 G_debug(3,
"gs_los_intersect");
217 bgn[
X] -= gs->x_trans;
218 bgn[
Y] -= gs->y_trans;
220 end[
X] -= gs->x_trans;
221 end[
Y] -= gs->y_trans;
260 G_debug(3,
" %d points to check", num);
266 usedx = (fabs(u_d[
X]) > fabs(u_d[
Y]));
268 incr = ((points[0][
X] - (los[
FROM][
X] - gs->x_trans)) / u_d[
X]);
271 incr = ((points[0][
Y] - (los[
FROM][
Y] - gs->y_trans)) / u_d[
Y]);
274 point[
X] = los[
FROM][
X] - gs->x_trans;
275 point[
Y] = los[
FROM][
Y] - gs->y_trans;
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;
296 if (a[Z] < points[0][Z]) {
311 for (i = 1; i < num; i++) {
313 incr = ((points[i][
X] - a1[
X]) / u_d[X]);
316 incr = ((points[i][
Y] - a1[
Y]) / u_d[Y]);
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]);
337 incr = ((a[
X] - b[
X]) / u_d[X]);
340 incr = ((a[
Y] - b[
Y]) / u_d[Y]);
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);
353 G_debug(3,
" line of sight error %d", ret);
385 int ph_num,
double *tresult,
int *pn)
387 double tnear, tfar, t, vn, vd;
388 int fnorm_num, bnorm_num;
396 vd = DOT3(dir, phdrn[ph_num]);
397 vn = DOT3(org, phdrn[ph_num]) + phdrn[ph_num][W];
471 float n,
s,
w, e,
b, t;
472 Point3 tlfront, brback;
478 tlfront[
X] = tlfront[
Y] = 0.0;
486 planes[0][
X] = planes[0][
Y] = 0.0;
488 planes[0][W] = -(DOT3(planes[0], tlfront));
491 planes[1][
X] = planes[1][
Y] = 0.0;
493 planes[1][W] = -(DOT3(planes[1], brback));
496 planes[2][
Y] = planes[2][Z] = 0.0;
498 planes[2][W] = -(DOT3(planes[2], tlfront));
501 planes[3][
Y] = planes[3][Z] = 0.0;
503 planes[3][W] = -(DOT3(planes[3], brback));
506 planes[4][
X] = planes[4][Z] = 0.0;
508 planes[4][W] = -(DOT3(planes[4], tlfront));
511 planes[5][
X] = planes[5][Z] = 0.0;
513 planes[5][W] = -(DOT3(planes[5], brback));
533 double dist, maxdist;
542 planes, num + 6, &dist, &retp);