GRASS GIS 8 Programmer's Manual 8.4.1(2025)-45ca3179ab
Loading...
Searching...
No Matches
gsdrape.c
Go to the documentation of this file.
1/*!
2 \file lib/ogsf/gsdrape.c
3
4 \brief OGSF library - functions to intersect line segments with edges of
5 surface polygons
6
7 GRASS OpenGL gsurf OGSF Library
8
9 For efficiency, intersections are found without respect to which
10 specific triangle edge is intersected, but on a broader sense with
11 the horizontal, vertical, and diagonal seams in the grid, then
12 the intersections are ordered. If quadstrips are used for drawing
13 rather than tmesh, triangulation is not consistent; for diagonal
14 intersections, the proper diagonal to intersect would need to be
15 determined according to the algorithm used by qstrip (look at nearby
16 normals). It may be faster to go ahead and find the intersections
17 with the other diagonals using the same methods, then at sorting
18 time determine which diagonal array to look at for each quad.
19 It would also require a mechanism for throwing out unused intersections
20 with the diagonals during the ordering phase.
21 Do intersections in 2D, fill line structure with 3D pts (maybe calling
22 routine will cache for redrawing). Get Z value by using linear interp
23 between corners.
24
25 - check for easy cases:
26 - single point
27 - colinear with horizontal or vertical edges
28 - colinear with diagonal edges of triangles
29 - calculate three arrays of ordered intersections:
30 - with vertical edges
31 - with horizontal edges
32 - with diagonal edges and interpolate Z, using simple linear interpolation.
33 - eliminate duplicate intersections (need only compare one coord for each)
34 - build ordered set of points.
35
36 Return static pointer to 3D set of points. Max number of intersections
37 will be rows + cols + diags, so should allocate this number to initialize.
38 Let calling routine worry about copying points for caching.
39
40 (C) 1999-2008 by the GRASS Development Team
41
42 This program is free software under the
43 GNU General Public License (>=v2).
44 Read the file COPYING that comes with GRASS
45 for details.
46
47 \author Bill Brown UI GMS Lab
48 \author Doxygenized by Martin Landa <landa.martin gmail.com> (May 2008)
49 */
50
51#include <stdlib.h>
52
53#include <grass/ogsf.h>
54#include <grass/glocale.h>
55
56#include "gsget.h"
57#include "rowcol.h"
58#include "math.h"
59
60#define DONT_INTERSECT 0
61#define DO_INTERSECT 1
62#define COLLINEAR 2
63
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)
67
68#define SAME_SIGNS(a, b) ((a >= 0 && b >= 0) || (a < 0 && b < 0))
69
70static int drape_line_init(int, int);
71static Point3 *_gsdrape_get_segments(geosurf *, float *, float *, int *);
72static float dist_squared_2d(float *, float *);
73
74/* array of points to be returned */
75static Point3 *I3d;
76
77/* make dependent on resolution? */
78static float EPSILON = 0.000001;
79
80/*vertical, horizontal, & diagonal intersections */
81static Point3 *Vi, *Hi, *Di;
82
83static typbuff *Ebuf; /* elevation buffer */
84static int Flat;
85
86/*!
87 \brief Initizalize
88
89 \param rows number of rows
90 \param cols number of columns
91
92 \return -1 on failure
93 \return 1 on success
94 */
95static int drape_line_init(int rows, int cols)
96{
97 /* use G_calloc() [-> G_fatal_error] instead of calloc ? */
98 if (NULL == (I3d = (Point3 *)calloc(2 * (rows + cols), sizeof(Point3)))) {
99 return (-1);
100 }
101
102 if (NULL == (Vi = (Point3 *)calloc(cols, sizeof(Point3)))) {
103 G_free(I3d);
104
105 return (-1);
106 }
107
108 if (NULL == (Hi = (Point3 *)calloc(rows, sizeof(Point3)))) {
109 G_free(I3d);
110 G_free(Vi);
111
112 return (-1);
113 }
114
115 if (NULL == (Di = (Point3 *)calloc(rows + cols, sizeof(Point3)))) {
116 G_free(I3d);
117 G_free(Vi);
118 G_free(Hi);
119
120 return (-1);
121 }
122
123 return (1);
124}
125
126/*!
127 \brief Get segments
128
129 \param gs surface (geosurf)
130 \param bgn begin point
131 \param end end point
132 \param num
133
134 \return pointer to Point3 struct
135 */
136static Point3 *_gsdrape_get_segments(geosurf *gs, float *bgn, float *end,
137 int *num)
138{
139 float f[3], l[3];
140 int vi, hi, di;
141 float dir[2], yres, xres;
142
143 xres = VXRES(gs);
144 yres = VYRES(gs);
145
146 vi = hi = di = 0;
147 GS_v2dir(bgn, end, dir);
148
149 if (dir[X]) {
150 vi = get_vert_intersects(gs, bgn, end, dir);
151 }
152
153 if (dir[Y]) {
154 hi = get_horz_intersects(gs, bgn, end, dir);
155 }
156
157 if (!((end[Y] - bgn[Y]) / (end[X] - bgn[X]) == yres / xres)) {
158 di = get_diag_intersects(gs, bgn, end, dir);
159 }
160
161 interp_first_last(gs, bgn, end, f, l);
162
163 *num = order_intersects(gs, f, l, vi, hi, di);
164 /* fills in return values, eliminates dupes (corners) */
165
166 G_debug(5, "_gsdrape_get_segments vi=%d, hi=%d, di=%d, num=%d", vi, hi, di,
167 *num);
168
169 return (I3d);
170}
171
172/*!
173 \brief Calculate 2D distance
174
175 \param p1 first point
176 \param p2 second point
177
178 \return distance
179 */
180static float dist_squared_2d(float *p1, float *p2)
181{
182 float dx, dy;
183
184 dx = p2[X] - p1[X];
185 dy = p2[Y] - p1[Y];
186
187 return (dx * dx + dy * dy);
188}
189
190/*!
191 \brief ADD
192
193 \param gs surface (geosurf)
194
195 \return -1 on failure
196 \return 1 on success
197 */
198int gsdrape_set_surface(geosurf *gs)
199{
200 static int first = 1;
201
202 if (first) {
203 first = 0;
204
205 if (0 > drape_line_init(gs->rows, gs->cols)) {
206 G_warning(_("Unable to process vector map - out of memory"));
207 Ebuf = NULL;
208
209 return (-1);
210 }
211 }
212
213 Ebuf = gs_get_att_typbuff(gs, ATT_TOPO, 0);
214
215 return (1);
216}
217
218/*!
219 \brief Check if segment intersect vector region
220
221 Clipping performed:
222 - bgn and end are replaced so that both points are within viewregion
223 - if seg intersects
224
225 \param gs surface (geosurf)
226 \param bgn begin point
227 \param end end point
228
229 \return 0 if segment doesn't intersect the viewregion, or intersects only at
230 corner \return otherwise returns 1
231 */
232int seg_intersect_vregion(geosurf *gs, float *bgn, float *end)
233{
234 float *replace, xl, yb, xr, yt, xi, yi;
235 int inside = 0;
236
237 xl = 0.0;
238 xr = VCOL2X(gs, VCOLS(gs));
239 yt = VROW2Y(gs, 0);
240 yb = VROW2Y(gs, VROWS(gs));
241
242 if (in_vregion(gs, bgn)) {
243 replace = end;
244 inside++;
245 }
246
247 if (in_vregion(gs, end)) {
248 replace = bgn;
249 inside++;
250 }
251
252 if (inside == 2) {
253 return (1);
254 }
255 else if (inside) {
256 /* one in & one out - replace gets first intersection */
257 if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yb, xl, yt, &xi,
258 &yi)) {
259 /* left */
260 }
261 else if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xr, yb, xr, yt,
262 &xi, &yi)) {
263 /* right */
264 }
265 else if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yb, xr, yb,
266 &xi, &yi)) {
267 /* bottom */
268 }
269 else if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yt, xr, yt,
270 &xi, &yi)) {
271 /* top */
272 }
273
274 replace[X] = xi;
275 replace[Y] = yi;
276 }
277 else {
278 /* both out - find 2 intersects & replace both */
279 float pt1[2], pt2[2];
280
281 replace = pt1;
282 if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yb, xl, yt, &xi,
283 &yi)) {
284 replace[X] = xi;
285 replace[Y] = yi;
286 replace = pt2;
287 inside++;
288 }
289
290 if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xr, yb, xr, yt, &xi,
291 &yi)) {
292 replace[X] = xi;
293 replace[Y] = yi;
294 replace = pt2;
295 inside++;
296 }
297
298 if (inside < 2) {
299 if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yb, xr, yb,
300 &xi, &yi)) {
301 replace[X] = xi;
302 replace[Y] = yi;
303 replace = pt2;
304 inside++;
305 }
306 }
307
308 if (inside < 2) {
309 if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yt, xr, yt,
310 &xi, &yi)) {
311 replace[X] = xi;
312 replace[Y] = yi;
313 inside++;
314 }
315 }
316
317 if (inside < 2) {
318 return (0); /* no intersect or only 1 point on corner */
319 }
320
321 /* compare dist of intersects to bgn - closest replaces bgn */
322 if (GS_P2distance(bgn, pt1) < GS_P2distance(bgn, pt2)) {
323 bgn[X] = pt1[X];
324 bgn[Y] = pt1[Y];
325 end[X] = pt2[X];
326 end[Y] = pt2[Y];
327 }
328 else {
329 bgn[X] = pt2[X];
330 bgn[Y] = pt2[Y];
331 end[X] = pt1[X];
332 end[Y] = pt1[Y];
333 }
334 }
335
336 return (1);
337}
338
339/*!
340 \brief ADD
341
342 \param gs surface (geosurf)
343 \param bgn begin point (x,y)
344 \param end end point (x,y)
345 \param num
346
347 \return pointer to Point3 struct
348 */
349Point3 *gsdrape_get_segments(geosurf *gs, float *bgn, float *end, int *num)
350{
352
353 if (!seg_intersect_vregion(gs, bgn, end)) {
354 *num = 0;
355
356 return (NULL);
357 }
358
359 if (CONST_ATT == gs_get_att_src(gs, ATT_TOPO)) {
360 /* will probably want a force_drape option to get all intersects */
361 I3d[0][X] = bgn[X];
362 I3d[0][Y] = bgn[Y];
363 I3d[0][Z] = gs->att[ATT_TOPO].constant;
364 I3d[1][X] = end[X];
365 I3d[1][Y] = end[Y];
366 I3d[1][Z] = gs->att[ATT_TOPO].constant;
367 *num = 2;
368
369 return (I3d);
370 }
371
372 if (bgn[X] == end[X] && bgn[Y] == end[Y]) {
373 float f[3], l[3];
374
375 interp_first_last(gs, bgn, end, f, l);
376 GS_v3eq(I3d[0], f);
377 GS_v3eq(I3d[1], l);
378
379 /* CHANGE (*num = 1) to reflect degenerate line ? */
380 *num = 2;
381
382 return (I3d);
383 }
384
385 Flat = 0;
386 return (_gsdrape_get_segments(gs, bgn, end, num));
387}
388
389/*!
390 \brief Get all segments
391
392 \param gs surface (geosurf)
393 \param bgn begin point
394 \param end end point
395 \param num
396
397 \return pointer to Point3 struct
398 */
399Point3 *gsdrape_get_allsegments(geosurf *gs, float *bgn, float *end, int *num)
400{
402
403 if (!seg_intersect_vregion(gs, bgn, end)) {
404 *num = 0;
405 return (NULL);
406 }
407
408 if (bgn[X] == end[X] && bgn[Y] == end[Y]) {
409 float f[3], l[3];
410
411 interp_first_last(gs, bgn, end, f, l);
412 GS_v3eq(I3d[0], f);
413 GS_v3eq(I3d[1], l);
414 *num = 2;
415
416 return (I3d);
417 }
418
419 if (CONST_ATT == gs_get_att_src(gs, ATT_TOPO)) {
420 Flat = 1;
421 }
422 else {
423 Flat = 0;
424 }
425
426 return (_gsdrape_get_segments(gs, bgn, end, num));
427}
428
429/*!
430 \brief ADD
431
432 \param gs surface (geosurf)
433 \param bgn begin point
434 \param end end point
435 \param f first
436 \param l last
437 */
438void interp_first_last(geosurf *gs, float *bgn, float *end, Point3 f, Point3 l)
439{
440 f[X] = bgn[X];
441 f[Y] = bgn[Y];
442
443 l[X] = end[X];
444 l[Y] = end[Y];
445
446 if (Flat) {
447 f[Z] = l[Z] = gs->att[ATT_TOPO].constant;
448 }
449 else {
450 viewcell_tri_interp(gs, Ebuf, f, 0);
451 viewcell_tri_interp(gs, Ebuf, l, 0);
452 }
453
454 return;
455}
456
457/*!
458 \brief ADD
459
460 \param gs surface (geosurf)
461 \param pt
462 */
463int _viewcell_tri_interp(geosurf *gs, Point3 pt)
464{
465 typbuff *buf;
466
467 buf = gs_get_att_typbuff(gs, ATT_TOPO, 0);
468
469 return (viewcell_tri_interp(gs, buf, pt, 0));
470}
471
472/*!
473 \brief ADD
474
475 In gsd_surf, tmesh draws polys like so:
476 <pre>
477 --------------
478 | /|
479 | / |
480 | / |
481 | / |
482 | / |
483 | / |
484 | / |
485 | / |
486 | / |
487 | / |
488 | / |
489 |/ |
490 --------------
491 </pre>
492
493 UNLESS the top right or bottom left point is masked, in which case a
494 single triangle with the opposite diagonal is drawn. This case is
495 not yet handled here & should only occur on edges.
496 pt has X & Y coordinates in it, we interpolate Z here
497
498 This could probably be much shorter, but not much faster.
499
500 \return 1 if point is in view region
501 \return otherwise 0 (if masked)
502 */
503int viewcell_tri_interp(geosurf *gs, typbuff *buf, Point3 pt, int check_mask)
504{
505 Point3 p1, p2, p3;
506 int offset, drow, dcol, vrow, vcol;
507 float xmax, ymin, ymax, alpha;
508
509 xmax = VCOL2X(gs, VCOLS(gs));
510 ymax = VROW2Y(gs, 0);
511 ymin = VROW2Y(gs, VROWS(gs));
512
513 if (check_mask) {
514 if (gs_point_is_masked(gs, pt)) {
515 return (0);
516 }
517 }
518
519 if (pt[X] < 0.0 || pt[Y] > ymax) {
520 /* outside on left or top */
521 return (0);
522 }
523
524 if (pt[Y] < ymin || pt[X] > xmax) {
525 /* outside on bottom or right */
526 return (0);
527 }
528
529 if (CONST_ATT == gs_get_att_src(gs, ATT_TOPO)) {
530 pt[Z] = gs->att[ATT_TOPO].constant;
531
532 return (1);
533 }
534 else if (MAP_ATT != gs_get_att_src(gs, ATT_TOPO)) {
535 return (0);
536 }
537
538 vrow = Y2VROW(gs, pt[Y]);
539 vcol = X2VCOL(gs, pt[X]);
540
541 if (vrow < VROWS(gs) && vcol < VCOLS(gs)) {
542 /*not on bottom or right edge */
543 if (pt[X] > 0.0 && pt[Y] < ymax) {
544 /* not on left or top edge */
545 p1[X] = VCOL2X(gs, vcol + 1);
546 p1[Y] = VROW2Y(gs, vrow);
547 drow = VROW2DROW(gs, vrow);
548 dcol = VCOL2DCOL(gs, vcol + 1);
549 offset = DRC2OFF(gs, drow, dcol);
550 GET_MAPATT(buf, offset, p1[Z]); /* top right */
551
552 p2[X] = VCOL2X(gs, vcol);
553 p2[Y] = VROW2Y(gs, vrow + 1);
554 drow = VROW2DROW(gs, vrow + 1);
555 dcol = VCOL2DCOL(gs, vcol);
556 offset = DRC2OFF(gs, drow, dcol);
557 GET_MAPATT(buf, offset, p2[Z]); /* bottom left */
558
559 if ((pt[X] - p2[X]) / VXRES(gs) > (pt[Y] - p2[Y]) / VYRES(gs)) {
560 /* lower triangle */
561 p3[X] = VCOL2X(gs, vcol + 1);
562 p3[Y] = VROW2Y(gs, vrow + 1);
563 drow = VROW2DROW(gs, vrow + 1);
564 dcol = VCOL2DCOL(gs, vcol + 1);
565 offset = DRC2OFF(gs, drow, dcol);
566 GET_MAPATT(buf, offset, p3[Z]); /* bottom right */
567 }
568 else {
569 /* upper triangle */
570 p3[X] = VCOL2X(gs, vcol);
571 p3[Y] = VROW2Y(gs, vrow);
572 drow = VROW2DROW(gs, vrow);
573 dcol = VCOL2DCOL(gs, vcol);
574 offset = DRC2OFF(gs, drow, dcol);
575 GET_MAPATT(buf, offset, p3[Z]); /* top left */
576 }
577
578 return (Point_on_plane(p1, p2, p3, pt));
579 }
580 else if (pt[X] == 0.0) {
581 /* on left edge */
582 if (pt[Y] < ymax) {
583 vrow = Y2VROW(gs, pt[Y]);
584 drow = VROW2DROW(gs, vrow);
585 offset = DRC2OFF(gs, drow, 0);
586 GET_MAPATT(buf, offset, p1[Z]);
587
588 drow = VROW2DROW(gs, vrow + 1);
589 offset = DRC2OFF(gs, drow, 0);
590 GET_MAPATT(buf, offset, p2[Z]);
591
592 alpha = (VROW2Y(gs, vrow) - pt[Y]) / VYRES(gs);
593 pt[Z] = LERP(alpha, p1[Z], p2[Z]);
594 }
595 else {
596 /* top left corner */
597 GET_MAPATT(buf, 0, pt[Z]);
598 }
599
600 return (1);
601 }
602 else if (pt[Y] == gs->yrange) {
603 /* on top edge, not a corner */
604 vcol = X2VCOL(gs, pt[X]);
605 dcol = VCOL2DCOL(gs, vcol);
606 GET_MAPATT(buf, dcol, p1[Z]);
607
608 dcol = VCOL2DCOL(gs, vcol + 1);
609 GET_MAPATT(buf, dcol, p2[Z]);
610
611 alpha = (pt[X] - VCOL2X(gs, vcol)) / VXRES(gs);
612 pt[Z] = LERP(alpha, p1[Z], p2[Z]);
613
614 return (1);
615 }
616 }
617 else if (vrow == VROWS(gs)) {
618 /* on bottom edge */
619 drow = VROW2DROW(gs, VROWS(gs));
620
621 if (pt[X] > 0.0 && pt[X] < xmax) {
622 /* not a corner */
623 vcol = X2VCOL(gs, pt[X]);
624 dcol = VCOL2DCOL(gs, vcol);
625 offset = DRC2OFF(gs, drow, dcol);
626 GET_MAPATT(buf, offset, p1[Z]);
627
628 dcol = VCOL2DCOL(gs, vcol + 1);
629 offset = DRC2OFF(gs, drow, dcol);
630 GET_MAPATT(buf, offset, p2[Z]);
631
632 alpha = (pt[X] - VCOL2X(gs, vcol)) / VXRES(gs);
633 pt[Z] = LERP(alpha, p1[Z], p2[Z]);
634
635 return (1);
636 }
637 else if (pt[X] == 0.0) {
638 /* bottom left corner */
639 offset = DRC2OFF(gs, drow, 0);
640 GET_MAPATT(buf, offset, pt[Z]);
641
642 return (1);
643 }
644 else {
645 /* bottom right corner */
646 dcol = VCOL2DCOL(gs, VCOLS(gs));
647 offset = DRC2OFF(gs, drow, dcol);
648 GET_MAPATT(buf, offset, pt[Z]);
649
650 return (1);
651 }
652 }
653 else {
654 /* on right edge, not bottom corner */
655 dcol = VCOL2DCOL(gs, VCOLS(gs));
656
657 if (pt[Y] < ymax) {
658 vrow = Y2VROW(gs, pt[Y]);
659 drow = VROW2DROW(gs, vrow);
660 offset = DRC2OFF(gs, drow, dcol);
661 GET_MAPATT(buf, offset, p1[Z]);
662
663 drow = VROW2DROW(gs, vrow + 1);
664 offset = DRC2OFF(gs, drow, dcol);
665 GET_MAPATT(buf, offset, p2[Z]);
666
667 alpha = (VROW2Y(gs, vrow) - pt[Y]) / VYRES(gs);
668 pt[Z] = LERP(alpha, p1[Z], p2[Z]);
669
670 return (1);
671 }
672 else {
673 /* top right corner */
674 GET_MAPATT(buf, dcol, pt[Z]);
675
676 return (1);
677 }
678 }
679
680 return (0);
681}
682
683/*!
684 \brief ADD
685
686 \param gs surface (geosurf)
687
688 \return 1
689 \return 0
690 */
691int in_vregion(geosurf *gs, float *pt)
692{
693 if (pt[X] >= 0.0 && pt[Y] <= gs->yrange) {
694 if (pt[X] <= VCOL2X(gs, VCOLS(gs))) {
695 return (pt[Y] >= VROW2Y(gs, VROWS(gs)));
696 }
697 }
698
699 return (0);
700}
701
702/*!
703 \brief ADD
704
705 After all the intersections between the segment and triangle
706 edges have been found, they are in three lists. (intersections
707 with vertical, horizontal, and diagonal triangle edges)
708
709 Each list is ordered in space from first to last segment points,
710 but now the lists need to be woven together. This routine
711 starts with the first point of the segment and then checks the
712 next point in each list to find the closest, eliminating duplicates
713 along the way and storing the result in I3d.
714
715 \param gs surface (geosurf)
716 \param first first point
717 \param last last point
718 \param vi
719 \param hi
720 \param di
721
722 \return
723 */
724int order_intersects(geosurf *gs, Point3 first, Point3 last, int vi, int hi,
725 int di)
726{
727 int num, i, found, cv, ch, cd, cnum;
728 float dv, dh, dd, big, cpoint[2];
729
730 cv = ch = cd = cnum = 0;
731 num = vi + hi + di;
732
733 cpoint[X] = first[X];
734 cpoint[Y] = first[Y];
735
736 if (in_vregion(gs, first)) {
737 I3d[cnum][X] = first[X];
738 I3d[cnum][Y] = first[Y];
739 I3d[cnum][Z] = first[Z];
740 cnum++;
741 }
742
743 /* TODO: big could still be less than first dist */
744 big = gs->yrange * gs->yrange + gs->xrange * gs->xrange; /*BIG distance */
745 dv = dh = dd = big;
746
747 for (i = 0; i < num; i = cv + ch + cd) {
748 if (cv < vi) {
749 dv = dist_squared_2d(Vi[cv], cpoint);
750
751 if (dv < EPSILON) {
752 cv++;
753 continue;
754 }
755 }
756 else {
757 dv = big;
758 }
759
760 if (ch < hi) {
761 dh = dist_squared_2d(Hi[ch], cpoint);
762
763 if (dh < EPSILON) {
764 ch++;
765 continue;
766 }
767 }
768 else {
769 dh = big;
770 }
771
772 if (cd < di) {
773 dd = dist_squared_2d(Di[cd], cpoint);
774
775 if (dd < EPSILON) {
776 cd++;
777 continue;
778 }
779 }
780 else {
781 dd = big;
782 }
783
784 found = 0;
785
786 if (cd < di) {
787 if (dd <= dv && dd <= dh) {
788 found = 1;
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];
792 cnum++;
793
794 if (EQUAL(dd, dv)) {
795 cv++;
796 }
797
798 if (EQUAL(dd, dh)) {
799 ch++;
800 }
801
802 cd++;
803 }
804 }
805
806 if (!found) {
807 if (cv < vi) {
808 if (dv <= dh) {
809 found = 1;
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];
813 cnum++;
814
815 if (EQUAL(dv, dh)) {
816 ch++;
817 }
818
819 cv++;
820 }
821 }
822 }
823
824 if (!found) {
825 if (ch < hi) {
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];
829 cnum++;
830 ch++;
831 }
832 }
833
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,
837 cd);
838 G_debug(5, "order_intersects(): dv = %f, dh = %f, dd = %f", dv, dh,
839 dd);
840
841 break;
842 }
843 }
844
845 if (EQUAL(last[X], cpoint[X]) && EQUAL(last[Y], cpoint[Y])) {
846 return (cnum);
847 }
848
849 if (in_vregion(gs, last)) {
850 /* TODO: check for last point on corner ? */
851 I3d[cnum][X] = last[X];
852 I3d[cnum][Y] = last[Y];
853 I3d[cnum][Z] = last[Z];
854 ++cnum;
855 }
856
857 return (cnum);
858}
859
860/*!
861 \brief ADD
862
863 \todo For consistency, need to decide how last row & last column are
864 displayed - would it look funny to always draw last row/col with
865 finer resolution if necessary, or would it be better to only show
866 full rows/cols?
867
868 Colinear already eliminated
869
870 \param gs surface (geosurf)
871 \param bgn begin point
872 \param end end point
873 \param dir direction
874
875 \return
876 */
877int get_vert_intersects(geosurf *gs, float *bgn, float *end, float *dir)
878{
879 int fcol, lcol, incr, hits, num, offset, drow1, drow2;
880 float xl, yb, xr, yt, z1, z2, alpha;
881 float yres, xi, yi;
882 int bgncol, endcol, cols, rows;
883
884 yres = VYRES(gs);
885 cols = VCOLS(gs);
886 rows = VROWS(gs);
887
888 bgncol = X2VCOL(gs, bgn[X]);
889 endcol = X2VCOL(gs, end[X]);
890
891 if (bgncol > cols && endcol > cols) {
892 return 0;
893 }
894
895 if (bgncol == endcol) {
896 return 0;
897 }
898
899 fcol = dir[X] > 0 ? bgncol + 1 : bgncol;
900 lcol = dir[X] > 0 ? endcol : endcol + 1;
901
902 /* assuming only showing FULL cols */
903 incr = lcol - fcol > 0 ? 1 : -1;
904
905 while (fcol > cols || fcol < 0) {
906 fcol += incr;
907 }
908
909 while (lcol > cols || lcol < 0) {
910 lcol -= incr;
911 }
912
913 num = abs(lcol - fcol) + 1;
914
915 yb = gs->yrange - (yres * rows) - EPSILON;
916 yt = gs->yrange + EPSILON;
917
918 for (hits = 0; hits < num; hits++) {
919 xl = xr = VCOL2X(gs, fcol);
920
921 if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yt, xr, yb, &xi,
922 &yi)) {
923 Vi[hits][X] = xi;
924 Vi[hits][Y] = yi;
925
926 /* find data rows */
927 if (Flat) {
928 Vi[hits][Z] = gs->att[ATT_TOPO].constant;
929 }
930 else {
931 drow1 = Y2VROW(gs, Vi[hits][Y]) * gs->y_mod;
932 drow2 = (1 + Y2VROW(gs, Vi[hits][Y])) * gs->y_mod;
933
934 if (drow2 >= gs->rows) {
935 drow2 = gs->rows - 1; /*bottom edge */
936 }
937
938 alpha = ((gs->yrange - drow1 * gs->yres) - Vi[hits][Y]) / yres;
939
940 offset = DRC2OFF(gs, drow1, fcol * gs->x_mod);
941 GET_MAPATT(Ebuf, offset, z1);
942 offset = DRC2OFF(gs, drow2, fcol * gs->x_mod);
943 GET_MAPATT(Ebuf, offset, z2);
944 Vi[hits][Z] = LERP(alpha, z1, z2);
945 }
946 }
947
948 /* if they don't intersect, something's wrong! */
949 /* should only happen on endpoint, so it will be added later */
950 else {
951 hits--;
952 num--;
953 }
954
955 fcol += incr;
956 }
957
958 return (hits);
959}
960
961/*!
962 \brief Get horizontal intersects
963
964 \param gs surface (geosurf)
965 \param bgn begin point
966 \param end end point
967 \param dir
968
969 \return number of intersects
970 */
971int get_horz_intersects(geosurf *gs, float *bgn, float *end, float *dir)
972{
973 int frow, lrow, incr, hits, num, offset, dcol1, dcol2;
974 float xl, yb, xr, yt, z1, z2, alpha;
975 float xres, xi, yi;
976 int bgnrow, endrow, rows, cols;
977
978 xres = VXRES(gs);
979 cols = VCOLS(gs);
980 rows = VROWS(gs);
981
982 bgnrow = Y2VROW(gs, bgn[Y]);
983 endrow = Y2VROW(gs, end[Y]);
984 if (bgnrow == endrow) {
985 return 0;
986 }
987
988 if (bgnrow > rows && endrow > rows) {
989 return 0;
990 }
991
992 frow = dir[Y] > 0 ? bgnrow : bgnrow + 1;
993 lrow = dir[Y] > 0 ? endrow + 1 : endrow;
994
995 /* assuming only showing FULL rows */
996 incr = lrow - frow > 0 ? 1 : -1;
997
998 while (frow > rows || frow < 0) {
999 frow += incr;
1000 }
1001
1002 while (lrow > rows || lrow < 0) {
1003 lrow -= incr;
1004 }
1005
1006 num = abs(lrow - frow) + 1;
1007
1008 xl = 0.0 - EPSILON;
1009 xr = xres * cols + EPSILON;
1010
1011 for (hits = 0; hits < num; hits++) {
1012 yb = yt = VROW2Y(gs, frow);
1013
1014 if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yt, xr, yb, &xi,
1015 &yi)) {
1016 Hi[hits][X] = xi;
1017 Hi[hits][Y] = yi;
1018
1019 /* find data cols */
1020 if (Flat) {
1021 Hi[hits][Z] = gs->att[ATT_TOPO].constant;
1022 }
1023 else {
1024 dcol1 = X2VCOL(gs, Hi[hits][X]) * gs->x_mod;
1025 dcol2 = (1 + X2VCOL(gs, Hi[hits][X])) * gs->x_mod;
1026
1027 if (dcol2 >= gs->cols) {
1028 dcol2 = gs->cols - 1; /* right edge */
1029 }
1030
1031 alpha = (Hi[hits][X] - (dcol1 * gs->xres)) / xres;
1032
1033 offset = DRC2OFF(gs, frow * gs->y_mod, dcol1);
1034 GET_MAPATT(Ebuf, offset, z1);
1035 offset = DRC2OFF(gs, frow * gs->y_mod, dcol2);
1036 GET_MAPATT(Ebuf, offset, z2);
1037 Hi[hits][Z] = LERP(alpha, z1, z2);
1038 }
1039 }
1040
1041 /* if they don't intersect, something's wrong! */
1042 /* should only happen on endpoint, so it will be added later */
1043 else {
1044 hits--;
1045 num--;
1046 }
1047
1048 frow += incr;
1049 }
1050
1051 return (hits);
1052}
1053
1054/*!
1055 \brief Get diagonal intersects
1056
1057 Colinear already eliminated
1058
1059 \param gs surface (geosurf)
1060 \param bgn begin point
1061 \param end end point
1062 \param dir ? (unused)
1063
1064 \return number of intersects
1065 */
1066int get_diag_intersects(geosurf *gs, float *bgn, float *end, float *dir UNUSED)
1067{
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;
1073 Point3 pt;
1074
1075 xres = VXRES(gs);
1076 yres = VYRES(gs);
1077 cols = VCOLS(gs);
1078 rows = VROWS(gs);
1079 diags = rows + cols; /* -1 ? */
1080
1081 /* determine upper/lower triangle for last */
1082 vrow = Y2VROW(gs, end[Y]);
1083 vcol = X2VCOL(gs, end[X]);
1084 pt[X] = VCOL2X(gs, vcol);
1085 pt[Y] = VROW2Y(gs, vrow + 1);
1086 lower = ((end[X] - pt[X]) / xres > (end[Y] - pt[Y]) / yres);
1087 ldig = lower ? vrow + vcol + 1 : vrow + vcol;
1088
1089 /* determine upper/lower triangle for first */
1090 vrow = Y2VROW(gs, bgn[Y]);
1091 vcol = X2VCOL(gs, bgn[X]);
1092 pt[X] = VCOL2X(gs, vcol);
1093 pt[Y] = VROW2Y(gs, vrow + 1);
1094 lower = ((bgn[X] - pt[X]) / xres > (bgn[Y] - pt[Y]) / yres);
1095 fdig = lower ? vrow + vcol + 1 : vrow + vcol;
1096
1097 /* adjust according to direction */
1098 if (ldig > fdig) {
1099 fdig++;
1100 }
1101
1102 if (fdig > ldig) {
1103 ldig++;
1104 }
1105
1106 incr = ldig - fdig > 0 ? 1 : -1;
1107
1108 while (fdig > diags || fdig < 0) {
1109 fdig += incr;
1110 }
1111
1112 while (ldig > diags || ldig < 0) {
1113 ldig -= incr;
1114 }
1115
1116 num = abs(ldig - fdig) + 1;
1117
1118 for (hits = 0; hits < num; hits++) {
1119 yb = gs->yrange - (yres * (fdig < rows ? fdig : rows)) - EPSILON;
1120 xl = VCOL2X(gs, (fdig < rows ? 0 : fdig - rows)) - EPSILON;
1121 yt = gs->yrange - (yres * (fdig < cols ? 0 : fdig - cols)) + EPSILON;
1122 xr = VCOL2X(gs, (fdig < cols ? fdig : cols)) + EPSILON;
1123
1124 if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yb, xr, yt, &xi,
1125 &yi)) {
1126 Di[hits][X] = xi;
1127 Di[hits][Y] = yi;
1128
1129 if (ISNODE(xi, xres)) {
1130 /* then it's also a ynode */
1131 num--;
1132 hits--;
1133 continue;
1134 }
1135
1136 /* find data rows */
1137 drow1 = Y2VROW(gs, Di[hits][Y]) * gs->y_mod;
1138 drow2 = (1 + Y2VROW(gs, Di[hits][Y])) * gs->y_mod;
1139
1140 if (drow2 >= gs->rows) {
1141 drow2 = gs->rows - 1; /* bottom edge */
1142 }
1143
1144 /* find data cols */
1145 if (Flat) {
1146 Di[hits][Z] = gs->att[ATT_TOPO].constant;
1147 }
1148 else {
1149 dcol1 = X2VCOL(gs, Di[hits][X]) * gs->x_mod;
1150 dcol2 = (1 + X2VCOL(gs, Di[hits][X])) * gs->x_mod;
1151
1152 if (dcol2 >= gs->cols) {
1153 dcol2 = gs->cols - 1; /* right edge */
1154 }
1155
1156 dx = DCOL2X(gs, dcol2) - Di[hits][X];
1157 dy = DROW2Y(gs, drow1) - Di[hits][Y];
1158 alpha =
1159 sqrt(dx * dx + dy * dy) / sqrt(xres * xres + yres * yres);
1160
1161 offset = DRC2OFF(gs, drow1, dcol2);
1162 GET_MAPATT(Ebuf, offset, z1);
1163 offset = DRC2OFF(gs, drow2, dcol1);
1164 GET_MAPATT(Ebuf, offset, z2);
1165 Di[hits][Z] = LERP(alpha, z1, z2);
1166 }
1167 }
1168
1169 /* if they don't intersect, something's wrong! */
1170 /* should only happen on endpoint, so it will be added later */
1171 else {
1172 hits--;
1173 num--;
1174 }
1175
1176 fdig += incr;
1177 }
1178
1179 return (hits);
1180}
1181
1182/*!
1183 \brief Line intersect
1184
1185 Author: Mukesh Prasad
1186 Modified for floating point: Bill Brown
1187
1188 This function computes whether two line segments,
1189 respectively joining the input points (x1,y1) -- (x2,y2)
1190 and the input points (x3,y3) -- (x4,y4) intersect.
1191 If the lines intersect, the output variables x, y are
1192 set to coordinates of the point of intersection.
1193
1194 \param x1,y1,x2,y2 coordinates of endpoints of one segment
1195 \param x3,y3,x4,y4 coordinates of endpoints of other segment
1196 \param[out] x,y coordinates of intersection point
1197
1198 \return 0 no intersection
1199 \return 1 intersect
1200 \return 2 collinear
1201 */
1202int segs_intersect(float x1, float y1, float x2, float y2, float x3, float y3,
1203 float x4, float y4, float *x, float *y)
1204{
1205 float a1, a2, b1, b2, c1, c2; /* Coefficients of line eqns. */
1206 float r1, r2, r3, r4; /* 'Sign' values */
1207 float denom, /* offset, */ num; /* Intermediate values */
1208
1209 /* Compute a1, b1, c1, where line joining points 1 and 2
1210 * is "a1 x + b1 y + c1 = 0".
1211 */
1212 a1 = y2 - y1;
1213 b1 = x1 - x2;
1214 c1 = x2 * y1 - x1 * y2;
1215
1216 /* Compute r3 and r4.
1217 */
1218 r3 = a1 * x3 + b1 * y3 + c1;
1219 r4 = a1 * x4 + b1 * y4 + c1;
1220
1221 /* Check signs of r3 and r4. If both point 3 and point 4 lie on
1222 * same side of line 1, the line segments do not intersect.
1223 */
1224
1225 if (!EQUAL(r3, 0.0) && !EQUAL(r4, 0.0) && SAME_SIGNS(r3, r4)) {
1226 return (DONT_INTERSECT);
1227 }
1228
1229 /* Compute a2, b2, c2 */
1230 a2 = y4 - y3;
1231 b2 = x3 - x4;
1232 c2 = x4 * y3 - x3 * y4;
1233
1234 /* Compute r1 and r2 */
1235 r1 = a2 * x1 + b2 * y1 + c2;
1236 r2 = a2 * x2 + b2 * y2 + c2;
1237
1238 /* Check signs of r1 and r2. If both point 1 and point 2 lie
1239 * on same side of second line segment, the line segments do
1240 * not intersect.
1241 */
1242
1243 if (!EQUAL(r1, 0.0) && !EQUAL(r2, 0.0) && SAME_SIGNS(r1, r2)) {
1244 return (DONT_INTERSECT);
1245 }
1246
1247 /* Line segments intersect: compute intersection point.
1248 */
1249 denom = a1 * b2 - a2 * b1;
1250
1251 if (denom == 0) {
1252 return (COLLINEAR);
1253 }
1254
1255 /* offset = denom < 0 ? -denom / 2 : denom / 2; */
1256
1257 /* The denom/2 is to get rounding instead of truncating. It
1258 * is added or subtracted to the numerator, depending upon the
1259 * sign of the numerator.
1260 */
1261 num = b1 * c2 - b2 * c1;
1262
1263 *x = num / denom;
1264
1265 num = a2 * c1 - a1 * c2;
1266 *y = num / denom;
1267
1268 return (DO_INTERSECT);
1269}
1270
1271/*!
1272 \brief Check if point is on plane
1273
1274 Plane defined by three points here; user fills in unk[X] & unk[Y]
1275
1276 \param p1,p2,p3 points defining plane
1277 \param unk point
1278
1279 \return 1 point on plane
1280 \return 0 point not on plane
1281 */
1282int Point_on_plane(Point3 p1, Point3 p2, Point3 p3, Point3 unk)
1283{
1284 float plane[4];
1285
1286 P3toPlane(p1, p2, p3, plane);
1287
1288 return (XY_intersect_plane(unk, plane));
1289}
1290
1291/*!
1292 \brief Check for intersection (point and plane)
1293
1294 Ax + By + Cz + D = 0, so z = (Ax + By + D) / -C
1295
1296 User fills in intersect[X] & intersect[Y]
1297
1298 \param[out] intersect intersect coordinates
1299 \param plane plane definition
1300
1301 \return 0 doesn't intersect
1302 \return 1 intesects
1303 */
1304int XY_intersect_plane(float *intersect, float *plane)
1305{
1306 float x, y;
1307
1308 if (!plane[Z]) {
1309 return (0); /* doesn't intersect */
1310 }
1311
1312 x = intersect[X];
1313 y = intersect[Y];
1314 intersect[Z] = (plane[X] * x + plane[Y] * y + plane[W]) / -plane[Z];
1315
1316 return (1);
1317}
1318
1319/*!
1320 \brief Define plane
1321
1322 \param p1,p2,p3 three point on plane
1323 \param[out] plane plane definition
1324
1325 \return 1
1326 */
1327int P3toPlane(Point3 p1, Point3 p2, Point3 p3, float *plane)
1328{
1329 Point3 v1, v2, norm;
1330
1331 v1[X] = p1[X] - p3[X];
1332 v1[Y] = p1[Y] - p3[Y];
1333 v1[Z] = p1[Z] - p3[Z];
1334
1335 v2[X] = p2[X] - p3[X];
1336 v2[Y] = p2[Y] - p3[Y];
1337 v2[Z] = p2[Z] - p3[Z];
1338
1339 V3Cross(v1, v2, norm);
1340
1341 plane[X] = norm[X];
1342 plane[Y] = norm[Y];
1343 plane[Z] = norm[Z];
1344 plane[W] = -p3[X] * norm[X] - p3[Y] * norm[Y] - p3[Z] * norm[Z];
1345
1346 return (1);
1347}
1348
1349/*!
1350 \brief Get cross product
1351
1352 \param a,b,c
1353
1354 \return cross product c = a cross b
1355 */
1356int V3Cross(Point3 a, Point3 b, Point3 c)
1357{
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]);
1361
1362 return (1);
1363}
void G_free(void *buf)
Free allocated memory.
Definition alloc.c:150
#define EPSILON
#define NULL
Definition ccmath.h:32
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition debug.c:66
double b
double l
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition gis/error.c:203
typbuff * gs_get_att_typbuff(geosurf *gs, int desc, int to_write)
Get attribute data buffer.
Definition gs.c:681
int gs_point_is_masked(geosurf *gs, float *pt)
Check if point is masked.
Definition gs.c:1314
int gs_get_att_src(geosurf *gs, int desc)
Get attribute source.
Definition gs.c:656
void GS_v2dir(float *v1, float *v2, float *v3)
Get a normalized direction from v1 to v2, store in v3 (2D)
Definition gs_util.c:382
void GS_v3eq(float *v1, float *v2)
Copy vector values.
Definition gs_util.c:178
float GS_P2distance(float *from, float *to)
Calculate distance in plane.
Definition gs_util.c:160
#define COLLINEAR
Definition gsdrape.c:62
int P3toPlane(Point3 p1, Point3 p2, Point3 p3, float *plane)
Define plane.
Definition gsdrape.c:1327
#define ISNODE(p, res)
Definition gsdrape.c:66
Point3 * gsdrape_get_segments(geosurf *gs, float *bgn, float *end, int *num)
ADD.
Definition gsdrape.c:349
int in_vregion(geosurf *gs, float *pt)
ADD.
Definition gsdrape.c:691
#define SAME_SIGNS(a, b)
Definition gsdrape.c:68
int order_intersects(geosurf *gs, Point3 first, Point3 last, int vi, int hi, int di)
ADD.
Definition gsdrape.c:724
int _viewcell_tri_interp(geosurf *gs, Point3 pt)
ADD.
Definition gsdrape.c:463
int gsdrape_set_surface(geosurf *gs)
ADD.
Definition gsdrape.c:198
#define EQUAL(a, b)
Definition gsdrape.c:65
void interp_first_last(geosurf *gs, float *bgn, float *end, Point3 f, Point3 l)
ADD.
Definition gsdrape.c:438
#define DONT_INTERSECT
Definition gsdrape.c:60
int V3Cross(Point3 a, Point3 b, Point3 c)
Get cross product.
Definition gsdrape.c:1356
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
int get_horz_intersects(geosurf *gs, float *bgn, float *end, float *dir)
Get horizontal intersects.
Definition gsdrape.c:971
#define LERP(a, l, h)
Definition gsdrape.c:64
int get_vert_intersects(geosurf *gs, float *bgn, float *end, float *dir)
ADD.
Definition gsdrape.c:877
int seg_intersect_vregion(geosurf *gs, float *bgn, float *end)
Check if segment intersect vector region.
Definition gsdrape.c:232
#define DO_INTERSECT
Definition gsdrape.c:61
Point3 * gsdrape_get_allsegments(geosurf *gs, float *bgn, float *end, int *num)
Get all segments.
Definition gsdrape.c:399
int get_diag_intersects(geosurf *gs, float *bgn, float *end, float *dir UNUSED)
Get diagonal intersects.
Definition gsdrape.c:1066
int XY_intersect_plane(float *intersect, float *plane)
Check for intersection (point and plane)
Definition gsdrape.c:1304
int Point_on_plane(Point3 p1, Point3 p2, Point3 p3, Point3 unk)
Check if point is on plane.
Definition gsdrape.c:1282
#define GET_MAPATT(buff, offset, att)
Definition gsget.h:29
#define VYRES(gs)
Definition rowcol.h:10
#define Y2VROW(gs, py)
Definition rowcol.h:27
#define VXRES(gs)
Definition rowcol.h:9
#define VCOL2X(gs, vcol)
Definition rowcol.h:40
#define VCOLS(gs)
Definition rowcol.h:14
#define VROWS(gs)
Definition rowcol.h:13
#define DRC2OFF(gs, drow, dcol)
Definition rowcol.h:17
#define DCOL2X(gs, dcol)
Definition rowcol.h:36
#define VROW2Y(gs, vrow)
Definition rowcol.h:39
#define VROW2DROW(gs, vrow)
Definition rowcol.h:31
#define X2VCOL(gs, px)
Definition rowcol.h:28
#define VCOL2DCOL(gs, vcol)
Definition rowcol.h:32
#define DROW2Y(gs, drow)
Definition rowcol.h:35
#define X(j)
#define x
#define Y(j)