GRASS Programmer's Manual  6.4.2(2012)
gsdrape.c
Go to the documentation of this file.
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 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines