GRASS Programmer's Manual
6.4.2(2012)
|
00001 00050 #include <stdlib.h> 00051 00052 #include <grass/gstypes.h> 00053 #include <grass/glocale.h> 00054 00055 #include "gsget.h" 00056 #include "rowcol.h" 00057 #include "math.h" 00058 00059 #define DONT_INTERSECT 0 00060 #define DO_INTERSECT 1 00061 #define COLLINEAR 2 00062 00063 #define LERP(a,l,h) ((l)+(((h)-(l))*(a))) 00064 #define EQUAL(a,b) (fabs((a)-(b))<EPSILON) 00065 #define ISNODE(p,res) (fmod((double)p,(double)res)<EPSILON) 00066 00067 #define SAME_SIGNS( a, b ) \ 00068 ((a >= 0 && b >= 0) || (a < 0 && b < 0)) 00069 00070 static int drape_line_init(int, int); 00071 static Point3 *_gsdrape_get_segments(geosurf *, float *, float *, int *); 00072 static float dist_squared_2d(float *, float *); 00073 00074 /* array of points to be returned */ 00075 static Point3 *I3d; 00076 00077 /* make dependent on resolution? */ 00078 static float EPSILON = 0.000001; 00079 00080 /*vertical, horizontal, & diagonal intersections */ 00081 static Point3 *Vi, *Hi, *Di; 00082 00083 static typbuff *Ebuf; /* elevation buffer */ 00084 static int Flat; 00085 00086 00096 static int drape_line_init(int rows, int cols) 00097 { 00098 /* use G_calloc() [-> G_fatal_error] instead of calloc ? */ 00099 if (NULL == (I3d = (Point3 *) calloc(2 * (rows + cols), sizeof(Point3)))) { 00100 return (-1); 00101 } 00102 00103 if (NULL == (Vi = (Point3 *) calloc(cols, sizeof(Point3)))) { 00104 G_free(I3d); 00105 00106 return (-1); 00107 } 00108 00109 if (NULL == (Hi = (Point3 *) calloc(rows, sizeof(Point3)))) { 00110 G_free(I3d); 00111 G_free(Vi); 00112 00113 return (-1); 00114 } 00115 00116 if (NULL == (Di = (Point3 *) calloc(rows + cols, sizeof(Point3)))) { 00117 G_free(I3d); 00118 G_free(Vi); 00119 G_free(Hi); 00120 00121 return (-1); 00122 } 00123 00124 return (1); 00125 } 00126 00137 static Point3 *_gsdrape_get_segments(geosurf * gs, float *bgn, float *end, 00138 int *num) 00139 { 00140 float f[3], l[3]; 00141 int vi, hi, di; 00142 float dir[2], yres, xres; 00143 00144 xres = VXRES(gs); 00145 yres = VYRES(gs); 00146 00147 vi = hi = di = 0; 00148 GS_v2dir(bgn, end, dir); 00149 00150 if (dir[X]) { 00151 vi = get_vert_intersects(gs, bgn, end, dir); 00152 } 00153 00154 if (dir[Y]) { 00155 hi = get_horz_intersects(gs, bgn, end, dir); 00156 } 00157 00158 if (!((end[Y] - bgn[Y]) / (end[X] - bgn[X]) == yres / xres)) { 00159 di = get_diag_intersects(gs, bgn, end, dir); 00160 } 00161 00162 interp_first_last(gs, bgn, end, f, l); 00163 00164 *num = order_intersects(gs, f, l, vi, hi, di); 00165 /* fills in return values, eliminates dupes (corners) */ 00166 00167 G_debug(5, "_gsdrape_get_segments vi=%d, hi=%d, di=%d, num=%d", 00168 vi, hi, di, *num); 00169 00170 return (I3d); 00171 } 00172 00181 static float dist_squared_2d(float *p1, float *p2) 00182 { 00183 float dx, dy; 00184 00185 dx = p2[X] - p1[X]; 00186 dy = p2[Y] - p1[Y]; 00187 00188 return (dx * dx + dy * dy); 00189 } 00190 00199 int gsdrape_set_surface(geosurf * gs) 00200 { 00201 static int first = 1; 00202 00203 if (first) { 00204 first = 0; 00205 00206 if (0 > drape_line_init(gs->rows, gs->cols)) { 00207 G_warning(_("Unable to process vector map - out of memory")); 00208 Ebuf = NULL; 00209 00210 return (-1); 00211 } 00212 } 00213 00214 Ebuf = gs_get_att_typbuff(gs, ATT_TOPO, 0); 00215 00216 return (1); 00217 } 00218 00233 int seg_intersect_vregion(geosurf * gs, float *bgn, float *end) 00234 { 00235 float *replace, xl, yb, xr, yt, xi, yi; 00236 int inside = 0; 00237 00238 xl = 0.0; 00239 xr = VCOL2X(gs, VCOLS(gs)); 00240 yt = VROW2Y(gs, 0); 00241 yb = VROW2Y(gs, VROWS(gs)); 00242 00243 if (in_vregion(gs, bgn)) { 00244 replace = end; 00245 inside++; 00246 } 00247 00248 if (in_vregion(gs, end)) { 00249 replace = bgn; 00250 inside++; 00251 } 00252 00253 if (inside == 2) { 00254 return (1); 00255 } 00256 else if (inside) { 00257 /* one in & one out - replace gets first intersection */ 00258 if (segs_intersect 00259 (bgn[X], bgn[Y], end[X], end[Y], xl, yb, xl, yt, &xi, &yi)) { 00260 /* left */ 00261 } 00262 else if (segs_intersect 00263 (bgn[X], bgn[Y], end[X], end[Y], xr, yb, xr, yt, &xi, &yi)) { 00264 /* right */ 00265 } 00266 else if (segs_intersect 00267 (bgn[X], bgn[Y], end[X], end[Y], xl, yb, xr, yb, &xi, &yi)) { 00268 /* bottom */ 00269 } 00270 else if (segs_intersect 00271 (bgn[X], bgn[Y], end[X], end[Y], xl, yt, xr, yt, &xi, &yi)) { 00272 /* top */ 00273 } 00274 00275 replace[X] = xi; 00276 replace[Y] = yi; 00277 } 00278 else { 00279 /* both out - find 2 intersects & replace both */ 00280 float pt1[2], pt2[2]; 00281 00282 replace = pt1; 00283 if (segs_intersect 00284 (bgn[X], bgn[Y], end[X], end[Y], xl, yb, xl, yt, &xi, &yi)) { 00285 replace[X] = xi; 00286 replace[Y] = yi; 00287 replace = pt2; 00288 inside++; 00289 } 00290 00291 if (segs_intersect 00292 (bgn[X], bgn[Y], end[X], end[Y], xr, yb, xr, yt, &xi, &yi)) { 00293 replace[X] = xi; 00294 replace[Y] = yi; 00295 replace = pt2; 00296 inside++; 00297 } 00298 00299 if (inside < 2) { 00300 if (segs_intersect 00301 (bgn[X], bgn[Y], end[X], end[Y], xl, yb, xr, yb, &xi, &yi)) { 00302 replace[X] = xi; 00303 replace[Y] = yi; 00304 replace = pt2; 00305 inside++; 00306 } 00307 } 00308 00309 if (inside < 2) { 00310 if (segs_intersect 00311 (bgn[X], bgn[Y], end[X], end[Y], xl, yt, xr, yt, &xi, &yi)) { 00312 replace[X] = xi; 00313 replace[Y] = yi; 00314 inside++; 00315 } 00316 } 00317 00318 if (inside < 2) { 00319 return (0); /* no intersect or only 1 point on corner */ 00320 } 00321 00322 /* compare dist of intersects to bgn - closest replaces bgn */ 00323 if (GS_P2distance(bgn, pt1) < GS_P2distance(bgn, pt2)) { 00324 bgn[X] = pt1[X]; 00325 bgn[Y] = pt1[Y]; 00326 end[X] = pt2[X]; 00327 end[Y] = pt2[Y]; 00328 } 00329 else { 00330 bgn[X] = pt2[X]; 00331 bgn[Y] = pt2[Y]; 00332 end[X] = pt1[X]; 00333 end[Y] = pt1[Y]; 00334 } 00335 } 00336 00337 return (1); 00338 } 00339 00350 Point3 *gsdrape_get_segments(geosurf * gs, float *bgn, float *end, int *num) 00351 { 00352 gsdrape_set_surface(gs); 00353 00354 if (!seg_intersect_vregion(gs, bgn, end)) { 00355 *num = 0; 00356 00357 return (NULL); 00358 } 00359 00360 if (CONST_ATT == gs_get_att_src(gs, ATT_TOPO)) { 00361 /* will probably want a force_drape option to get all intersects */ 00362 I3d[0][X] = bgn[X]; 00363 I3d[0][Y] = bgn[Y]; 00364 I3d[0][Z] = gs->att[ATT_TOPO].constant; 00365 I3d[1][X] = end[X]; 00366 I3d[1][Y] = end[Y]; 00367 I3d[1][Z] = gs->att[ATT_TOPO].constant; 00368 *num = 2; 00369 00370 return (I3d); 00371 } 00372 00373 if (bgn[X] == end[X] && bgn[Y] == end[Y]) { 00374 float f[3], l[3]; 00375 00376 interp_first_last(gs, bgn, end, f, l); 00377 GS_v3eq(I3d[0], f); 00378 GS_v3eq(I3d[1], l); 00379 00380 /* CHANGE (*num = 1) to reflect degenerate line ? */ 00381 *num = 2; 00382 00383 return (I3d); 00384 } 00385 00386 Flat = 0; 00387 return (_gsdrape_get_segments(gs, bgn, end, num)); 00388 } 00389 00390 00401 Point3 *gsdrape_get_allsegments(geosurf * gs, float *bgn, float *end, 00402 int *num) 00403 { 00404 gsdrape_set_surface(gs); 00405 00406 if (!seg_intersect_vregion(gs, bgn, end)) { 00407 *num = 0; 00408 return (NULL); 00409 } 00410 00411 if (bgn[X] == end[X] && bgn[Y] == end[Y]) { 00412 float f[3], l[3]; 00413 00414 interp_first_last(gs, bgn, end, f, l); 00415 GS_v3eq(I3d[0], f); 00416 GS_v3eq(I3d[1], l); 00417 *num = 2; 00418 00419 return (I3d); 00420 } 00421 00422 if (CONST_ATT == gs_get_att_src(gs, ATT_TOPO)) { 00423 Flat = 1; 00424 } 00425 else { 00426 Flat = 0; 00427 } 00428 00429 return (_gsdrape_get_segments(gs, bgn, end, num)); 00430 } 00431 00441 void interp_first_last(geosurf * gs, float *bgn, float *end, Point3 f, 00442 Point3 l) 00443 { 00444 f[X] = bgn[X]; 00445 f[Y] = bgn[Y]; 00446 00447 l[X] = end[X]; 00448 l[Y] = end[Y]; 00449 00450 if (Flat) { 00451 f[Z] = l[Z] = gs->att[ATT_TOPO].constant; 00452 } 00453 else { 00454 viewcell_tri_interp(gs, Ebuf, f, 0); 00455 viewcell_tri_interp(gs, Ebuf, l, 0); 00456 } 00457 00458 return; 00459 } 00460 00467 int _viewcell_tri_interp(geosurf * gs, Point3 pt) 00468 { 00469 typbuff *buf; 00470 00471 buf = gs_get_att_typbuff(gs, ATT_TOPO, 0); 00472 00473 return (viewcell_tri_interp(gs, buf, pt, 0)); 00474 } 00475 00507 int viewcell_tri_interp(geosurf * gs, typbuff * buf, Point3 pt, 00508 int check_mask) 00509 { 00510 Point3 p1, p2, p3; 00511 int offset, drow, dcol, vrow, vcol; 00512 float xmax, ymin, ymax, alpha; 00513 00514 xmax = VCOL2X(gs, VCOLS(gs)); 00515 ymax = VROW2Y(gs, 0); 00516 ymin = VROW2Y(gs, VROWS(gs)); 00517 00518 if (check_mask) { 00519 if (gs_point_is_masked(gs, pt)) { 00520 return (0); 00521 } 00522 } 00523 00524 if (pt[X] < 0.0 || pt[Y] > ymax) { 00525 /* outside on left or top */ 00526 return (0); 00527 } 00528 00529 if (pt[Y] < ymin || pt[X] > xmax) { 00530 /* outside on bottom or right */ 00531 return (0); 00532 } 00533 00534 if (CONST_ATT == gs_get_att_src(gs, ATT_TOPO)) { 00535 pt[Z] = gs->att[ATT_TOPO].constant; 00536 00537 return (1); 00538 } 00539 else if (MAP_ATT != gs_get_att_src(gs, ATT_TOPO)) { 00540 return (0); 00541 } 00542 00543 vrow = Y2VROW(gs, pt[Y]); 00544 vcol = X2VCOL(gs, pt[X]); 00545 00546 if (vrow < VROWS(gs) && vcol < VCOLS(gs)) { 00547 /*not on bottom or right edge */ 00548 if (pt[X] > 0.0 && pt[Y] < ymax) { 00549 /* not on left or top edge */ 00550 p1[X] = VCOL2X(gs, vcol + 1); 00551 p1[Y] = VROW2Y(gs, vrow); 00552 drow = VROW2DROW(gs, vrow); 00553 dcol = VCOL2DCOL(gs, vcol + 1); 00554 offset = DRC2OFF(gs, drow, dcol); 00555 GET_MAPATT(buf, offset, p1[Z]); /* top right */ 00556 00557 p2[X] = VCOL2X(gs, vcol); 00558 p2[Y] = VROW2Y(gs, vrow + 1); 00559 drow = VROW2DROW(gs, vrow + 1); 00560 dcol = VCOL2DCOL(gs, vcol); 00561 offset = DRC2OFF(gs, drow, dcol); 00562 GET_MAPATT(buf, offset, p2[Z]); /* bottom left */ 00563 00564 if ((pt[X] - p2[X]) / VXRES(gs) > (pt[Y] - p2[Y]) / VYRES(gs)) { 00565 /* lower triangle */ 00566 p3[X] = VCOL2X(gs, vcol + 1); 00567 p3[Y] = VROW2Y(gs, vrow + 1); 00568 drow = VROW2DROW(gs, vrow + 1); 00569 dcol = VCOL2DCOL(gs, vcol + 1); 00570 offset = DRC2OFF(gs, drow, dcol); 00571 GET_MAPATT(buf, offset, p3[Z]); /* bottom right */ 00572 } 00573 else { 00574 /* upper triangle */ 00575 p3[X] = VCOL2X(gs, vcol); 00576 p3[Y] = VROW2Y(gs, vrow); 00577 drow = VROW2DROW(gs, vrow); 00578 dcol = VCOL2DCOL(gs, vcol); 00579 offset = DRC2OFF(gs, drow, dcol); 00580 GET_MAPATT(buf, offset, p3[Z]); /* top left */ 00581 } 00582 00583 return (Point_on_plane(p1, p2, p3, pt)); 00584 } 00585 else if (pt[X] == 0.0) { 00586 /* on left edge */ 00587 if (pt[Y] < ymax) { 00588 vrow = Y2VROW(gs, pt[Y]); 00589 drow = VROW2DROW(gs, vrow); 00590 offset = DRC2OFF(gs, drow, 0); 00591 GET_MAPATT(buf, offset, p1[Z]); 00592 00593 drow = VROW2DROW(gs, vrow + 1); 00594 offset = DRC2OFF(gs, drow, 0); 00595 GET_MAPATT(buf, offset, p2[Z]); 00596 00597 alpha = (VROW2Y(gs, vrow) - pt[Y]) / VYRES(gs); 00598 pt[Z] = LERP(alpha, p1[Z], p2[Z]); 00599 } 00600 else { 00601 /* top left corner */ 00602 GET_MAPATT(buf, 0, pt[Z]); 00603 } 00604 00605 return (1); 00606 } 00607 else if (pt[Y] == gs->yrange) { 00608 /* on top edge, not a corner */ 00609 vcol = X2VCOL(gs, pt[X]); 00610 dcol = VCOL2DCOL(gs, vcol); 00611 GET_MAPATT(buf, dcol, p1[Z]); 00612 00613 dcol = VCOL2DCOL(gs, vcol + 1); 00614 GET_MAPATT(buf, dcol, p2[Z]); 00615 00616 alpha = (pt[X] - VCOL2X(gs, vcol)) / VXRES(gs); 00617 pt[Z] = LERP(alpha, p1[Z], p2[Z]); 00618 00619 return (1); 00620 } 00621 } 00622 else if (vrow == VROWS(gs)) { 00623 /* on bottom edge */ 00624 drow = VROW2DROW(gs, VROWS(gs)); 00625 00626 if (pt[X] > 0.0 && pt[X] < xmax) { 00627 /* not a corner */ 00628 vcol = X2VCOL(gs, pt[X]); 00629 dcol = VCOL2DCOL(gs, vcol); 00630 offset = DRC2OFF(gs, drow, dcol); 00631 GET_MAPATT(buf, offset, p1[Z]); 00632 00633 dcol = VCOL2DCOL(gs, vcol + 1); 00634 offset = DRC2OFF(gs, drow, dcol); 00635 GET_MAPATT(buf, offset, p2[Z]); 00636 00637 alpha = (pt[X] - VCOL2X(gs, vcol)) / VXRES(gs); 00638 pt[Z] = LERP(alpha, p1[Z], p2[Z]); 00639 00640 return (1); 00641 } 00642 else if (pt[X] == 0.0) { 00643 /* bottom left corner */ 00644 offset = DRC2OFF(gs, drow, 0); 00645 GET_MAPATT(buf, offset, pt[Z]); 00646 00647 return (1); 00648 } 00649 else { 00650 /* bottom right corner */ 00651 dcol = VCOL2DCOL(gs, VCOLS(gs)); 00652 offset = DRC2OFF(gs, drow, dcol); 00653 GET_MAPATT(buf, offset, pt[Z]); 00654 00655 return (1); 00656 } 00657 } 00658 else { 00659 /* on right edge, not bottom corner */ 00660 dcol = VCOL2DCOL(gs, VCOLS(gs)); 00661 00662 if (pt[Y] < ymax) { 00663 vrow = Y2VROW(gs, pt[Y]); 00664 drow = VROW2DROW(gs, vrow); 00665 offset = DRC2OFF(gs, drow, dcol); 00666 GET_MAPATT(buf, offset, p1[Z]); 00667 00668 drow = VROW2DROW(gs, vrow + 1); 00669 offset = DRC2OFF(gs, drow, dcol); 00670 GET_MAPATT(buf, offset, p2[Z]); 00671 00672 alpha = (VROW2Y(gs, vrow) - pt[Y]) / VYRES(gs); 00673 pt[Z] = LERP(alpha, p1[Z], p2[Z]); 00674 00675 return (1); 00676 } 00677 else { 00678 /* top right corner */ 00679 GET_MAPATT(buf, dcol, pt[Z]); 00680 00681 return (1); 00682 } 00683 } 00684 00685 return (0); 00686 } 00687 00696 int in_vregion(geosurf * gs, float *pt) 00697 { 00698 if (pt[X] >= 0.0 && pt[Y] <= gs->yrange) { 00699 if (pt[X] <= VCOL2X(gs, VCOLS(gs))) { 00700 return (pt[Y] >= VROW2Y(gs, VROWS(gs))); 00701 } 00702 } 00703 00704 return (0); 00705 } 00706 00729 int order_intersects(geosurf * gs, Point3 first, Point3 last, int vi, int hi, 00730 int di) 00731 { 00732 int num, i, found, cv, ch, cd, cnum; 00733 float dv, dh, dd, big, cpoint[2]; 00734 00735 cv = ch = cd = cnum = 0; 00736 num = vi + hi + di; 00737 00738 cpoint[X] = first[X]; 00739 cpoint[Y] = first[Y]; 00740 00741 if (in_vregion(gs, first)) { 00742 I3d[cnum][X] = first[X]; 00743 I3d[cnum][Y] = first[Y]; 00744 I3d[cnum][Z] = first[Z]; 00745 cnum++; 00746 } 00747 00748 /* TODO: big could still be less than first dist */ 00749 big = gs->yrange * gs->yrange + gs->xrange * gs->xrange; /*BIG distance */ 00750 dv = dh = dd = big; 00751 00752 for (i = 0; i < num; i = cv + ch + cd) { 00753 if (cv < vi) { 00754 dv = dist_squared_2d(Vi[cv], cpoint); 00755 00756 if (dv < EPSILON) { 00757 cv++; 00758 continue; 00759 } 00760 } 00761 else { 00762 dv = big; 00763 } 00764 00765 if (ch < hi) { 00766 dh = dist_squared_2d(Hi[ch], cpoint); 00767 00768 if (dh < EPSILON) { 00769 ch++; 00770 continue; 00771 } 00772 } 00773 else { 00774 dh = big; 00775 } 00776 00777 if (cd < di) { 00778 dd = dist_squared_2d(Di[cd], cpoint); 00779 00780 if (dd < EPSILON) { 00781 cd++; 00782 continue; 00783 } 00784 } 00785 else { 00786 dd = big; 00787 } 00788 00789 found = 0; 00790 00791 if (cd < di) { 00792 if (dd <= dv && dd <= dh) { 00793 found = 1; 00794 cpoint[X] = I3d[cnum][X] = Di[cd][X]; 00795 cpoint[Y] = I3d[cnum][Y] = Di[cd][Y]; 00796 I3d[cnum][Z] = Di[cd][Z]; 00797 cnum++; 00798 00799 if (EQUAL(dd, dv)) { 00800 cv++; 00801 } 00802 00803 if (EQUAL(dd, dh)) { 00804 ch++; 00805 } 00806 00807 cd++; 00808 } 00809 } 00810 00811 if (!found) { 00812 if (cv < vi) { 00813 if (dv <= dh) { 00814 found = 1; 00815 cpoint[X] = I3d[cnum][X] = Vi[cv][X]; 00816 cpoint[Y] = I3d[cnum][Y] = Vi[cv][Y]; 00817 I3d[cnum][Z] = Vi[cv][Z]; 00818 cnum++; 00819 00820 if (EQUAL(dv, dh)) { 00821 ch++; 00822 } 00823 00824 cv++; 00825 } 00826 } 00827 } 00828 00829 if (!found) { 00830 if (ch < hi) { 00831 cpoint[X] = I3d[cnum][X] = Hi[ch][X]; 00832 cpoint[Y] = I3d[cnum][Y] = Hi[ch][Y]; 00833 I3d[cnum][Z] = Hi[ch][Z]; 00834 cnum++; 00835 ch++; 00836 } 00837 } 00838 00839 if (i == cv + ch + cd) { 00840 G_debug(5, "order_intersects(): stuck on %d", cnum); 00841 G_debug(5, "order_intersects(): cv = %d, ch = %d, cd = %d", cv, 00842 ch, cd); 00843 G_debug(5, "order_intersects(): dv = %f, dh = %f, dd = %f", dv, 00844 dh, dd); 00845 00846 break; 00847 } 00848 } 00849 00850 if (EQUAL(last[X], cpoint[X]) && EQUAL(last[Y], cpoint[Y])) { 00851 return (cnum); 00852 } 00853 00854 if (in_vregion(gs, last)) { 00855 /* TODO: check for last point on corner ? */ 00856 I3d[cnum][X] = last[X]; 00857 I3d[cnum][Y] = last[Y]; 00858 I3d[cnum][Z] = last[Z]; 00859 ++cnum; 00860 } 00861 00862 return (cnum); 00863 } 00864 00882 int get_vert_intersects(geosurf * gs, float *bgn, float *end, float *dir) 00883 { 00884 int fcol, lcol, incr, hits, num, offset, drow1, drow2; 00885 float xl, yb, xr, yt, z1, z2, alpha; 00886 float xres, yres, xi, yi; 00887 int bgncol, endcol, cols, rows; 00888 00889 xres = VXRES(gs); 00890 yres = VYRES(gs); 00891 cols = VCOLS(gs); 00892 rows = VROWS(gs); 00893 00894 bgncol = X2VCOL(gs, bgn[X]); 00895 endcol = X2VCOL(gs, end[X]); 00896 00897 if (bgncol > cols && endcol > cols) { 00898 return 0; 00899 } 00900 00901 if (bgncol == endcol) { 00902 return 0; 00903 } 00904 00905 fcol = dir[X] > 0 ? bgncol + 1 : bgncol; 00906 lcol = dir[X] > 0 ? endcol : endcol + 1; 00907 00908 /* assuming only showing FULL cols */ 00909 incr = lcol - fcol > 0 ? 1 : -1; 00910 00911 while (fcol > cols || fcol < 0) { 00912 fcol += incr; 00913 } 00914 00915 while (lcol > cols || lcol < 0) { 00916 lcol -= incr; 00917 } 00918 00919 num = abs(lcol - fcol) + 1; 00920 00921 yb = gs->yrange - (yres * rows) - EPSILON; 00922 yt = gs->yrange + EPSILON; 00923 00924 for (hits = 0; hits < num; hits++) { 00925 xl = xr = VCOL2X(gs, fcol); 00926 00927 if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yt, xr, yb, 00928 &xi, &yi)) { 00929 Vi[hits][X] = xi; 00930 Vi[hits][Y] = yi; 00931 00932 /* find data rows */ 00933 if (Flat) { 00934 Vi[hits][Z] = gs->att[ATT_TOPO].constant; 00935 } 00936 else { 00937 drow1 = Y2VROW(gs, Vi[hits][Y]) * gs->y_mod; 00938 drow2 = (1 + Y2VROW(gs, Vi[hits][Y])) * gs->y_mod; 00939 00940 if (drow2 >= gs->rows) { 00941 drow2 = gs->rows - 1; /*bottom edge */ 00942 } 00943 00944 alpha = 00945 ((gs->yrange - drow1 * gs->yres) - Vi[hits][Y]) / yres; 00946 00947 offset = DRC2OFF(gs, drow1, fcol * gs->x_mod); 00948 GET_MAPATT(Ebuf, offset, z1); 00949 offset = DRC2OFF(gs, drow2, fcol * gs->x_mod); 00950 GET_MAPATT(Ebuf, offset, z2); 00951 Vi[hits][Z] = LERP(alpha, z1, z2); 00952 } 00953 } 00954 00955 /* if they don't intersect, something's wrong! */ 00956 /* should only happen on endpoint, so it will be added later */ 00957 else { 00958 hits--; 00959 num--; 00960 } 00961 00962 fcol += incr; 00963 } 00964 00965 return (hits); 00966 } 00967 00978 int get_horz_intersects(geosurf * gs, float *bgn, float *end, float *dir) 00979 { 00980 int frow, lrow, incr, hits, num, offset, dcol1, dcol2; 00981 float xl, yb, xr, yt, z1, z2, alpha; 00982 float xres, yres, xi, yi; 00983 int bgnrow, endrow, rows, cols; 00984 00985 xres = VXRES(gs); 00986 yres = VYRES(gs); 00987 cols = VCOLS(gs); 00988 rows = VROWS(gs); 00989 00990 bgnrow = Y2VROW(gs, bgn[Y]); 00991 endrow = Y2VROW(gs, end[Y]); 00992 if (bgnrow == endrow) { 00993 return 0; 00994 } 00995 00996 if (bgnrow > rows && endrow > rows) { 00997 return 0; 00998 } 00999 01000 frow = dir[Y] > 0 ? bgnrow : bgnrow + 1; 01001 lrow = dir[Y] > 0 ? endrow + 1 : endrow; 01002 01003 /* assuming only showing FULL rows */ 01004 incr = lrow - frow > 0 ? 1 : -1; 01005 01006 while (frow > rows || frow < 0) { 01007 frow += incr; 01008 } 01009 01010 while (lrow > rows || lrow < 0) { 01011 lrow -= incr; 01012 } 01013 01014 num = abs(lrow - frow) + 1; 01015 01016 xl = 0.0 - EPSILON; 01017 xr = xres * cols + EPSILON; 01018 01019 for (hits = 0; hits < num; hits++) { 01020 yb = yt = VROW2Y(gs, frow); 01021 01022 if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yt, xr, yb, 01023 &xi, &yi)) { 01024 Hi[hits][X] = xi; 01025 Hi[hits][Y] = yi; 01026 01027 /* find data cols */ 01028 if (Flat) { 01029 Hi[hits][Z] = gs->att[ATT_TOPO].constant; 01030 } 01031 else { 01032 dcol1 = X2VCOL(gs, Hi[hits][X]) * gs->x_mod; 01033 dcol2 = (1 + X2VCOL(gs, Hi[hits][X])) * gs->x_mod; 01034 01035 if (dcol2 >= gs->cols) { 01036 dcol2 = gs->cols - 1; /* right edge */ 01037 } 01038 01039 alpha = (Hi[hits][X] - (dcol1 * gs->xres)) / xres; 01040 01041 offset = DRC2OFF(gs, frow * gs->y_mod, dcol1); 01042 GET_MAPATT(Ebuf, offset, z1); 01043 offset = DRC2OFF(gs, frow * gs->y_mod, dcol2); 01044 GET_MAPATT(Ebuf, offset, z2); 01045 Hi[hits][Z] = LERP(alpha, z1, z2); 01046 } 01047 } 01048 01049 /* if they don't intersect, something's wrong! */ 01050 /* should only happen on endpoint, so it will be added later */ 01051 else { 01052 hits--; 01053 num--; 01054 } 01055 01056 frow += incr; 01057 } 01058 01059 return (hits); 01060 } 01061 01074 int get_diag_intersects(geosurf * gs, float *bgn, float *end, float *dir) 01075 { 01076 int fdig, ldig, incr, hits, num, offset; 01077 int vrow, vcol, drow1, drow2, dcol1, dcol2; 01078 float xl, yb, xr, yt, z1, z2, alpha; 01079 float xres, yres, xi, yi, dx, dy; 01080 int diags, cols, rows, lower; 01081 Point3 pt; 01082 01083 xres = VXRES(gs); 01084 yres = VYRES(gs); 01085 cols = VCOLS(gs); 01086 rows = VROWS(gs); 01087 diags = rows + cols; /* -1 ? */ 01088 01089 /* determine upper/lower triangle for last */ 01090 vrow = Y2VROW(gs, end[Y]); 01091 vcol = X2VCOL(gs, end[X]); 01092 pt[X] = VCOL2X(gs, vcol); 01093 pt[Y] = VROW2Y(gs, vrow + 1); 01094 lower = ((end[X] - pt[X]) / xres > (end[Y] - pt[Y]) / yres); 01095 ldig = lower ? vrow + vcol + 1 : vrow + vcol; 01096 01097 /* determine upper/lower triangle for first */ 01098 vrow = Y2VROW(gs, bgn[Y]); 01099 vcol = X2VCOL(gs, bgn[X]); 01100 pt[X] = VCOL2X(gs, vcol); 01101 pt[Y] = VROW2Y(gs, vrow + 1); 01102 lower = ((bgn[X] - pt[X]) / xres > (bgn[Y] - pt[Y]) / yres); 01103 fdig = lower ? vrow + vcol + 1 : vrow + vcol; 01104 01105 /* adjust according to direction */ 01106 if (ldig > fdig) { 01107 fdig++; 01108 } 01109 01110 if (fdig > ldig) { 01111 ldig++; 01112 } 01113 01114 incr = ldig - fdig > 0 ? 1 : -1; 01115 01116 while (fdig > diags || fdig < 0) { 01117 fdig += incr; 01118 } 01119 01120 while (ldig > diags || ldig < 0) { 01121 ldig -= incr; 01122 } 01123 01124 num = abs(ldig - fdig) + 1; 01125 01126 for (hits = 0; hits < num; hits++) { 01127 yb = gs->yrange - (yres * (fdig < rows ? fdig : rows)) - EPSILON; 01128 xl = VCOL2X(gs, (fdig < rows ? 0 : fdig - rows)) - EPSILON; 01129 yt = gs->yrange - (yres * (fdig < cols ? 0 : fdig - cols)) + EPSILON; 01130 xr = VCOL2X(gs, (fdig < cols ? fdig : cols)) + EPSILON; 01131 01132 if (segs_intersect(bgn[X], bgn[Y], end[X], end[Y], xl, yb, xr, yt, 01133 &xi, &yi)) { 01134 Di[hits][X] = xi; 01135 Di[hits][Y] = yi; 01136 01137 if (ISNODE(xi, xres)) { 01138 /* then it's also a ynode */ 01139 num--; 01140 hits--; 01141 continue; 01142 } 01143 01144 /* find data rows */ 01145 drow1 = Y2VROW(gs, Di[hits][Y]) * gs->y_mod; 01146 drow2 = (1 + Y2VROW(gs, Di[hits][Y])) * gs->y_mod; 01147 01148 if (drow2 >= gs->rows) { 01149 drow2 = gs->rows - 1; /* bottom edge */ 01150 } 01151 01152 /* find data cols */ 01153 if (Flat) { 01154 Di[hits][Z] = gs->att[ATT_TOPO].constant; 01155 } 01156 else { 01157 dcol1 = X2VCOL(gs, Di[hits][X]) * gs->x_mod; 01158 dcol2 = (1 + X2VCOL(gs, Di[hits][X])) * gs->x_mod; 01159 01160 if (dcol2 >= gs->cols) { 01161 dcol2 = gs->cols - 1; /* right edge */ 01162 } 01163 01164 dx = DCOL2X(gs, dcol2) - Di[hits][X]; 01165 dy = DROW2Y(gs, drow1) - Di[hits][Y]; 01166 alpha = 01167 sqrt(dx * dx + dy * dy) / sqrt(xres * xres + yres * yres); 01168 01169 offset = DRC2OFF(gs, drow1, dcol2); 01170 GET_MAPATT(Ebuf, offset, z1); 01171 offset = DRC2OFF(gs, drow2, dcol1); 01172 GET_MAPATT(Ebuf, offset, z2); 01173 Di[hits][Z] = LERP(alpha, z1, z2); 01174 } 01175 } 01176 01177 /* if they don't intersect, something's wrong! */ 01178 /* should only happen on endpoint, so it will be added later */ 01179 else { 01180 hits--; 01181 num--; 01182 } 01183 01184 fdig += incr; 01185 } 01186 01187 return (hits); 01188 } 01189 01210 int segs_intersect(float x1, float y1, float x2, float y2, float x3, float y3, 01211 float x4, float y4, float *x, float *y) 01212 { 01213 float a1, a2, b1, b2, c1, c2; /* Coefficients of line eqns. */ 01214 float r1, r2, r3, r4; /* 'Sign' values */ 01215 float denom, offset, num; /* Intermediate values */ 01216 01217 /* Compute a1, b1, c1, where line joining points 1 and 2 01218 * is "a1 x + b1 y + c1 = 0". 01219 */ 01220 a1 = y2 - y1; 01221 b1 = x1 - x2; 01222 c1 = x2 * y1 - x1 * y2; 01223 01224 /* Compute r3 and r4. 01225 */ 01226 r3 = a1 * x3 + b1 * y3 + c1; 01227 r4 = a1 * x4 + b1 * y4 + c1; 01228 01229 /* Check signs of r3 and r4. If both point 3 and point 4 lie on 01230 * same side of line 1, the line segments do not intersect. 01231 */ 01232 01233 if (!EQUAL(r3, 0.0) && !EQUAL(r4, 0.0) && SAME_SIGNS(r3, r4)) { 01234 return (DONT_INTERSECT); 01235 } 01236 01237 /* Compute a2, b2, c2 */ 01238 a2 = y4 - y3; 01239 b2 = x3 - x4; 01240 c2 = x4 * y3 - x3 * y4; 01241 01242 /* Compute r1 and r2 */ 01243 r1 = a2 * x1 + b2 * y1 + c2; 01244 r2 = a2 * x2 + b2 * y2 + c2; 01245 01246 /* Check signs of r1 and r2. If both point 1 and point 2 lie 01247 * on same side of second line segment, the line segments do 01248 * not intersect. 01249 */ 01250 01251 if (!EQUAL(r1, 0.0) && !EQUAL(r2, 0.0) && SAME_SIGNS(r1, r2)) { 01252 return (DONT_INTERSECT); 01253 } 01254 01255 /* Line segments intersect: compute intersection point. 01256 */ 01257 denom = a1 * b2 - a2 * b1; 01258 01259 if (denom == 0) { 01260 return (COLLINEAR); 01261 } 01262 01263 offset = denom < 0 ? -denom / 2 : denom / 2; 01264 01265 /* The denom/2 is to get rounding instead of truncating. It 01266 * is added or subtracted to the numerator, depending upon the 01267 * sign of the numerator. 01268 */ 01269 num = b1 * c2 - b2 * c1; 01270 01271 *x = num / denom; 01272 01273 num = a2 * c1 - a1 * c2; 01274 *y = num / denom; 01275 01276 return (DO_INTERSECT); 01277 } 01278 01290 int Point_on_plane(Point3 p1, Point3 p2, Point3 p3, Point3 unk) 01291 { 01292 float plane[4]; 01293 01294 P3toPlane(p1, p2, p3, plane); 01295 01296 return (XY_intersect_plane(unk, plane)); 01297 } 01298 01312 int XY_intersect_plane(float *intersect, float *plane) 01313 { 01314 float x, y; 01315 01316 if (!plane[Z]) { 01317 return (0); /* doesn't intersect */ 01318 } 01319 01320 x = intersect[X]; 01321 y = intersect[Y]; 01322 intersect[Z] = (plane[X] * x + plane[Y] * y + plane[W]) / -plane[Z]; 01323 01324 return (1); 01325 } 01326 01335 int P3toPlane(Point3 p1, Point3 p2, Point3 p3, float *plane) 01336 { 01337 Point3 v1, v2, norm; 01338 01339 v1[X] = p1[X] - p3[X]; 01340 v1[Y] = p1[Y] - p3[Y]; 01341 v1[Z] = p1[Z] - p3[Z]; 01342 01343 v2[X] = p2[X] - p3[X]; 01344 v2[Y] = p2[Y] - p3[Y]; 01345 v2[Z] = p2[Z] - p3[Z]; 01346 01347 V3Cross(v1, v2, norm); 01348 01349 plane[X] = norm[X]; 01350 plane[Y] = norm[Y]; 01351 plane[Z] = norm[Z]; 01352 plane[W] = -p3[X] * norm[X] - p3[Y] * norm[Y] - p3[Z] * norm[Z]; 01353 01354 return (1); 01355 } 01356 01357 01365 int V3Cross(Point3 a, Point3 b, Point3 c) 01366 { 01367 c[X] = (a[Y] * b[Z]) - (a[Z] * b[Y]); 01368 c[Y] = (a[Z] * b[X]) - (a[X] * b[Z]); 01369 c[Z] = (a[X] * b[Y]) - (a[Y] * b[X]); 01370 01371 return (1); 01372 }