GRASS 8 Programmer's Manual 8.6.0dev(2026)-5a31d549cc
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 Initialize
88
89 \param[in] rows number of rows
90 \param[in] 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[in] bgn begin point
131 \param[in] end end point
132 \param[out] 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[in] p1 first point
176 \param[in] 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 */
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[in,out] bgn begin point
227 \param[in,out] end end point
228
229 \return 0 if segment doesn't intersect the viewregion, or intersects only at
230 corner
231 \return otherwise returns 1
232 */
233int seg_intersect_vregion(geosurf *gs, float *bgn, float *end)
234{
235 float *replace, xl, yb, xr, yt, xi, yi;
236 int inside = 0;
237
238 xl = 0.0;
239 xr = VCOL2X(gs, VCOLS(gs));
240 yt = VROW2Y(gs, 0);
241 yb = VROW2Y(gs, VROWS(gs));
242
243 if (in_vregion(gs, bgn)) {
244 replace = end;
245 inside++;
246 }
247
248 if (in_vregion(gs, end)) {
249 replace = bgn;
250 inside++;
251 }
252
253 if (inside == 2) {
254 return (1);
255 }
256 else if (inside) {
257 /* one in & one out - replace gets first intersection */
258 if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yb, xl, yt, &xi,
259 &yi)) {
260 /* left */
261 }
262 else if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xr, yb, xr, yt,
263 &xi, &yi)) {
264 /* right */
265 }
266 else if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yb, xr, yb,
267 &xi, &yi)) {
268 /* bottom */
269 }
270 else if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yt, xr, yt,
271 &xi, &yi)) {
272 /* top */
273 }
274
275 replace[X] = xi;
276 replace[Y] = yi;
277 }
278 else {
279 /* both out - find 2 intersects & replace both */
280 float pt1[2], pt2[2];
281
282 replace = pt1;
283 if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yb, xl, yt, &xi,
284 &yi)) {
285 replace[X] = xi;
286 replace[Y] = yi;
287 replace = pt2;
288 inside++;
289 }
290
291 if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xr, yb, xr, yt, &xi,
292 &yi)) {
293 replace[X] = xi;
294 replace[Y] = yi;
295 replace = pt2;
296 inside++;
297 }
298
299 if (inside < 2) {
300 if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yb, xr, yb,
301 &xi, &yi)) {
302 replace[X] = xi;
303 replace[Y] = yi;
304 replace = pt2;
305 inside++;
306 }
307 }
308
309 if (inside < 2) {
310 if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yt, xr, yt,
311 &xi, &yi)) {
312 replace[X] = xi;
313 replace[Y] = yi;
314 inside++;
315 }
316 }
317
318 if (inside < 2) {
319 return (0); /* no intersect or only 1 point on corner */
320 }
321
322 /* compare dist of intersects to bgn - closest replaces bgn */
324 bgn[X] = pt1[X];
325 bgn[Y] = pt1[Y];
326 end[X] = pt2[X];
327 end[Y] = pt2[Y];
328 }
329 else {
330 bgn[X] = pt2[X];
331 bgn[Y] = pt2[Y];
332 end[X] = pt1[X];
333 end[Y] = pt1[Y];
334 }
335 }
336
337 return (1);
338}
339
340/*!
341 \brief ADD
342
343 \param gs surface (geosurf)
344 \param[in,out] bgn begin point (x,y)
345 \param[in,out] end end point (x,y)
346 \param[out] num
347
348 \return pointer to Point3 struct
349 */
350Point3 *gsdrape_get_segments(geosurf *gs, float *bgn, float *end, int *num)
351{
353
354 if (!seg_intersect_vregion(gs, bgn, end)) {
355 *num = 0;
356
357 return (NULL);
358 }
359
361 /* will probably want a force_drape option to get all intersects */
362 I3d[0][X] = bgn[X];
363 I3d[0][Y] = bgn[Y];
364 I3d[0][Z] = gs->att[ATT_TOPO].constant;
365 I3d[1][X] = end[X];
366 I3d[1][Y] = end[Y];
367 I3d[1][Z] = gs->att[ATT_TOPO].constant;
368 *num = 2;
369
370 return (I3d);
371 }
372
373 if (bgn[X] == end[X] && bgn[Y] == end[Y]) {
374 float f[3], l[3];
375
376 interp_first_last(gs, bgn, end, f, l);
377 GS_v3eq(I3d[0], f);
378 GS_v3eq(I3d[1], l);
379
380 /* CHANGE (*num = 1) to reflect degenerate line ? */
381 *num = 2;
382
383 return (I3d);
384 }
385
386 Flat = 0;
387 return (_gsdrape_get_segments(gs, bgn, end, num));
388}
389
390/*!
391 \brief Get all segments
392
393 \param gs surface (geosurf)
394 \param[in,out] bgn begin point
395 \param[in,out] end end point
396 \param[out] num
397
398 \return pointer to Point3 struct
399 */
400Point3 *gsdrape_get_allsegments(geosurf *gs, float *bgn, float *end, int *num)
401{
403
404 if (!seg_intersect_vregion(gs, bgn, end)) {
405 *num = 0;
406 return (NULL);
407 }
408
409 if (bgn[X] == end[X] && bgn[Y] == end[Y]) {
410 float f[3], l[3];
411
412 interp_first_last(gs, bgn, end, f, l);
413 GS_v3eq(I3d[0], f);
414 GS_v3eq(I3d[1], l);
415 *num = 2;
416
417 return (I3d);
418 }
419
421 Flat = 1;
422 }
423 else {
424 Flat = 0;
425 }
426
427 return (_gsdrape_get_segments(gs, bgn, end, num));
428}
429
430/*!
431 \brief ADD
432
433 \param gs surface (geosurf)
434 \param[in] bgn begin point
435 \param[in] end end point
436 \param[out] f first
437 \param[out] l last
438 */
439void interp_first_last(geosurf *gs, float *bgn, float *end, Point3 f, Point3 l)
440{
441 f[X] = bgn[X];
442 f[Y] = bgn[Y];
443
444 l[X] = end[X];
445 l[Y] = end[Y];
446
447 if (Flat) {
448 f[Z] = l[Z] = gs->att[ATT_TOPO].constant;
449 }
450 else {
451 viewcell_tri_interp(gs, Ebuf, f, 0);
452 viewcell_tri_interp(gs, Ebuf, l, 0);
453 }
454
455 return;
456}
457
458/*!
459 \brief ADD
460
461 \param gs surface (geosurf)
462 \param[in,out] pt
463 */
465{
466 typbuff *buf;
467
468 buf = gs_get_att_typbuff(gs, ATT_TOPO, 0);
469
470 return (viewcell_tri_interp(gs, buf, pt, 0));
471}
472
473/*!
474 \brief ADD
475
476 In gsd_surf, tmesh draws polys like so:
477 <pre>
478 --------------
479 | /|
480 | / |
481 | / |
482 | / |
483 | / |
484 | / |
485 | / |
486 | / |
487 | / |
488 | / |
489 | / |
490 |/ |
491 --------------
492 </pre>
493
494 UNLESS the top right or bottom left point is masked, in which case a
495 single triangle with the opposite diagonal is drawn. This case is
496 not yet handled here & should only occur on edges.
497 pt has X & Y coordinates in it, we interpolate Z here
498
499 This could probably be much shorter, but not much faster.
500
501 \param gs
502 \param buf
503 \param[in,out] pt
504 \param[in] check_mask
505
506 \return 1 if point is in view region
507 \return otherwise 0 (if masked)
508 */
510{
511 Point3 p1, p2, p3;
512 int offset, drow, dcol, vrow, vcol;
513 float xmax, ymin, ymax, alpha;
514
515 xmax = VCOL2X(gs, VCOLS(gs));
516 ymax = VROW2Y(gs, 0);
517 ymin = VROW2Y(gs, VROWS(gs));
518
519 if (check_mask) {
520 if (gs_point_is_masked(gs, pt)) {
521 return (0);
522 }
523 }
524
525 if (pt[X] < 0.0 || pt[Y] > ymax) {
526 /* outside on left or top */
527 return (0);
528 }
529
530 if (pt[Y] < ymin || pt[X] > xmax) {
531 /* outside on bottom or right */
532 return (0);
533 }
534
536 pt[Z] = gs->att[ATT_TOPO].constant;
537
538 return (1);
539 }
540 else if (MAP_ATT != gs_get_att_src(gs, ATT_TOPO)) {
541 return (0);
542 }
543
544 vrow = Y2VROW(gs, pt[Y]);
545 vcol = X2VCOL(gs, pt[X]);
546
547 if (vrow < VROWS(gs) && vcol < VCOLS(gs)) {
548 /*not on bottom or right edge */
549 if (pt[X] > 0.0 && pt[Y] < ymax) {
550 /* not on left or top edge */
551 p1[X] = VCOL2X(gs, vcol + 1);
552 p1[Y] = VROW2Y(gs, vrow);
553 drow = VROW2DROW(gs, vrow);
554 dcol = VCOL2DCOL(gs, vcol + 1);
555 offset = DRC2OFF(gs, drow, dcol);
556 GET_MAPATT(buf, offset, p1[Z]); /* top right */
557
558 p2[X] = VCOL2X(gs, vcol);
559 p2[Y] = VROW2Y(gs, vrow + 1);
560 drow = VROW2DROW(gs, vrow + 1);
561 dcol = VCOL2DCOL(gs, vcol);
562 offset = DRC2OFF(gs, drow, dcol);
563 GET_MAPATT(buf, offset, p2[Z]); /* bottom left */
564
565 if ((pt[X] - p2[X]) / VXRES(gs) > (pt[Y] - p2[Y]) / VYRES(gs)) {
566 /* lower triangle */
567 p3[X] = VCOL2X(gs, vcol + 1);
568 p3[Y] = VROW2Y(gs, vrow + 1);
569 drow = VROW2DROW(gs, vrow + 1);
570 dcol = VCOL2DCOL(gs, vcol + 1);
571 offset = DRC2OFF(gs, drow, dcol);
572 GET_MAPATT(buf, offset, p3[Z]); /* bottom right */
573 }
574 else {
575 /* upper triangle */
576 p3[X] = VCOL2X(gs, vcol);
577 p3[Y] = VROW2Y(gs, vrow);
578 drow = VROW2DROW(gs, vrow);
579 dcol = VCOL2DCOL(gs, vcol);
580 offset = DRC2OFF(gs, drow, dcol);
581 GET_MAPATT(buf, offset, p3[Z]); /* top left */
582 }
583
584 return (Point_on_plane(p1, p2, p3, pt));
585 }
586 else if (pt[X] == 0.0) {
587 /* on left edge */
588 if (pt[Y] < ymax) {
589 vrow = Y2VROW(gs, pt[Y]);
590 drow = VROW2DROW(gs, vrow);
591 offset = DRC2OFF(gs, drow, 0);
592 GET_MAPATT(buf, offset, p1[Z]);
593
594 drow = VROW2DROW(gs, vrow + 1);
595 offset = DRC2OFF(gs, drow, 0);
596 GET_MAPATT(buf, offset, p2[Z]);
597
598 alpha = (VROW2Y(gs, vrow) - pt[Y]) / VYRES(gs);
599 pt[Z] = LERP(alpha, p1[Z], p2[Z]);
600 }
601 else {
602 /* top left corner */
603 GET_MAPATT(buf, 0, pt[Z]);
604 }
605
606 return (1);
607 }
608 else if (pt[Y] == gs->yrange) {
609 /* on top edge, not a corner */
610 vcol = X2VCOL(gs, pt[X]);
611 dcol = VCOL2DCOL(gs, vcol);
612 GET_MAPATT(buf, dcol, p1[Z]);
613
614 dcol = VCOL2DCOL(gs, vcol + 1);
615 GET_MAPATT(buf, dcol, p2[Z]);
616
617 alpha = (pt[X] - VCOL2X(gs, vcol)) / VXRES(gs);
618 pt[Z] = LERP(alpha, p1[Z], p2[Z]);
619
620 return (1);
621 }
622 }
623 else if (vrow == VROWS(gs)) {
624 /* on bottom edge */
625 drow = VROW2DROW(gs, VROWS(gs));
626
627 if (pt[X] > 0.0 && pt[X] < xmax) {
628 /* not a corner */
629 vcol = X2VCOL(gs, pt[X]);
630 dcol = VCOL2DCOL(gs, vcol);
631 offset = DRC2OFF(gs, drow, dcol);
632 GET_MAPATT(buf, offset, p1[Z]);
633
634 dcol = VCOL2DCOL(gs, vcol + 1);
635 offset = DRC2OFF(gs, drow, dcol);
636 GET_MAPATT(buf, offset, p2[Z]);
637
638 alpha = (pt[X] - VCOL2X(gs, vcol)) / VXRES(gs);
639 pt[Z] = LERP(alpha, p1[Z], p2[Z]);
640
641 return (1);
642 }
643 else if (pt[X] == 0.0) {
644 /* bottom left corner */
645 offset = DRC2OFF(gs, drow, 0);
646 GET_MAPATT(buf, offset, pt[Z]);
647
648 return (1);
649 }
650 else {
651 /* bottom right corner */
652 dcol = VCOL2DCOL(gs, VCOLS(gs));
653 offset = DRC2OFF(gs, drow, dcol);
654 GET_MAPATT(buf, offset, pt[Z]);
655
656 return (1);
657 }
658 }
659 else {
660 /* on right edge, not bottom corner */
661 dcol = VCOL2DCOL(gs, VCOLS(gs));
662
663 if (pt[Y] < ymax) {
664 vrow = Y2VROW(gs, pt[Y]);
665 drow = VROW2DROW(gs, vrow);
666 offset = DRC2OFF(gs, drow, dcol);
667 GET_MAPATT(buf, offset, p1[Z]);
668
669 drow = VROW2DROW(gs, vrow + 1);
670 offset = DRC2OFF(gs, drow, dcol);
671 GET_MAPATT(buf, offset, p2[Z]);
672
673 alpha = (VROW2Y(gs, vrow) - pt[Y]) / VYRES(gs);
674 pt[Z] = LERP(alpha, p1[Z], p2[Z]);
675
676 return (1);
677 }
678 else {
679 /* top right corner */
680 GET_MAPATT(buf, dcol, pt[Z]);
681
682 return (1);
683 }
684 }
685
686 return (0);
687}
688
689/*!
690 \brief ADD
691
692 \param gs surface (geosurf)
693 \param[in] pt
694
695 \return 1
696 \return 0
697 */
698int in_vregion(geosurf *gs, float *pt)
699{
700 if (pt[X] >= 0.0 && pt[Y] <= gs->yrange) {
701 if (pt[X] <= VCOL2X(gs, VCOLS(gs))) {
702 return (pt[Y] >= VROW2Y(gs, VROWS(gs)));
703 }
704 }
705
706 return (0);
707}
708
709/*!
710 \brief ADD
711
712 After all the intersections between the segment and triangle
713 edges have been found, they are in three lists. (intersections
714 with vertical, horizontal, and diagonal triangle edges)
715
716 Each list is ordered in space from first to last segment points,
717 but now the lists need to be woven together. This routine
718 starts with the first point of the segment and then checks the
719 next point in each list to find the closest, eliminating duplicates
720 along the way and storing the result in I3d.
721
722 \param gs surface (geosurf)
723 \param[in] first first point
724 \param[in] last last point
725 \param[in] vi
726 \param[in] hi
727 \param[in] di
728
729 \return
730 */
731int order_intersects(geosurf *gs, Point3 first, Point3 last, int vi, int hi,
732 int di)
733{
734 int num, i, found, cv, ch, cd, cnum;
735 float dv, dh, dd, big, cpoint[2];
736
737 cv = ch = cd = cnum = 0;
738 num = vi + hi + di;
739
740 cpoint[X] = first[X];
741 cpoint[Y] = first[Y];
742
743 if (in_vregion(gs, first)) {
744 I3d[cnum][X] = first[X];
745 I3d[cnum][Y] = first[Y];
746 I3d[cnum][Z] = first[Z];
747 cnum++;
748 }
749
750 /* TODO: big could still be less than first dist */
751 big = gs->yrange * gs->yrange + gs->xrange * gs->xrange; /*BIG distance */
752 dv = dh = dd = big;
753
754 for (i = 0; i < num; i = cv + ch + cd) {
755 if (cv < vi) {
756 dv = dist_squared_2d(Vi[cv], cpoint);
757
758 if (dv < EPSILON) {
759 cv++;
760 continue;
761 }
762 }
763 else {
764 dv = big;
765 }
766
767 if (ch < hi) {
768 dh = dist_squared_2d(Hi[ch], cpoint);
769
770 if (dh < EPSILON) {
771 ch++;
772 continue;
773 }
774 }
775 else {
776 dh = big;
777 }
778
779 if (cd < di) {
780 dd = dist_squared_2d(Di[cd], cpoint);
781
782 if (dd < EPSILON) {
783 cd++;
784 continue;
785 }
786 }
787 else {
788 dd = big;
789 }
790
791 found = 0;
792
793 if (cd < di) {
794 if (dd <= dv && dd <= dh) {
795 found = 1;
796 cpoint[X] = I3d[cnum][X] = Di[cd][X];
797 cpoint[Y] = I3d[cnum][Y] = Di[cd][Y];
798 I3d[cnum][Z] = Di[cd][Z];
799 cnum++;
800
801 if (EQUAL(dd, dv)) {
802 cv++;
803 }
804
805 if (EQUAL(dd, dh)) {
806 ch++;
807 }
808
809 cd++;
810 }
811 }
812
813 if (!found) {
814 if (cv < vi) {
815 if (dv <= dh) {
816 found = 1;
817 cpoint[X] = I3d[cnum][X] = Vi[cv][X];
818 cpoint[Y] = I3d[cnum][Y] = Vi[cv][Y];
819 I3d[cnum][Z] = Vi[cv][Z];
820 cnum++;
821
822 if (EQUAL(dv, dh)) {
823 ch++;
824 }
825
826 cv++;
827 }
828 }
829 }
830
831 if (!found) {
832 if (ch < hi) {
833 cpoint[X] = I3d[cnum][X] = Hi[ch][X];
834 cpoint[Y] = I3d[cnum][Y] = Hi[ch][Y];
835 I3d[cnum][Z] = Hi[ch][Z];
836 cnum++;
837 ch++;
838 }
839 }
840
841 if (i == cv + ch + cd) {
842 G_debug(5, "order_intersects(): stuck on %d", cnum);
843 G_debug(5, "order_intersects(): cv = %d, ch = %d, cd = %d", cv, ch,
844 cd);
845 G_debug(5, "order_intersects(): dv = %f, dh = %f, dd = %f", dv, dh,
846 dd);
847
848 break;
849 }
850 }
851
852 if (EQUAL(last[X], cpoint[X]) && EQUAL(last[Y], cpoint[Y])) {
853 return (cnum);
854 }
855
856 if (in_vregion(gs, last)) {
857 /* TODO: check for last point on corner ? */
858 I3d[cnum][X] = last[X];
859 I3d[cnum][Y] = last[Y];
860 I3d[cnum][Z] = last[Z];
861 ++cnum;
862 }
863
864 return (cnum);
865}
866
867/*!
868 \brief ADD
869
870 \todo For consistency, need to decide how last row & last column are
871 displayed - would it look funny to always draw last row/col with
872 finer resolution if necessary, or would it be better to only show
873 full rows/cols?
874
875 Colinear already eliminated
876
877 \param gs surface (geosurf)
878 \param[in] bgn begin point
879 \param[in] end end point
880 \param[in] dir direction
881
882 \return
883 */
884int get_vert_intersects(geosurf *gs, float *bgn, float *end, float *dir)
885{
886 int fcol, lcol, incr, hits, num, offset, drow1, drow2;
887 float xl, yb, xr, yt, z1, z2, alpha;
888 float yres, xi, yi;
889 int bgncol, endcol, cols, rows;
890
891 yres = VYRES(gs);
892 cols = VCOLS(gs);
893 rows = VROWS(gs);
894
895 bgncol = X2VCOL(gs, bgn[X]);
896 endcol = X2VCOL(gs, end[X]);
897
898 if (bgncol > cols && endcol > cols) {
899 return 0;
900 }
901
902 if (bgncol == endcol) {
903 return 0;
904 }
905
906 fcol = dir[X] > 0 ? bgncol + 1 : bgncol;
907 lcol = dir[X] > 0 ? endcol : endcol + 1;
908
909 /* assuming only showing FULL cols */
910 incr = lcol - fcol > 0 ? 1 : -1;
911
912 while (fcol > cols || fcol < 0) {
913 fcol += incr;
914 }
915
916 while (lcol > cols || lcol < 0) {
917 lcol -= incr;
918 }
919
920 num = abs(lcol - fcol) + 1;
921
922 yb = gs->yrange - (yres * rows) - EPSILON;
923 yt = gs->yrange + EPSILON;
924
925 for (hits = 0; hits < num; hits++) {
926 xl = xr = VCOL2X(gs, fcol);
927
928 if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yt, xr, yb, &xi,
929 &yi)) {
930 Vi[hits][X] = xi;
931 Vi[hits][Y] = yi;
932
933 /* find data rows */
934 if (Flat) {
935 Vi[hits][Z] = gs->att[ATT_TOPO].constant;
936 }
937 else {
938 drow1 = Y2VROW(gs, Vi[hits][Y]) * gs->y_mod;
939 drow2 = (1 + Y2VROW(gs, Vi[hits][Y])) * gs->y_mod;
940
941 if (drow2 >= gs->rows) {
942 drow2 = gs->rows - 1; /*bottom edge */
943 }
944
945 alpha = ((gs->yrange - drow1 * gs->yres) - Vi[hits][Y]) / yres;
946
947 offset = DRC2OFF(gs, drow1, fcol * gs->x_mod);
948 GET_MAPATT(Ebuf, offset, z1);
949 offset = DRC2OFF(gs, drow2, fcol * gs->x_mod);
950 GET_MAPATT(Ebuf, offset, z2);
951 Vi[hits][Z] = LERP(alpha, z1, z2);
952 }
953 }
954
955 /* if they don't intersect, something's wrong! */
956 /* should only happen on endpoint, so it will be added later */
957 else {
958 hits--;
959 num--;
960 }
961
962 fcol += incr;
963 }
964
965 return (hits);
966}
967
968/*!
969 \brief Get horizontal intersects
970
971 \param gs surface (geosurf)
972 \param[in] bgn begin point
973 \param[in] end end point
974 \param[in] dir
975
976 \return number of intersects
977 */
978int get_horz_intersects(geosurf *gs, float *bgn, float *end, float *dir)
979{
980 int frow, lrow, incr, hits, num, offset, dcol1, dcol2;
981 float xl, yb, xr, yt, z1, z2, alpha;
982 float xres, xi, yi;
983 int bgnrow, endrow, rows, cols;
984
985 xres = VXRES(gs);
986 cols = VCOLS(gs);
987 rows = VROWS(gs);
988
989 bgnrow = Y2VROW(gs, bgn[Y]);
990 endrow = Y2VROW(gs, end[Y]);
991 if (bgnrow == endrow) {
992 return 0;
993 }
994
995 if (bgnrow > rows && endrow > rows) {
996 return 0;
997 }
998
999 frow = dir[Y] > 0 ? bgnrow : bgnrow + 1;
1000 lrow = dir[Y] > 0 ? endrow + 1 : endrow;
1001
1002 /* assuming only showing FULL rows */
1003 incr = lrow - frow > 0 ? 1 : -1;
1004
1005 while (frow > rows || frow < 0) {
1006 frow += incr;
1007 }
1008
1009 while (lrow > rows || lrow < 0) {
1010 lrow -= incr;
1011 }
1012
1013 num = abs(lrow - frow) + 1;
1014
1015 xl = 0.0 - EPSILON;
1016 xr = xres * cols + EPSILON;
1017
1018 for (hits = 0; hits < num; hits++) {
1019 yb = yt = VROW2Y(gs, frow);
1020
1021 if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yt, xr, yb, &xi,
1022 &yi)) {
1023 Hi[hits][X] = xi;
1024 Hi[hits][Y] = yi;
1025
1026 /* find data cols */
1027 if (Flat) {
1028 Hi[hits][Z] = gs->att[ATT_TOPO].constant;
1029 }
1030 else {
1031 dcol1 = X2VCOL(gs, Hi[hits][X]) * gs->x_mod;
1032 dcol2 = (1 + X2VCOL(gs, Hi[hits][X])) * gs->x_mod;
1033
1034 if (dcol2 >= gs->cols) {
1035 dcol2 = gs->cols - 1; /* right edge */
1036 }
1037
1038 alpha = (Hi[hits][X] - (dcol1 * gs->xres)) / xres;
1039
1040 offset = DRC2OFF(gs, frow * gs->y_mod, dcol1);
1041 GET_MAPATT(Ebuf, offset, z1);
1042 offset = DRC2OFF(gs, frow * gs->y_mod, dcol2);
1043 GET_MAPATT(Ebuf, offset, z2);
1044 Hi[hits][Z] = LERP(alpha, z1, z2);
1045 }
1046 }
1047
1048 /* if they don't intersect, something's wrong! */
1049 /* should only happen on endpoint, so it will be added later */
1050 else {
1051 hits--;
1052 num--;
1053 }
1054
1055 frow += incr;
1056 }
1057
1058 return (hits);
1059}
1060
1061/*!
1062 \brief Get diagonal intersects
1063
1064 Colinear already eliminated
1065
1066 \param gs surface (geosurf)
1067 \param[in] bgn begin point
1068 \param[in] end end point
1069 \param dir ? (unused)
1070
1071 \return number of intersects
1072 */
1073int get_diag_intersects(geosurf *gs, float *bgn, float *end,
1074 float *dir G_UNUSED)
1075{
1076 int fdig, ldig, incr, hits, num, offset;
1077 int vrow, vcol, drow1, drow2, dcol1, dcol2;
1078 float xl, yb, xr, yt, z1, z2, alpha;
1079 float xres, yres, xi, yi, dx, dy;
1080 int diags, cols, rows, lower;
1081 Point3 pt;
1082
1083 xres = VXRES(gs);
1084 yres = VYRES(gs);
1085 cols = VCOLS(gs);
1086 rows = VROWS(gs);
1087 diags = rows + cols; /* -1 ? */
1088
1089 /* determine upper/lower triangle for last */
1090 vrow = Y2VROW(gs, end[Y]);
1091 vcol = X2VCOL(gs, end[X]);
1092 pt[X] = VCOL2X(gs, vcol);
1093 pt[Y] = VROW2Y(gs, vrow + 1);
1094 lower = ((end[X] - pt[X]) / xres > (end[Y] - pt[Y]) / yres);
1095 ldig = lower ? vrow + vcol + 1 : vrow + vcol;
1096
1097 /* determine upper/lower triangle for first */
1098 vrow = Y2VROW(gs, bgn[Y]);
1099 vcol = X2VCOL(gs, bgn[X]);
1100 pt[X] = VCOL2X(gs, vcol);
1101 pt[Y] = VROW2Y(gs, vrow + 1);
1102 lower = ((bgn[X] - pt[X]) / xres > (bgn[Y] - pt[Y]) / yres);
1103 fdig = lower ? vrow + vcol + 1 : vrow + vcol;
1104
1105 /* adjust according to direction */
1106 if (ldig > fdig) {
1107 fdig++;
1108 }
1109
1110 if (fdig > ldig) {
1111 ldig++;
1112 }
1113
1114 incr = ldig - fdig > 0 ? 1 : -1;
1115
1116 while (fdig > diags || fdig < 0) {
1117 fdig += incr;
1118 }
1119
1120 while (ldig > diags || ldig < 0) {
1121 ldig -= incr;
1122 }
1123
1124 num = abs(ldig - fdig) + 1;
1125
1126 for (hits = 0; hits < num; hits++) {
1127 yb = gs->yrange - (yres * (fdig < rows ? fdig : rows)) - EPSILON;
1128 xl = VCOL2X(gs, (fdig < rows ? 0 : fdig - rows)) - EPSILON;
1129 yt = gs->yrange - (yres * (fdig < cols ? 0 : fdig - cols)) + EPSILON;
1130 xr = VCOL2X(gs, (fdig < cols ? fdig : cols)) + EPSILON;
1131
1132 if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yb, xr, yt, &xi,
1133 &yi)) {
1134 Di[hits][X] = xi;
1135 Di[hits][Y] = yi;
1136
1137 if (ISNODE(xi, xres)) {
1138 /* then it's also a ynode */
1139 num--;
1140 hits--;
1141 continue;
1142 }
1143
1144 /* find data rows */
1145 drow1 = Y2VROW(gs, Di[hits][Y]) * gs->y_mod;
1146 drow2 = (1 + Y2VROW(gs, Di[hits][Y])) * gs->y_mod;
1147
1148 if (drow2 >= gs->rows) {
1149 drow2 = gs->rows - 1; /* bottom edge */
1150 }
1151
1152 /* find data cols */
1153 if (Flat) {
1154 Di[hits][Z] = gs->att[ATT_TOPO].constant;
1155 }
1156 else {
1157 dcol1 = X2VCOL(gs, Di[hits][X]) * gs->x_mod;
1158 dcol2 = (1 + X2VCOL(gs, Di[hits][X])) * gs->x_mod;
1159
1160 if (dcol2 >= gs->cols) {
1161 dcol2 = gs->cols - 1; /* right edge */
1162 }
1163
1164 dx = DCOL2X(gs, dcol2) - Di[hits][X];
1165 dy = DROW2Y(gs, drow1) - Di[hits][Y];
1166 alpha =
1167 sqrt(dx * dx + dy * dy) / sqrt(xres * xres + yres * yres);
1168
1169 offset = DRC2OFF(gs, drow1, dcol2);
1170 GET_MAPATT(Ebuf, offset, z1);
1171 offset = DRC2OFF(gs, drow2, dcol1);
1172 GET_MAPATT(Ebuf, offset, z2);
1173 Di[hits][Z] = LERP(alpha, z1, z2);
1174 }
1175 }
1176
1177 /* if they don't intersect, something's wrong! */
1178 /* should only happen on endpoint, so it will be added later */
1179 else {
1180 hits--;
1181 num--;
1182 }
1183
1184 fdig += incr;
1185 }
1186
1187 return (hits);
1188}
1189
1190/*!
1191 \brief Line intersect
1192
1193 Author: Mukesh Prasad
1194 Modified for floating point: Bill Brown
1195
1196 This function computes whether two line segments,
1197 respectively joining the input points (x1,y1) -- (x2,y2)
1198 and the input points (x3,y3) -- (x4,y4) intersect.
1199 If the lines intersect, the output variables x, y are
1200 set to coordinates of the point of intersection.
1201
1202 \param[in] x1,y1,x2,y2 coordinates of endpoints of one segment
1203 \param[in] x3,y3,x4,y4 coordinates of endpoints of other segment
1204 \param[out] x,y coordinates of intersection point
1205
1206 \return 0 no intersection
1207 \return 1 intersect
1208 \return 2 collinear
1209 */
1210int segs_intersect(float x1, float y1, float x2, float y2, float x3, float y3,
1211 float x4, float y4, float *x, float *y)
1212{
1213 float a1, a2, b1, b2, c1, c2; /* Coefficients of line eqns. */
1214 float r1, r2, r3, r4; /* 'Sign' values */
1215 float denom, /* offset, */ num; /* Intermediate values */
1216
1217 /* Compute a1, b1, c1, where line joining points 1 and 2
1218 * is "a1 x + b1 y + c1 = 0".
1219 */
1220 a1 = y2 - y1;
1221 b1 = x1 - x2;
1222 c1 = x2 * y1 - x1 * y2;
1223
1224 /* Compute r3 and r4.
1225 */
1226 r3 = a1 * x3 + b1 * y3 + c1;
1227 r4 = a1 * x4 + b1 * y4 + c1;
1228
1229 /* Check signs of r3 and r4. If both point 3 and point 4 lie on
1230 * same side of line 1, the line segments do not intersect.
1231 */
1232
1233 if (!EQUAL(r3, 0.0) && !EQUAL(r4, 0.0) && SAME_SIGNS(r3, r4)) {
1234 return (DONT_INTERSECT);
1235 }
1236
1237 /* Compute a2, b2, c2 */
1238 a2 = y4 - y3;
1239 b2 = x3 - x4;
1240 c2 = x4 * y3 - x3 * y4;
1241
1242 /* Compute r1 and r2 */
1243 r1 = a2 * x1 + b2 * y1 + c2;
1244 r2 = a2 * x2 + b2 * y2 + c2;
1245
1246 /* Check signs of r1 and r2. If both point 1 and point 2 lie
1247 * on same side of second line segment, the line segments do
1248 * not intersect.
1249 */
1250
1251 if (!EQUAL(r1, 0.0) && !EQUAL(r2, 0.0) && SAME_SIGNS(r1, r2)) {
1252 return (DONT_INTERSECT);
1253 }
1254
1255 /* Line segments intersect: compute intersection point.
1256 */
1257 denom = a1 * b2 - a2 * b1;
1258
1259 if (denom == 0) {
1260 return (COLLINEAR);
1261 }
1262
1263 /* offset = denom < 0 ? -denom / 2 : denom / 2; */
1264
1265 /* The denom/2 is to get rounding instead of truncating. It
1266 * is added or subtracted to the numerator, depending upon the
1267 * sign of the numerator.
1268 */
1269 num = b1 * c2 - b2 * c1;
1270
1271 *x = num / denom;
1272
1273 num = a2 * c1 - a1 * c2;
1274 *y = num / denom;
1275
1276 return (DO_INTERSECT);
1277}
1278
1279/*!
1280 \brief Check if point is on plane
1281
1282 Plane defined by three points here; user fills in unk[X] & unk[Y]
1283
1284 \param[in,out] p1,p2,p3 points defining plane
1285 \param[in,out] unk point
1286
1287 \return 1 point on plane
1288 \return 0 point not on plane
1289 */
1291{
1292 float plane[4];
1293
1294 P3toPlane(p1, p2, p3, plane);
1295
1296 return (XY_intersect_plane(unk, plane));
1297}
1298
1299/*!
1300 \brief Check for intersection (point and plane)
1301
1302 Ax + By + Cz + D = 0, so z = (Ax + By + D) / -C
1303
1304 User fills in intersect[X] & intersect[Y]
1305
1306 \param[in,out] intersect intersect coordinates
1307 \param[in] plane plane definition
1308
1309 \return 0 doesn't intersect
1310 \return 1 intesects
1311 */
1312int XY_intersect_plane(float *intersect, float *plane)
1313{
1314 float x, y;
1315
1316 if (!plane[Z]) {
1317 return (0); /* doesn't intersect */
1318 }
1319
1320 x = intersect[X];
1321 y = intersect[Y];
1322 intersect[Z] = (plane[X] * x + plane[Y] * y + plane[W]) / -plane[Z];
1323
1324 return (1);
1325}
1326
1327/*!
1328 \brief Define plane
1329
1330 \param[in] p1,p2,p3 three point on plane
1331 \param[out] plane plane definition
1332
1333 \return 1
1334 */
1335int P3toPlane(Point3 p1, Point3 p2, Point3 p3, float *plane)
1336{
1337 Point3 v1, v2, norm;
1338
1339 v1[X] = p1[X] - p3[X];
1340 v1[Y] = p1[Y] - p3[Y];
1341 v1[Z] = p1[Z] - p3[Z];
1342
1343 v2[X] = p2[X] - p3[X];
1344 v2[Y] = p2[Y] - p3[Y];
1345 v2[Z] = p2[Z] - p3[Z];
1346
1347 V3Cross(v1, v2, norm);
1348
1349 plane[X] = norm[X];
1350 plane[Y] = norm[Y];
1351 plane[Z] = norm[Z];
1352 plane[W] = -p3[X] * norm[X] - p3[Y] * norm[Y] - p3[Z] * norm[Z];
1353
1354 return (1);
1355}
1356
1357/*!
1358 \brief Get cross product
1359
1360 \param[in] a,b
1361 \param[out] c
1362
1363 \return cross product c = a cross b
1364 */
1366{
1367 c[X] = (a[Y] * b[Z]) - (a[Z] * b[Y]);
1368 c[Y] = (a[Z] * b[X]) - (a[X] * b[Z]);
1369 c[Z] = (a[X] * b[Y]) - (a[Y] * b[X]);
1370
1371 return (1);
1372}
#define EPSILON
#define NULL
Definition ccmath.h:32
void G_free(void *)
Free allocated memory.
Definition gis/alloc.c:147
void G_warning(const char *,...) __attribute__((format(printf
int G_debug(int, const char *,...) __attribute__((format(printf
int gs_point_is_masked(geosurf *, float *)
Check if point is masked.
Definition gs.c:1314
float GS_P2distance(float *, float *)
Calculate distance in plane.
Definition gs_util.c:160
int gs_get_att_src(geosurf *, int)
Get attribute source.
Definition gs.c:656
void GS_v3eq(float *, float *)
Copy vector values.
Definition gs_util.c:178
typbuff * gs_get_att_typbuff(geosurf *, int, int)
Get attribute data buffer.
Definition gs.c:681
void GS_v2dir(float *, float *, float *)
Get a normalized direction from v1 to v2, store in v3 (2D)
Definition gs_util.c:382
#define G_UNUSED
A macro for an attribute, if attached to a variable, indicating that the variable is not used.
Definition gis.h:46
#define _(str)
Definition glocale.h:10
#define COLLINEAR
Definition gsdrape.c:62
int P3toPlane(Point3 p1, Point3 p2, Point3 p3, float *plane)
Define plane.
Definition gsdrape.c:1335
#define ISNODE(p, res)
Definition gsdrape.c:66
Point3 * gsdrape_get_segments(geosurf *gs, float *bgn, float *end, int *num)
ADD.
Definition gsdrape.c:350
int in_vregion(geosurf *gs, float *pt)
ADD.
Definition gsdrape.c:698
#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:731
int _viewcell_tri_interp(geosurf *gs, Point3 pt)
ADD.
Definition gsdrape.c:464
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:439
#define DONT_INTERSECT
Definition gsdrape.c:60
int V3Cross(Point3 a, Point3 b, Point3 c)
Get cross product.
Definition gsdrape.c:1365
int viewcell_tri_interp(geosurf *gs, typbuff *buf, Point3 pt, int check_mask)
ADD.
Definition gsdrape.c:509
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
int get_horz_intersects(geosurf *gs, float *bgn, float *end, float *dir)
Get horizontal intersects.
Definition gsdrape.c:978
#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:884
int get_diag_intersects(geosurf *gs, float *bgn, float *end, float *dir)
Get diagonal intersects.
Definition gsdrape.c:1073
int seg_intersect_vregion(geosurf *gs, float *bgn, float *end)
Check if segment intersect vector region.
Definition gsdrape.c:233
#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:400
int XY_intersect_plane(float *intersect, float *plane)
Check for intersection (point and plane)
Definition gsdrape.c:1312
int Point_on_plane(Point3 p1, Point3 p2, Point3 p3, Point3 unk)
Check if point is on plane.
Definition gsdrape.c:1290
#define GET_MAPATT(buff, offset, att)
Definition gsget.h:29
OGSF header file (structures)
#define X
Definition ogsf.h:140
#define ATT_TOPO
Definition ogsf.h:75
float Point3[3]
Definition ogsf.h:205
#define Z
Definition ogsf.h:142
#define W
Definition ogsf.h:143
#define Y
Definition ogsf.h:141
#define MAP_ATT
Definition ogsf.h:85
#define CONST_ATT
Definition ogsf.h:86
double b
Definition r_raster.c:39
double l
Definition r_raster.c:39
#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
Definition ogsf.h:266
Definition clip.h:6
#define x