Vlib/poly.c

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

Generated on Sun Apr 6 17:32:44 2008 for GRASS by  doxygen 1.5.5