1.1 --- a/latlon-v0005.c	Fri Sep 09 23:44:28 2016 +0200
     1.2 +++ b/latlon-v0005.c	Mon Sep 12 18:00:52 2016 +0200
     1.3 @@ -489,7 +489,13 @@
     1.4  }
     1.5  
     1.6  /* check if point is inside cluster */
     1.7 -static bool pgl_point_in_cluster(pgl_point *point, pgl_cluster *cluster) {
     1.8 +/* (if point is on perimeter, then true is returned if and only if
     1.9 +   strict == false) */
    1.10 +static bool pgl_point_in_cluster(
    1.11 +  pgl_point *point,
    1.12 +  pgl_cluster *cluster,
    1.13 +  bool strict
    1.14 +) {
    1.15    int i, j, k;  /* i: entry, j: point in entry, k: next point in entry */
    1.16    int entrytype;         /* type of entry */
    1.17    int npoints;           /* number of points in entry */
    1.18 @@ -504,8 +510,11 @@
    1.19    int counter = 0;       /* counter for intersections east of point */
    1.20    /* iterate over all entries */
    1.21    for (i=0; i<cluster->nentries; i++) {
    1.22 -    /* get properties of entry */
    1.23 +    /* get type of entry */
    1.24      entrytype = cluster->entries[i].entrytype;
    1.25 +    /* skip all entries but polygons if perimeters are excluded */
    1.26 +    if (strict && entrytype != PGL_ENTRY_POLYGON) continue;
    1.27 +    /* get points of entry */
    1.28      npoints = cluster->entries[i].npoints;
    1.29      points = PGL_ENTRY_POINTS(cluster, i);
    1.30      /* determine east/west orientation of first point of entry and calculate
    1.31 @@ -521,8 +530,8 @@
    1.32      else if (lon_dir > 0 && lon0 < lon_break) lon0 = pgl_round(lon0 + 360);
    1.33      /* iterate over all edges and vertices */
    1.34      for (j=0; j<npoints; j++) {
    1.35 -      /* return true if point is on vertex of polygon */
    1.36 -      if (pgl_point_cmp(point, &(points[j])) == 0) return true;
    1.37 +      /* return if point is on vertex of polygon */
    1.38 +      if (pgl_point_cmp(point, &(points[j])) == 0) return !strict;
    1.39        /* calculate index of next vertex */
    1.40        k = (j+1) % npoints;
    1.41        /* skip last edge unless entry is (closed) outline or polygon */
    1.42 @@ -549,17 +558,17 @@
    1.43        /* consider longitude wrap-around for next vertex */
    1.44        if      (lon_dir < 0 && lon2 > lon_break) lon2 = pgl_round(lon2 - 360);
    1.45        else if (lon_dir > 0 && lon2 < lon_break) lon2 = pgl_round(lon2 + 360);
    1.46 -      /* return true if point is on horizontal (west to east) edge of polygon */
    1.47 +      /* return if point is on horizontal (west to east) edge of polygon */
    1.48        if (
    1.49          lat0 == lat1 && lat0 == lat2 &&
    1.50          ( (lon0 >= lon1 && lon0 <= lon2) || (lon0 >= lon2 && lon0 <= lon1) )
    1.51 -      ) return true;
    1.52 +      ) return !strict;
    1.53        /* check if edge crosses east/west line of point */
    1.54        if ((lat1 < lat0 && lat2 >= lat0) || (lat2 < lat0 && lat1 >= lat0)) {
    1.55          /* calculate longitude of intersection */
    1.56          lon = (lon1 * (lat2-lat0) + lon2 * (lat0-lat1)) / (lat2-lat1);
    1.57 -        /* return true if intersection goes (approximately) through point */
    1.58 -        if (pgl_round(lon) == lon0) return true;
    1.59 +        /* return if intersection goes (approximately) through point */
    1.60 +        if (pgl_round(lon) == lon0) return !strict;
    1.61          /* count intersection if east of point and entry is polygon*/
    1.62          if (entrytype == PGL_ENTRY_POLYGON && lon > lon0) counter++;
    1.63        }
    1.64 @@ -569,8 +578,9 @@
    1.65    return counter & 1;
    1.66  }
    1.67  
    1.68 -/* check if all points of the second cluster are inside the first cluster */
    1.69 -static inline bool pgl_all_cluster_points_in_cluster(
    1.70 +/* check if all points of the second cluster are strictly inside the first
    1.71 +   cluster */
    1.72 +static inline bool pgl_all_cluster_points_strictly_in_cluster(
    1.73    pgl_cluster *outer, pgl_cluster *inner
    1.74  ) {
    1.75    int i, j;           /* i: entry, j: point in entry */
    1.76 @@ -584,7 +594,7 @@
    1.77      /* iterate over all points in entry of "inner" cluster */
    1.78      for (j=0; j<npoints; j++) {
    1.79        /* return false if one point of inner cluster is not in outer cluster */
    1.80 -      if (!pgl_point_in_cluster(points+j, outer)) return false;
    1.81 +      if (!pgl_point_in_cluster(points+j, outer, true)) return false;
    1.82      }
    1.83    }
    1.84    /* otherwise return true */
    1.85 @@ -606,36 +616,33 @@
    1.86      /* iterate over all points in entry of "inner" cluster */
    1.87      for (j=0; j<npoints; j++) {
    1.88        /* return true if one point of inner cluster is in outer cluster */
    1.89 -      if (pgl_point_in_cluster(points+j, outer)) return true;
    1.90 +      if (pgl_point_in_cluster(points+j, outer, false)) return true;
    1.91      }
    1.92    }
    1.93    /* otherwise return false */
    1.94    return false;
    1.95  }
    1.96  
    1.97 -/* check if line segment crosses line */
    1.98 -/* returns -1 if yes, 1 if no, and 0 in corner cases */
    1.99 -/* NOTE: each line (segment) must have a length greater than zero */
   1.100 -static inline double pgl_lseg_crosses_line(
   1.101 +/* check if line segment strictly crosses line (not just touching) */
   1.102 +static inline bool pgl_lseg_crosses_line(
   1.103    double seg_x1,  double seg_y1,  double seg_x2,  double seg_y2,
   1.104 -  double line_x1, double line_y1, double line_x2, double line_y2,
   1.105 -  bool strict
   1.106 +  double line_x1, double line_y1, double line_x2, double line_y2
   1.107  ) {
   1.108 -  double value = (
   1.109 -    (seg_x1-line_x1) * (line_y2-line_y1) -
   1.110 -    (seg_y1-line_y1) * (line_x2-line_x1)
   1.111 -  ) * (
   1.112 -    (seg_x2-line_x1) * (line_y2-line_y1) -
   1.113 -    (seg_y2-line_y1) * (line_x2-line_x1)
   1.114 -  );
   1.115 -  if (strict) return value < 0;
   1.116 -  else return value <= 0;
   1.117 +  return (
   1.118 +    (
   1.119 +      (seg_x1-line_x1) * (line_y2-line_y1) -
   1.120 +      (seg_y1-line_y1) * (line_x2-line_x1)
   1.121 +    ) * (
   1.122 +      (seg_x2-line_x1) * (line_y2-line_y1) -
   1.123 +      (seg_y2-line_y1) * (line_x2-line_x1)
   1.124 +    )
   1.125 +  ) < 0;
   1.126  }
   1.127  
   1.128 -/* check if paths and outlines of two clusters overlap */
   1.129 -/* (set strict to true to disregard corner cases) */
   1.130 +/* check if paths and outlines of two clusters strictly overlap (not just
   1.131 +   touching) */
   1.132  static bool pgl_outlines_overlap(
   1.133 -  pgl_cluster *cluster1, pgl_cluster *cluster2, bool strict
   1.134 +  pgl_cluster *cluster1, pgl_cluster *cluster2
   1.135  ) {
   1.136    int i1, j1, k1;  /* i: entry, j: point in entry, k: next point in entry */
   1.137    int i2, j2, k2;
   1.138 @@ -758,12 +765,10 @@
   1.139            if (
   1.140              pgl_lseg_crosses_line(
   1.141                lat11, lon11, lat12, lon12,
   1.142 -              lat21, lon21, lat22, lon22,
   1.143 -              strict
   1.144 +              lat21, lon21, lat22, lon22
   1.145              ) && pgl_lseg_crosses_line(
   1.146                lat21, lon21, lat22, lon22,
   1.147 -              lat11, lon11, lat12, lon12,
   1.148 -              strict
   1.149 +              lat11, lon11, lat12, lon12
   1.150              )
   1.151            ) {
   1.152              return true;
   1.153 @@ -778,8 +783,9 @@
   1.154  
   1.155  /* check if second cluster is completely contained in first cluster */
   1.156  static bool pgl_cluster_in_cluster(pgl_cluster *outer, pgl_cluster *inner) {
   1.157 -  if (!pgl_all_cluster_points_in_cluster(outer, inner)) return false;
   1.158 -  if (pgl_outlines_overlap(outer, inner, true)) return false;
   1.159 +  if (!pgl_all_cluster_points_strictly_in_cluster(outer, inner)) return false;
   1.160 +  if (pgl_any_cluster_points_in_cluster(inner, outer)) return false;
   1.161 +  if (pgl_outlines_overlap(outer, inner)) return false;
   1.162    return true;
   1.163  }
   1.164  
   1.165 @@ -789,7 +795,7 @@
   1.166  ) {
   1.167    if (pgl_any_cluster_points_in_cluster(cluster1, cluster2)) return true;
   1.168    if (pgl_any_cluster_points_in_cluster(cluster2, cluster1)) return true;
   1.169 -  if (pgl_outlines_overlap(cluster1, cluster2, false)) return true;
   1.170 +  if (pgl_outlines_overlap(cluster1, cluster2)) return true;
   1.171    return false;
   1.172  }
   1.173  
   1.174 @@ -812,7 +818,7 @@
   1.175    double dist;           /* distance calculated in one step */
   1.176    double min_dist = INFINITY;   /* minimum distance */
   1.177    /* distance is zero if point is contained in cluster */
   1.178 -  if (pgl_point_in_cluster(point, cluster)) return 0;
   1.179 +  if (pgl_point_in_cluster(point, cluster, false)) return 0;
   1.180    /* iterate over all entries */
   1.181    for (i=0; i<cluster->nentries; i++) {
   1.182      /* get properties of entry */
   1.183 @@ -2343,7 +2349,7 @@
   1.184        cluster->bounding.center.lat, cluster->bounding.center.lon
   1.185      ) > cluster->bounding.radius
   1.186    ) retval = false;
   1.187 -  else retval = pgl_point_in_cluster(point, cluster);
   1.188 +  else retval = pgl_point_in_cluster(point, cluster, false);
   1.189    PG_FREE_IF_COPY(cluster, 1);
   1.190    PG_RETURN_BOOL(retval);
   1.191  }