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

Generated on Mon Jan 1 19:49:16 2007 for GRASS by  doxygen 1.5.1