53#include <grass/ogsf.h>
54#include <grass/glocale.h>
60#define DONT_INTERSECT 0
64#define LERP(a, l, h) ((l) + (((h) - (l)) * (a)))
65#define EQUAL(a, b) (fabs((a) - (b)) < EPSILON)
66#define ISNODE(p, res) (fmod((double)p, (double)res) < EPSILON)
68#define SAME_SIGNS(a, b) ((a >= 0 && b >= 0) || (a < 0 && b < 0))
70static int drape_line_init(
int,
int);
71static Point3 *_gsdrape_get_segments(geosurf *,
float *,
float *,
int *);
72static float dist_squared_2d(
float *,
float *);
81static Point3 *Vi, *Hi, *Di;
95static int drape_line_init(
int rows,
int cols)
98 if (
NULL == (I3d = (Point3 *)calloc(2 * (rows + cols),
sizeof(Point3)))) {
102 if (
NULL == (Vi = (Point3 *)calloc(cols,
sizeof(Point3)))) {
108 if (
NULL == (Hi = (Point3 *)calloc(rows,
sizeof(Point3)))) {
115 if (
NULL == (Di = (Point3 *)calloc(rows + cols,
sizeof(Point3)))) {
136static Point3 *_gsdrape_get_segments(geosurf *gs,
float *bgn,
float *end,
141 float dir[2], yres, xres;
157 if (!((end[
Y] - bgn[
Y]) / (end[
X] - bgn[
X]) == yres / xres)) {
166 G_debug(5,
"_gsdrape_get_segments vi=%d, hi=%d, di=%d, num=%d", vi, hi, di,
180static float dist_squared_2d(
float *p1,
float *p2)
187 return (dx * dx + dy * dy);
200 static int first = 1;
205 if (0 > drape_line_init(gs->rows, gs->cols)) {
206 G_warning(_(
"Unable to process vector map - out of memory"));
234 float *replace, xl, yb, xr, yt, xi, yi;
279 float pt1[2], pt2[2];
363 I3d[0][Z] = gs->att[ATT_TOPO].constant;
366 I3d[1][Z] = gs->att[ATT_TOPO].constant;
372 if (bgn[
X] == end[
X] && bgn[
Y] == end[
Y]) {
386 return (_gsdrape_get_segments(gs, bgn, end, num));
408 if (bgn[
X] == end[
X] && bgn[
Y] == end[
Y]) {
426 return (_gsdrape_get_segments(gs, bgn, end, num));
447 f[Z] =
l[Z] = gs->att[ATT_TOPO].constant;
506 int offset, drow, dcol, vrow, vcol;
507 float xmax, ymin, ymax, alpha;
519 if (pt[
X] < 0.0 || pt[
Y] > ymax) {
524 if (pt[
Y] < ymin || pt[
X] > xmax) {
530 pt[Z] = gs->att[ATT_TOPO].constant;
543 if (pt[
X] > 0.0 && pt[
Y] < ymax) {
549 offset =
DRC2OFF(gs, drow, dcol);
556 offset =
DRC2OFF(gs, drow, dcol);
565 offset =
DRC2OFF(gs, drow, dcol);
574 offset =
DRC2OFF(gs, drow, dcol);
580 else if (pt[
X] == 0.0) {
593 pt[Z] =
LERP(alpha, p1[Z], p2[Z]);
602 else if (pt[
Y] == gs->yrange) {
612 pt[Z] =
LERP(alpha, p1[Z], p2[Z]);
617 else if (vrow ==
VROWS(gs)) {
621 if (pt[
X] > 0.0 && pt[
X] < xmax) {
625 offset =
DRC2OFF(gs, drow, dcol);
629 offset =
DRC2OFF(gs, drow, dcol);
633 pt[Z] =
LERP(alpha, p1[Z], p2[Z]);
637 else if (pt[
X] == 0.0) {
647 offset =
DRC2OFF(gs, drow, dcol);
660 offset =
DRC2OFF(gs, drow, dcol);
664 offset =
DRC2OFF(gs, drow, dcol);
668 pt[Z] =
LERP(alpha, p1[Z], p2[Z]);
693 if (pt[
X] >= 0.0 && pt[
Y] <= gs->yrange) {
727 int num, i, found, cv, ch, cd, cnum;
728 float dv, dh, dd, big, cpoint[2];
730 cv = ch = cd = cnum = 0;
733 cpoint[
X] = first[
X];
734 cpoint[
Y] = first[
Y];
737 I3d[cnum][
X] = first[
X];
738 I3d[cnum][
Y] = first[
Y];
739 I3d[cnum][Z] = first[Z];
744 big = gs->yrange * gs->yrange + gs->xrange * gs->xrange;
747 for (i = 0; i < num; i = cv + ch + cd) {
749 dv = dist_squared_2d(Vi[cv], cpoint);
761 dh = dist_squared_2d(Hi[ch], cpoint);
773 dd = dist_squared_2d(Di[cd], cpoint);
787 if (dd <= dv && dd <= dh) {
789 cpoint[
X] = I3d[cnum][
X] = Di[cd][
X];
790 cpoint[
Y] = I3d[cnum][
Y] = Di[cd][
Y];
791 I3d[cnum][Z] = Di[cd][Z];
810 cpoint[
X] = I3d[cnum][
X] = Vi[cv][
X];
811 cpoint[
Y] = I3d[cnum][
Y] = Vi[cv][
Y];
812 I3d[cnum][Z] = Vi[cv][Z];
826 cpoint[
X] = I3d[cnum][
X] = Hi[ch][
X];
827 cpoint[
Y] = I3d[cnum][
Y] = Hi[ch][
Y];
828 I3d[cnum][Z] = Hi[ch][Z];
834 if (i == cv + ch + cd) {
835 G_debug(5,
"order_intersects(): stuck on %d", cnum);
836 G_debug(5,
"order_intersects(): cv = %d, ch = %d, cd = %d", cv, ch,
838 G_debug(5,
"order_intersects(): dv = %f, dh = %f, dd = %f", dv, dh,
851 I3d[cnum][
X] = last[
X];
852 I3d[cnum][
Y] = last[
Y];
853 I3d[cnum][Z] = last[Z];
879 int fcol, lcol, incr, hits, num, offset, drow1, drow2;
880 float xl, yb, xr, yt, z1, z2, alpha;
882 int bgncol, endcol, cols, rows;
891 if (bgncol > cols && endcol > cols) {
895 if (bgncol == endcol) {
899 fcol = dir[
X] > 0 ? bgncol + 1 : bgncol;
900 lcol = dir[
X] > 0 ? endcol : endcol + 1;
903 incr = lcol - fcol > 0 ? 1 : -1;
905 while (fcol > cols || fcol < 0) {
909 while (lcol > cols || lcol < 0) {
913 num = abs(lcol - fcol) + 1;
915 yb = gs->yrange - (yres * rows) -
EPSILON;
918 for (hits = 0; hits < num; hits++) {
919 xl = xr =
VCOL2X(gs, fcol);
928 Vi[hits][Z] = gs->att[ATT_TOPO].constant;
931 drow1 =
Y2VROW(gs, Vi[hits][
Y]) * gs->y_mod;
932 drow2 = (1 +
Y2VROW(gs, Vi[hits][
Y])) * gs->y_mod;
934 if (drow2 >= gs->rows) {
935 drow2 = gs->rows - 1;
938 alpha = ((gs->yrange - drow1 * gs->yres) - Vi[hits][
Y]) / yres;
940 offset =
DRC2OFF(gs, drow1, fcol * gs->x_mod);
942 offset =
DRC2OFF(gs, drow2, fcol * gs->x_mod);
944 Vi[hits][Z] =
LERP(alpha, z1, z2);
973 int frow, lrow, incr, hits, num, offset, dcol1, dcol2;
974 float xl, yb, xr, yt, z1, z2, alpha;
976 int bgnrow, endrow, rows, cols;
984 if (bgnrow == endrow) {
988 if (bgnrow > rows && endrow > rows) {
992 frow = dir[
Y] > 0 ? bgnrow : bgnrow + 1;
993 lrow = dir[
Y] > 0 ? endrow + 1 : endrow;
996 incr = lrow - frow > 0 ? 1 : -1;
998 while (frow > rows || frow < 0) {
1002 while (lrow > rows || lrow < 0) {
1006 num = abs(lrow - frow) + 1;
1011 for (hits = 0; hits < num; hits++) {
1012 yb = yt =
VROW2Y(gs, frow);
1021 Hi[hits][Z] = gs->att[ATT_TOPO].constant;
1024 dcol1 =
X2VCOL(gs, Hi[hits][
X]) * gs->x_mod;
1025 dcol2 = (1 +
X2VCOL(gs, Hi[hits][
X])) * gs->x_mod;
1027 if (dcol2 >= gs->cols) {
1028 dcol2 = gs->cols - 1;
1031 alpha = (Hi[hits][
X] - (dcol1 * gs->xres)) / xres;
1033 offset =
DRC2OFF(gs, frow * gs->y_mod, dcol1);
1035 offset =
DRC2OFF(gs, frow * gs->y_mod, dcol2);
1037 Hi[hits][Z] =
LERP(alpha, z1, z2);
1068 int fdig, ldig, incr, hits, num, offset;
1069 int vrow, vcol, drow1, drow2, dcol1, dcol2;
1070 float xl, yb, xr, yt, z1, z2, alpha;
1071 float xres, yres, xi, yi, dx, dy;
1072 int diags, cols, rows, lower;
1079 diags = rows + cols;
1086 lower = ((end[
X] - pt[
X]) / xres > (end[
Y] - pt[
Y]) / yres);
1087 ldig = lower ? vrow + vcol + 1 : vrow + vcol;
1094 lower = ((bgn[
X] - pt[
X]) / xres > (bgn[
Y] - pt[
Y]) / yres);
1095 fdig = lower ? vrow + vcol + 1 : vrow + vcol;
1106 incr = ldig - fdig > 0 ? 1 : -1;
1108 while (fdig > diags || fdig < 0) {
1112 while (ldig > diags || ldig < 0) {
1116 num = abs(ldig - fdig) + 1;
1118 for (hits = 0; hits < num; hits++) {
1119 yb = gs->yrange - (yres * (fdig < rows ? fdig : rows)) -
EPSILON;
1121 yt = gs->yrange - (yres * (fdig < cols ? 0 : fdig - cols)) +
EPSILON;
1137 drow1 =
Y2VROW(gs, Di[hits][
Y]) * gs->y_mod;
1138 drow2 = (1 +
Y2VROW(gs, Di[hits][
Y])) * gs->y_mod;
1140 if (drow2 >= gs->rows) {
1141 drow2 = gs->rows - 1;
1146 Di[hits][Z] = gs->att[ATT_TOPO].constant;
1149 dcol1 =
X2VCOL(gs, Di[hits][
X]) * gs->x_mod;
1150 dcol2 = (1 +
X2VCOL(gs, Di[hits][
X])) * gs->x_mod;
1152 if (dcol2 >= gs->cols) {
1153 dcol2 = gs->cols - 1;
1156 dx =
DCOL2X(gs, dcol2) - Di[hits][
X];
1157 dy =
DROW2Y(gs, drow1) - Di[hits][
Y];
1159 sqrt(dx * dx + dy * dy) / sqrt(xres * xres + yres * yres);
1161 offset =
DRC2OFF(gs, drow1, dcol2);
1163 offset =
DRC2OFF(gs, drow2, dcol1);
1165 Di[hits][Z] =
LERP(alpha, z1, z2);
1203 float x4,
float y4,
float *
x,
float *y)
1205 float a1, a2, b1, b2, c1, c2;
1206 float r1, r2, r3, r4;
1214 c1 = x2 * y1 - x1 * y2;
1218 r3 = a1 * x3 + b1 * y3 + c1;
1219 r4 = a1 * x4 + b1 * y4 + c1;
1232 c2 = x4 * y3 - x3 * y4;
1235 r1 = a2 * x1 + b2 * y1 + c2;
1236 r2 = a2 * x2 + b2 * y2 + c2;
1249 denom = a1 * b2 - a2 * b1;
1261 num = b1 * c2 - b2 * c1;
1265 num = a2 * c1 - a1 * c2;
1314 intersect[Z] = (plane[
X] *
x + plane[
Y] * y + plane[W]) / -plane[Z];
1329 Point3 v1, v2, norm;
1331 v1[
X] = p1[
X] - p3[
X];
1332 v1[
Y] = p1[
Y] - p3[
Y];
1333 v1[Z] = p1[Z] - p3[Z];
1335 v2[
X] = p2[
X] - p3[
X];
1336 v2[
Y] = p2[
Y] - p3[
Y];
1337 v2[Z] = p2[Z] - p3[Z];
1344 plane[W] = -p3[
X] * norm[
X] - p3[
Y] * norm[
Y] - p3[Z] * norm[Z];
1358 c[
X] = (a[
Y] *
b[Z]) - (a[Z] *
b[
Y]);
1359 c[
Y] = (a[Z] *
b[
X]) - (a[
X] *
b[Z]);
1360 c[Z] = (a[
X] *
b[
Y]) - (a[
Y] *
b[
X]);
void G_free(void *buf)
Free allocated memory.
int G_debug(int level, const char *msg,...)
Print debugging message.
void G_warning(const char *msg,...)
Print a warning message to stderr.
typbuff * gs_get_att_typbuff(geosurf *gs, int desc, int to_write)
Get attribute data buffer.
int gs_point_is_masked(geosurf *gs, float *pt)
Check if point is masked.
int gs_get_att_src(geosurf *gs, int desc)
Get attribute source.
void GS_v2dir(float *v1, float *v2, float *v3)
Get a normalized direction from v1 to v2, store in v3 (2D)
void GS_v3eq(float *v1, float *v2)
Copy vector values.
float GS_P2distance(float *from, float *to)
Calculate distance in plane.
int P3toPlane(Point3 p1, Point3 p2, Point3 p3, float *plane)
Define plane.
Point3 * gsdrape_get_segments(geosurf *gs, float *bgn, float *end, int *num)
ADD.
int in_vregion(geosurf *gs, float *pt)
ADD.
int order_intersects(geosurf *gs, Point3 first, Point3 last, int vi, int hi, int di)
ADD.
int _viewcell_tri_interp(geosurf *gs, Point3 pt)
ADD.
int gsdrape_set_surface(geosurf *gs)
ADD.
void interp_first_last(geosurf *gs, float *bgn, float *end, Point3 f, Point3 l)
ADD.
int V3Cross(Point3 a, Point3 b, Point3 c)
Get cross product.
int viewcell_tri_interp(geosurf *gs, typbuff *buf, Point3 pt, int check_mask)
ADD.
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.
int get_horz_intersects(geosurf *gs, float *bgn, float *end, float *dir)
Get horizontal intersects.
int get_vert_intersects(geosurf *gs, float *bgn, float *end, float *dir)
ADD.
int seg_intersect_vregion(geosurf *gs, float *bgn, float *end)
Check if segment intersect vector region.
Point3 * gsdrape_get_allsegments(geosurf *gs, float *bgn, float *end, int *num)
Get all segments.
int get_diag_intersects(geosurf *gs, float *bgn, float *end, float *dir UNUSED)
Get diagonal intersects.
int XY_intersect_plane(float *intersect, float *plane)
Check for intersection (point and plane)
int Point_on_plane(Point3 p1, Point3 p2, Point3 p3, Point3 unk)
Check if point is on plane.
#define GET_MAPATT(buff, offset, att)
#define DRC2OFF(gs, drow, dcol)
#define VROW2DROW(gs, vrow)
#define VCOL2DCOL(gs, vcol)