1 /*!
2 \file lib/ogsf/gsdrape.c
3
4 \brief OGSF library - functions to intersect line segments with edges of surface polygons
5
6 GRASS OpenGL gsurf OGSF Library
7
8 For efficiency, intersections are found without respect to which
9 specific triangle edge is intersected, but on a broader sense with
10 the horizontal, vertical, and diagonal seams in the grid, then
11 the intersections are ordered. If quadstrips are used for drawing
12 rather than tmesh, triangulation is not consistent; for diagonal
13 intersections, the proper diagonal to intersect would need to be
14 determined according to the algorithm used by qstrip (look at nearby
15 normals). It may be faster to go ahead and find the intersections
16 with the other diagonals using the same methods, then at sorting
17 time determine which diagonal array to look at for each quad.
18 It would also require a mechanism for throwing out unused intersections
19 with the diagonals during the ordering phase.
20 Do intersections in 2D, fill line structure with 3D pts (maybe calling
21 routine will cache for redrawing). Get Z value by using linear interp
22 between corners.
23
24 - check for easy cases:
25 - single point
26 - colinear with horizontal or vertical edges
27 - colinear with diagonal edges of triangles
28 - calculate three arrays of ordered intersections:
29 - with vertical edges
30 - with horizontal edges
31 - with diagonal edges and interpolate Z, using simple linear interpolation.
32 - eliminate duplicate intersections (need only compare one coord for each)
33 - build ordered set of points.
34
35 Return static pointer to 3D set of points. Max number of intersections
36 will be rows + cols + diags, so should allocate this number to initialize.
37 Let calling routine worry about copying points for caching.
38
39 (C) 1999-2008 by the GRASS Development Team
40
41 This program is free software under the
42 GNU General Public License (>=v2).
43 Read the file COPYING that comes with GRASS
44 for details.
45
46 \author Bill Brown UI GMS Lab
47 \author Doxygenized by Martin Landa <landa.martin gmail.com> (May 2008)
48 */
49
50 #include <stdlib.h>
51
52 #include <grass/ogsf.h>
53 #include <grass/glocale.h>
54
55 #include "gsget.h"
56 #include "rowcol.h"
57 #include "math.h"
58
59 #define DONT_INTERSECT 0
60 #define DO_INTERSECT 1
61 #define COLLINEAR 2
62
63 #define LERP(a,l,h) ((l)+(((h)-(l))*(a)))
64 #define EQUAL(a,b) (fabs((a)-(b))<EPSILON)
65 #define ISNODE(p,res) (fmod((double)p,(double)res)<EPSILON)
66
67 #define SAME_SIGNS( a, b ) \
68 ((a >= 0 && b >= 0) || (a < 0 && b < 0))
69
70 static int drape_line_init(int, int);
71 static Point3 *_gsdrape_get_segments(geosurf *, float *, float *, int *);
72 static float dist_squared_2d(float *, float *);
73
74 /* array of points to be returned */
75 static Point3 *I3d;
76
77 /* make dependent on resolution? */
78 static float EPSILON = 0.000001;
79
80 /*vertical, horizontal, & diagonal intersections */
81 static Point3 *Vi, *Hi, *Di;
82
83 static typbuff *Ebuf; /* elevation buffer */
84 static int Flat;
85
86
87 /*!
88 \brief Initizalize
89
90 \param rows number of rows
91 \param cols number of columns
92
93 \return -1 on failure
94 \return 1 on success
95 */
drape_line_init(int rows,int cols)96 static int drape_line_init(int rows, int cols)
97 {
98 /* use G_calloc() [-> G_fatal_error] instead of calloc ? */
99 if (NULL == (I3d = (Point3 *) calloc(2 * (rows + cols), sizeof(Point3)))) {
100 return (-1);
101 }
102
103 if (NULL == (Vi = (Point3 *) calloc(cols, sizeof(Point3)))) {
104 G_free(I3d);
105
106 return (-1);
107 }
108
109 if (NULL == (Hi = (Point3 *) calloc(rows, sizeof(Point3)))) {
110 G_free(I3d);
111 G_free(Vi);
112
113 return (-1);
114 }
115
116 if (NULL == (Di = (Point3 *) calloc(rows + cols, sizeof(Point3)))) {
117 G_free(I3d);
118 G_free(Vi);
119 G_free(Hi);
120
121 return (-1);
122 }
123
124 return (1);
125 }
126
127 /*!
128 \brief Get segments
129
130 \param gs surface (geosurf)
131 \param bgn begin point
132 \param end end point
133 \param num
134
135 \return pointer to Point3 struct
136 */
_gsdrape_get_segments(geosurf * gs,float * bgn,float * end,int * num)137 static Point3 *_gsdrape_get_segments(geosurf * gs, float *bgn, float *end,
138 int *num)
139 {
140 float f[3], l[3];
141 int vi, hi, di;
142 float dir[2], yres, xres;
143
144 xres = VXRES(gs);
145 yres = VYRES(gs);
146
147 vi = hi = di = 0;
148 GS_v2dir(bgn, end, dir);
149
150 if (dir[X]) {
151 vi = get_vert_intersects(gs, bgn, end, dir);
152 }
153
154 if (dir[Y]) {
155 hi = get_horz_intersects(gs, bgn, end, dir);
156 }
157
158 if (!((end[Y] - bgn[Y]) / (end[X] - bgn[X]) == yres / xres)) {
159 di = get_diag_intersects(gs, bgn, end, dir);
160 }
161
162 interp_first_last(gs, bgn, end, f, l);
163
164 *num = order_intersects(gs, f, l, vi, hi, di);
165 /* fills in return values, eliminates dupes (corners) */
166
167 G_debug(5, "_gsdrape_get_segments vi=%d, hi=%d, di=%d, num=%d",
168 vi, hi, di, *num);
169
170 return (I3d);
171 }
172
173 /*!
174 \brief Calculate 2D distance
175
176 \param p1 first point
177 \param p2 second point
178
179 \return distance
180 */
dist_squared_2d(float * p1,float * p2)181 static float dist_squared_2d(float *p1, float *p2)
182 {
183 float dx, dy;
184
185 dx = p2[X] - p1[X];
186 dy = p2[Y] - p1[Y];
187
188 return (dx * dx + dy * dy);
189 }
190
191 /*!
192 \brief ADD
193
194 \param gs surface (geosurf)
195
196 \return -1 on failure
197 \return 1 on success
198 */
gsdrape_set_surface(geosurf * gs)199 int gsdrape_set_surface(geosurf * gs)
200 {
201 static int first = 1;
202
203 if (first) {
204 first = 0;
205
206 if (0 > drape_line_init(gs->rows, gs->cols)) {
207 G_warning(_("Unable to process vector map - out of memory"));
208 Ebuf = NULL;
209
210 return (-1);
211 }
212 }
213
214 Ebuf = gs_get_att_typbuff(gs, ATT_TOPO, 0);
215
216 return (1);
217 }
218
219 /*!
220 \brief Check if segment intersect vector region
221
222 Clipping performed:
223 - bgn and end are replaced so that both points are within viewregion
224 - if seg intersects
225
226 \param gs surface (geosurf)
227 \param bgn begin point
228 \param end end point
229
230 \return 0 if segment doesn't intersect the viewregion, or intersects only at corner
231 \return otherwise returns 1
232 */
seg_intersect_vregion(geosurf * gs,float * bgn,float * end)233 int 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
259 (bgn[X], bgn[Y], end[X], end[Y], xl, yb, xl, yt, &xi, &yi)) {
260 /* left */
261 }
262 else if (segs_intersect
263 (bgn[X], bgn[Y], end[X], end[Y], xr, yb, xr, yt, &xi, &yi)) {
264 /* right */
265 }
266 else if (segs_intersect
267 (bgn[X], bgn[Y], end[X], end[Y], xl, yb, xr, yb, &xi, &yi)) {
268 /* bottom */
269 }
270 else if (segs_intersect
271 (bgn[X], bgn[Y], end[X], end[Y], xl, yt, xr, yt, &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
284 (bgn[X], bgn[Y], end[X], end[Y], xl, yb, xl, yt, &xi, &yi)) {
285 replace[X] = xi;
286 replace[Y] = yi;
287 replace = pt2;
288 inside++;
289 }
290
291 if (segs_intersect
292 (bgn[X], bgn[Y], end[X], end[Y], xr, yb, xr, yt, &xi, &yi)) {
293 replace[X] = xi;
294 replace[Y] = yi;
295 replace = pt2;
296 inside++;
297 }
298
299 if (inside < 2) {
300 if (segs_intersect
301 (bgn[X], bgn[Y], end[X], end[Y], xl, yb, xr, yb, &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
311 (bgn[X], bgn[Y], end[X], end[Y], xl, yt, xr, yt, &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 */
323 if (GS_P2distance(bgn, pt1) < GS_P2distance(bgn, pt2)) {
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 bgn begin point (x,y)
345 \param end end point (x,y)
346 \param num
347
348 \return pointer to Point3 struct
349 */
gsdrape_get_segments(geosurf * gs,float * bgn,float * end,int * num)350 Point3 *gsdrape_get_segments(geosurf * gs, float *bgn, float *end, int *num)
351 {
352 gsdrape_set_surface(gs);
353
354 if (!seg_intersect_vregion(gs, bgn, end)) {
355 *num = 0;
356
357 return (NULL);
358 }
359
360 if (CONST_ATT == gs_get_att_src(gs, ATT_TOPO)) {
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 /*!
392 \brief Get all segments
393
394 \param gs surface (geosurf)
395 \param bgn begin point
396 \param end end point
397 \param num
398
399 \return pointer to Point3 struct
400 */
gsdrape_get_allsegments(geosurf * gs,float * bgn,float * end,int * num)401 Point3 *gsdrape_get_allsegments(geosurf * gs, float *bgn, float *end,
402 int *num)
403 {
404 gsdrape_set_surface(gs);
405
406 if (!seg_intersect_vregion(gs, bgn, end)) {
407 *num = 0;
408 return (NULL);
409 }
410
411 if (bgn[X] == end[X] && bgn[Y] == end[Y]) {
412 float f[3], l[3];
413
414 interp_first_last(gs, bgn, end, f, l);
415 GS_v3eq(I3d[0], f);
416 GS_v3eq(I3d[1], l);
417 *num = 2;
418
419 return (I3d);
420 }
421
422 if (CONST_ATT == gs_get_att_src(gs, ATT_TOPO)) {
423 Flat = 1;
424 }
425 else {
426 Flat = 0;
427 }
428
429 return (_gsdrape_get_segments(gs, bgn, end, num));
430 }
431
432 /*!
433 \brief ADD
434
435 \param gs surface (geosurf)
436 \param bgn begin point
437 \param end end point
438 \param f first
439 \param l last
440 */
interp_first_last(geosurf * gs,float * bgn,float * end,Point3 f,Point3 l)441 void interp_first_last(geosurf * gs, float *bgn, float *end, Point3 f,
442 Point3 l)
443 {
444 f[X] = bgn[X];
445 f[Y] = bgn[Y];
446
447 l[X] = end[X];
448 l[Y] = end[Y];
449
450 if (Flat) {
451 f[Z] = l[Z] = gs->att[ATT_TOPO].constant;
452 }
453 else {
454 viewcell_tri_interp(gs, Ebuf, f, 0);
455 viewcell_tri_interp(gs, Ebuf, l, 0);
456 }
457
458 return;
459 }
460
461 /*!
462 \brief ADD
463
464 \param gs surface (geosurf)
465 \param pt
466 */
_viewcell_tri_interp(geosurf * gs,Point3 pt)467 int _viewcell_tri_interp(geosurf * gs, Point3 pt)
468 {
469 typbuff *buf;
470
471 buf = gs_get_att_typbuff(gs, ATT_TOPO, 0);
472
473 return (viewcell_tri_interp(gs, buf, pt, 0));
474 }
475
476 /*!
477 \brief ADD
478
479 In gsd_surf, tmesh draws polys like so:
480 <pre>
481 --------------
482 | /|
483 | / |
484 | / |
485 | / |
486 | / |
487 | / |
488 | / |
489 | / |
490 | / |
491 | / |
492 | / |
493 |/ |
494 --------------
495 </pre>
496
497 UNLESS the top right or bottom left point is masked, in which case a
498 single triangle with the opposite diagonal is drawn. This case is
499 not yet handled here & should only occur on edges.
500 pt has X & Y coordinates in it, we interpolate Z here
501
502 This could probably be much shorter, but not much faster.
503
504 \return 1 if point is in view region
505 \return otherwise 0 (if masked)
506 */
viewcell_tri_interp(geosurf * gs,typbuff * buf,Point3 pt,int check_mask)507 int viewcell_tri_interp(geosurf * gs, typbuff * buf, Point3 pt,
508 int check_mask)
509 {
510 Point3 p1, p2, p3;
511 int offset, drow, dcol, vrow, vcol;
512 float xmax, ymin, ymax, alpha;
513
514 xmax = VCOL2X(gs, VCOLS(gs));
515 ymax = VROW2Y(gs, 0);
516 ymin = VROW2Y(gs, VROWS(gs));
517
518 if (check_mask) {
519 if (gs_point_is_masked(gs, pt)) {
520 return (0);
521 }
522 }
523
524 if (pt[X] < 0.0 || pt[Y] > ymax) {
525 /* outside on left or top */
526 return (0);
527 }
528
529 if (pt[Y] < ymin || pt[X] > xmax) {
530 /* outside on bottom or right */
531 return (0);
532 }
533
534 if (CONST_ATT == gs_get_att_src(gs, ATT_TOPO)) {
535 pt[Z] = gs->att[ATT_TOPO].constant;
536
537 return (1);
538 }
539 else if (MAP_ATT != gs_get_att_src(gs, ATT_TOPO)) {
540 return (0);
541 }
542
543 vrow = Y2VROW(gs, pt[Y]);
544 vcol = X2VCOL(gs, pt[X]);
545
546 if (vrow < VROWS(gs) && vcol < VCOLS(gs)) {
547 /*not on bottom or right edge */
548 if (pt[X] > 0.0 && pt[Y] < ymax) {
549 /* not on left or top edge */
550 p1[X] = VCOL2X(gs, vcol + 1);
551 p1[Y] = VROW2Y(gs, vrow);
552 drow = VROW2DROW(gs, vrow);
553 dcol = VCOL2DCOL(gs, vcol + 1);
554 offset = DRC2OFF(gs, drow, dcol);
555 GET_MAPATT(buf, offset, p1[Z]); /* top right */
556
557 p2[X] = VCOL2X(gs, vcol);
558 p2[Y] = VROW2Y(gs, vrow + 1);
559 drow = VROW2DROW(gs, vrow + 1);
560 dcol = VCOL2DCOL(gs, vcol);
561 offset = DRC2OFF(gs, drow, dcol);
562 GET_MAPATT(buf, offset, p2[Z]); /* bottom left */
563
564 if ((pt[X] - p2[X]) / VXRES(gs) > (pt[Y] - p2[Y]) / VYRES(gs)) {
565 /* lower triangle */
566 p3[X] = VCOL2X(gs, vcol + 1);
567 p3[Y] = VROW2Y(gs, vrow + 1);
568 drow = VROW2DROW(gs, vrow + 1);
569 dcol = VCOL2DCOL(gs, vcol + 1);
570 offset = DRC2OFF(gs, drow, dcol);
571 GET_MAPATT(buf, offset, p3[Z]); /* bottom right */
572 }
573 else {
574 /* upper triangle */
575 p3[X] = VCOL2X(gs, vcol);
576 p3[Y] = VROW2Y(gs, vrow);
577 drow = VROW2DROW(gs, vrow);
578 dcol = VCOL2DCOL(gs, vcol);
579 offset = DRC2OFF(gs, drow, dcol);
580 GET_MAPATT(buf, offset, p3[Z]); /* top left */
581 }
582
583 return (Point_on_plane(p1, p2, p3, pt));
584 }
585 else if (pt[X] == 0.0) {
586 /* on left edge */
587 if (pt[Y] < ymax) {
588 vrow = Y2VROW(gs, pt[Y]);
589 drow = VROW2DROW(gs, vrow);
590 offset = DRC2OFF(gs, drow, 0);
591 GET_MAPATT(buf, offset, p1[Z]);
592
593 drow = VROW2DROW(gs, vrow + 1);
594 offset = DRC2OFF(gs, drow, 0);
595 GET_MAPATT(buf, offset, p2[Z]);
596
597 alpha = (VROW2Y(gs, vrow) - pt[Y]) / VYRES(gs);
598 pt[Z] = LERP(alpha, p1[Z], p2[Z]);
599 }
600 else {
601 /* top left corner */
602 GET_MAPATT(buf, 0, pt[Z]);
603 }
604
605 return (1);
606 }
607 else if (pt[Y] == gs->yrange) {
608 /* on top edge, not a corner */
609 vcol = X2VCOL(gs, pt[X]);
610 dcol = VCOL2DCOL(gs, vcol);
611 GET_MAPATT(buf, dcol, p1[Z]);
612
613 dcol = VCOL2DCOL(gs, vcol + 1);
614 GET_MAPATT(buf, dcol, p2[Z]);
615
616 alpha = (pt[X] - VCOL2X(gs, vcol)) / VXRES(gs);
617 pt[Z] = LERP(alpha, p1[Z], p2[Z]);
618
619 return (1);
620 }
621 }
622 else if (vrow == VROWS(gs)) {
623 /* on bottom edge */
624 drow = VROW2DROW(gs, VROWS(gs));
625
626 if (pt[X] > 0.0 && pt[X] < xmax) {
627 /* not a corner */
628 vcol = X2VCOL(gs, pt[X]);
629 dcol = VCOL2DCOL(gs, vcol);
630 offset = DRC2OFF(gs, drow, dcol);
631 GET_MAPATT(buf, offset, p1[Z]);
632
633 dcol = VCOL2DCOL(gs, vcol + 1);
634 offset = DRC2OFF(gs, drow, dcol);
635 GET_MAPATT(buf, offset, p2[Z]);
636
637 alpha = (pt[X] - VCOL2X(gs, vcol)) / VXRES(gs);
638 pt[Z] = LERP(alpha, p1[Z], p2[Z]);
639
640 return (1);
641 }
642 else if (pt[X] == 0.0) {
643 /* bottom left corner */
644 offset = DRC2OFF(gs, drow, 0);
645 GET_MAPATT(buf, offset, pt[Z]);
646
647 return (1);
648 }
649 else {
650 /* bottom right corner */
651 dcol = VCOL2DCOL(gs, VCOLS(gs));
652 offset = DRC2OFF(gs, drow, dcol);
653 GET_MAPATT(buf, offset, pt[Z]);
654
655 return (1);
656 }
657 }
658 else {
659 /* on right edge, not bottom corner */
660 dcol = VCOL2DCOL(gs, VCOLS(gs));
661
662 if (pt[Y] < ymax) {
663 vrow = Y2VROW(gs, pt[Y]);
664 drow = VROW2DROW(gs, vrow);
665 offset = DRC2OFF(gs, drow, dcol);
666 GET_MAPATT(buf, offset, p1[Z]);
667
668 drow = VROW2DROW(gs, vrow + 1);
669 offset = DRC2OFF(gs, drow, dcol);
670 GET_MAPATT(buf, offset, p2[Z]);
671
672 alpha = (VROW2Y(gs, vrow) - pt[Y]) / VYRES(gs);
673 pt[Z] = LERP(alpha, p1[Z], p2[Z]);
674
675 return (1);
676 }
677 else {
678 /* top right corner */
679 GET_MAPATT(buf, dcol, pt[Z]);
680
681 return (1);
682 }
683 }
684
685 return (0);
686 }
687
688 /*!
689 \brief ADD
690
691 \param gs surface (geosurf)
692
693 \return 1
694 \return 0
695 */
in_vregion(geosurf * gs,float * pt)696 int in_vregion(geosurf * gs, float *pt)
697 {
698 if (pt[X] >= 0.0 && pt[Y] <= gs->yrange) {
699 if (pt[X] <= VCOL2X(gs, VCOLS(gs))) {
700 return (pt[Y] >= VROW2Y(gs, VROWS(gs)));
701 }
702 }
703
704 return (0);
705 }
706
707 /*!
708 \brief ADD
709
710 After all the intersections between the segment and triangle
711 edges have been found, they are in three lists. (intersections
712 with vertical, horizontal, and diagonal triangle edges)
713
714 Each list is ordered in space from first to last segment points,
715 but now the lists need to be woven together. This routine
716 starts with the first point of the segment and then checks the
717 next point in each list to find the closest, eliminating duplicates
718 along the way and storing the result in I3d.
719
720 \param gs surface (geosurf)
721 \param first first point
722 \param last last point
723 \param vi
724 \param hi
725 \param di
726
727 \return
728 */
order_intersects(geosurf * gs,Point3 first,Point3 last,int vi,int hi,int di)729 int order_intersects(geosurf * gs, Point3 first, Point3 last, int vi, int hi,
730 int di)
731 {
732 int num, i, found, cv, ch, cd, cnum;
733 float dv, dh, dd, big, cpoint[2];
734
735 cv = ch = cd = cnum = 0;
736 num = vi + hi + di;
737
738 cpoint[X] = first[X];
739 cpoint[Y] = first[Y];
740
741 if (in_vregion(gs, first)) {
742 I3d[cnum][X] = first[X];
743 I3d[cnum][Y] = first[Y];
744 I3d[cnum][Z] = first[Z];
745 cnum++;
746 }
747
748 /* TODO: big could still be less than first dist */
749 big = gs->yrange * gs->yrange + gs->xrange * gs->xrange; /*BIG distance */
750 dv = dh = dd = big;
751
752 for (i = 0; i < num; i = cv + ch + cd) {
753 if (cv < vi) {
754 dv = dist_squared_2d(Vi[cv], cpoint);
755
756 if (dv < EPSILON) {
757 cv++;
758 continue;
759 }
760 }
761 else {
762 dv = big;
763 }
764
765 if (ch < hi) {
766 dh = dist_squared_2d(Hi[ch], cpoint);
767
768 if (dh < EPSILON) {
769 ch++;
770 continue;
771 }
772 }
773 else {
774 dh = big;
775 }
776
777 if (cd < di) {
778 dd = dist_squared_2d(Di[cd], cpoint);
779
780 if (dd < EPSILON) {
781 cd++;
782 continue;
783 }
784 }
785 else {
786 dd = big;
787 }
788
789 found = 0;
790
791 if (cd < di) {
792 if (dd <= dv && dd <= dh) {
793 found = 1;
794 cpoint[X] = I3d[cnum][X] = Di[cd][X];
795 cpoint[Y] = I3d[cnum][Y] = Di[cd][Y];
796 I3d[cnum][Z] = Di[cd][Z];
797 cnum++;
798
799 if (EQUAL(dd, dv)) {
800 cv++;
801 }
802
803 if (EQUAL(dd, dh)) {
804 ch++;
805 }
806
807 cd++;
808 }
809 }
810
811 if (!found) {
812 if (cv < vi) {
813 if (dv <= dh) {
814 found = 1;
815 cpoint[X] = I3d[cnum][X] = Vi[cv][X];
816 cpoint[Y] = I3d[cnum][Y] = Vi[cv][Y];
817 I3d[cnum][Z] = Vi[cv][Z];
818 cnum++;
819
820 if (EQUAL(dv, dh)) {
821 ch++;
822 }
823
824 cv++;
825 }
826 }
827 }
828
829 if (!found) {
830 if (ch < hi) {
831 cpoint[X] = I3d[cnum][X] = Hi[ch][X];
832 cpoint[Y] = I3d[cnum][Y] = Hi[ch][Y];
833 I3d[cnum][Z] = Hi[ch][Z];
834 cnum++;
835 ch++;
836 }
837 }
838
839 if (i == cv + ch + cd) {
840 G_debug(5, "order_intersects(): stuck on %d", cnum);
841 G_debug(5, "order_intersects(): cv = %d, ch = %d, cd = %d", cv,
842 ch, cd);
843 G_debug(5, "order_intersects(): dv = %f, dh = %f, dd = %f", dv,
844 dh, dd);
845
846 break;
847 }
848 }
849
850 if (EQUAL(last[X], cpoint[X]) && EQUAL(last[Y], cpoint[Y])) {
851 return (cnum);
852 }
853
854 if (in_vregion(gs, last)) {
855 /* TODO: check for last point on corner ? */
856 I3d[cnum][X] = last[X];
857 I3d[cnum][Y] = last[Y];
858 I3d[cnum][Z] = last[Z];
859 ++cnum;
860 }
861
862 return (cnum);
863 }
864
865 /*!
866 \brief ADD
867
868 \todo For consistancy, need to decide how last row & last column are
869 displayed - would it look funny to always draw last row/col with
870 finer resolution if necessary, or would it be better to only show
871 full rows/cols?
872
873 Colinear already eliminated
874
875 \param gs surface (geosurf)
876 \param bgn begin point
877 \param end end point
878 \param dir direction
879
880 \return
881 */
get_vert_intersects(geosurf * gs,float * bgn,float * end,float * dir)882 int get_vert_intersects(geosurf * gs, float *bgn, float *end, float *dir)
883 {
884 int fcol, lcol, incr, hits, num, offset, drow1, drow2;
885 float xl, yb, xr, yt, z1, z2, alpha;
886 float xres, yres, xi, yi;
887 int bgncol, endcol, cols, rows;
888
889 xres = VXRES(gs);
890 yres = VYRES(gs);
891 cols = VCOLS(gs);
892 rows = VROWS(gs);
893
894 bgncol = X2VCOL(gs, bgn[X]);
895 endcol = X2VCOL(gs, end[X]);
896
897 if (bgncol > cols && endcol > cols) {
898 return 0;
899 }
900
901 if (bgncol == endcol) {
902 return 0;
903 }
904
905 fcol = dir[X] > 0 ? bgncol + 1 : bgncol;
906 lcol = dir[X] > 0 ? endcol : endcol + 1;
907
908 /* assuming only showing FULL cols */
909 incr = lcol - fcol > 0 ? 1 : -1;
910
911 while (fcol > cols || fcol < 0) {
912 fcol += incr;
913 }
914
915 while (lcol > cols || lcol < 0) {
916 lcol -= incr;
917 }
918
919 num = abs(lcol - fcol) + 1;
920
921 yb = gs->yrange - (yres * rows) - EPSILON;
922 yt = gs->yrange + EPSILON;
923
924 for (hits = 0; hits < num; hits++) {
925 xl = xr = VCOL2X(gs, fcol);
926
927 if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yt, xr, yb,
928 &xi, &yi)) {
929 Vi[hits][X] = xi;
930 Vi[hits][Y] = yi;
931
932 /* find data rows */
933 if (Flat) {
934 Vi[hits][Z] = gs->att[ATT_TOPO].constant;
935 }
936 else {
937 drow1 = Y2VROW(gs, Vi[hits][Y]) * gs->y_mod;
938 drow2 = (1 + Y2VROW(gs, Vi[hits][Y])) * gs->y_mod;
939
940 if (drow2 >= gs->rows) {
941 drow2 = gs->rows - 1; /*bottom edge */
942 }
943
944 alpha =
945 ((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 bgn begin point
973 \param end end point
974 \param dir
975
976 \return number of intersects
977 */
get_horz_intersects(geosurf * gs,float * bgn,float * end,float * dir)978 int 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, yres, xi, yi;
983 int bgnrow, endrow, rows, cols;
984
985 xres = VXRES(gs);
986 yres = VYRES(gs);
987 cols = VCOLS(gs);
988 rows = VROWS(gs);
989
990 bgnrow = Y2VROW(gs, bgn[Y]);
991 endrow = Y2VROW(gs, end[Y]);
992 if (bgnrow == endrow) {
993 return 0;
994 }
995
996 if (bgnrow > rows && endrow > rows) {
997 return 0;
998 }
999
1000 frow = dir[Y] > 0 ? bgnrow : bgnrow + 1;
1001 lrow = dir[Y] > 0 ? endrow + 1 : endrow;
1002
1003 /* assuming only showing FULL rows */
1004 incr = lrow - frow > 0 ? 1 : -1;
1005
1006 while (frow > rows || frow < 0) {
1007 frow += incr;
1008 }
1009
1010 while (lrow > rows || lrow < 0) {
1011 lrow -= incr;
1012 }
1013
1014 num = abs(lrow - frow) + 1;
1015
1016 xl = 0.0 - EPSILON;
1017 xr = xres * cols + EPSILON;
1018
1019 for (hits = 0; hits < num; hits++) {
1020 yb = yt = VROW2Y(gs, frow);
1021
1022 if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yt, xr, yb,
1023 &xi, &yi)) {
1024 Hi[hits][X] = xi;
1025 Hi[hits][Y] = yi;
1026
1027 /* find data cols */
1028 if (Flat) {
1029 Hi[hits][Z] = gs->att[ATT_TOPO].constant;
1030 }
1031 else {
1032 dcol1 = X2VCOL(gs, Hi[hits][X]) * gs->x_mod;
1033 dcol2 = (1 + X2VCOL(gs, Hi[hits][X])) * gs->x_mod;
1034
1035 if (dcol2 >= gs->cols) {
1036 dcol2 = gs->cols - 1; /* right edge */
1037 }
1038
1039 alpha = (Hi[hits][X] - (dcol1 * gs->xres)) / xres;
1040
1041 offset = DRC2OFF(gs, frow * gs->y_mod, dcol1);
1042 GET_MAPATT(Ebuf, offset, z1);
1043 offset = DRC2OFF(gs, frow * gs->y_mod, dcol2);
1044 GET_MAPATT(Ebuf, offset, z2);
1045 Hi[hits][Z] = LERP(alpha, z1, z2);
1046 }
1047 }
1048
1049 /* if they don't intersect, something's wrong! */
1050 /* should only happen on endpoint, so it will be added later */
1051 else {
1052 hits--;
1053 num--;
1054 }
1055
1056 frow += incr;
1057 }
1058
1059 return (hits);
1060 }
1061
1062 /*!
1063 \brief Get diagonal intersects
1064
1065 Colinear already eliminated
1066
1067 \param gs surface (geosurf)
1068 \param bgn begin point
1069 \param end end point
1070 \param dir ? (unused)
1071
1072 \return number of intersects
1073 */
get_diag_intersects(geosurf * gs,float * bgn,float * end,float * dir)1074 int get_diag_intersects(geosurf * gs, float *bgn, float *end, float *dir)
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,
1133 &xi, &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 x1,y1,x2,y2 coordinates of endpoints of one segment
1203 \param 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 */
segs_intersect(float x1,float y1,float x2,float y2,float x3,float y3,float x4,float y4,float * x,float * y)1210 int 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 p1,p2,p3 points defining plane
1285 \param unk point
1286
1287 \return 1 point on plane
1288 \return 0 point not on plane
1289 */
Point_on_plane(Point3 p1,Point3 p2,Point3 p3,Point3 unk)1290 int Point_on_plane(Point3 p1, Point3 p2, Point3 p3, Point3 unk)
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[out] intersect intersect coordinates
1307 \param plane plane definition
1308
1309 \return 0 doesn't intersect
1310 \return 1 intesects
1311 */
XY_intersect_plane(float * intersect,float * plane)1312 int 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 p1,p2,p3 three point on plane
1331 \param[out] plane plane definition
1332
1333 \return 1
1334 */
P3toPlane(Point3 p1,Point3 p2,Point3 p3,float * plane)1335 int 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 /*!
1359 \brief Get cross product
1360
1361 \param a,b,c
1362
1363 \return cross product c = a cross b
1364 */
V3Cross(Point3 a,Point3 b,Point3 c)1365 int V3Cross(Point3 a, Point3 b, Point3 c)
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 }
1373