1.1 --- a/latlon-v0004.c	Sat Sep 03 16:00:22 2016 +0200
     1.2 +++ b/latlon-v0004.c	Fri Sep 09 19:22:30 2016 +0200
     1.3 @@ -502,14 +502,6 @@
     1.4    double lat2, lon2;     /* latitude and (adjusted) longitude of next vertex */
     1.5    double lon;            /* longitude of intersection */
     1.6    int counter = 0;       /* counter for intersections east of point */
     1.7 -  /* points outside bounding circle are always assumed to be non-overlapping */
     1.8 -  /* (necessary for consistent table and index scans) */
     1.9 -  if (
    1.10 -    pgl_distance(
    1.11 -      point->lat, point->lon,
    1.12 -      cluster->bounding.center.lat, cluster->bounding.center.lon
    1.13 -    ) > cluster->bounding.radius
    1.14 -  ) return false;
    1.15    /* iterate over all entries */
    1.16    for (i=0; i<cluster->nentries; i++) {
    1.17      /* get properties of entry */
    1.18 @@ -539,14 +531,22 @@
    1.19          entrytype != PGL_ENTRY_OUTLINE &&
    1.20          entrytype != PGL_ENTRY_POLYGON
    1.21        ) continue;
    1.22 -      /* get latitude and longitude values of edge */
    1.23 -      lat1 = points[j].lat;
    1.24 +      /* use previously calculated values for lat1 and lon1 if possible */
    1.25 +      if (j) {
    1.26 +        lat1 = lat2;
    1.27 +        lon1 = lon2;
    1.28 +      } else {
    1.29 +        /* otherwise get latitude and longitude values of first vertex */
    1.30 +        lat1 = points[0].lat;
    1.31 +        lon1 = points[0].lon;
    1.32 +        /* and consider longitude wrap-around for first vertex */
    1.33 +        if      (lon_dir < 0 && lon1 > lon_break) lon1 = pgl_round(lon1 - 360);
    1.34 +        else if (lon_dir > 0 && lon1 < lon_break) lon1 = pgl_round(lon1 + 360);
    1.35 +      }
    1.36 +      /* get latitude and longitude of next vertex */
    1.37        lat2 = points[k].lat;
    1.38 -      lon1 = points[j].lon;
    1.39        lon2 = points[k].lon;
    1.40 -      /* consider longitude wrap-around for edge */
    1.41 -      if      (lon_dir < 0 && lon1 > lon_break) lon1 = pgl_round(lon1 - 360);
    1.42 -      else if (lon_dir > 0 && lon1 < lon_break) lon1 = pgl_round(lon1 + 360);
    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 @@ -569,6 +569,231 @@
    1.48    return counter & 1;
    1.49  }
    1.50  
    1.51 +/* check if all points of the second cluster are inside the first cluster */
    1.52 +static inline bool pgl_all_cluster_points_in_cluster(
    1.53 +  pgl_cluster *outer, pgl_cluster *inner
    1.54 +) {
    1.55 +  int i, j;           /* i: entry, j: point in entry */
    1.56 +  int npoints;        /* number of points in entry */
    1.57 +  pgl_point *points;  /* array of points in entry */
    1.58 +  /* iterate over all entries of "inner" cluster */
    1.59 +  for (i=0; i<inner->nentries; i++) {
    1.60 +    /* get properties of entry */
    1.61 +    npoints = inner->entries[i].npoints;
    1.62 +    points = PGL_ENTRY_POINTS(inner, i);
    1.63 +    /* iterate over all points in entry of "inner" cluster */
    1.64 +    for (j=0; j<npoints; j++) {
    1.65 +      /* return false if one point of inner cluster is not in outer cluster */
    1.66 +      if (!pgl_point_in_cluster(points+j, outer)) return false;
    1.67 +    }
    1.68 +  }
    1.69 +  /* otherwise return true */
    1.70 +  return true;
    1.71 +}
    1.72 +
    1.73 +/* check if any point the second cluster is inside the first cluster */
    1.74 +static inline bool pgl_any_cluster_points_in_cluster(
    1.75 +  pgl_cluster *outer, pgl_cluster *inner
    1.76 +) {
    1.77 +  int i, j;           /* i: entry, j: point in entry */
    1.78 +  int npoints;        /* number of points in entry */
    1.79 +  pgl_point *points;  /* array of points in entry */
    1.80 +  /* iterate over all entries of "inner" cluster */
    1.81 +  for (i=0; i<inner->nentries; i++) {
    1.82 +    /* get properties of entry */
    1.83 +    npoints = inner->entries[i].npoints;
    1.84 +    points = PGL_ENTRY_POINTS(inner, i);
    1.85 +    /* iterate over all points in entry of "inner" cluster */
    1.86 +    for (j=0; j<npoints; j++) {
    1.87 +      /* return true if one point of inner cluster is in outer cluster */
    1.88 +      if (pgl_point_in_cluster(points+j, outer)) return true;
    1.89 +    }
    1.90 +  }
    1.91 +  /* otherwise return false */
    1.92 +  return false;
    1.93 +}
    1.94 +
    1.95 +/* check if line segment crosses line */
    1.96 +/* returns -1 if yes, 1 if no, and 0 in corner cases */
    1.97 +/* NOTE: each line (segment) must have a length greater than zero */
    1.98 +static inline double pgl_lseg_crosses_line(
    1.99 +  double seg_x1,  double seg_y1,  double seg_x2,  double seg_y2,
   1.100 +  double line_x1, double line_y1, double line_x2, double line_y2,
   1.101 +  bool strict
   1.102 +) {
   1.103 +  double value = (
   1.104 +    (seg_x1-line_x1) * (line_y2-line_y1) -
   1.105 +    (seg_y1-line_y1) * (line_x2-line_x1)
   1.106 +  ) * (
   1.107 +    (seg_x2-line_x1) * (line_y2-line_y1) -
   1.108 +    (seg_y2-line_y1) * (line_x2-line_x1)
   1.109 +  );
   1.110 +  if (strict) return value < 0;
   1.111 +  else return value <= 0;
   1.112 +}
   1.113 +
   1.114 +/* check if paths and outlines of two clusters overlap */
   1.115 +/* (set strict to true to disregard corner cases) */
   1.116 +static bool pgl_outlines_overlap(
   1.117 +  pgl_cluster *cluster1, pgl_cluster *cluster2, bool strict
   1.118 +) {
   1.119 +  int i1, j1, k1;  /* i: entry, j: point in entry, k: next point in entry */
   1.120 +  int i2, j2, k2;
   1.121 +  int entrytype1, entrytype2;     /* type of entry */
   1.122 +  int npoints1, npoints2;         /* number of points in entry */
   1.123 +  pgl_point *points1;             /* array of points in entry of cluster1 */
   1.124 +  pgl_point *points2;             /* array of points in entry of cluster2 */
   1.125 +  int lon_dir1, lon_dir2;         /* first vertex west (-1) or east (+1) */
   1.126 +  double lon_break1, lon_break2;  /* antipodal longitude of first vertex */
   1.127 +  double lat11, lon11;  /* latitude and (adjusted) longitude of vertex */
   1.128 +  double lat12, lon12;  /* latitude and (adjusted) longitude of next vertex */
   1.129 +  double lat21, lon21;  /* latitude and (adjusted) longitudes for cluster2 */
   1.130 +  double lat22, lon22;
   1.131 +  double wrapvalue;     /* temporary helper value to adjust wrap-around */  
   1.132 +  /* iterate over all entries of cluster1 */
   1.133 +  for (i1=0; i1<cluster1->nentries; i1++) {
   1.134 +    /* get properties of entry in cluster1 and skip points */
   1.135 +    npoints1 = cluster1->entries[i1].npoints;
   1.136 +    if (npoints1 < 2) continue;
   1.137 +    entrytype1 = cluster1->entries[i1].entrytype;
   1.138 +    points1 = PGL_ENTRY_POINTS(cluster1, i1);
   1.139 +    /* determine east/west orientation of first point and calculate antipodal
   1.140 +       longitude */
   1.141 +    lon_break1 = points1[0].lon;
   1.142 +    if (lon_break1 < 0) {
   1.143 +      lon_dir1   = -1;
   1.144 +      lon_break1 = pgl_round(lon_break1 + 180);
   1.145 +    } else if (lon_break1 > 0) {
   1.146 +      lon_dir1   = 1;
   1.147 +      lon_break1 = pgl_round(lon_break1 - 180);
   1.148 +    } else lon_dir1 = 0;
   1.149 +    /* iterate over all edges and vertices in cluster1 */
   1.150 +    for (j1=0; j1<npoints1; j1++) {
   1.151 +      /* calculate index of next vertex */
   1.152 +      k1 = (j1+1) % npoints1;
   1.153 +      /* skip last edge unless entry is (closed) outline or polygon */
   1.154 +      if (
   1.155 +        k1 == 0 &&
   1.156 +        entrytype1 != PGL_ENTRY_OUTLINE &&
   1.157 +        entrytype1 != PGL_ENTRY_POLYGON
   1.158 +      ) continue;
   1.159 +      /* use previously calculated values for lat1 and lon1 if possible */
   1.160 +      if (j1) {
   1.161 +        lat11 = lat12;
   1.162 +        lon11 = lon12;
   1.163 +      } else {
   1.164 +        /* otherwise get latitude and longitude values of first vertex */
   1.165 +        lat11 = points1[0].lat;
   1.166 +        lon11 = points1[0].lon;
   1.167 +        /* and consider longitude wrap-around for first vertex */
   1.168 +        if      (lon_dir1<0 && lon11>lon_break1) lon11 = pgl_round(lon11-360);
   1.169 +        else if (lon_dir1>0 && lon11<lon_break1) lon11 = pgl_round(lon11+360);
   1.170 +      }
   1.171 +      /* get latitude and longitude of next vertex */
   1.172 +      lat12 = points1[k1].lat;
   1.173 +      lon12 = points1[k1].lon;
   1.174 +      /* consider longitude wrap-around for next vertex */
   1.175 +      if      (lon_dir1<0 && lon12>lon_break1) lon12 = pgl_round(lon12-360);
   1.176 +      else if (lon_dir1>0 && lon12<lon_break1) lon12 = pgl_round(lon12+360);
   1.177 +      /* skip degenerated edges */
   1.178 +      if (lat11 == lat12 && lon11 == lon12) continue;
   1.179 +      /* iterate over all entries of cluster2 */
   1.180 +      for (i2=0; i2<cluster2->nentries; i2++) {
   1.181 +        /* get points and number of points of entry in cluster2 */
   1.182 +        npoints2 = cluster2->entries[i2].npoints;
   1.183 +        if (npoints2 < 2) continue;
   1.184 +        entrytype2 = cluster2->entries[i2].entrytype;
   1.185 +        points2 = PGL_ENTRY_POINTS(cluster2, i2);
   1.186 +        /* determine east/west orientation of first point and calculate antipodal
   1.187 +           longitude */
   1.188 +        lon_break2 = points2[0].lon;
   1.189 +        if (lon_break2 < 0) {
   1.190 +          lon_dir2   = -1;
   1.191 +          lon_break2 = pgl_round(lon_break2 + 180);
   1.192 +        } else if (lon_break2 > 0) {
   1.193 +          lon_dir2   = 1;
   1.194 +          lon_break2 = pgl_round(lon_break2 - 180);
   1.195 +        } else lon_dir2 = 0;
   1.196 +        /* iterate over all edges and vertices in cluster2 */
   1.197 +        for (j2=0; j2<npoints2; j2++) {
   1.198 +          /* calculate index of next vertex */
   1.199 +          k2 = (j2+1) % npoints2;
   1.200 +          /* skip last edge unless entry is (closed) outline or polygon */
   1.201 +          if (
   1.202 +            k2 == 0 &&
   1.203 +            entrytype2 != PGL_ENTRY_OUTLINE &&
   1.204 +            entrytype2 != PGL_ENTRY_POLYGON
   1.205 +          ) continue;
   1.206 +          /* use previously calculated values for lat1 and lon1 if possible */
   1.207 +          if (j2) {
   1.208 +            lat21 = lat22;
   1.209 +            lon21 = lon22;
   1.210 +          } else {
   1.211 +            /* otherwise get latitude and longitude values of first vertex */
   1.212 +            lat21 = points2[0].lat;
   1.213 +            lon21 = points2[0].lon;
   1.214 +            /* and consider longitude wrap-around for first vertex */
   1.215 +            if      (lon_dir2<0 && lon21>lon_break2) lon21 = pgl_round(lon21-360);
   1.216 +            else if (lon_dir2>0 && lon21<lon_break2) lon21 = pgl_round(lon21+360);
   1.217 +          }
   1.218 +          /* get latitude and longitude of next vertex */
   1.219 +          lat22 = points2[k2].lat;
   1.220 +          lon22 = points2[k2].lon;
   1.221 +          /* consider longitude wrap-around for next vertex */
   1.222 +          if      (lon_dir2<0 && lon22>lon_break2) lon22 = pgl_round(lon22-360);
   1.223 +          else if (lon_dir2>0 && lon22<lon_break2) lon22 = pgl_round(lon22+360);
   1.224 +          /* skip degenerated edges */
   1.225 +          if (lat21 == lat22 && lon21 == lon22) continue;
   1.226 +          /* perform another wrap-around where necessary */
   1.227 +          /* TODO: improve performance of whole wrap-around mechanism */
   1.228 +          wrapvalue = (lon21 + lon22) - (lon11 + lon12);
   1.229 +          if (wrapvalue > 360) {
   1.230 +            lon21 = pgl_round(lon21 - 360);
   1.231 +            lon22 = pgl_round(lon22 - 360);
   1.232 +          } else if (wrapvalue < -360) {
   1.233 +            lon21 = pgl_round(lon21 + 360);
   1.234 +            lon22 = pgl_round(lon22 + 360);
   1.235 +          }
   1.236 +          /* return true if segments overlap */
   1.237 +          if (
   1.238 +            pgl_lseg_crosses_line(
   1.239 +              lat11, lon11, lat12, lon12,
   1.240 +              lat21, lon21, lat22, lon22,
   1.241 +              strict
   1.242 +            ) && pgl_lseg_crosses_line(
   1.243 +              lat21, lon21, lat22, lon22,
   1.244 +              lat11, lon11, lat12, lon12,
   1.245 +              strict
   1.246 +            )
   1.247 +          ) {
   1.248 +            return true;
   1.249 +          }
   1.250 +        }
   1.251 +      }
   1.252 +    }
   1.253 +  }
   1.254 +  /* otherwise return false */
   1.255 +  return false;
   1.256 +}
   1.257 +
   1.258 +/* check if second cluster is completely contained in first cluster */
   1.259 +static bool pgl_cluster_in_cluster(pgl_cluster *outer, pgl_cluster *inner) {
   1.260 +  if (!pgl_all_cluster_points_in_cluster(outer, inner)) return false;
   1.261 +  if (pgl_outlines_overlap(outer, inner, true)) return false;
   1.262 +  return true;
   1.263 +}
   1.264 +
   1.265 +/* check if two clusters overlap */
   1.266 +static bool pgl_clusters_overlap(
   1.267 +  pgl_cluster *cluster1, pgl_cluster *cluster2
   1.268 +) {
   1.269 +  if (pgl_any_cluster_points_in_cluster(cluster1, cluster2)) return true;
   1.270 +  if (pgl_any_cluster_points_in_cluster(cluster2, cluster1)) return true;
   1.271 +  if (pgl_outlines_overlap(cluster1, cluster2, false)) return true;
   1.272 +  return false;
   1.273 +}
   1.274 +
   1.275 +
   1.276  /* calculate (approximate) distance between point and cluster */
   1.277  static double pgl_point_cluster_distance(pgl_point *point, pgl_cluster *cluster) {
   1.278    int i, j, k;  /* i: entry, j: point in entry, k: next point in entry */
   1.279 @@ -622,12 +847,18 @@
   1.280      else if (lon_dir > 0 && lon0 < lon_break) lon0 += 360;
   1.281      /* iterate over all edges and vertices */
   1.282      for (j=0; j<npoints; j++) {
   1.283 -      /* get latitude and longitude values of current point */
   1.284 -      lat1 = points[j].lat;
   1.285 -      lon1 = points[j].lon;
   1.286 -      /* consider longitude wrap-around for current point */
   1.287 -      if      (lon_dir < 0 && lon1 > lon_break) lon1 -= 360;
   1.288 -      else if (lon_dir > 0 && lon1 < lon_break) lon1 += 360;
   1.289 +      /* use previously calculated values for lat1 and lon1 if possible */
   1.290 +      if (j) {
   1.291 +        lat1 = lat2;
   1.292 +        lon1 = lon2;
   1.293 +      } else {
   1.294 +        /* otherwise get latitude and longitude values of first vertex */
   1.295 +        lat1 = points[0].lat;
   1.296 +        lon1 = points[0].lon;
   1.297 +        /* and consider longitude wrap-around for first vertex */
   1.298 +        if      (lon_dir < 0 && lon1 > lon_break) lon1 -= 360;
   1.299 +        else if (lon_dir > 0 && lon1 < lon_break) lon1 += 360;
   1.300 +      }
   1.301        /* calculate distance to vertex */
   1.302        dist = pgl_distance(lat0, lon0, lat1, lon1);
   1.303        /* store calculated distance if smallest */
   1.304 @@ -640,10 +871,10 @@
   1.305          entrytype != PGL_ENTRY_OUTLINE &&
   1.306          entrytype != PGL_ENTRY_POLYGON
   1.307        ) continue;
   1.308 -      /* get latitude and longitude values of next point */
   1.309 +      /* get latitude and longitude of next vertex */
   1.310        lat2 = points[k].lat;
   1.311        lon2 = points[k].lon;
   1.312 -      /* consider longitude wrap-around for next point */
   1.313 +      /* consider longitude wrap-around for next vertex */
   1.314        if      (lon_dir < 0 && lon2 > lon_break) lon2 -= 360;
   1.315        else if (lon_dir > 0 && lon2 < lon_break) lon2 += 360;
   1.316        /* go to next vertex and edge if edge is degenerated */
   1.317 @@ -669,65 +900,73 @@
   1.318    return min_dist;
   1.319  }
   1.320  
   1.321 +/* calculate (approximate) distance between two clusters */
   1.322 +static double pgl_cluster_distance(pgl_cluster *cluster1, pgl_cluster *cluster2) {
   1.323 +  int i, j;                    /* i: entry, j: point in entry */
   1.324 +  int npoints;                 /* number of points in entry */
   1.325 +  pgl_point *points;           /* array of points in entry */
   1.326 +  double dist;                 /* distance calculated in one step */
   1.327 +  double min_dist = INFINITY;  /* minimum distance */
   1.328 +  /* consider distance from each point in one cluster to the whole other */
   1.329 +  for (i=0; i<cluster1->nentries; i++) {
   1.330 +    npoints = cluster1->entries[i].npoints;
   1.331 +    points = PGL_ENTRY_POINTS(cluster1, i);
   1.332 +    for (j=0; j<npoints; j++) {
   1.333 +      dist = pgl_point_cluster_distance(points+j, cluster2);
   1.334 +      if (dist == 0) return dist;
   1.335 +      if (dist < min_dist) min_dist = dist;
   1.336 +    }
   1.337 +  }
   1.338 +  /* consider distance from each point in other cluster to the first cluster */
   1.339 +  for (i=0; i<cluster2->nentries; i++) {
   1.340 +    npoints = cluster2->entries[i].npoints;
   1.341 +    points = PGL_ENTRY_POINTS(cluster2, i);
   1.342 +    for (j=0; j<npoints; j++) {
   1.343 +      dist = pgl_point_cluster_distance(points+j, cluster1);
   1.344 +      if (dist == 0) return dist;
   1.345 +      if (dist < min_dist) min_dist = dist;
   1.346 +    }
   1.347 +  }
   1.348 +  return min_dist;
   1.349 +}
   1.350 +
   1.351  /* estimator function for distance between box and point */
   1.352 -/* allowed to return smaller values than actually correct */
   1.353 +/* always returns a smaller value than actually correct or zero */
   1.354  static double pgl_estimate_point_box_distance(pgl_point *point, pgl_box *box) {
   1.355 -  double dlon;  /* longitude range of box (delta longitude) */
   1.356 -  double h;     /* half of distance along meridian */
   1.357 -  double d;     /* distance between both southern or both northern points */
   1.358 -  double cur_dist;  /* calculated distance */
   1.359 -  double min_dist;  /* minimum distance calculated */
   1.360 -  /* return infinity if bounding box is empty */
   1.361 +  double dlon;      /* longitude range of box (delta longitude) */
   1.362 +  double distance;  /* return value */
   1.363 +  /* return infinity if box is empty */
   1.364    if (box->lat_min > box->lat_max) return INFINITY;
   1.365 -  /* return zero if point is inside bounding box */
   1.366 +  /* return zero if point is inside box */
   1.367    if (pgl_point_in_box(point, box)) return 0;
   1.368    /* calculate delta longitude */
   1.369    dlon = box->lon_max - box->lon_min;
   1.370    if (dlon < 0) dlon += 360;  /* 180th meridian crossed */
   1.371 -  /* if delta longitude is greater than 180 degrees, perform safe fall-back */
   1.372 -  if (dlon > 180) return 0;
   1.373 -  /* calculate half of distance along meridian */
   1.374 -  h = pgl_distance(box->lat_min, 0, box->lat_max, 0) / 2;
   1.375 -  /* calculate full distance between southern points */
   1.376 -  d = pgl_distance(box->lat_min, 0, box->lat_min, dlon);
   1.377 -  /* calculate maximum of full distance and half distance */
   1.378 -  if (h > d) d = h;
   1.379 -  /* calculate distance from point to first southern vertex and substract
   1.380 -     maximum error */
   1.381 -  min_dist = pgl_distance(
   1.382 -    point->lat, point->lon, box->lat_min, box->lon_min
   1.383 -  ) - d;
   1.384 -  /* return zero if estimated distance is smaller than zero */
   1.385 -  if (min_dist <= 0) return 0;
   1.386 -  /* repeat procedure with second southern vertex */
   1.387 -  cur_dist = pgl_distance(
   1.388 -    point->lat, point->lon, box->lat_min, box->lon_max
   1.389 -  ) - d;
   1.390 -  if (cur_dist <= 0) return 0;
   1.391 -  if (cur_dist < min_dist) min_dist = cur_dist;
   1.392 -  /* calculate full distance between northern points */
   1.393 -  d = pgl_distance(box->lat_max, 0, box->lat_max, dlon);
   1.394 -  /* calculate maximum of full distance and half distance */
   1.395 -  if (h > d) d = h;
   1.396 -  /* repeat procedure with northern vertices */
   1.397 -  cur_dist = pgl_distance(
   1.398 -    point->lat, point->lon, box->lat_max, box->lon_max
   1.399 -  ) - d;
   1.400 -  if (cur_dist <= 0) return 0;
   1.401 -  if (cur_dist < min_dist) min_dist = cur_dist;
   1.402 -  cur_dist = pgl_distance(
   1.403 -    point->lat, point->lon, box->lat_max, box->lon_min
   1.404 -  ) - d;
   1.405 -  if (cur_dist <= 0) return 0;
   1.406 -  if (cur_dist < min_dist) min_dist = cur_dist;
   1.407 -  /* return smallest value (unless already returned zero) */
   1.408 -  return min_dist;
   1.409 +  /* if delta longitude is greater than 150 degrees, perform safe fall-back */
   1.410 +  if (dlon > 150) return 0;
   1.411 +  /* calculate lower limit for distance (formula below requires dlon <= 150) */
   1.412 +  /* TODO: provide better estimation function to improve performance */
   1.413 +  distance = (
   1.414 +    (1.0-1e-14) * pgl_distance(
   1.415 +      point->lat,
   1.416 +      point->lon,
   1.417 +      (box->lat_min + box->lat_max) / 2,
   1.418 +      box->lon_min + dlon/2
   1.419 +    ) - pgl_distance(
   1.420 +      box->lat_min, box->lon_min,
   1.421 +      box->lat_max, box->lon_max
   1.422 +    )
   1.423 +  );
   1.424 +  /* truncate negative results to zero */
   1.425 +  if (distance <= 0) distance = 0;
   1.426 +  /* return result */
   1.427 +  return distance;
   1.428  }
   1.429  
   1.430  
   1.431 -/*----------------------------*
   1.432 - *  fractal geographic index  *
   1.433 - *----------------------------*/
   1.434 +/*-------------------------------------------------*
   1.435 + *  geographic index based on space-filling curve  *
   1.436 + *-------------------------------------------------*/
   1.437  
   1.438  /* number of bytes used for geographic (center) position in keys */
   1.439  #define PGL_KEY_LATLON_BYTELEN 7
   1.440 @@ -735,9 +974,6 @@
   1.441  /* maximum reference value for logarithmic size of geographic objects */
   1.442  #define PGL_AREAKEY_REFOBJSIZE (PGL_DIAMETER/3.0)  /* can be tweaked */
   1.443  
   1.444 -/* safety margin to avoid floating point errors in distance estimation */
   1.445 -#define PGL_FPE_SAFETY (1.0+1e-14)  /* slightly greater than 1.0 */
   1.446 -
   1.447  /* pointer to index key (either pgl_pointkey or pgl_areakey) */
   1.448  typedef unsigned char *pgl_keyptr;
   1.449  
   1.450 @@ -1094,7 +1330,7 @@
   1.451  }
   1.452  
   1.453  /* estimator function for distance between point and index key */
   1.454 -/* allowed to return smaller values than actually correct */
   1.455 +/* always returns a smaller value than actually correct or zero */
   1.456  static double pgl_estimate_key_distance(pgl_keyptr key, pgl_point *point) {
   1.457    pgl_box box;  /* center(!) bounding box of area index key */
   1.458    /* calculate center(!) bounding box and maximum radius of objects covered
   1.459 @@ -1103,11 +1339,7 @@
   1.460    /* calculate estimated distance between bounding box of center point of
   1.461       indexed object and point passed as second argument, then substract maximum
   1.462       radius of objects covered by index key */
   1.463 -  /* (use PGL_FPE_SAFETY factor to cope with minor floating point errors) */
   1.464 -  distance = (
   1.465 -    pgl_estimate_point_box_distance(point, &box) / PGL_FPE_SAFETY -
   1.466 -    distance * PGL_FPE_SAFETY
   1.467 -  );
   1.468 +  distance = pgl_estimate_point_box_distance(point, &box) - distance;
   1.469    /* truncate negative results to zero */
   1.470    if (distance <= 0) distance = 0;
   1.471    /* return result */
   1.472 @@ -2102,7 +2334,16 @@
   1.473  Datum pgl_epoint_ecluster_overlap(PG_FUNCTION_ARGS) {
   1.474    pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
   1.475    pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
   1.476 -  bool retval = pgl_point_in_cluster(point, cluster);
   1.477 +  bool retval;
   1.478 +  /* points outside bounding circle are always assumed to be non-overlapping
   1.479 +     (necessary for consistent table and index scans) */
   1.480 +  if (
   1.481 +    pgl_distance(
   1.482 +      point->lat, point->lon,
   1.483 +      cluster->bounding.center.lat, cluster->bounding.center.lon
   1.484 +    ) > cluster->bounding.radius
   1.485 +  ) retval = false;
   1.486 +  else retval = pgl_point_in_cluster(point, cluster);
   1.487    PG_FREE_IF_COPY(cluster, 1);
   1.488    PG_RETURN_BOOL(retval);
   1.489  }
   1.490 @@ -2189,6 +2430,27 @@
   1.491    PG_RETURN_BOOL(retval);
   1.492  }
   1.493  
   1.494 +/* check if two clusters overlap (overlap operator "&&") in SQL */
   1.495 +PG_FUNCTION_INFO_V1(pgl_ecluster_overlap);
   1.496 +Datum pgl_ecluster_overlap(PG_FUNCTION_ARGS) {
   1.497 +  pgl_cluster *cluster1 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
   1.498 +  pgl_cluster *cluster2 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
   1.499 +  bool retval;
   1.500 +  /* clusters with non-touching bounding circles are always assumed to be
   1.501 +     non-overlapping (improves performance and is necessary for consistent
   1.502 +     table and index scans) */
   1.503 +  if (
   1.504 +    pgl_distance(
   1.505 +      cluster1->bounding.center.lat, cluster1->bounding.center.lon,
   1.506 +      cluster2->bounding.center.lat, cluster2->bounding.center.lon
   1.507 +    ) > cluster1->bounding.radius + cluster2->bounding.radius
   1.508 +  ) retval = false;
   1.509 +  else retval = pgl_clusters_overlap(cluster1, cluster2);
   1.510 +  PG_FREE_IF_COPY(cluster1, 0);
   1.511 +  PG_FREE_IF_COPY(cluster2, 1);
   1.512 +  PG_RETURN_BOOL(retval);
   1.513 +}
   1.514 +
   1.515  /* check if two clusters may overlap (lossy overlap operator "&&+") in SQL */
   1.516  PG_FUNCTION_INFO_V1(pgl_ecluster_may_overlap);
   1.517  Datum pgl_ecluster_may_overlap(PG_FUNCTION_ARGS) {
   1.518 @@ -2203,6 +2465,27 @@
   1.519    PG_RETURN_BOOL(retval);
   1.520  }
   1.521  
   1.522 +/* check if second cluster is in first cluster (cont. operator "@>) in SQL */
   1.523 +PG_FUNCTION_INFO_V1(pgl_ecluster_contains);
   1.524 +Datum pgl_ecluster_contains(PG_FUNCTION_ARGS) {
   1.525 +  pgl_cluster *outer = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
   1.526 +  pgl_cluster *inner = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
   1.527 +  bool retval;
   1.528 +  /* clusters with non-touching bounding circles are always assumed to be
   1.529 +     non-overlapping (improves performance and is necessary for consistent
   1.530 +     table and index scans) */
   1.531 +  if (
   1.532 +    pgl_distance(
   1.533 +      outer->bounding.center.lat, outer->bounding.center.lon,
   1.534 +      inner->bounding.center.lat, inner->bounding.center.lon
   1.535 +    ) > outer->bounding.radius + inner->bounding.radius
   1.536 +  ) retval = false;
   1.537 +  else retval = pgl_cluster_in_cluster(outer, inner);
   1.538 +  PG_FREE_IF_COPY(outer, 0);
   1.539 +  PG_FREE_IF_COPY(inner, 1);
   1.540 +  PG_RETURN_BOOL(retval);
   1.541 +}
   1.542 +
   1.543  /* calculate distance between two points ("<->" operator) in SQL */
   1.544  PG_FUNCTION_INFO_V1(pgl_epoint_distance);
   1.545  Datum pgl_epoint_distance(PG_FUNCTION_ARGS) {
   1.546 @@ -2258,6 +2541,17 @@
   1.547    PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance);
   1.548  }
   1.549  
   1.550 +/* calculate distance between two clusters ("<->" operator) in SQL */
   1.551 +PG_FUNCTION_INFO_V1(pgl_ecluster_distance);
   1.552 +Datum pgl_ecluster_distance(PG_FUNCTION_ARGS) {
   1.553 +  pgl_cluster *cluster1 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
   1.554 +  pgl_cluster *cluster2 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
   1.555 +  double retval = pgl_cluster_distance(cluster1, cluster2);
   1.556 +  PG_FREE_IF_COPY(cluster1, 0);
   1.557 +  PG_FREE_IF_COPY(cluster2, 1);
   1.558 +  PG_RETURN_FLOAT8(retval);
   1.559 +}
   1.560 +
   1.561  
   1.562  /*-----------------------------------------------------------*
   1.563   *  B-tree comparison operators and index support functions  *