Vlib/poly.c

Go to the documentation of this file.
00001 
00021 #include <math.h>
00022 #include <stdlib.h>
00023 #include <grass/Vect.h>
00024 #include <grass/gis.h>
00025 #include <grass/linkm.h>
00026 #include <grass/glocale.h>
00027 
00028 struct Slink
00029 {
00030     double x;
00031     struct Slink *next;
00032 };
00033 
00034 
00035 /* function prototypes */
00036 static int comp_double(double *, double *);
00037 static int V__within(double, double, double);
00038 int Vect__intersect_line_with_poly();
00039 static void destroy_links(struct Slink *);
00040 static int Vect__divide_and_conquer(struct Slink *, struct line_pnts *,
00041                                     struct link_head *, double *, double *,
00042                                     int);
00043 
00044 
00060 int
00061 Vect_get_point_in_area(struct Map_info *Map, int area, double *X, double *Y)
00062 {
00063     static struct line_pnts *Points;
00064     static struct line_pnts **IPoints;
00065     static int first_time = 1;
00066     static int isl_allocated = 0;
00067     int i, n_isles;
00068 
00069     G_debug(3, "Vect_get_point_in_area()");
00070 
00071     if (first_time) {
00072         Points = Vect_new_line_struct();
00073         IPoints = NULL;
00074         first_time = 0;
00075     }
00076     n_isles = Vect_get_area_num_isles(Map, area);
00077     if (n_isles > isl_allocated) {
00078         IPoints = (struct line_pnts **)
00079             G_realloc(IPoints, (1 + n_isles) * sizeof(struct line_pnts *));
00080         for (i = isl_allocated; i < n_isles; i++)
00081             IPoints[i] = Vect_new_line_struct();
00082         isl_allocated = n_isles;
00083     }
00084 
00085     if (0 > Vect_get_area_points(Map, area, Points))
00086         return -1;
00087 
00088     for (i = 0; i < n_isles; i++) {
00089         IPoints[i]->alloc_points = 0;
00090         if (0 >
00091             Vect_get_isle_points(Map, Vect_get_area_isle(Map, area, i),
00092                                  IPoints[i]))
00093             return -1;
00094     }
00095     return (Vect_get_point_in_poly_isl(Points, IPoints, n_isles, X, Y));
00096 
00097     return -1;
00098 }
00099 
00100 static int comp_double(double *i, double *j)
00101 {
00102     if (*i < *j)
00103         return -1;
00104 
00105     if (*i > *j)
00106         return 1;
00107 
00108     return 0;
00109 }
00110 
00111 static int V__within(double a, double x, double b)
00112 {
00113     double tmp;
00114 
00115     if (a > b) {
00116         tmp = a;
00117         a = b;
00118         b = tmp;
00119     }
00120 
00121     return (x >= a && x <= b);
00122 }
00123 
00124 /*
00125    \brief Intersects line with polygon
00126 
00127    For each intersection of a polygon w/ a line, stuff the X value in
00128    the Inter Points array.  I used line_pnts, just cuz the memory
00129    management was already there.  I am getting real tired of managing
00130    realloc stuff.  Assumes that no vertex of polygon lies on Y This is
00131    taken care of by functions calling this function
00132 
00133    \param Points line
00134    \param y ?
00135    \param Iter intersection ?
00136 
00137    \return 0 on success
00138    \return -1 on error 
00139  */
00140 int
00141 Vect__intersect_line_with_poly(struct line_pnts *Points,
00142                                double y, struct line_pnts *Inter)
00143 {
00144     int i;
00145     double a, b, c, d, x;
00146     double perc;
00147 
00148     for (i = 1; i < Points->n_points; i++) {
00149         a = Points->y[i - 1];
00150         b = Points->y[i];
00151 
00152         c = Points->x[i - 1];
00153         d = Points->x[i];
00154 
00155         if (V__within(a, y, b)) {
00156             if (a == b)
00157                 continue;
00158 
00159             perc = (y - a) / (b - a);
00160             x = perc * (d - c) + c;     /* interp X */
00161 
00162             if (0 > Vect_append_point(Inter, x, y, 0))
00163                 return -1;
00164         }
00165     }
00166     return 0;
00167 }
00168 
00180 int Vect_get_point_in_poly(struct line_pnts *Points, double *X, double *Y)
00181 {
00182     double cent_x, cent_y;
00183     struct Slink *Head;
00184     static struct link_head *Token;
00185     struct Slink *tmp;
00186     static int first_time = 1;
00187     register int i;
00188     double x_max, x_min;
00189     int ret;
00190 
00191     /* get centroid */
00192     Vect_find_poly_centroid(Points, &cent_x, &cent_y);
00193     /* is it w/in poly? */
00194     if (Vect_point_in_poly(cent_x, cent_y, Points) == 1) {
00195         *X = cent_x;
00196         *Y = cent_y;
00197         return 0;
00198     }
00199 
00200     /* guess we have to do it the hard way... */
00201     /* get min and max x values */
00202     x_max = x_min = Points->x[0];
00203     for (i = 0; i < Points->n_points; i++) {
00204         if (x_min > Points->x[i])
00205             x_min = Points->x[i];
00206         if (x_max < Points->x[i])
00207             x_max = Points->x[i];
00208     }
00209 
00210 
00211     /* init the linked list */
00212     if (first_time) {
00213         /* will never call link_cleanup ()  */
00214         link_exit_on_error(1);  /* kill program if out of memory */
00215         Token = (struct link_head *)link_init(sizeof(struct Slink));
00216         first_time = 0;
00217     }
00218 
00219     Head = (struct Slink *)link_new(Token);
00220     tmp = (struct Slink *)link_new(Token);
00221 
00222     Head->next = tmp;
00223     tmp->next = NULL;
00224 
00225     Head->x = x_min;
00226     tmp->x = x_max;
00227 
00228     *Y = cent_y;                /* pick line segment (x_min, cent_y) - (x_max, cent_y) */
00229     ret = Vect__divide_and_conquer(Head, Points, Token, X, Y, 10);
00230 
00231     destroy_links(Head);
00232 
00233     if (ret < 0) {
00234         G_warning("Vect_get_point_in_poly(): %s",
00235                   _("Unable to find point in polygon"));
00236         return -1;
00237     }
00238 
00239     G_debug(3, "Found point in %d iterations", 10 - ret);
00240 
00241     return 0;
00242 }
00243 
00244 
00245 /*
00246    \brief Provide a breadth first binary division of real space along line segment.
00247 
00248    Looking for a point w/in the polygon.
00249 
00250    This routine walks along the list of points on line segment
00251    and divides each pair in half. It sticks that new point right into
00252    the list, and then checks to see if it is inside the poly. 
00253 
00254    After going through the whole list, it calls itself.  The list 
00255    now has a whole extra set of points to divide again.
00256 
00257    \param Head 
00258    \param Points
00259    \param Token
00260    \param X,Y
00261    \param levels
00262 
00263    \return # levels it took
00264    \return -1 if exceeded # of levels
00265  */
00266 static int
00267 Vect__divide_and_conquer(struct Slink *Head,
00268                          struct line_pnts *Points,
00269                          struct link_head *Token,
00270                          double *X, double *Y, int levels)
00271 {
00272     struct Slink *A, *B, *C;
00273 
00274     G_debug(3, "Vect__divide_and_conquer(): LEVEL %d", levels);
00275     A = Head;
00276     B = Head->next;
00277 
00278     do {
00279         C = (struct Slink *)link_new(Token);
00280         A->next = C;
00281         C->next = B;
00282 
00283         C->x = (A->x + B->x) / 2.;
00284 
00285         if (Vect_point_in_poly(C->x, *Y, Points) == 1) {
00286             *X = C->x;
00287             return levels;
00288         }
00289 
00290         A = B;
00291         B = B->next;
00292     }
00293     while (B != NULL);
00294 
00295     /*
00296      **  If it got through the entire loop and still no hits,
00297      **   then lets go a level deeper and divide again.
00298      */
00299 
00300     if (levels <= 0)
00301         return -1;
00302 
00303     return Vect__divide_and_conquer(Head, Points, Token, X, Y, --levels);
00304 }
00305 
00306 static void destroy_links(struct Slink *Head)
00307 {
00308     struct Slink *p, *tmp;
00309 
00310     p = Head;
00311 
00312     while (p != NULL) {
00313         tmp = p->next;
00314         link_dispose((struct link_head *)Head, (VOID_T *) p);
00315         p = tmp;
00316     }
00317 }
00318 
00319 
00329 int
00330 Vect_find_poly_centroid(struct line_pnts *points,
00331                         double *cent_x, double *cent_y)
00332 {
00333     int i;
00334     double *xptr1, *yptr1;
00335     double *xptr2, *yptr2;
00336     double cent_weight_x, cent_weight_y;
00337     double len, tot_len;
00338 
00339     tot_len = 0.0;
00340     cent_weight_x = 0.0;
00341     cent_weight_y = 0.0;
00342 
00343     xptr1 = points->x;
00344     yptr1 = points->y;
00345     xptr2 = points->x + 1;
00346     yptr2 = points->y + 1;
00347 
00348     for (i = 1; i < points->n_points; i++) {
00349         len = hypot(*xptr1 - *xptr2, *yptr1 - *yptr2);
00350         cent_weight_x += len * ((*xptr1 + *xptr2) / 2.);
00351         cent_weight_y += len * ((*yptr1 + *yptr2) / 2.);
00352         tot_len += len;
00353         xptr1++;
00354         xptr2++;
00355         yptr1++;
00356         yptr2++;
00357     }
00358 
00359     if (tot_len == 0.0)
00360         return -1;
00361 
00362     *cent_x = cent_weight_x / tot_len;
00363     *cent_y = cent_weight_y / tot_len;
00364 
00365     return 0;
00366 }
00367 
00368 
00369 /*
00370  ** returns true if point is in any of islands /w in area
00371  ** returns 0 if not
00372  ** returns -1 on error
00373  */
00374 /*
00375    int 
00376    Vect_point_in_islands (
00377    struct Map_info *Map,
00378    int area,
00379    double cent_x, double cent_y)
00380    {
00381    P_AREA *Area;
00382    static struct line_pnts *TPoints;
00383    static int first_time = 1;
00384    int isle;
00385 
00386    if (first_time == 1)
00387    {
00388    TPoints = Vect_new_line_struct ();
00389    first_time = 0;
00390    }
00391 
00392    Area = &(Map->plus.Area[area]);
00393 
00394    for (isle = 0; isle < Area->n_isles; isle++)
00395    {
00396    if (0 > Vect_get_isle_points (Map, Area->isles[isle], TPoints))
00397    return -1;
00398 
00399    if ( Vect_point_in_poly (cent_x, cent_y, TPoints) == 1 )
00400    return 1;
00401    }
00402 
00403    return 0;
00404    }
00405  */
00406 
00423 int
00424 Vect_get_point_in_poly_isl(struct line_pnts *Points,
00425                            struct line_pnts **IPoints, int n_isles,
00426                            double *att_x, double *att_y)
00427 {
00428     static struct line_pnts *Intersects;
00429     static int first_time = 1;
00430     double cent_x, cent_y;
00431     register int i, j;
00432     double max, hi_y, lo_y;
00433     int maxpos;
00434     int point_in_sles = 0;
00435     double diff;
00436 
00437     G_debug(3, "Vect_get_point_in_poly_isl(): n_isles = %d", n_isles);
00438 
00439     if (first_time) {
00440         Intersects = Vect_new_line_struct();
00441         first_time = 0;
00442     }
00443 
00444     if (Points->n_points < 3) { /* test */
00445         if (Points->n_points > 0) {
00446             *att_x = Points->x[0];
00447             *att_y = Points->y[0];
00448             return 0;
00449         }
00450         return -1;
00451     }
00452 
00453     /* get centroid */
00454     Vect_find_poly_centroid(Points, &cent_x, &cent_y);
00455     /* is it w/in poly? */
00456     if (Vect_point_in_poly(cent_x, cent_y, Points) == 1)
00457         /* if the point is iside the polygon */
00458     {
00459         for (i = 0; i < n_isles; i++) {
00460             if (Vect_point_in_poly(cent_x, cent_y, IPoints[i]) >= 1) {
00461                 point_in_sles = 1;
00462                 break;
00463             }
00464         }
00465         if (!point_in_sles) {
00466             *att_x = cent_x;
00467             *att_y = cent_y;
00468             return 0;
00469         }
00470     }
00471     /* guess we have to do it the hard way... */
00472 
00473     /* first find att_y close to cent_y so that no points lie on the line */
00474     /* find the point closest to line from below, and point close to line
00475        from above and take average of their y-coordinates */
00476 
00477     /* first initializing lo_y,hi_y to be any 2 pnts on either side of cent_y */
00478     hi_y = cent_y - 1;
00479     lo_y = cent_y + 1;
00480     for (i = 0; i < Points->n_points; i++) {
00481         if ((lo_y < cent_y) && (hi_y >= cent_y))
00482             break;              /* already initialized */
00483         if (Points->y[i] < cent_y)
00484             lo_y = Points->y[i];
00485         if (Points->y[i] >= cent_y)
00486             hi_y = Points->y[i];
00487     }
00488     /* first going throught boundary points */
00489     for (i = 0; i < Points->n_points; i++) {
00490         if ((Points->y[i] < cent_y) &&
00491             ((cent_y - Points->y[i]) < (cent_y - lo_y)))
00492             lo_y = Points->y[i];
00493         if ((Points->y[i] >= cent_y) &&
00494             ((Points->y[i] - cent_y) < (hi_y - cent_y)))
00495             hi_y = Points->y[i];
00496     }
00497     for (i = 0; i < n_isles; i++)
00498         for (j = 0; j < IPoints[i]->n_points; j++) {
00499             if ((IPoints[i]->y[j] < cent_y) &&
00500                 ((cent_y - IPoints[i]->y[j]) < (cent_y - lo_y)))
00501                 lo_y = IPoints[i]->y[j];
00502 
00503             if ((IPoints[i]->y[j] >= cent_y) &&
00504                 ((IPoints[i]->y[j] - cent_y) < (hi_y - cent_y)))
00505                 hi_y = IPoints[i]->y[j];
00506         }
00507 
00508     if (lo_y == hi_y)
00509         return (-1);            /* area is empty */
00510     else
00511         *att_y = (hi_y + lo_y) / 2.0;
00512 
00513     Intersects->n_points = 0;
00514     Vect__intersect_line_with_poly(Points, *att_y, Intersects);
00515 
00516     /* add in intersections w/ holes */
00517     for (i = 0; i < n_isles; i++) {
00518         if (0 >
00519             Vect__intersect_line_with_poly(IPoints[i], *att_y, Intersects))
00520             return -1;
00521     }
00522 
00523     if (Intersects->n_points < 2)       /* test */
00524         return -1;
00525 
00526     qsort(Intersects->x, (size_t) Intersects->n_points, sizeof(double),
00527           (void *)comp_double);
00528 
00529     max = 0;
00530     maxpos = 0;
00531 
00532     /* find area of MAX distance */
00533     for (i = 0; i < Intersects->n_points; i += 2) {
00534         diff = Intersects->x[i + 1] - Intersects->x[i];
00535 
00536         if (diff > max) {
00537             max = diff;
00538             maxpos = i;
00539         }
00540     }
00541     if (max == 0.0)             /* area was empty: example ((x1,y1), (x2,y2), (x1,y1)) */
00542         return -1;
00543 
00544     *att_x = (Intersects->x[maxpos] + Intersects->x[maxpos + 1]) / 2.;
00545 
00546     return 0;
00547 }
00548 
00549 
00550 /* Intersect segments of Points with ray from point X,Y to the right.
00551  * Returns: -1 point exactly on segment
00552  *          number of intersections
00553  */
00554 static int segments_x_ray(double X, double Y, struct line_pnts *Points)
00555 {
00556     double x1, x2, y1, y2;
00557     double x_inter;
00558     int n_intersects;
00559     int n;
00560 
00561     G_debug(3, "segments_x_ray(): x = %f y = %f n_points = %d", X, Y,
00562             Points->n_points);
00563 
00564     /* Follow the ray from X,Y along positive x and find number of intersections.
00565      * Coordinates exactly on ray are considered to be slightly above. */
00566 
00567     n_intersects = 0;
00568     for (n = 0; n < Points->n_points - 1; n++) {
00569         x1 = Points->x[n];
00570         y1 = Points->y[n];
00571         x2 = Points->x[n + 1];
00572         y2 = Points->y[n + 1];
00573 
00574         G_debug(3, "X = %f Y = %f x1 = %f y1 = %f x2 = %f y2 = %f", X, Y, x1,
00575                 y1, x2, y2);
00576 
00577         /* I know, it should be possible to do that with less conditions, but it should be 
00578          * enough readable also! */
00579 
00580         /* segment left from X -> no intersection */
00581         if (x1 < X && x2 < X)
00582             continue;
00583 
00584         /* point on vertex */
00585         if ((x1 == X && y1 == Y) || (x2 == X && y2 == Y))
00586             return -1;
00587 
00588         /* on vertical boundary */
00589         if ((x1 == x2 && x1 == X) &&
00590             ((y1 <= Y && y2 >= Y) || (y1 >= Y && y2 <= Y)))
00591             return -1;
00592 
00593         /* on horizontal boundary */
00594         if ((y1 == y2 && y1 == Y) &&
00595             ((x1 <= X && x2 >= X) || (x1 >= X && x2 <= X)))
00596             return -1;
00597 
00598         /* segment on ray (X is not important) */
00599         if (y1 == Y && y2 == Y)
00600             continue;
00601 
00602         /* segment above (X is not important) */
00603         if (y1 > Y && y2 > Y)
00604             continue;
00605 
00606         /* segment below (X is not important) */
00607         if (y1 < Y && y2 < Y)
00608             continue;
00609 
00610         /* one end on Y second above (X is not important) */
00611         if ((y1 == Y && y2 > Y) || (y2 == Y && y1 > Y))
00612             continue;
00613 
00614         /* For following cases we know that at least one of x1 and x2 is  >= X */
00615 
00616         /* one end of segment on Y second below Y */
00617         if (y1 == Y && y2 < Y) {
00618             if (x1 >= X)        /* x of the end on the ray is >= X */
00619                 n_intersects++;
00620             continue;
00621         }
00622         if (y2 == Y && y1 < Y) {
00623             if (x2 >= X)
00624                 n_intersects++;
00625             continue;
00626         }
00627 
00628         /* one end of segment above Y second below Y */
00629         if ((y1 < Y && y2 > Y) || (y1 > Y && y2 < Y)) {
00630             if (x1 >= X && x2 >= X) {
00631                 n_intersects++;
00632                 continue;
00633             }
00634 
00635             /* now either x1 < X && x2 > X or x1 > X && x2 < X -> calculate intersection */
00636             x_inter = dig_x_intersect(x1, x2, y1, y2, Y);
00637             G_debug(3, "x_inter = %f", x_inter);
00638             if (x_inter == X)
00639                 return 1;
00640             else if (x_inter > X)
00641                 n_intersects++;
00642 
00643             continue;           /* would not be necessary, just to check, see below */
00644         }
00645         /* should not be reached (one condition is not necessary, but it is may be better readable
00646          * and it is a check) */
00647         G_warning
00648             ("segments_x_ray() %s: X = %f Y = %f x1 = %f y1 = %f x2 = %f y2 = %f",
00649              _("conditions failed"), X, Y, x1, y1, x2, y2);
00650     }
00651 
00652     return n_intersects;
00653 }
00654 
00665 int Vect_point_in_poly(double X, double Y, struct line_pnts *Points)
00666 {
00667     int n_intersects;
00668 
00669     G_debug(3, "Vect_point_in_poly(): x = %f y = %f n_points = %d", X, Y,
00670             Points->n_points);
00671 
00672     n_intersects = segments_x_ray(X, Y, Points);
00673 
00674     if (n_intersects == -1)
00675         return 2;
00676 
00677     if (n_intersects % 2)
00678         return 1;
00679     else
00680         return 0;
00681 }
00682 
00694 int
00695 Vect_point_in_area_outer_ring(double X, double Y, struct Map_info *Map,
00696                               int area)
00697 {
00698     static int first = 1;
00699     int n_intersects, inter;
00700     int i, line;
00701     static struct line_pnts *Points;
00702     struct Plus_head *Plus;
00703     P_LINE *Line;
00704     P_AREA *Area;
00705 
00706     G_debug(3, "Vect_point_in_area_outer_ring(): x = %f y = %f area = %d", X,
00707             Y, area);
00708 
00709     if (first == 1) {
00710         Points = Vect_new_line_struct();
00711         first = 0;
00712     }
00713 
00714     Plus = &(Map->plus);
00715     Area = Plus->Area[area];
00716 
00717     /* First it must be in box */
00718     if (X < Area->W || X > Area->E || Y > Area->N || Y < Area->S)
00719         return 0;
00720 
00721     n_intersects = 0;
00722     for (i = 0; i < Area->n_lines; i++) {
00723         line = abs(Area->lines[i]);
00724         G_debug(3, "  line[%d] = %d", i, line);
00725 
00726         Line = Plus->Line[line];
00727 
00728         /* dont check lines that obviously do not intersect with test ray */
00729         if ((Line->N < Y) || (Line->S > Y) || (Line->E < X))
00730             continue;
00731 
00732         Vect_read_line(Map, Points, NULL, line);
00733 
00734         inter = segments_x_ray(X, Y, Points);
00735         G_debug(3, "  inter = %d", inter);
00736 
00737         if (inter == -1)
00738             return 2;
00739         n_intersects += inter;
00740         G_debug(3, "  n_intersects = %d", n_intersects);
00741     }
00742 
00743     if (n_intersects % 2)
00744         return 1;
00745     else
00746         return 0;
00747 }
00748 
00760 int Vect_point_in_island(double X, double Y, struct Map_info *Map, int isle)
00761 {
00762     static int first = 1;
00763     int n_intersects, inter;
00764     int i, line;
00765     static struct line_pnts *Points;
00766     struct Plus_head *Plus;
00767     P_LINE *Line;
00768     P_ISLE *Isle;
00769 
00770     G_debug(3, "Vect_point_in_island(): x = %f y = %f isle = %d", X, Y, isle);
00771 
00772     if (first == 1) {
00773         Points = Vect_new_line_struct();
00774         first = 0;
00775     }
00776 
00777     Plus = &(Map->plus);
00778     Isle = Plus->Isle[isle];
00779 
00780     if (X < Isle->W || X > Isle->E || Y > Isle->N || Y < Isle->S)
00781         return 0;
00782 
00783     n_intersects = 0;
00784     for (i = 0; i < Isle->n_lines; i++) {
00785         line = abs(Isle->lines[i]);
00786 
00787         Line = Plus->Line[line];
00788 
00789         /* dont check lines that obviously do not intersect with test ray */
00790         if ((Line->N < Y) || (Line->S > Y) || (Line->E < X))
00791             continue;
00792 
00793         Vect_read_line(Map, Points, NULL, line);
00794 
00795         inter = segments_x_ray(X, Y, Points);
00796         if (inter == -1)
00797             return 2;
00798         n_intersects += inter;
00799     }
00800 
00801     if (n_intersects % 2)
00802         return 1;
00803     else
00804         return 0;
00805 }

Generated on Sat Oct 24 03:25:20 2009 for GRASS Programmer's Manual by  doxygen 1.6.1