| rev |
line source |
|
jbe@0
|
1
|
|
jbe@0
|
2 /*-------------*
|
|
jbe@0
|
3 * C prelude *
|
|
jbe@0
|
4 *-------------*/
|
|
jbe@0
|
5
|
|
jbe@0
|
6 #include "postgres.h"
|
|
jbe@0
|
7 #include "fmgr.h"
|
|
jbe@0
|
8 #include "libpq/pqformat.h"
|
|
jbe@0
|
9 #include "access/gist.h"
|
|
jbe@0
|
10 #include "access/stratnum.h"
|
|
jbe@0
|
11 #include "utils/array.h"
|
|
jbe@51
|
12 #include <limits.h>
|
|
jbe@0
|
13 #include <math.h>
|
|
jbe@0
|
14
|
|
jbe@0
|
15 #ifdef PG_MODULE_MAGIC
|
|
jbe@0
|
16 PG_MODULE_MAGIC;
|
|
jbe@0
|
17 #endif
|
|
jbe@0
|
18
|
|
jbe@0
|
19 #if INT_MAX < 2147483647
|
|
jbe@0
|
20 #error Expected int type to be at least 32 bit wide
|
|
jbe@0
|
21 #endif
|
|
jbe@0
|
22
|
|
jbe@0
|
23
|
|
jbe@0
|
24 /*---------------------------------*
|
|
jbe@0
|
25 * distance calculation on earth *
|
|
jbe@0
|
26 * (using WGS-84 spheroid) *
|
|
jbe@0
|
27 *---------------------------------*/
|
|
jbe@0
|
28
|
|
jbe@0
|
29 /* WGS-84 spheroid with following parameters:
|
|
jbe@0
|
30 semi-major axis a = 6378137
|
|
jbe@0
|
31 semi-minor axis b = a * (1 - 1/298.257223563)
|
|
jbe@0
|
32 estimated diameter = 2 * (2*a+b)/3
|
|
jbe@0
|
33 */
|
|
jbe@0
|
34 #define PGL_SPHEROID_A 6378137.0 /* semi major axis */
|
|
jbe@0
|
35 #define PGL_SPHEROID_F (1.0/298.257223563) /* flattening */
|
|
jbe@0
|
36 #define PGL_SPHEROID_B (PGL_SPHEROID_A * (1.0-PGL_SPHEROID_F))
|
|
jbe@0
|
37 #define PGL_EPS2 ( ( PGL_SPHEROID_A * PGL_SPHEROID_A - \
|
|
jbe@0
|
38 PGL_SPHEROID_B * PGL_SPHEROID_B ) / \
|
|
jbe@0
|
39 ( PGL_SPHEROID_A * PGL_SPHEROID_A ) )
|
|
jbe@0
|
40 #define PGL_SUBEPS2 (1.0-PGL_EPS2)
|
|
jbe@42
|
41 #define PGL_RADIUS ((2.0*PGL_SPHEROID_A + PGL_SPHEROID_B) / 3.0)
|
|
jbe@42
|
42 #define PGL_DIAMETER (2.0 * PGL_RADIUS)
|
|
jbe@0
|
43 #define PGL_SCALE (PGL_SPHEROID_A / PGL_DIAMETER) /* semi-major ref. */
|
|
jbe@42
|
44 #define PGL_MAXDIST (PGL_RADIUS * M_PI) /* maximum distance */
|
|
jbe@42
|
45 #define PGL_FADELIMIT (PGL_MAXDIST / 3.0) /* 1/6 circumference */
|
|
jbe@0
|
46
|
|
jbe@0
|
47 /* calculate distance between two points on earth (given in degrees) */
|
|
jbe@0
|
48 static inline double pgl_distance(
|
|
jbe@0
|
49 double lat1, double lon1, double lat2, double lon2
|
|
jbe@0
|
50 ) {
|
|
jbe@0
|
51 float8 lat1cos, lat1sin, lat2cos, lat2sin, lon2cos, lon2sin;
|
|
jbe@0
|
52 float8 nphi1, nphi2, x1, z1, x2, y2, z2, g, s, t;
|
|
jbe@0
|
53 /* normalize delta longitude (lon2 > 0 && lon1 = 0) */
|
|
jbe@0
|
54 /* lon1 = 0 (not used anymore) */
|
|
jbe@0
|
55 lon2 = fabs(lon2-lon1);
|
|
jbe@0
|
56 /* convert to radians (first divide, then multiply) */
|
|
jbe@0
|
57 lat1 = (lat1 / 180.0) * M_PI;
|
|
jbe@0
|
58 lat2 = (lat2 / 180.0) * M_PI;
|
|
jbe@0
|
59 lon2 = (lon2 / 180.0) * M_PI;
|
|
jbe@0
|
60 /* make lat2 >= lat1 to ensure reversal-symmetry despite floating point
|
|
jbe@0
|
61 operations (lon2 >= lon1 is already ensured in a previous step) */
|
|
jbe@0
|
62 if (lat2 < lat1) { float8 swap = lat1; lat1 = lat2; lat2 = swap; }
|
|
jbe@0
|
63 /* calculate 3d coordinates on scaled ellipsoid which has an average diameter
|
|
jbe@0
|
64 of 1.0 */
|
|
jbe@0
|
65 lat1cos = cos(lat1); lat1sin = sin(lat1);
|
|
jbe@0
|
66 lat2cos = cos(lat2); lat2sin = sin(lat2);
|
|
jbe@0
|
67 lon2cos = cos(lon2); lon2sin = sin(lon2);
|
|
jbe@0
|
68 nphi1 = PGL_SCALE / sqrt(1 - PGL_EPS2 * lat1sin * lat1sin);
|
|
jbe@0
|
69 nphi2 = PGL_SCALE / sqrt(1 - PGL_EPS2 * lat2sin * lat2sin);
|
|
jbe@0
|
70 x1 = nphi1 * lat1cos;
|
|
jbe@0
|
71 z1 = nphi1 * PGL_SUBEPS2 * lat1sin;
|
|
jbe@0
|
72 x2 = nphi2 * lat2cos * lon2cos;
|
|
jbe@0
|
73 y2 = nphi2 * lat2cos * lon2sin;
|
|
jbe@0
|
74 z2 = nphi2 * PGL_SUBEPS2 * lat2sin;
|
|
jbe@0
|
75 /* calculate tunnel distance through scaled (diameter 1.0) ellipsoid */
|
|
jbe@0
|
76 g = sqrt((x2-x1)*(x2-x1) + y2*y2 + (z2-z1)*(z2-z1));
|
|
jbe@0
|
77 /* convert tunnel distance through scaled ellipsoid to approximated surface
|
|
jbe@0
|
78 distance on original ellipsoid */
|
|
jbe@0
|
79 if (g > 1.0) g = 1.0;
|
|
jbe@0
|
80 s = PGL_DIAMETER * asin(g);
|
|
jbe@0
|
81 /* return result only if small enough to be precise (less than 1/3 of
|
|
jbe@0
|
82 maximum possible distance) */
|
|
jbe@0
|
83 if (s <= PGL_FADELIMIT) return s;
|
|
jbe@0
|
84 /* calculate tunnel distance to antipodal point through scaled ellipsoid */
|
|
jbe@3
|
85 g = sqrt((x2+x1)*(x2+x1) + y2*y2 + (z2+z1)*(z2+z1));
|
|
jbe@0
|
86 /* convert tunnel distance to antipodal point through scaled ellipsoid to
|
|
jbe@0
|
87 approximated surface distance to antipodal point on original ellipsoid */
|
|
jbe@0
|
88 if (g > 1.0) g = 1.0;
|
|
jbe@0
|
89 t = PGL_DIAMETER * asin(g);
|
|
jbe@0
|
90 /* surface distance between original points can now be approximated by
|
|
jbe@0
|
91 substracting antipodal distance from maximum possible distance;
|
|
jbe@0
|
92 return result only if small enough (less than 1/3 of maximum possible
|
|
jbe@0
|
93 distance) */
|
|
jbe@0
|
94 if (t <= PGL_FADELIMIT) return PGL_MAXDIST-t;
|
|
jbe@0
|
95 /* otherwise crossfade direct and antipodal result to ensure monotonicity */
|
|
jbe@0
|
96 return (
|
|
jbe@0
|
97 (s * (t-PGL_FADELIMIT) + (PGL_MAXDIST-t) * (s-PGL_FADELIMIT)) /
|
|
jbe@0
|
98 (s + t - 2*PGL_FADELIMIT)
|
|
jbe@0
|
99 );
|
|
jbe@0
|
100 }
|
|
jbe@0
|
101
|
|
jbe@0
|
102 /* finite distance that can not be reached on earth */
|
|
jbe@0
|
103 #define PGL_ULTRA_DISTANCE (3 * PGL_MAXDIST)
|
|
jbe@0
|
104
|
|
jbe@0
|
105
|
|
jbe@0
|
106 /*--------------------------------*
|
|
jbe@0
|
107 * simple geographic data types *
|
|
jbe@0
|
108 *--------------------------------*/
|
|
jbe@0
|
109
|
|
jbe@0
|
110 /* point on earth given by latitude and longitude in degrees */
|
|
jbe@0
|
111 /* (type "epoint" in SQL) */
|
|
jbe@0
|
112 typedef struct {
|
|
jbe@0
|
113 double lat; /* between -90 and 90 (both inclusive) */
|
|
jbe@0
|
114 double lon; /* between -180 and 180 (both inclusive) */
|
|
jbe@0
|
115 } pgl_point;
|
|
jbe@0
|
116
|
|
jbe@0
|
117 /* box delimited by two parallels and two meridians (all in degrees) */
|
|
jbe@0
|
118 /* (type "ebox" in SQL) */
|
|
jbe@0
|
119 typedef struct {
|
|
jbe@0
|
120 double lat_min; /* between -90 and 90 (both inclusive) */
|
|
jbe@0
|
121 double lat_max; /* between -90 and 90 (both inclusive) */
|
|
jbe@0
|
122 double lon_min; /* between -180 and 180 (both inclusive) */
|
|
jbe@0
|
123 double lon_max; /* between -180 and 180 (both inclusive) */
|
|
jbe@0
|
124 /* if lat_min > lat_max, then box is empty */
|
|
jbe@0
|
125 /* if lon_min > lon_max, then 180th meridian is crossed */
|
|
jbe@0
|
126 } pgl_box;
|
|
jbe@0
|
127
|
|
jbe@0
|
128 /* circle on earth surface (for radial searches with fixed radius) */
|
|
jbe@0
|
129 /* (type "ecircle" in SQL) */
|
|
jbe@0
|
130 typedef struct {
|
|
jbe@0
|
131 pgl_point center;
|
|
jbe@0
|
132 double radius; /* positive (including +0 but excluding -0), or -INFINITY */
|
|
jbe@0
|
133 /* A negative radius (i.e. -INFINITY) denotes nothing (i.e. no point),
|
|
jbe@0
|
134 zero radius (0) denotes a single point,
|
|
jbe@0
|
135 a finite radius (0 < radius < INFINITY) denotes a filled circle, and
|
|
jbe@0
|
136 a radius of INFINITY is valid and means complete coverage of earth. */
|
|
jbe@0
|
137 } pgl_circle;
|
|
jbe@0
|
138
|
|
jbe@0
|
139
|
|
jbe@0
|
140 /*----------------------------------*
|
|
jbe@0
|
141 * geographic "cluster" data type *
|
|
jbe@0
|
142 *----------------------------------*/
|
|
jbe@0
|
143
|
|
jbe@0
|
144 /* A cluster is a collection of points, paths, outlines, and polygons. If two
|
|
jbe@0
|
145 polygons in a cluster overlap, the area covered by both polygons does not
|
|
jbe@0
|
146 belong to the cluster. This way, a cluster can be used to describe complex
|
|
jbe@0
|
147 shapes like polygons with holes. Outlines are non-filled polygons. Paths are
|
|
jbe@0
|
148 open by default (i.e. the last point in the list is not connected with the
|
|
jbe@0
|
149 first point in the list). Note that each outline or polygon in a cluster
|
|
jbe@0
|
150 must cover a longitude range of less than 180 degrees to avoid ambiguities.
|
|
jbe@0
|
151 Areas which are larger may be split into multiple polygons. */
|
|
jbe@0
|
152
|
|
jbe@0
|
153 /* maximum number of points in a cluster */
|
|
jbe@0
|
154 /* (limited to avoid integer overflows, e.g. when allocating memory) */
|
|
jbe@0
|
155 #define PGL_CLUSTER_MAXPOINTS 16777216
|
|
jbe@0
|
156
|
|
jbe@0
|
157 /* types of cluster entries */
|
|
jbe@0
|
158 #define PGL_ENTRY_POINT 1 /* a point */
|
|
jbe@0
|
159 #define PGL_ENTRY_PATH 2 /* a path from first point to last point */
|
|
jbe@0
|
160 #define PGL_ENTRY_OUTLINE 3 /* a non-filled polygon with given vertices */
|
|
jbe@0
|
161 #define PGL_ENTRY_POLYGON 4 /* a filled polygon with given vertices */
|
|
jbe@0
|
162
|
|
jbe@0
|
163 /* Entries of a cluster are described by two different structs: pgl_newentry
|
|
jbe@0
|
164 and pgl_entry. The first is used only during construction of a cluster, the
|
|
jbe@0
|
165 second is used in all other cases (e.g. when reading clusters from the
|
|
jbe@0
|
166 database, performing operations, etc). */
|
|
jbe@0
|
167
|
|
jbe@0
|
168 /* entry for new geographic cluster during construction of that cluster */
|
|
jbe@0
|
169 typedef struct {
|
|
jbe@0
|
170 int32_t entrytype;
|
|
jbe@0
|
171 int32_t npoints;
|
|
jbe@0
|
172 pgl_point *points; /* pointer to an array of points (pgl_point) */
|
|
jbe@0
|
173 } pgl_newentry;
|
|
jbe@0
|
174
|
|
jbe@0
|
175 /* entry of geographic cluster */
|
|
jbe@0
|
176 typedef struct {
|
|
jbe@0
|
177 int32_t entrytype; /* type of entry: point, path, outline, polygon */
|
|
jbe@0
|
178 int32_t npoints; /* number of stored points (set to 1 for point entry) */
|
|
jbe@0
|
179 int32_t offset; /* offset of pgl_point array from cluster base address */
|
|
jbe@0
|
180 /* use macro PGL_ENTRY_POINTS to obtain a pointer to the array of points */
|
|
jbe@0
|
181 } pgl_entry;
|
|
jbe@0
|
182
|
|
jbe@0
|
183 /* geographic cluster which is a collection of points, (open) paths, polygons,
|
|
jbe@0
|
184 and outlines (non-filled polygons) */
|
|
jbe@0
|
185 typedef struct {
|
|
jbe@0
|
186 char header[VARHDRSZ]; /* PostgreSQL header for variable size data types */
|
|
jbe@0
|
187 int32_t nentries; /* number of stored points */
|
|
jbe@0
|
188 pgl_circle bounding; /* bounding circle */
|
|
jbe@0
|
189 /* Note: bounding circle ensures alignment of pgl_cluster for points */
|
|
jbe@0
|
190 pgl_entry entries[FLEXIBLE_ARRAY_MEMBER]; /* var-length data */
|
|
jbe@0
|
191 } pgl_cluster;
|
|
jbe@0
|
192
|
|
jbe@0
|
193 /* macro to determine memory alignment of points */
|
|
jbe@0
|
194 /* (needed to store pgl_point array after entries in pgl_cluster) */
|
|
jbe@0
|
195 typedef struct { char dummy; pgl_point aligned; } pgl_point_alignment;
|
|
jbe@0
|
196 #define PGL_POINT_ALIGNMENT offsetof(pgl_point_alignment, aligned)
|
|
jbe@0
|
197
|
|
jbe@0
|
198 /* macro to extract a pointer to the array of points of a cluster entry */
|
|
jbe@0
|
199 #define PGL_ENTRY_POINTS(cluster, idx) \
|
|
jbe@0
|
200 ((pgl_point *)(((intptr_t)cluster)+(cluster)->entries[idx].offset))
|
|
jbe@0
|
201
|
|
jbe@0
|
202 /* convert pgl_newentry array to pgl_cluster */
|
|
jbe@42
|
203 /* NOTE: requires pgl_finalize_cluster to be called to finalize result */
|
|
jbe@0
|
204 static pgl_cluster *pgl_new_cluster(int nentries, pgl_newentry *entries) {
|
|
jbe@0
|
205 int i; /* index of current entry */
|
|
jbe@0
|
206 int npoints = 0; /* number of points in whole cluster */
|
|
jbe@0
|
207 int entry_npoints; /* number of points in current entry */
|
|
jbe@0
|
208 int points_offset = PGL_POINT_ALIGNMENT * (
|
|
jbe@0
|
209 ( offsetof(pgl_cluster, entries) +
|
|
jbe@0
|
210 nentries * sizeof(pgl_entry) +
|
|
jbe@0
|
211 PGL_POINT_ALIGNMENT - 1
|
|
jbe@0
|
212 ) / PGL_POINT_ALIGNMENT
|
|
jbe@0
|
213 ); /* offset of pgl_point array from base address (considering alignment) */
|
|
jbe@0
|
214 pgl_cluster *cluster; /* new cluster to be returned */
|
|
jbe@0
|
215 /* determine total number of points */
|
|
jbe@0
|
216 for (i=0; i<nentries; i++) npoints += entries[i].npoints;
|
|
jbe@0
|
217 /* allocate memory for cluster (including entries and points) */
|
|
jbe@0
|
218 cluster = palloc(points_offset + npoints * sizeof(pgl_point));
|
|
jbe@0
|
219 /* re-count total number of points to determine offset for each entry */
|
|
jbe@0
|
220 npoints = 0;
|
|
jbe@0
|
221 /* copy entries and points */
|
|
jbe@0
|
222 for (i=0; i<nentries; i++) {
|
|
jbe@0
|
223 /* determine number of points in entry */
|
|
jbe@0
|
224 entry_npoints = entries[i].npoints;
|
|
jbe@0
|
225 /* copy entry */
|
|
jbe@0
|
226 cluster->entries[i].entrytype = entries[i].entrytype;
|
|
jbe@0
|
227 cluster->entries[i].npoints = entry_npoints;
|
|
jbe@0
|
228 /* calculate offset (in bytes) of pgl_point array */
|
|
jbe@0
|
229 cluster->entries[i].offset = points_offset + npoints * sizeof(pgl_point);
|
|
jbe@0
|
230 /* copy points */
|
|
jbe@0
|
231 memcpy(
|
|
jbe@0
|
232 PGL_ENTRY_POINTS(cluster, i),
|
|
jbe@0
|
233 entries[i].points,
|
|
jbe@0
|
234 entry_npoints * sizeof(pgl_point)
|
|
jbe@0
|
235 );
|
|
jbe@0
|
236 /* update total number of points processed */
|
|
jbe@0
|
237 npoints += entry_npoints;
|
|
jbe@0
|
238 }
|
|
jbe@0
|
239 /* set number of entries in cluster */
|
|
jbe@0
|
240 cluster->nentries = nentries;
|
|
jbe@0
|
241 /* set PostgreSQL header for variable sized data */
|
|
jbe@0
|
242 SET_VARSIZE(cluster, points_offset + npoints * sizeof(pgl_point));
|
|
jbe@0
|
243 /* return newly created cluster */
|
|
jbe@0
|
244 return cluster;
|
|
jbe@0
|
245 }
|
|
jbe@0
|
246
|
|
jbe@0
|
247
|
|
jbe@46
|
248 /*----------------------------------------------*
|
|
jbe@46
|
249 * Geographic point with integer sample count *
|
|
jbe@46
|
250 * (needed for fair distance calculation) *
|
|
jbe@46
|
251 *----------------------------------------------*/
|
|
jbe@46
|
252
|
|
jbe@46
|
253 typedef struct {
|
|
jbe@46
|
254 pgl_point point; /* NOTE: point first to allow C cast to pgl_point */
|
|
jbe@46
|
255 int32 samples;
|
|
jbe@46
|
256 } pgl_point_sc;
|
|
jbe@46
|
257
|
|
jbe@46
|
258
|
|
jbe@0
|
259 /*----------------------------------------*
|
|
jbe@0
|
260 * C functions on geographic data types *
|
|
jbe@0
|
261 *----------------------------------------*/
|
|
jbe@0
|
262
|
|
jbe@0
|
263 /* round latitude or longitude to 12 digits after decimal point */
|
|
jbe@0
|
264 static inline double pgl_round(double val) {
|
|
jbe@0
|
265 return round(val * 1e12) / 1e12;
|
|
jbe@0
|
266 }
|
|
jbe@0
|
267
|
|
jbe@70
|
268 /* normalize longitude to be between -180 and 180 */
|
|
jbe@70
|
269 static inline double pgl_normalize(double lon, bool warn) {
|
|
jbe@70
|
270 if (lon < -180) {
|
|
jbe@70
|
271 if (warn) {
|
|
jbe@70
|
272 ereport(NOTICE, (errmsg("longitude west of 180th meridian normalized")));
|
|
jbe@70
|
273 }
|
|
jbe@70
|
274 lon += 360 - trunc(lon / 360) * 360;
|
|
jbe@70
|
275 } else if (lon > 180) {
|
|
jbe@70
|
276 if (warn) {
|
|
jbe@70
|
277 ereport(NOTICE, (errmsg("longitude east of 180th meridian normalized")));
|
|
jbe@70
|
278 }
|
|
jbe@70
|
279 lon -= 360 + trunc(lon / 360) * 360;
|
|
jbe@70
|
280 }
|
|
jbe@70
|
281 return lon;
|
|
jbe@70
|
282 }
|
|
jbe@70
|
283
|
|
jbe@0
|
284 /* compare two points */
|
|
jbe@0
|
285 /* (equality when same point on earth is described, otherwise an arbitrary
|
|
jbe@0
|
286 linear order) */
|
|
jbe@0
|
287 static int pgl_point_cmp(pgl_point *point1, pgl_point *point2) {
|
|
jbe@0
|
288 double lon1, lon2; /* modified longitudes for special cases */
|
|
jbe@0
|
289 /* use latitude as first ordering criterion */
|
|
jbe@0
|
290 if (point1->lat < point2->lat) return -1;
|
|
jbe@0
|
291 if (point1->lat > point2->lat) return 1;
|
|
jbe@0
|
292 /* determine modified longitudes (considering special case of poles and
|
|
jbe@0
|
293 180th meridian which can be described as W180 or E180) */
|
|
jbe@0
|
294 if (point1->lat == -90 || point1->lat == 90) lon1 = 0;
|
|
jbe@0
|
295 else if (point1->lon == 180) lon1 = -180;
|
|
jbe@0
|
296 else lon1 = point1->lon;
|
|
jbe@0
|
297 if (point2->lat == -90 || point2->lat == 90) lon2 = 0;
|
|
jbe@0
|
298 else if (point2->lon == 180) lon2 = -180;
|
|
jbe@0
|
299 else lon2 = point2->lon;
|
|
jbe@0
|
300 /* use (modified) longitude as secondary ordering criterion */
|
|
jbe@0
|
301 if (lon1 < lon2) return -1;
|
|
jbe@0
|
302 if (lon1 > lon2) return 1;
|
|
jbe@0
|
303 /* no difference found, points are equal */
|
|
jbe@0
|
304 return 0;
|
|
jbe@0
|
305 }
|
|
jbe@0
|
306
|
|
jbe@0
|
307 /* compare two boxes */
|
|
jbe@0
|
308 /* (equality when same box on earth is described, otherwise an arbitrary linear
|
|
jbe@0
|
309 order) */
|
|
jbe@0
|
310 static int pgl_box_cmp(pgl_box *box1, pgl_box *box2) {
|
|
jbe@0
|
311 /* two empty boxes are equal, and an empty box is always considered "less
|
|
jbe@0
|
312 than" a non-empty box */
|
|
jbe@0
|
313 if (box1->lat_min> box1->lat_max && box2->lat_min<=box2->lat_max) return -1;
|
|
jbe@0
|
314 if (box1->lat_min> box1->lat_max && box2->lat_min> box2->lat_max) return 0;
|
|
jbe@0
|
315 if (box1->lat_min<=box1->lat_max && box2->lat_min> box2->lat_max) return 1;
|
|
jbe@0
|
316 /* use southern border as first ordering criterion */
|
|
jbe@0
|
317 if (box1->lat_min < box2->lat_min) return -1;
|
|
jbe@0
|
318 if (box1->lat_min > box2->lat_min) return 1;
|
|
jbe@0
|
319 /* use northern border as second ordering criterion */
|
|
jbe@0
|
320 if (box1->lat_max < box2->lat_max) return -1;
|
|
jbe@0
|
321 if (box1->lat_max > box2->lat_max) return 1;
|
|
jbe@0
|
322 /* use western border as third ordering criterion */
|
|
jbe@0
|
323 if (box1->lon_min < box2->lon_min) return -1;
|
|
jbe@0
|
324 if (box1->lon_min > box2->lon_min) return 1;
|
|
jbe@0
|
325 /* use eastern border as fourth ordering criterion */
|
|
jbe@0
|
326 if (box1->lon_max < box2->lon_max) return -1;
|
|
jbe@0
|
327 if (box1->lon_max > box2->lon_max) return 1;
|
|
jbe@0
|
328 /* no difference found, boxes are equal */
|
|
jbe@0
|
329 return 0;
|
|
jbe@0
|
330 }
|
|
jbe@0
|
331
|
|
jbe@0
|
332 /* compare two circles */
|
|
jbe@0
|
333 /* (equality when same circle on earth is described, otherwise an arbitrary
|
|
jbe@0
|
334 linear order) */
|
|
jbe@0
|
335 static int pgl_circle_cmp(pgl_circle *circle1, pgl_circle *circle2) {
|
|
jbe@0
|
336 /* two circles with same infinite radius (positive or negative infinity) are
|
|
jbe@0
|
337 considered equal independently of center point */
|
|
jbe@0
|
338 if (
|
|
jbe@0
|
339 !isfinite(circle1->radius) && !isfinite(circle2->radius) &&
|
|
jbe@0
|
340 circle1->radius == circle2->radius
|
|
jbe@0
|
341 ) return 0;
|
|
jbe@0
|
342 /* use radius as first ordering criterion */
|
|
jbe@0
|
343 if (circle1->radius < circle2->radius) return -1;
|
|
jbe@0
|
344 if (circle1->radius > circle2->radius) return 1;
|
|
jbe@0
|
345 /* use center point as secondary ordering criterion */
|
|
jbe@0
|
346 return pgl_point_cmp(&(circle1->center), &(circle2->center));
|
|
jbe@0
|
347 }
|
|
jbe@0
|
348
|
|
jbe@0
|
349 /* set box to empty box*/
|
|
jbe@0
|
350 static void pgl_box_set_empty(pgl_box *box) {
|
|
jbe@0
|
351 box->lat_min = INFINITY;
|
|
jbe@0
|
352 box->lat_max = -INFINITY;
|
|
jbe@0
|
353 box->lon_min = 0;
|
|
jbe@0
|
354 box->lon_max = 0;
|
|
jbe@0
|
355 }
|
|
jbe@0
|
356
|
|
jbe@0
|
357 /* check if point is inside a box */
|
|
jbe@0
|
358 static bool pgl_point_in_box(pgl_point *point, pgl_box *box) {
|
|
jbe@0
|
359 return (
|
|
jbe@0
|
360 point->lat >= box->lat_min && point->lat <= box->lat_max && (
|
|
jbe@0
|
361 (box->lon_min > box->lon_max) ? (
|
|
jbe@0
|
362 /* box crosses 180th meridian */
|
|
jbe@0
|
363 point->lon >= box->lon_min || point->lon <= box->lon_max
|
|
jbe@0
|
364 ) : (
|
|
jbe@0
|
365 /* box does not cross the 180th meridian */
|
|
jbe@0
|
366 point->lon >= box->lon_min && point->lon <= box->lon_max
|
|
jbe@0
|
367 )
|
|
jbe@0
|
368 )
|
|
jbe@0
|
369 );
|
|
jbe@0
|
370 }
|
|
jbe@0
|
371
|
|
jbe@0
|
372 /* check if two boxes overlap */
|
|
jbe@0
|
373 static bool pgl_boxes_overlap(pgl_box *box1, pgl_box *box2) {
|
|
jbe@0
|
374 return (
|
|
jbe@0
|
375 box2->lat_max >= box2->lat_min && /* ensure box2 is not empty */
|
|
jbe@0
|
376 ( box2->lat_min >= box1->lat_min || box2->lat_max >= box1->lat_min ) &&
|
|
jbe@0
|
377 ( box2->lat_min <= box1->lat_max || box2->lat_max <= box1->lat_max ) && (
|
|
jbe@0
|
378 (
|
|
jbe@0
|
379 /* check if one and only one box crosses the 180th meridian */
|
|
jbe@0
|
380 ((box1->lon_min > box1->lon_max) ? 1 : 0) ^
|
|
jbe@0
|
381 ((box2->lon_min > box2->lon_max) ? 1 : 0)
|
|
jbe@0
|
382 ) ? (
|
|
jbe@0
|
383 /* exactly one box crosses the 180th meridian */
|
|
jbe@0
|
384 box2->lon_min >= box1->lon_min || box2->lon_max >= box1->lon_min ||
|
|
jbe@0
|
385 box2->lon_min <= box1->lon_max || box2->lon_max <= box1->lon_max
|
|
jbe@0
|
386 ) : (
|
|
jbe@0
|
387 /* no box or both boxes cross the 180th meridian */
|
|
jbe@0
|
388 (
|
|
jbe@0
|
389 (box2->lon_min >= box1->lon_min || box2->lon_max >= box1->lon_min) &&
|
|
jbe@0
|
390 (box2->lon_min <= box1->lon_max || box2->lon_max <= box1->lon_max)
|
|
jbe@0
|
391 ) ||
|
|
jbe@0
|
392 /* handle W180 == E180 */
|
|
jbe@0
|
393 ( box1->lon_min == -180 && box2->lon_max == 180 ) ||
|
|
jbe@0
|
394 ( box2->lon_min == -180 && box1->lon_max == 180 )
|
|
jbe@0
|
395 )
|
|
jbe@0
|
396 )
|
|
jbe@0
|
397 );
|
|
jbe@0
|
398 }
|
|
jbe@0
|
399
|
|
jbe@0
|
400 /* check unambiguousness of east/west orientation of cluster entries and set
|
|
jbe@0
|
401 bounding circle of cluster */
|
|
jbe@0
|
402 static bool pgl_finalize_cluster(pgl_cluster *cluster) {
|
|
jbe@0
|
403 int i, j; /* i: index of entry, j: index of point in entry */
|
|
jbe@0
|
404 int npoints; /* number of points in entry */
|
|
jbe@0
|
405 int total_npoints = 0; /* total number of points in cluster */
|
|
jbe@0
|
406 pgl_point *points; /* points in entry */
|
|
jbe@0
|
407 int lon_dir; /* first point of entry west (-1) or east (+1) */
|
|
jbe@0
|
408 double lon_break = 0; /* antipodal longitude of first point in entry */
|
|
jbe@0
|
409 double lon_min, lon_max; /* covered longitude range of entry */
|
|
jbe@0
|
410 double value; /* temporary variable */
|
|
jbe@0
|
411 /* reset bounding circle center to empty circle at 0/0 coordinates */
|
|
jbe@0
|
412 cluster->bounding.center.lat = 0;
|
|
jbe@0
|
413 cluster->bounding.center.lon = 0;
|
|
jbe@0
|
414 cluster->bounding.radius = -INFINITY;
|
|
jbe@0
|
415 /* if cluster is not empty */
|
|
jbe@0
|
416 if (cluster->nentries != 0) {
|
|
jbe@0
|
417 /* iterate over all cluster entries and ensure they each cover a longitude
|
|
jbe@0
|
418 range less than 180 degrees */
|
|
jbe@0
|
419 for (i=0; i<cluster->nentries; i++) {
|
|
jbe@0
|
420 /* get properties of entry */
|
|
jbe@0
|
421 npoints = cluster->entries[i].npoints;
|
|
jbe@0
|
422 points = PGL_ENTRY_POINTS(cluster, i);
|
|
jbe@0
|
423 /* get longitude of first point of entry */
|
|
jbe@0
|
424 value = points[0].lon;
|
|
jbe@0
|
425 /* initialize lon_min and lon_max with longitude of first point */
|
|
jbe@0
|
426 lon_min = value;
|
|
jbe@0
|
427 lon_max = value;
|
|
jbe@0
|
428 /* determine east/west orientation of first point and calculate antipodal
|
|
jbe@0
|
429 longitude (Note: rounding required here) */
|
|
jbe@0
|
430 if (value < 0) { lon_dir = -1; lon_break = pgl_round(value + 180); }
|
|
jbe@0
|
431 else if (value > 0) { lon_dir = 1; lon_break = pgl_round(value - 180); }
|
|
jbe@0
|
432 else lon_dir = 0;
|
|
jbe@0
|
433 /* iterate over all other points in entry */
|
|
jbe@0
|
434 for (j=1; j<npoints; j++) {
|
|
jbe@0
|
435 /* consider longitude wrap-around */
|
|
jbe@0
|
436 value = points[j].lon;
|
|
jbe@0
|
437 if (lon_dir<0 && value>lon_break) value = pgl_round(value - 360);
|
|
jbe@0
|
438 else if (lon_dir>0 && value<lon_break) value = pgl_round(value + 360);
|
|
jbe@0
|
439 /* update lon_min and lon_max */
|
|
jbe@0
|
440 if (value < lon_min) lon_min = value;
|
|
jbe@0
|
441 else if (value > lon_max) lon_max = value;
|
|
jbe@0
|
442 /* return false if 180 degrees or more are covered */
|
|
jbe@0
|
443 if (lon_max - lon_min >= 180) return false;
|
|
jbe@0
|
444 }
|
|
jbe@0
|
445 }
|
|
jbe@0
|
446 /* iterate over all points of all entries and calculate arbitrary center
|
|
jbe@0
|
447 point for bounding circle (best if center point minimizes the radius,
|
|
jbe@0
|
448 but some error is allowed here) */
|
|
jbe@0
|
449 for (i=0; i<cluster->nentries; i++) {
|
|
jbe@0
|
450 /* get properties of entry */
|
|
jbe@0
|
451 npoints = cluster->entries[i].npoints;
|
|
jbe@0
|
452 points = PGL_ENTRY_POINTS(cluster, i);
|
|
jbe@0
|
453 /* check if first entry */
|
|
jbe@0
|
454 if (i==0) {
|
|
jbe@0
|
455 /* get longitude of first point of first entry in whole cluster */
|
|
jbe@0
|
456 value = points[0].lon;
|
|
jbe@0
|
457 /* initialize lon_min and lon_max with longitude of first point of
|
|
jbe@0
|
458 first entry in whole cluster (used to determine if whole cluster
|
|
jbe@0
|
459 covers a longitude range of 180 degrees or more) */
|
|
jbe@0
|
460 lon_min = value;
|
|
jbe@0
|
461 lon_max = value;
|
|
jbe@0
|
462 /* determine east/west orientation of first point and calculate
|
|
jbe@0
|
463 antipodal longitude (Note: rounding not necessary here) */
|
|
jbe@0
|
464 if (value < 0) { lon_dir = -1; lon_break = value + 180; }
|
|
jbe@0
|
465 else if (value > 0) { lon_dir = 1; lon_break = value - 180; }
|
|
jbe@0
|
466 else lon_dir = 0;
|
|
jbe@0
|
467 }
|
|
jbe@0
|
468 /* iterate over all points in entry */
|
|
jbe@0
|
469 for (j=0; j<npoints; j++) {
|
|
jbe@0
|
470 /* longitude wrap-around (Note: rounding not necessary here) */
|
|
jbe@0
|
471 value = points[j].lon;
|
|
jbe@0
|
472 if (lon_dir < 0 && value > lon_break) value -= 360;
|
|
jbe@0
|
473 else if (lon_dir > 0 && value < lon_break) value += 360;
|
|
jbe@0
|
474 if (value < lon_min) lon_min = value;
|
|
jbe@0
|
475 else if (value > lon_max) lon_max = value;
|
|
jbe@46
|
476 /* set bounding circle to cover whole earth if 180 degrees or more are
|
|
jbe@46
|
477 covered */
|
|
jbe@0
|
478 if (lon_max - lon_min >= 180) {
|
|
jbe@0
|
479 cluster->bounding.center.lat = 0;
|
|
jbe@0
|
480 cluster->bounding.center.lon = 0;
|
|
jbe@0
|
481 cluster->bounding.radius = INFINITY;
|
|
jbe@0
|
482 return true;
|
|
jbe@0
|
483 }
|
|
jbe@0
|
484 /* add point to bounding circle center (for average calculation) */
|
|
jbe@0
|
485 cluster->bounding.center.lat += points[j].lat;
|
|
jbe@0
|
486 cluster->bounding.center.lon += value;
|
|
jbe@0
|
487 }
|
|
jbe@0
|
488 /* count total number of points */
|
|
jbe@0
|
489 total_npoints += npoints;
|
|
jbe@0
|
490 }
|
|
jbe@0
|
491 /* determine average latitude and longitude of cluster */
|
|
jbe@0
|
492 cluster->bounding.center.lat /= total_npoints;
|
|
jbe@0
|
493 cluster->bounding.center.lon /= total_npoints;
|
|
jbe@0
|
494 /* normalize longitude of center of cluster bounding circle */
|
|
jbe@0
|
495 if (cluster->bounding.center.lon < -180) {
|
|
jbe@0
|
496 cluster->bounding.center.lon += 360;
|
|
jbe@0
|
497 }
|
|
jbe@0
|
498 else if (cluster->bounding.center.lon > 180) {
|
|
jbe@0
|
499 cluster->bounding.center.lon -= 360;
|
|
jbe@0
|
500 }
|
|
jbe@0
|
501 /* round bounding circle center (useful if it is used by other functions) */
|
|
jbe@0
|
502 cluster->bounding.center.lat = pgl_round(cluster->bounding.center.lat);
|
|
jbe@0
|
503 cluster->bounding.center.lon = pgl_round(cluster->bounding.center.lon);
|
|
jbe@0
|
504 /* calculate radius of bounding circle */
|
|
jbe@0
|
505 for (i=0; i<cluster->nentries; i++) {
|
|
jbe@0
|
506 npoints = cluster->entries[i].npoints;
|
|
jbe@0
|
507 points = PGL_ENTRY_POINTS(cluster, i);
|
|
jbe@0
|
508 for (j=0; j<npoints; j++) {
|
|
jbe@0
|
509 value = pgl_distance(
|
|
jbe@0
|
510 cluster->bounding.center.lat, cluster->bounding.center.lon,
|
|
jbe@0
|
511 points[j].lat, points[j].lon
|
|
jbe@0
|
512 );
|
|
jbe@0
|
513 if (value > cluster->bounding.radius) cluster->bounding.radius = value;
|
|
jbe@0
|
514 }
|
|
jbe@0
|
515 }
|
|
jbe@0
|
516 }
|
|
jbe@0
|
517 /* return true (east/west orientation is unambiguous) */
|
|
jbe@0
|
518 return true;
|
|
jbe@0
|
519 }
|
|
jbe@0
|
520
|
|
jbe@0
|
521 /* check if point is inside cluster */
|
|
jbe@20
|
522 /* (if point is on perimeter, then true is returned if and only if
|
|
jbe@20
|
523 strict == false) */
|
|
jbe@20
|
524 static bool pgl_point_in_cluster(
|
|
jbe@20
|
525 pgl_point *point,
|
|
jbe@20
|
526 pgl_cluster *cluster,
|
|
jbe@20
|
527 bool strict
|
|
jbe@20
|
528 ) {
|
|
jbe@0
|
529 int i, j, k; /* i: entry, j: point in entry, k: next point in entry */
|
|
jbe@0
|
530 int entrytype; /* type of entry */
|
|
jbe@0
|
531 int npoints; /* number of points in entry */
|
|
jbe@0
|
532 pgl_point *points; /* array of points in entry */
|
|
jbe@0
|
533 int lon_dir = 0; /* first vertex west (-1) or east (+1) */
|
|
jbe@0
|
534 double lon_break = 0; /* antipodal longitude of first vertex */
|
|
jbe@0
|
535 double lat0 = point->lat; /* latitude of point */
|
|
jbe@0
|
536 double lon0; /* (adjusted) longitude of point */
|
|
jbe@0
|
537 double lat1, lon1; /* latitude and (adjusted) longitude of vertex */
|
|
jbe@0
|
538 double lat2, lon2; /* latitude and (adjusted) longitude of next vertex */
|
|
jbe@0
|
539 double lon; /* longitude of intersection */
|
|
jbe@0
|
540 int counter = 0; /* counter for intersections east of point */
|
|
jbe@0
|
541 /* iterate over all entries */
|
|
jbe@0
|
542 for (i=0; i<cluster->nentries; i++) {
|
|
jbe@20
|
543 /* get type of entry */
|
|
jbe@0
|
544 entrytype = cluster->entries[i].entrytype;
|
|
jbe@20
|
545 /* skip all entries but polygons if perimeters are excluded */
|
|
jbe@20
|
546 if (strict && entrytype != PGL_ENTRY_POLYGON) continue;
|
|
jbe@20
|
547 /* get points of entry */
|
|
jbe@0
|
548 npoints = cluster->entries[i].npoints;
|
|
jbe@0
|
549 points = PGL_ENTRY_POINTS(cluster, i);
|
|
jbe@0
|
550 /* determine east/west orientation of first point of entry and calculate
|
|
jbe@0
|
551 antipodal longitude */
|
|
jbe@0
|
552 lon_break = points[0].lon;
|
|
jbe@0
|
553 if (lon_break < 0) { lon_dir = -1; lon_break += 180; }
|
|
jbe@0
|
554 else if (lon_break > 0) { lon_dir = 1; lon_break -= 180; }
|
|
jbe@0
|
555 else lon_dir = 0;
|
|
jbe@0
|
556 /* get longitude of point */
|
|
jbe@0
|
557 lon0 = point->lon;
|
|
jbe@0
|
558 /* consider longitude wrap-around for point */
|
|
jbe@0
|
559 if (lon_dir < 0 && lon0 > lon_break) lon0 = pgl_round(lon0 - 360);
|
|
jbe@0
|
560 else if (lon_dir > 0 && lon0 < lon_break) lon0 = pgl_round(lon0 + 360);
|
|
jbe@0
|
561 /* iterate over all edges and vertices */
|
|
jbe@0
|
562 for (j=0; j<npoints; j++) {
|
|
jbe@20
|
563 /* return if point is on vertex of polygon */
|
|
jbe@20
|
564 if (pgl_point_cmp(point, &(points[j])) == 0) return !strict;
|
|
jbe@0
|
565 /* calculate index of next vertex */
|
|
jbe@0
|
566 k = (j+1) % npoints;
|
|
jbe@0
|
567 /* skip last edge unless entry is (closed) outline or polygon */
|
|
jbe@0
|
568 if (
|
|
jbe@0
|
569 k == 0 &&
|
|
jbe@0
|
570 entrytype != PGL_ENTRY_OUTLINE &&
|
|
jbe@0
|
571 entrytype != PGL_ENTRY_POLYGON
|
|
jbe@0
|
572 ) continue;
|
|
jbe@16
|
573 /* use previously calculated values for lat1 and lon1 if possible */
|
|
jbe@16
|
574 if (j) {
|
|
jbe@16
|
575 lat1 = lat2;
|
|
jbe@16
|
576 lon1 = lon2;
|
|
jbe@16
|
577 } else {
|
|
jbe@16
|
578 /* otherwise get latitude and longitude values of first vertex */
|
|
jbe@16
|
579 lat1 = points[0].lat;
|
|
jbe@16
|
580 lon1 = points[0].lon;
|
|
jbe@16
|
581 /* and consider longitude wrap-around for first vertex */
|
|
jbe@16
|
582 if (lon_dir < 0 && lon1 > lon_break) lon1 = pgl_round(lon1 - 360);
|
|
jbe@16
|
583 else if (lon_dir > 0 && lon1 < lon_break) lon1 = pgl_round(lon1 + 360);
|
|
jbe@16
|
584 }
|
|
jbe@16
|
585 /* get latitude and longitude of next vertex */
|
|
jbe@0
|
586 lat2 = points[k].lat;
|
|
jbe@0
|
587 lon2 = points[k].lon;
|
|
jbe@16
|
588 /* consider longitude wrap-around for next vertex */
|
|
jbe@0
|
589 if (lon_dir < 0 && lon2 > lon_break) lon2 = pgl_round(lon2 - 360);
|
|
jbe@0
|
590 else if (lon_dir > 0 && lon2 < lon_break) lon2 = pgl_round(lon2 + 360);
|
|
jbe@20
|
591 /* return if point is on horizontal (west to east) edge of polygon */
|
|
jbe@0
|
592 if (
|
|
jbe@0
|
593 lat0 == lat1 && lat0 == lat2 &&
|
|
jbe@0
|
594 ( (lon0 >= lon1 && lon0 <= lon2) || (lon0 >= lon2 && lon0 <= lon1) )
|
|
jbe@20
|
595 ) return !strict;
|
|
jbe@0
|
596 /* check if edge crosses east/west line of point */
|
|
jbe@0
|
597 if ((lat1 < lat0 && lat2 >= lat0) || (lat2 < lat0 && lat1 >= lat0)) {
|
|
jbe@0
|
598 /* calculate longitude of intersection */
|
|
jbe@0
|
599 lon = (lon1 * (lat2-lat0) + lon2 * (lat0-lat1)) / (lat2-lat1);
|
|
jbe@20
|
600 /* return if intersection goes (approximately) through point */
|
|
jbe@20
|
601 if (pgl_round(lon) == lon0) return !strict;
|
|
jbe@0
|
602 /* count intersection if east of point and entry is polygon*/
|
|
jbe@0
|
603 if (entrytype == PGL_ENTRY_POLYGON && lon > lon0) counter++;
|
|
jbe@0
|
604 }
|
|
jbe@0
|
605 }
|
|
jbe@0
|
606 }
|
|
jbe@0
|
607 /* return true if number of intersections is odd */
|
|
jbe@0
|
608 return counter & 1;
|
|
jbe@0
|
609 }
|
|
jbe@0
|
610
|
|
jbe@20
|
611 /* check if all points of the second cluster are strictly inside the first
|
|
jbe@20
|
612 cluster */
|
|
jbe@20
|
613 static inline bool pgl_all_cluster_points_strictly_in_cluster(
|
|
jbe@16
|
614 pgl_cluster *outer, pgl_cluster *inner
|
|
jbe@16
|
615 ) {
|
|
jbe@16
|
616 int i, j; /* i: entry, j: point in entry */
|
|
jbe@16
|
617 int npoints; /* number of points in entry */
|
|
jbe@16
|
618 pgl_point *points; /* array of points in entry */
|
|
jbe@16
|
619 /* iterate over all entries of "inner" cluster */
|
|
jbe@16
|
620 for (i=0; i<inner->nentries; i++) {
|
|
jbe@16
|
621 /* get properties of entry */
|
|
jbe@16
|
622 npoints = inner->entries[i].npoints;
|
|
jbe@16
|
623 points = PGL_ENTRY_POINTS(inner, i);
|
|
jbe@16
|
624 /* iterate over all points in entry of "inner" cluster */
|
|
jbe@16
|
625 for (j=0; j<npoints; j++) {
|
|
jbe@16
|
626 /* return false if one point of inner cluster is not in outer cluster */
|
|
jbe@20
|
627 if (!pgl_point_in_cluster(points+j, outer, true)) return false;
|
|
jbe@16
|
628 }
|
|
jbe@16
|
629 }
|
|
jbe@16
|
630 /* otherwise return true */
|
|
jbe@16
|
631 return true;
|
|
jbe@16
|
632 }
|
|
jbe@16
|
633
|
|
jbe@16
|
634 /* check if any point the second cluster is inside the first cluster */
|
|
jbe@16
|
635 static inline bool pgl_any_cluster_points_in_cluster(
|
|
jbe@16
|
636 pgl_cluster *outer, pgl_cluster *inner
|
|
jbe@16
|
637 ) {
|
|
jbe@16
|
638 int i, j; /* i: entry, j: point in entry */
|
|
jbe@16
|
639 int npoints; /* number of points in entry */
|
|
jbe@16
|
640 pgl_point *points; /* array of points in entry */
|
|
jbe@16
|
641 /* iterate over all entries of "inner" cluster */
|
|
jbe@16
|
642 for (i=0; i<inner->nentries; i++) {
|
|
jbe@16
|
643 /* get properties of entry */
|
|
jbe@16
|
644 npoints = inner->entries[i].npoints;
|
|
jbe@16
|
645 points = PGL_ENTRY_POINTS(inner, i);
|
|
jbe@16
|
646 /* iterate over all points in entry of "inner" cluster */
|
|
jbe@16
|
647 for (j=0; j<npoints; j++) {
|
|
jbe@16
|
648 /* return true if one point of inner cluster is in outer cluster */
|
|
jbe@20
|
649 if (pgl_point_in_cluster(points+j, outer, false)) return true;
|
|
jbe@16
|
650 }
|
|
jbe@16
|
651 }
|
|
jbe@16
|
652 /* otherwise return false */
|
|
jbe@16
|
653 return false;
|
|
jbe@16
|
654 }
|
|
jbe@16
|
655
|
|
jbe@20
|
656 /* check if line segment strictly crosses line (not just touching) */
|
|
jbe@20
|
657 static inline bool pgl_lseg_crosses_line(
|
|
jbe@16
|
658 double seg_x1, double seg_y1, double seg_x2, double seg_y2,
|
|
jbe@20
|
659 double line_x1, double line_y1, double line_x2, double line_y2
|
|
jbe@16
|
660 ) {
|
|
jbe@20
|
661 return (
|
|
jbe@20
|
662 (
|
|
jbe@20
|
663 (seg_x1-line_x1) * (line_y2-line_y1) -
|
|
jbe@20
|
664 (seg_y1-line_y1) * (line_x2-line_x1)
|
|
jbe@20
|
665 ) * (
|
|
jbe@20
|
666 (seg_x2-line_x1) * (line_y2-line_y1) -
|
|
jbe@20
|
667 (seg_y2-line_y1) * (line_x2-line_x1)
|
|
jbe@20
|
668 )
|
|
jbe@20
|
669 ) < 0;
|
|
jbe@16
|
670 }
|
|
jbe@16
|
671
|
|
jbe@20
|
672 /* check if paths and outlines of two clusters strictly overlap (not just
|
|
jbe@20
|
673 touching) */
|
|
jbe@16
|
674 static bool pgl_outlines_overlap(
|
|
jbe@20
|
675 pgl_cluster *cluster1, pgl_cluster *cluster2
|
|
jbe@16
|
676 ) {
|
|
jbe@16
|
677 int i1, j1, k1; /* i: entry, j: point in entry, k: next point in entry */
|
|
jbe@16
|
678 int i2, j2, k2;
|
|
jbe@16
|
679 int entrytype1, entrytype2; /* type of entry */
|
|
jbe@16
|
680 int npoints1, npoints2; /* number of points in entry */
|
|
jbe@16
|
681 pgl_point *points1; /* array of points in entry of cluster1 */
|
|
jbe@16
|
682 pgl_point *points2; /* array of points in entry of cluster2 */
|
|
jbe@16
|
683 int lon_dir1, lon_dir2; /* first vertex west (-1) or east (+1) */
|
|
jbe@16
|
684 double lon_break1, lon_break2; /* antipodal longitude of first vertex */
|
|
jbe@16
|
685 double lat11, lon11; /* latitude and (adjusted) longitude of vertex */
|
|
jbe@16
|
686 double lat12, lon12; /* latitude and (adjusted) longitude of next vertex */
|
|
jbe@16
|
687 double lat21, lon21; /* latitude and (adjusted) longitudes for cluster2 */
|
|
jbe@16
|
688 double lat22, lon22;
|
|
jbe@16
|
689 double wrapvalue; /* temporary helper value to adjust wrap-around */
|
|
jbe@16
|
690 /* iterate over all entries of cluster1 */
|
|
jbe@16
|
691 for (i1=0; i1<cluster1->nentries; i1++) {
|
|
jbe@16
|
692 /* get properties of entry in cluster1 and skip points */
|
|
jbe@16
|
693 npoints1 = cluster1->entries[i1].npoints;
|
|
jbe@16
|
694 if (npoints1 < 2) continue;
|
|
jbe@16
|
695 entrytype1 = cluster1->entries[i1].entrytype;
|
|
jbe@16
|
696 points1 = PGL_ENTRY_POINTS(cluster1, i1);
|
|
jbe@16
|
697 /* determine east/west orientation of first point and calculate antipodal
|
|
jbe@16
|
698 longitude */
|
|
jbe@16
|
699 lon_break1 = points1[0].lon;
|
|
jbe@16
|
700 if (lon_break1 < 0) {
|
|
jbe@16
|
701 lon_dir1 = -1;
|
|
jbe@16
|
702 lon_break1 = pgl_round(lon_break1 + 180);
|
|
jbe@16
|
703 } else if (lon_break1 > 0) {
|
|
jbe@16
|
704 lon_dir1 = 1;
|
|
jbe@16
|
705 lon_break1 = pgl_round(lon_break1 - 180);
|
|
jbe@16
|
706 } else lon_dir1 = 0;
|
|
jbe@16
|
707 /* iterate over all edges and vertices in cluster1 */
|
|
jbe@16
|
708 for (j1=0; j1<npoints1; j1++) {
|
|
jbe@16
|
709 /* calculate index of next vertex */
|
|
jbe@16
|
710 k1 = (j1+1) % npoints1;
|
|
jbe@16
|
711 /* skip last edge unless entry is (closed) outline or polygon */
|
|
jbe@16
|
712 if (
|
|
jbe@16
|
713 k1 == 0 &&
|
|
jbe@16
|
714 entrytype1 != PGL_ENTRY_OUTLINE &&
|
|
jbe@16
|
715 entrytype1 != PGL_ENTRY_POLYGON
|
|
jbe@16
|
716 ) continue;
|
|
jbe@16
|
717 /* use previously calculated values for lat1 and lon1 if possible */
|
|
jbe@16
|
718 if (j1) {
|
|
jbe@16
|
719 lat11 = lat12;
|
|
jbe@16
|
720 lon11 = lon12;
|
|
jbe@16
|
721 } else {
|
|
jbe@16
|
722 /* otherwise get latitude and longitude values of first vertex */
|
|
jbe@16
|
723 lat11 = points1[0].lat;
|
|
jbe@16
|
724 lon11 = points1[0].lon;
|
|
jbe@16
|
725 /* and consider longitude wrap-around for first vertex */
|
|
jbe@16
|
726 if (lon_dir1<0 && lon11>lon_break1) lon11 = pgl_round(lon11-360);
|
|
jbe@16
|
727 else if (lon_dir1>0 && lon11<lon_break1) lon11 = pgl_round(lon11+360);
|
|
jbe@16
|
728 }
|
|
jbe@16
|
729 /* get latitude and longitude of next vertex */
|
|
jbe@16
|
730 lat12 = points1[k1].lat;
|
|
jbe@16
|
731 lon12 = points1[k1].lon;
|
|
jbe@16
|
732 /* consider longitude wrap-around for next vertex */
|
|
jbe@16
|
733 if (lon_dir1<0 && lon12>lon_break1) lon12 = pgl_round(lon12-360);
|
|
jbe@16
|
734 else if (lon_dir1>0 && lon12<lon_break1) lon12 = pgl_round(lon12+360);
|
|
jbe@16
|
735 /* skip degenerated edges */
|
|
jbe@16
|
736 if (lat11 == lat12 && lon11 == lon12) continue;
|
|
jbe@16
|
737 /* iterate over all entries of cluster2 */
|
|
jbe@16
|
738 for (i2=0; i2<cluster2->nentries; i2++) {
|
|
jbe@16
|
739 /* get points and number of points of entry in cluster2 */
|
|
jbe@16
|
740 npoints2 = cluster2->entries[i2].npoints;
|
|
jbe@16
|
741 if (npoints2 < 2) continue;
|
|
jbe@16
|
742 entrytype2 = cluster2->entries[i2].entrytype;
|
|
jbe@16
|
743 points2 = PGL_ENTRY_POINTS(cluster2, i2);
|
|
jbe@16
|
744 /* determine east/west orientation of first point and calculate antipodal
|
|
jbe@16
|
745 longitude */
|
|
jbe@16
|
746 lon_break2 = points2[0].lon;
|
|
jbe@16
|
747 if (lon_break2 < 0) {
|
|
jbe@16
|
748 lon_dir2 = -1;
|
|
jbe@16
|
749 lon_break2 = pgl_round(lon_break2 + 180);
|
|
jbe@16
|
750 } else if (lon_break2 > 0) {
|
|
jbe@16
|
751 lon_dir2 = 1;
|
|
jbe@16
|
752 lon_break2 = pgl_round(lon_break2 - 180);
|
|
jbe@16
|
753 } else lon_dir2 = 0;
|
|
jbe@16
|
754 /* iterate over all edges and vertices in cluster2 */
|
|
jbe@16
|
755 for (j2=0; j2<npoints2; j2++) {
|
|
jbe@16
|
756 /* calculate index of next vertex */
|
|
jbe@16
|
757 k2 = (j2+1) % npoints2;
|
|
jbe@16
|
758 /* skip last edge unless entry is (closed) outline or polygon */
|
|
jbe@16
|
759 if (
|
|
jbe@16
|
760 k2 == 0 &&
|
|
jbe@16
|
761 entrytype2 != PGL_ENTRY_OUTLINE &&
|
|
jbe@16
|
762 entrytype2 != PGL_ENTRY_POLYGON
|
|
jbe@16
|
763 ) continue;
|
|
jbe@16
|
764 /* use previously calculated values for lat1 and lon1 if possible */
|
|
jbe@16
|
765 if (j2) {
|
|
jbe@16
|
766 lat21 = lat22;
|
|
jbe@16
|
767 lon21 = lon22;
|
|
jbe@16
|
768 } else {
|
|
jbe@16
|
769 /* otherwise get latitude and longitude values of first vertex */
|
|
jbe@16
|
770 lat21 = points2[0].lat;
|
|
jbe@16
|
771 lon21 = points2[0].lon;
|
|
jbe@16
|
772 /* and consider longitude wrap-around for first vertex */
|
|
jbe@16
|
773 if (lon_dir2<0 && lon21>lon_break2) lon21 = pgl_round(lon21-360);
|
|
jbe@16
|
774 else if (lon_dir2>0 && lon21<lon_break2) lon21 = pgl_round(lon21+360);
|
|
jbe@16
|
775 }
|
|
jbe@16
|
776 /* get latitude and longitude of next vertex */
|
|
jbe@16
|
777 lat22 = points2[k2].lat;
|
|
jbe@16
|
778 lon22 = points2[k2].lon;
|
|
jbe@16
|
779 /* consider longitude wrap-around for next vertex */
|
|
jbe@16
|
780 if (lon_dir2<0 && lon22>lon_break2) lon22 = pgl_round(lon22-360);
|
|
jbe@16
|
781 else if (lon_dir2>0 && lon22<lon_break2) lon22 = pgl_round(lon22+360);
|
|
jbe@16
|
782 /* skip degenerated edges */
|
|
jbe@16
|
783 if (lat21 == lat22 && lon21 == lon22) continue;
|
|
jbe@16
|
784 /* perform another wrap-around where necessary */
|
|
jbe@16
|
785 /* TODO: improve performance of whole wrap-around mechanism */
|
|
jbe@16
|
786 wrapvalue = (lon21 + lon22) - (lon11 + lon12);
|
|
jbe@16
|
787 if (wrapvalue > 360) {
|
|
jbe@16
|
788 lon21 = pgl_round(lon21 - 360);
|
|
jbe@16
|
789 lon22 = pgl_round(lon22 - 360);
|
|
jbe@16
|
790 } else if (wrapvalue < -360) {
|
|
jbe@16
|
791 lon21 = pgl_round(lon21 + 360);
|
|
jbe@16
|
792 lon22 = pgl_round(lon22 + 360);
|
|
jbe@16
|
793 }
|
|
jbe@16
|
794 /* return true if segments overlap */
|
|
jbe@16
|
795 if (
|
|
jbe@16
|
796 pgl_lseg_crosses_line(
|
|
jbe@16
|
797 lat11, lon11, lat12, lon12,
|
|
jbe@20
|
798 lat21, lon21, lat22, lon22
|
|
jbe@16
|
799 ) && pgl_lseg_crosses_line(
|
|
jbe@16
|
800 lat21, lon21, lat22, lon22,
|
|
jbe@20
|
801 lat11, lon11, lat12, lon12
|
|
jbe@16
|
802 )
|
|
jbe@16
|
803 ) {
|
|
jbe@16
|
804 return true;
|
|
jbe@16
|
805 }
|
|
jbe@16
|
806 }
|
|
jbe@16
|
807 }
|
|
jbe@16
|
808 }
|
|
jbe@16
|
809 }
|
|
jbe@16
|
810 /* otherwise return false */
|
|
jbe@16
|
811 return false;
|
|
jbe@16
|
812 }
|
|
jbe@16
|
813
|
|
jbe@16
|
814 /* check if second cluster is completely contained in first cluster */
|
|
jbe@16
|
815 static bool pgl_cluster_in_cluster(pgl_cluster *outer, pgl_cluster *inner) {
|
|
jbe@20
|
816 if (!pgl_all_cluster_points_strictly_in_cluster(outer, inner)) return false;
|
|
jbe@20
|
817 if (pgl_any_cluster_points_in_cluster(inner, outer)) return false;
|
|
jbe@20
|
818 if (pgl_outlines_overlap(outer, inner)) return false;
|
|
jbe@16
|
819 return true;
|
|
jbe@16
|
820 }
|
|
jbe@16
|
821
|
|
jbe@16
|
822 /* check if two clusters overlap */
|
|
jbe@16
|
823 static bool pgl_clusters_overlap(
|
|
jbe@16
|
824 pgl_cluster *cluster1, pgl_cluster *cluster2
|
|
jbe@16
|
825 ) {
|
|
jbe@16
|
826 if (pgl_any_cluster_points_in_cluster(cluster1, cluster2)) return true;
|
|
jbe@16
|
827 if (pgl_any_cluster_points_in_cluster(cluster2, cluster1)) return true;
|
|
jbe@20
|
828 if (pgl_outlines_overlap(cluster1, cluster2)) return true;
|
|
jbe@16
|
829 return false;
|
|
jbe@16
|
830 }
|
|
jbe@16
|
831
|
|
jbe@0
|
832 /* calculate (approximate) distance between point and cluster */
|
|
jbe@0
|
833 static double pgl_point_cluster_distance(pgl_point *point, pgl_cluster *cluster) {
|
|
jbe@24
|
834 double comp; /* square of compression of meridians */
|
|
jbe@0
|
835 int i, j, k; /* i: entry, j: point in entry, k: next point in entry */
|
|
jbe@0
|
836 int entrytype; /* type of entry */
|
|
jbe@0
|
837 int npoints; /* number of points in entry */
|
|
jbe@0
|
838 pgl_point *points; /* array of points in entry */
|
|
jbe@0
|
839 int lon_dir = 0; /* first vertex west (-1) or east (+1) */
|
|
jbe@0
|
840 double lon_break = 0; /* antipodal longitude of first vertex */
|
|
jbe@0
|
841 double lon_min = 0; /* minimum (adjusted) longitude of entry vertices */
|
|
jbe@0
|
842 double lon_max = 0; /* maximum (adjusted) longitude of entry vertices */
|
|
jbe@0
|
843 double lat0 = point->lat; /* latitude of point */
|
|
jbe@0
|
844 double lon0; /* (adjusted) longitude of point */
|
|
jbe@0
|
845 double lat1, lon1; /* latitude and (adjusted) longitude of vertex */
|
|
jbe@0
|
846 double lat2, lon2; /* latitude and (adjusted) longitude of next vertex */
|
|
jbe@0
|
847 double s; /* scalar for vector calculations */
|
|
jbe@0
|
848 double dist; /* distance calculated in one step */
|
|
jbe@0
|
849 double min_dist = INFINITY; /* minimum distance */
|
|
jbe@0
|
850 /* distance is zero if point is contained in cluster */
|
|
jbe@20
|
851 if (pgl_point_in_cluster(point, cluster, false)) return 0;
|
|
jbe@30
|
852 /* calculate approximate square compression of meridians */
|
|
jbe@24
|
853 comp = cos((lat0 / 180.0) * M_PI);
|
|
jbe@24
|
854 comp *= comp;
|
|
jbe@30
|
855 /* calculate exact square compression of meridians */
|
|
jbe@30
|
856 comp *= (
|
|
jbe@30
|
857 (1.0 - PGL_EPS2 * (1.0-comp)) *
|
|
jbe@30
|
858 (1.0 - PGL_EPS2 * (1.0-comp)) /
|
|
jbe@30
|
859 (PGL_SUBEPS2 * PGL_SUBEPS2)
|
|
jbe@30
|
860 );
|
|
jbe@0
|
861 /* iterate over all entries */
|
|
jbe@0
|
862 for (i=0; i<cluster->nentries; i++) {
|
|
jbe@0
|
863 /* get properties of entry */
|
|
jbe@0
|
864 entrytype = cluster->entries[i].entrytype;
|
|
jbe@0
|
865 npoints = cluster->entries[i].npoints;
|
|
jbe@0
|
866 points = PGL_ENTRY_POINTS(cluster, i);
|
|
jbe@0
|
867 /* determine east/west orientation of first point of entry and calculate
|
|
jbe@0
|
868 antipodal longitude */
|
|
jbe@0
|
869 lon_break = points[0].lon;
|
|
jbe@0
|
870 if (lon_break < 0) { lon_dir = -1; lon_break += 180; }
|
|
jbe@0
|
871 else if (lon_break > 0) { lon_dir = 1; lon_break -= 180; }
|
|
jbe@0
|
872 else lon_dir = 0;
|
|
jbe@0
|
873 /* determine covered longitude range */
|
|
jbe@0
|
874 for (j=0; j<npoints; j++) {
|
|
jbe@0
|
875 /* get longitude of vertex */
|
|
jbe@0
|
876 lon1 = points[j].lon;
|
|
jbe@0
|
877 /* adjust longitude to fix potential wrap-around */
|
|
jbe@0
|
878 if (lon_dir < 0 && lon1 > lon_break) lon1 -= 360;
|
|
jbe@0
|
879 else if (lon_dir > 0 && lon1 < lon_break) lon1 += 360;
|
|
jbe@0
|
880 /* update minimum and maximum longitude of polygon */
|
|
jbe@0
|
881 if (j == 0 || lon1 < lon_min) lon_min = lon1;
|
|
jbe@0
|
882 if (j == 0 || lon1 > lon_max) lon_max = lon1;
|
|
jbe@0
|
883 }
|
|
jbe@0
|
884 /* adjust longitude wrap-around according to full longitude range */
|
|
jbe@0
|
885 lon_break = (lon_max + lon_min) / 2;
|
|
jbe@0
|
886 if (lon_break < 0) { lon_dir = -1; lon_break += 180; }
|
|
jbe@0
|
887 else if (lon_break > 0) { lon_dir = 1; lon_break -= 180; }
|
|
jbe@0
|
888 /* get longitude of point */
|
|
jbe@0
|
889 lon0 = point->lon;
|
|
jbe@0
|
890 /* consider longitude wrap-around for point */
|
|
jbe@0
|
891 if (lon_dir < 0 && lon0 > lon_break) lon0 -= 360;
|
|
jbe@0
|
892 else if (lon_dir > 0 && lon0 < lon_break) lon0 += 360;
|
|
jbe@0
|
893 /* iterate over all edges and vertices */
|
|
jbe@0
|
894 for (j=0; j<npoints; j++) {
|
|
jbe@16
|
895 /* use previously calculated values for lat1 and lon1 if possible */
|
|
jbe@16
|
896 if (j) {
|
|
jbe@16
|
897 lat1 = lat2;
|
|
jbe@16
|
898 lon1 = lon2;
|
|
jbe@16
|
899 } else {
|
|
jbe@16
|
900 /* otherwise get latitude and longitude values of first vertex */
|
|
jbe@16
|
901 lat1 = points[0].lat;
|
|
jbe@16
|
902 lon1 = points[0].lon;
|
|
jbe@16
|
903 /* and consider longitude wrap-around for first vertex */
|
|
jbe@16
|
904 if (lon_dir < 0 && lon1 > lon_break) lon1 -= 360;
|
|
jbe@16
|
905 else if (lon_dir > 0 && lon1 < lon_break) lon1 += 360;
|
|
jbe@16
|
906 }
|
|
jbe@0
|
907 /* calculate distance to vertex */
|
|
jbe@0
|
908 dist = pgl_distance(lat0, lon0, lat1, lon1);
|
|
jbe@0
|
909 /* store calculated distance if smallest */
|
|
jbe@0
|
910 if (dist < min_dist) min_dist = dist;
|
|
jbe@0
|
911 /* calculate index of next vertex */
|
|
jbe@0
|
912 k = (j+1) % npoints;
|
|
jbe@0
|
913 /* skip last edge unless entry is (closed) outline or polygon */
|
|
jbe@0
|
914 if (
|
|
jbe@0
|
915 k == 0 &&
|
|
jbe@0
|
916 entrytype != PGL_ENTRY_OUTLINE &&
|
|
jbe@0
|
917 entrytype != PGL_ENTRY_POLYGON
|
|
jbe@0
|
918 ) continue;
|
|
jbe@16
|
919 /* get latitude and longitude of next vertex */
|
|
jbe@0
|
920 lat2 = points[k].lat;
|
|
jbe@0
|
921 lon2 = points[k].lon;
|
|
jbe@16
|
922 /* consider longitude wrap-around for next vertex */
|
|
jbe@0
|
923 if (lon_dir < 0 && lon2 > lon_break) lon2 -= 360;
|
|
jbe@0
|
924 else if (lon_dir > 0 && lon2 < lon_break) lon2 += 360;
|
|
jbe@0
|
925 /* go to next vertex and edge if edge is degenerated */
|
|
jbe@0
|
926 if (lat1 == lat2 && lon1 == lon2) continue;
|
|
jbe@0
|
927 /* otherwise test if point can be projected onto edge of polygon */
|
|
jbe@0
|
928 s = (
|
|
jbe@24
|
929 ((lat0-lat1) * (lat2-lat1) + comp * (lon0-lon1) * (lon2-lon1)) /
|
|
jbe@24
|
930 ((lat2-lat1) * (lat2-lat1) + comp * (lon2-lon1) * (lon2-lon1))
|
|
jbe@0
|
931 );
|
|
jbe@0
|
932 /* go to next vertex and edge if point cannot be projected */
|
|
jbe@0
|
933 if (!(s > 0 && s < 1)) continue;
|
|
jbe@0
|
934 /* calculate distance from original point to projected point */
|
|
jbe@0
|
935 dist = pgl_distance(
|
|
jbe@0
|
936 lat0, lon0,
|
|
jbe@0
|
937 lat1 + s * (lat2-lat1),
|
|
jbe@0
|
938 lon1 + s * (lon2-lon1)
|
|
jbe@0
|
939 );
|
|
jbe@0
|
940 /* store calculated distance if smallest */
|
|
jbe@0
|
941 if (dist < min_dist) min_dist = dist;
|
|
jbe@0
|
942 }
|
|
jbe@0
|
943 }
|
|
jbe@0
|
944 /* return minimum distance */
|
|
jbe@0
|
945 return min_dist;
|
|
jbe@0
|
946 }
|
|
jbe@0
|
947
|
|
jbe@16
|
948 /* calculate (approximate) distance between two clusters */
|
|
jbe@16
|
949 static double pgl_cluster_distance(pgl_cluster *cluster1, pgl_cluster *cluster2) {
|
|
jbe@16
|
950 int i, j; /* i: entry, j: point in entry */
|
|
jbe@16
|
951 int npoints; /* number of points in entry */
|
|
jbe@16
|
952 pgl_point *points; /* array of points in entry */
|
|
jbe@16
|
953 double dist; /* distance calculated in one step */
|
|
jbe@16
|
954 double min_dist = INFINITY; /* minimum distance */
|
|
jbe@16
|
955 /* consider distance from each point in one cluster to the whole other */
|
|
jbe@16
|
956 for (i=0; i<cluster1->nentries; i++) {
|
|
jbe@16
|
957 npoints = cluster1->entries[i].npoints;
|
|
jbe@16
|
958 points = PGL_ENTRY_POINTS(cluster1, i);
|
|
jbe@16
|
959 for (j=0; j<npoints; j++) {
|
|
jbe@16
|
960 dist = pgl_point_cluster_distance(points+j, cluster2);
|
|
jbe@16
|
961 if (dist == 0) return dist;
|
|
jbe@16
|
962 if (dist < min_dist) min_dist = dist;
|
|
jbe@16
|
963 }
|
|
jbe@16
|
964 }
|
|
jbe@16
|
965 /* consider distance from each point in other cluster to the first cluster */
|
|
jbe@16
|
966 for (i=0; i<cluster2->nentries; i++) {
|
|
jbe@16
|
967 npoints = cluster2->entries[i].npoints;
|
|
jbe@16
|
968 points = PGL_ENTRY_POINTS(cluster2, i);
|
|
jbe@16
|
969 for (j=0; j<npoints; j++) {
|
|
jbe@16
|
970 dist = pgl_point_cluster_distance(points+j, cluster1);
|
|
jbe@16
|
971 if (dist == 0) return dist;
|
|
jbe@16
|
972 if (dist < min_dist) min_dist = dist;
|
|
jbe@16
|
973 }
|
|
jbe@16
|
974 }
|
|
jbe@16
|
975 return min_dist;
|
|
jbe@16
|
976 }
|
|
jbe@16
|
977
|
|
jbe@0
|
978 /* estimator function for distance between box and point */
|
|
jbe@16
|
979 /* always returns a smaller value than actually correct or zero */
|
|
jbe@0
|
980 static double pgl_estimate_point_box_distance(pgl_point *point, pgl_box *box) {
|
|
jbe@16
|
981 double dlon; /* longitude range of box (delta longitude) */
|
|
jbe@16
|
982 double distance; /* return value */
|
|
jbe@16
|
983 /* return infinity if box is empty */
|
|
jbe@0
|
984 if (box->lat_min > box->lat_max) return INFINITY;
|
|
jbe@16
|
985 /* return zero if point is inside box */
|
|
jbe@0
|
986 if (pgl_point_in_box(point, box)) return 0;
|
|
jbe@0
|
987 /* calculate delta longitude */
|
|
jbe@0
|
988 dlon = box->lon_max - box->lon_min;
|
|
jbe@0
|
989 if (dlon < 0) dlon += 360; /* 180th meridian crossed */
|
|
jbe@16
|
990 /* if delta longitude is greater than 150 degrees, perform safe fall-back */
|
|
jbe@16
|
991 if (dlon > 150) return 0;
|
|
jbe@16
|
992 /* calculate lower limit for distance (formula below requires dlon <= 150) */
|
|
jbe@16
|
993 /* TODO: provide better estimation function to improve performance */
|
|
jbe@16
|
994 distance = (
|
|
jbe@46
|
995 (1.0-PGL_SPHEROID_F) * /* safety margin due to flattening and approx. */
|
|
jbe@46
|
996 pgl_distance(
|
|
jbe@16
|
997 point->lat,
|
|
jbe@16
|
998 point->lon,
|
|
jbe@16
|
999 (box->lat_min + box->lat_max) / 2,
|
|
jbe@16
|
1000 box->lon_min + dlon/2
|
|
jbe@16
|
1001 )
|
|
jbe@46
|
1002 ) - pgl_distance(
|
|
jbe@46
|
1003 box->lat_min, box->lon_min,
|
|
jbe@46
|
1004 box->lat_max, box->lon_max
|
|
jbe@16
|
1005 );
|
|
jbe@16
|
1006 /* truncate negative results to zero */
|
|
jbe@16
|
1007 if (distance <= 0) distance = 0;
|
|
jbe@16
|
1008 /* return result */
|
|
jbe@16
|
1009 return distance;
|
|
jbe@0
|
1010 }
|
|
jbe@0
|
1011
|
|
jbe@0
|
1012
|
|
jbe@45
|
1013 /*------------------------------------------------------------*
|
|
jbe@45
|
1014 * Functions using numerical integration (Monte Carlo like) *
|
|
jbe@45
|
1015 *------------------------------------------------------------*/
|
|
jbe@42
|
1016
|
|
jbe@42
|
1017 /* half of (spherical) earth's surface area */
|
|
jbe@42
|
1018 #define PGL_HALF_SURFACE (PGL_RADIUS * PGL_DIAMETER * M_PI)
|
|
jbe@42
|
1019
|
|
jbe@42
|
1020 /* golden angle in radians */
|
|
jbe@42
|
1021 #define PGL_GOLDEN_ANGLE (M_PI * (sqrt(5) - 1.0))
|
|
jbe@42
|
1022
|
|
jbe@42
|
1023 /* create a list of sample points covering a bounding circle
|
|
jbe@42
|
1024 and return covered area */
|
|
jbe@42
|
1025 static double pgl_sample_points(
|
|
jbe@42
|
1026 pgl_point *center, /* center of bounding circle */
|
|
jbe@42
|
1027 double radius, /* radius of bounding circle */
|
|
jbe@46
|
1028 int samples, /* number of sample points (MUST be positive!) */
|
|
jbe@42
|
1029 pgl_point *result /* pointer to result array */
|
|
jbe@42
|
1030 ) {
|
|
jbe@42
|
1031 double double_share = 2.0; /* double of covered share of earth's surface */
|
|
jbe@42
|
1032 double double_share_div_samples; /* double_share divided by sample count */
|
|
jbe@42
|
1033 int i;
|
|
jbe@42
|
1034 double t; /* parameter of spiral laid on (spherical) earth's surface */
|
|
jbe@42
|
1035 double x, y, z; /* normalized coordinates of point on non-rotated spiral */
|
|
jbe@42
|
1036 double sin_phi; /* sine of sph. coordinate of point of non-rotated spiral */
|
|
jbe@42
|
1037 double lambda; /* other sph. coordinate of point of non-rotated spiral */
|
|
jbe@42
|
1038 double rot = (0.5 - center->lat / 180.0) * M_PI; /* needed rot. (in rad) */
|
|
jbe@42
|
1039 double cos_rot = cos(rot); /* cosine of rotation by latitude */
|
|
jbe@42
|
1040 double sin_rot = sin(rot); /* sine of rotation by latitude */
|
|
jbe@42
|
1041 double x_rot, z_rot; /* normalized coordinates of point on rotated spiral */
|
|
jbe@42
|
1042 double center_lon = center->lon; /* second rotation in degree */
|
|
jbe@42
|
1043 /* add safety margin to bounding circle because of spherical approximation */
|
|
jbe@42
|
1044 radius *= PGL_SPHEROID_A / PGL_RADIUS;
|
|
jbe@42
|
1045 /* if whole earth is covered, use initialized value, otherwise calculate
|
|
jbe@42
|
1046 share of covered area (multiplied by 2) */
|
|
jbe@42
|
1047 if (radius < PGL_MAXDIST) double_share = 1.0 - cos(radius / PGL_RADIUS);
|
|
jbe@42
|
1048 /* divide double_share by sample count for later calculations */
|
|
jbe@42
|
1049 double_share_div_samples = double_share / samples;
|
|
jbe@42
|
1050 /* generate sample points */
|
|
jbe@42
|
1051 for (i=0; i<samples; i++) {
|
|
jbe@42
|
1052 /* use an offset of 1/2 to avoid too dense clustering at spiral center */
|
|
jbe@42
|
1053 t = 0.5 + i;
|
|
jbe@42
|
1054 /* calculate normalized coordinates of point on non-rotated spiral */
|
|
jbe@42
|
1055 z = 1.0 - double_share_div_samples * t;
|
|
jbe@42
|
1056 sin_phi = sqrt(1.0 - z*z);
|
|
jbe@42
|
1057 lambda = t * PGL_GOLDEN_ANGLE;
|
|
jbe@42
|
1058 x = sin_phi * cos(lambda);
|
|
jbe@42
|
1059 y = sin_phi * sin(lambda);
|
|
jbe@42
|
1060 /* rotate spiral by latitude value of bounding circle */
|
|
jbe@42
|
1061 x_rot = cos_rot * x + sin_rot * z;
|
|
jbe@42
|
1062 z_rot = cos_rot * z - sin_rot * x;
|
|
jbe@42
|
1063 /* set resulting sample point in result array */
|
|
jbe@42
|
1064 /* (while performing second rotation by bounding circle longitude) */
|
|
jbe@42
|
1065 result[i].lat = 180.0 * (atan(z_rot / fabs(x_rot)) / M_PI);
|
|
jbe@42
|
1066 result[i].lon = center_lon + 180.0 * (atan2(y, x_rot) / M_PI);
|
|
jbe@42
|
1067 }
|
|
jbe@42
|
1068 /* return covered area */
|
|
jbe@42
|
1069 return PGL_HALF_SURFACE * double_share;
|
|
jbe@42
|
1070 }
|
|
jbe@42
|
1071
|
|
jbe@42
|
1072 /* fair distance between point and cluster (see README file for explanation) */
|
|
jbe@46
|
1073 /* NOTE: sample count passed as third argument MUST be positive */
|
|
jbe@42
|
1074 static double pgl_fair_distance(
|
|
jbe@42
|
1075 pgl_point *point, pgl_cluster *cluster, int samples
|
|
jbe@42
|
1076 ) {
|
|
jbe@42
|
1077 double distance; /* shortest distance from point to cluster */
|
|
jbe@46
|
1078 pgl_point *points; /* sample points for numerical integration */
|
|
jbe@42
|
1079 double area; /* area covered by sample points */
|
|
jbe@42
|
1080 int i;
|
|
jbe@46
|
1081 int inner = 0; /* number of sample points within cluster */
|
|
jbe@46
|
1082 int outer = 0; /* number of sample points outside cluster but
|
|
jbe@46
|
1083 within cluster enlarged by distance */
|
|
jbe@46
|
1084 double result;
|
|
jbe@42
|
1085 /* calculate shortest distance from point to cluster */
|
|
jbe@42
|
1086 distance = pgl_point_cluster_distance(point, cluster);
|
|
jbe@42
|
1087 /* if cluster consists of a single point or has no bounding circle with
|
|
jbe@42
|
1088 positive radius, simply return distance */
|
|
jbe@42
|
1089 if (
|
|
jbe@42
|
1090 (cluster->nentries==1 && cluster->entries[0].entrytype==PGL_ENTRY_POINT) ||
|
|
jbe@42
|
1091 !(cluster->bounding.radius > 0)
|
|
jbe@42
|
1092 ) return distance;
|
|
jbe@42
|
1093 /* if cluster consists of two points which are twice as far apart, return
|
|
jbe@42
|
1094 distance between point and cluster multiplied by square root of two */
|
|
jbe@42
|
1095 if (
|
|
jbe@42
|
1096 cluster->nentries == 2 &&
|
|
jbe@42
|
1097 cluster->entries[0].entrytype == PGL_ENTRY_POINT &&
|
|
jbe@42
|
1098 cluster->entries[1].entrytype == PGL_ENTRY_POINT &&
|
|
jbe@42
|
1099 pgl_distance(
|
|
jbe@42
|
1100 PGL_ENTRY_POINTS(cluster, 0)[0].lat,
|
|
jbe@42
|
1101 PGL_ENTRY_POINTS(cluster, 0)[0].lon,
|
|
jbe@42
|
1102 PGL_ENTRY_POINTS(cluster, 1)[0].lat,
|
|
jbe@42
|
1103 PGL_ENTRY_POINTS(cluster, 1)[0].lon
|
|
jbe@42
|
1104 ) >= 2.0 * distance
|
|
jbe@42
|
1105 ) {
|
|
jbe@42
|
1106 return distance * M_SQRT2;
|
|
jbe@42
|
1107 }
|
|
jbe@46
|
1108 /* otherwise create sample points for numerical integration and determine
|
|
jbe@46
|
1109 area covered by sample points */
|
|
jbe@42
|
1110 points = palloc(samples * sizeof(pgl_point));
|
|
jbe@42
|
1111 area = pgl_sample_points(
|
|
jbe@42
|
1112 &cluster->bounding.center,
|
|
jbe@42
|
1113 cluster->bounding.radius + distance, /* pad bounding circle by distance */
|
|
jbe@42
|
1114 samples,
|
|
jbe@42
|
1115 points
|
|
jbe@42
|
1116 );
|
|
jbe@46
|
1117 /* perform numerical integration */
|
|
jbe@42
|
1118 if (distance > 0) {
|
|
jbe@46
|
1119 /* point (that was passed as argument) is outside cluster */
|
|
jbe@42
|
1120 for (i=0; i<samples; i++) {
|
|
jbe@46
|
1121 /* count sample points within cluster */
|
|
jbe@46
|
1122 if (pgl_point_in_cluster(points+i, cluster, true)) inner++;
|
|
jbe@46
|
1123 /* count sample points outside of cluster but within cluster enlarged by
|
|
jbe@46
|
1124 distance between point (that was passed as argument) and cluster */
|
|
jbe@42
|
1125 else if (
|
|
jbe@42
|
1126 pgl_point_cluster_distance(points+i, cluster) < distance
|
|
jbe@46
|
1127 ) outer++;
|
|
jbe@42
|
1128 }
|
|
jbe@42
|
1129 } else {
|
|
jbe@46
|
1130 /* if point is within cluster, just count sample points within cluster */
|
|
jbe@42
|
1131 for (i=0; i<samples; i++) {
|
|
jbe@46
|
1132 if (pgl_point_in_cluster(points+i, cluster, true)) inner++;
|
|
jbe@42
|
1133 }
|
|
jbe@42
|
1134 }
|
|
jbe@46
|
1135 /* release memory for sample points needed for numerical integration */
|
|
jbe@42
|
1136 pfree(points);
|
|
jbe@46
|
1137 /* if enlargement was less than doubling the area, then combine inner and
|
|
jbe@46
|
1138 outer sample point counts with different weighting */
|
|
jbe@46
|
1139 /* (ensures fairness in such a way that the integral of the squared result
|
|
jbe@46
|
1140 over all possible point parameters is independent of the cluster) */
|
|
jbe@46
|
1141 if (outer < inner) result = (2*inner + 4*outer) / 3.0;
|
|
jbe@46
|
1142 /* otherwise weigh inner and outer points the same */
|
|
jbe@46
|
1143 else result = inner + outer;
|
|
jbe@46
|
1144 /* convert area into distance (i.e. radius of a circle with the same area) */
|
|
jbe@46
|
1145 result = sqrt(area * (result / samples) / M_PI);
|
|
jbe@42
|
1146 /* return result only if it is greater than the distance between point and
|
|
jbe@42
|
1147 cluster to avoid unexpected results because of errors due to limited
|
|
jbe@42
|
1148 precision */
|
|
jbe@42
|
1149 if (result > distance) return result;
|
|
jbe@42
|
1150 /* otherwise return distance between point and cluster */
|
|
jbe@42
|
1151 else return distance;
|
|
jbe@42
|
1152 }
|
|
jbe@42
|
1153
|
|
jbe@42
|
1154
|
|
jbe@16
|
1155 /*-------------------------------------------------*
|
|
jbe@16
|
1156 * geographic index based on space-filling curve *
|
|
jbe@16
|
1157 *-------------------------------------------------*/
|
|
jbe@0
|
1158
|
|
jbe@0
|
1159 /* number of bytes used for geographic (center) position in keys */
|
|
jbe@0
|
1160 #define PGL_KEY_LATLON_BYTELEN 7
|
|
jbe@0
|
1161
|
|
jbe@0
|
1162 /* maximum reference value for logarithmic size of geographic objects */
|
|
jbe@0
|
1163 #define PGL_AREAKEY_REFOBJSIZE (PGL_DIAMETER/3.0) /* can be tweaked */
|
|
jbe@0
|
1164
|
|
jbe@0
|
1165 /* pointer to index key (either pgl_pointkey or pgl_areakey) */
|
|
jbe@0
|
1166 typedef unsigned char *pgl_keyptr;
|
|
jbe@0
|
1167
|
|
jbe@0
|
1168 /* index key for points (objects with zero area) on the spheroid */
|
|
jbe@0
|
1169 /* bit 0..55: interspersed bits of latitude and longitude,
|
|
jbe@0
|
1170 bit 56..57: always zero,
|
|
jbe@0
|
1171 bit 58..63: node depth in hypothetic (full) tree from 0 to 56 (incl.) */
|
|
jbe@0
|
1172 typedef unsigned char pgl_pointkey[PGL_KEY_LATLON_BYTELEN+1];
|
|
jbe@0
|
1173
|
|
jbe@0
|
1174 /* index key for geographic objects on spheroid with area greater than zero */
|
|
jbe@0
|
1175 /* bit 0..55: interspersed bits of latitude and longitude of center point,
|
|
jbe@0
|
1176 bit 56: always set to 1,
|
|
jbe@0
|
1177 bit 57..63: node depth in hypothetic (full) tree from 0 to (2*56)+1 (incl.),
|
|
jbe@0
|
1178 bit 64..71: logarithmic object size from 0 to 56+1 = 57 (incl.), but set to
|
|
jbe@0
|
1179 PGL_KEY_OBJSIZE_EMPTY (with interspersed bits = 0 and node depth
|
|
jbe@0
|
1180 = 113) for empty objects, and set to PGL_KEY_OBJSIZE_UNIVERSAL
|
|
jbe@0
|
1181 (with interspersed bits = 0 and node depth = 0) for keys which
|
|
jbe@0
|
1182 cover both empty and non-empty objects */
|
|
jbe@0
|
1183
|
|
jbe@0
|
1184 typedef unsigned char pgl_areakey[PGL_KEY_LATLON_BYTELEN+2];
|
|
jbe@0
|
1185
|
|
jbe@0
|
1186 /* helper macros for reading/writing index keys */
|
|
jbe@0
|
1187 #define PGL_KEY_NODEDEPTH_OFFSET PGL_KEY_LATLON_BYTELEN
|
|
jbe@0
|
1188 #define PGL_KEY_OBJSIZE_OFFSET (PGL_KEY_NODEDEPTH_OFFSET+1)
|
|
jbe@0
|
1189 #define PGL_POINTKEY_MAXDEPTH (PGL_KEY_LATLON_BYTELEN*8)
|
|
jbe@0
|
1190 #define PGL_AREAKEY_MAXDEPTH (2*PGL_POINTKEY_MAXDEPTH+1)
|
|
jbe@0
|
1191 #define PGL_AREAKEY_MAXOBJSIZE (PGL_POINTKEY_MAXDEPTH+1)
|
|
jbe@0
|
1192 #define PGL_AREAKEY_TYPEMASK 0x80
|
|
jbe@0
|
1193 #define PGL_KEY_LATLONBIT(key, n) ((key)[(n)/8] & (0x80 >> ((n)%8)))
|
|
jbe@0
|
1194 #define PGL_KEY_LATLONBIT_DIFF(key1, key2, n) \
|
|
jbe@0
|
1195 ( PGL_KEY_LATLONBIT(key1, n) ^ \
|
|
jbe@0
|
1196 PGL_KEY_LATLONBIT(key2, n) )
|
|
jbe@0
|
1197 #define PGL_KEY_IS_AREAKEY(key) ((key)[PGL_KEY_NODEDEPTH_OFFSET] & \
|
|
jbe@0
|
1198 PGL_AREAKEY_TYPEMASK)
|
|
jbe@0
|
1199 #define PGL_KEY_NODEDEPTH(key) ((key)[PGL_KEY_NODEDEPTH_OFFSET] & \
|
|
jbe@0
|
1200 (PGL_AREAKEY_TYPEMASK-1))
|
|
jbe@0
|
1201 #define PGL_KEY_OBJSIZE(key) ((key)[PGL_KEY_OBJSIZE_OFFSET])
|
|
jbe@0
|
1202 #define PGL_KEY_OBJSIZE_EMPTY 126
|
|
jbe@0
|
1203 #define PGL_KEY_OBJSIZE_UNIVERSAL 127
|
|
jbe@0
|
1204 #define PGL_KEY_IS_EMPTY(key) ( PGL_KEY_IS_AREAKEY(key) && \
|
|
jbe@0
|
1205 (key)[PGL_KEY_OBJSIZE_OFFSET] == \
|
|
jbe@0
|
1206 PGL_KEY_OBJSIZE_EMPTY )
|
|
jbe@0
|
1207 #define PGL_KEY_IS_UNIVERSAL(key) ( PGL_KEY_IS_AREAKEY(key) && \
|
|
jbe@0
|
1208 (key)[PGL_KEY_OBJSIZE_OFFSET] == \
|
|
jbe@0
|
1209 PGL_KEY_OBJSIZE_UNIVERSAL )
|
|
jbe@0
|
1210
|
|
jbe@0
|
1211 /* set area key to match empty objects only */
|
|
jbe@0
|
1212 static void pgl_key_set_empty(pgl_keyptr key) {
|
|
jbe@0
|
1213 memset(key, 0, sizeof(pgl_areakey));
|
|
jbe@0
|
1214 /* Note: setting node depth to maximum is required for picksplit function */
|
|
jbe@0
|
1215 key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | PGL_AREAKEY_MAXDEPTH;
|
|
jbe@0
|
1216 key[PGL_KEY_OBJSIZE_OFFSET] = PGL_KEY_OBJSIZE_EMPTY;
|
|
jbe@0
|
1217 }
|
|
jbe@0
|
1218
|
|
jbe@0
|
1219 /* set area key to match any object (including empty objects) */
|
|
jbe@0
|
1220 static void pgl_key_set_universal(pgl_keyptr key) {
|
|
jbe@0
|
1221 memset(key, 0, sizeof(pgl_areakey));
|
|
jbe@0
|
1222 key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK;
|
|
jbe@0
|
1223 key[PGL_KEY_OBJSIZE_OFFSET] = PGL_KEY_OBJSIZE_UNIVERSAL;
|
|
jbe@0
|
1224 }
|
|
jbe@0
|
1225
|
|
jbe@0
|
1226 /* convert a point on earth into a max-depth key to be used in index */
|
|
jbe@0
|
1227 static void pgl_point_to_key(pgl_point *point, pgl_keyptr key) {
|
|
jbe@0
|
1228 double lat = point->lat;
|
|
jbe@0
|
1229 double lon = point->lon;
|
|
jbe@0
|
1230 int i;
|
|
jbe@0
|
1231 /* clear latitude and longitude bits */
|
|
jbe@0
|
1232 memset(key, 0, PGL_KEY_LATLON_BYTELEN);
|
|
jbe@0
|
1233 /* set node depth to maximum and type bit to zero */
|
|
jbe@0
|
1234 key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_POINTKEY_MAXDEPTH;
|
|
jbe@0
|
1235 /* iterate over all latitude/longitude bit pairs */
|
|
jbe@0
|
1236 for (i=0; i<PGL_POINTKEY_MAXDEPTH/2; i++) {
|
|
jbe@0
|
1237 /* determine latitude bit */
|
|
jbe@0
|
1238 if (lat >= 0) {
|
|
jbe@0
|
1239 key[i/4] |= 0x80 >> (2*(i%4));
|
|
jbe@0
|
1240 lat *= 2; lat -= 90;
|
|
jbe@0
|
1241 } else {
|
|
jbe@0
|
1242 lat *= 2; lat += 90;
|
|
jbe@0
|
1243 }
|
|
jbe@0
|
1244 /* determine longitude bit */
|
|
jbe@0
|
1245 if (lon >= 0) {
|
|
jbe@0
|
1246 key[i/4] |= 0x80 >> (2*(i%4)+1);
|
|
jbe@0
|
1247 lon *= 2; lon -= 180;
|
|
jbe@0
|
1248 } else {
|
|
jbe@0
|
1249 lon *= 2; lon += 180;
|
|
jbe@0
|
1250 }
|
|
jbe@0
|
1251 }
|
|
jbe@0
|
1252 }
|
|
jbe@0
|
1253
|
|
jbe@0
|
1254 /* convert a circle on earth into a max-depth key to be used in an index */
|
|
jbe@0
|
1255 static void pgl_circle_to_key(pgl_circle *circle, pgl_keyptr key) {
|
|
jbe@0
|
1256 /* handle special case of empty circle */
|
|
jbe@0
|
1257 if (circle->radius < 0) {
|
|
jbe@0
|
1258 pgl_key_set_empty(key);
|
|
jbe@0
|
1259 return;
|
|
jbe@0
|
1260 }
|
|
jbe@0
|
1261 /* perform same action as for point keys */
|
|
jbe@0
|
1262 pgl_point_to_key(&(circle->center), key);
|
|
jbe@0
|
1263 /* but overwrite type and node depth to fit area index key */
|
|
jbe@0
|
1264 key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | PGL_AREAKEY_MAXDEPTH;
|
|
jbe@0
|
1265 /* check if radius is greater than (or equal to) reference size */
|
|
jbe@0
|
1266 /* (treat equal values as greater values for numerical safety) */
|
|
jbe@0
|
1267 if (circle->radius >= PGL_AREAKEY_REFOBJSIZE) {
|
|
jbe@0
|
1268 /* if yes, set logarithmic size to zero */
|
|
jbe@0
|
1269 key[PGL_KEY_OBJSIZE_OFFSET] = 0;
|
|
jbe@0
|
1270 } else {
|
|
jbe@0
|
1271 /* otherwise, determine logarithmic size iteratively */
|
|
jbe@0
|
1272 /* (one step is equivalent to a factor of sqrt(2)) */
|
|
jbe@0
|
1273 double reference = PGL_AREAKEY_REFOBJSIZE / M_SQRT2;
|
|
jbe@0
|
1274 int objsize = 1;
|
|
jbe@0
|
1275 while (objsize < PGL_AREAKEY_MAXOBJSIZE) {
|
|
jbe@0
|
1276 /* stop when radius is greater than (or equal to) adjusted reference */
|
|
jbe@0
|
1277 /* (treat equal values as greater values for numerical safety) */
|
|
jbe@0
|
1278 if (circle->radius >= reference) break;
|
|
jbe@0
|
1279 reference /= M_SQRT2;
|
|
jbe@0
|
1280 objsize++;
|
|
jbe@0
|
1281 }
|
|
jbe@0
|
1282 /* set logarithmic size to determined value */
|
|
jbe@0
|
1283 key[PGL_KEY_OBJSIZE_OFFSET] = objsize;
|
|
jbe@0
|
1284 }
|
|
jbe@0
|
1285 }
|
|
jbe@0
|
1286
|
|
jbe@0
|
1287 /* check if one key is subkey of another key or vice versa */
|
|
jbe@0
|
1288 static bool pgl_keys_overlap(pgl_keyptr key1, pgl_keyptr key2) {
|
|
jbe@0
|
1289 int i; /* key bit offset (includes both lat/lon and log. obj. size bits) */
|
|
jbe@0
|
1290 /* determine smallest depth */
|
|
jbe@0
|
1291 int depth1 = PGL_KEY_NODEDEPTH(key1);
|
|
jbe@0
|
1292 int depth2 = PGL_KEY_NODEDEPTH(key2);
|
|
jbe@0
|
1293 int depth = (depth1 < depth2) ? depth1 : depth2;
|
|
jbe@0
|
1294 /* check if keys are area keys (assuming that both keys have same type) */
|
|
jbe@0
|
1295 if (PGL_KEY_IS_AREAKEY(key1)) {
|
|
jbe@0
|
1296 int j = 0; /* bit offset for logarithmic object size bits */
|
|
jbe@0
|
1297 int k = 0; /* bit offset for latitude and longitude */
|
|
jbe@0
|
1298 /* fetch logarithmic object size information */
|
|
jbe@0
|
1299 int objsize1 = PGL_KEY_OBJSIZE(key1);
|
|
jbe@0
|
1300 int objsize2 = PGL_KEY_OBJSIZE(key2);
|
|
jbe@0
|
1301 /* handle special cases for empty objects (universal and empty keys) */
|
|
jbe@0
|
1302 if (
|
|
jbe@0
|
1303 objsize1 == PGL_KEY_OBJSIZE_UNIVERSAL ||
|
|
jbe@0
|
1304 objsize2 == PGL_KEY_OBJSIZE_UNIVERSAL
|
|
jbe@0
|
1305 ) return true;
|
|
jbe@0
|
1306 if (
|
|
jbe@0
|
1307 objsize1 == PGL_KEY_OBJSIZE_EMPTY ||
|
|
jbe@0
|
1308 objsize2 == PGL_KEY_OBJSIZE_EMPTY
|
|
jbe@0
|
1309 ) return objsize1 == objsize2;
|
|
jbe@0
|
1310 /* iterate through key bits */
|
|
jbe@0
|
1311 for (i=0; i<depth; i++) {
|
|
jbe@0
|
1312 /* every second bit is a bit describing the object size */
|
|
jbe@0
|
1313 if (i%2 == 0) {
|
|
jbe@0
|
1314 /* check if object size bit is different in both keys (objsize1 and
|
|
jbe@0
|
1315 objsize2 describe the minimum index when object size bit is set) */
|
|
jbe@0
|
1316 if (
|
|
jbe@0
|
1317 (objsize1 <= j && objsize2 > j) ||
|
|
jbe@0
|
1318 (objsize2 <= j && objsize1 > j)
|
|
jbe@0
|
1319 ) {
|
|
jbe@0
|
1320 /* bit differs, therefore keys are in separate branches */
|
|
jbe@0
|
1321 return false;
|
|
jbe@0
|
1322 }
|
|
jbe@0
|
1323 /* increase bit counter for object size bits */
|
|
jbe@0
|
1324 j++;
|
|
jbe@0
|
1325 }
|
|
jbe@0
|
1326 /* all other bits describe latitude and longitude */
|
|
jbe@0
|
1327 else {
|
|
jbe@0
|
1328 /* check if bit differs in both keys */
|
|
jbe@0
|
1329 if (PGL_KEY_LATLONBIT_DIFF(key1, key2, k)) {
|
|
jbe@0
|
1330 /* bit differs, therefore keys are in separate branches */
|
|
jbe@0
|
1331 return false;
|
|
jbe@0
|
1332 }
|
|
jbe@0
|
1333 /* increase bit counter for latitude/longitude bits */
|
|
jbe@0
|
1334 k++;
|
|
jbe@0
|
1335 }
|
|
jbe@0
|
1336 }
|
|
jbe@0
|
1337 }
|
|
jbe@0
|
1338 /* if not, keys are point keys */
|
|
jbe@0
|
1339 else {
|
|
jbe@0
|
1340 /* iterate through key bits */
|
|
jbe@0
|
1341 for (i=0; i<depth; i++) {
|
|
jbe@0
|
1342 /* check if bit differs in both keys */
|
|
jbe@0
|
1343 if (PGL_KEY_LATLONBIT_DIFF(key1, key2, i)) {
|
|
jbe@0
|
1344 /* bit differs, therefore keys are in separate branches */
|
|
jbe@0
|
1345 return false;
|
|
jbe@0
|
1346 }
|
|
jbe@0
|
1347 }
|
|
jbe@0
|
1348 }
|
|
jbe@0
|
1349 /* return true because keys are in the same branch */
|
|
jbe@0
|
1350 return true;
|
|
jbe@0
|
1351 }
|
|
jbe@0
|
1352
|
|
jbe@0
|
1353 /* combine two keys into new key which covers both original keys */
|
|
jbe@0
|
1354 /* (result stored in first argument) */
|
|
jbe@0
|
1355 static void pgl_unite_keys(pgl_keyptr dst, pgl_keyptr src) {
|
|
jbe@0
|
1356 int i; /* key bit offset (includes both lat/lon and log. obj. size bits) */
|
|
jbe@0
|
1357 /* determine smallest depth */
|
|
jbe@0
|
1358 int depth1 = PGL_KEY_NODEDEPTH(dst);
|
|
jbe@0
|
1359 int depth2 = PGL_KEY_NODEDEPTH(src);
|
|
jbe@0
|
1360 int depth = (depth1 < depth2) ? depth1 : depth2;
|
|
jbe@0
|
1361 /* check if keys are area keys (assuming that both keys have same type) */
|
|
jbe@0
|
1362 if (PGL_KEY_IS_AREAKEY(dst)) {
|
|
jbe@0
|
1363 pgl_areakey dstbuf = { 0, }; /* destination buffer (cleared) */
|
|
jbe@0
|
1364 int j = 0; /* bit offset for logarithmic object size bits */
|
|
jbe@0
|
1365 int k = 0; /* bit offset for latitude and longitude */
|
|
jbe@0
|
1366 /* fetch logarithmic object size information */
|
|
jbe@0
|
1367 int objsize1 = PGL_KEY_OBJSIZE(dst);
|
|
jbe@0
|
1368 int objsize2 = PGL_KEY_OBJSIZE(src);
|
|
jbe@0
|
1369 /* handle special cases for empty objects (universal and empty keys) */
|
|
jbe@0
|
1370 if (
|
|
jbe@0
|
1371 objsize1 > PGL_AREAKEY_MAXOBJSIZE ||
|
|
jbe@0
|
1372 objsize2 > PGL_AREAKEY_MAXOBJSIZE
|
|
jbe@0
|
1373 ) {
|
|
jbe@0
|
1374 if (
|
|
jbe@0
|
1375 objsize1 == PGL_KEY_OBJSIZE_EMPTY &&
|
|
jbe@0
|
1376 objsize2 == PGL_KEY_OBJSIZE_EMPTY
|
|
jbe@0
|
1377 ) pgl_key_set_empty(dst);
|
|
jbe@0
|
1378 else pgl_key_set_universal(dst);
|
|
jbe@0
|
1379 return;
|
|
jbe@0
|
1380 }
|
|
jbe@0
|
1381 /* iterate through key bits */
|
|
jbe@0
|
1382 for (i=0; i<depth; i++) {
|
|
jbe@0
|
1383 /* every second bit is a bit describing the object size */
|
|
jbe@0
|
1384 if (i%2 == 0) {
|
|
jbe@0
|
1385 /* increase bit counter for object size bits first */
|
|
jbe@0
|
1386 /* (handy when setting objsize variable) */
|
|
jbe@0
|
1387 j++;
|
|
jbe@0
|
1388 /* check if object size bit is set in neither key */
|
|
jbe@0
|
1389 if (objsize1 >= j && objsize2 >= j) {
|
|
jbe@0
|
1390 /* set objsize in destination buffer to indicate that size bit is
|
|
jbe@0
|
1391 unset in destination buffer at the current bit position */
|
|
jbe@0
|
1392 dstbuf[PGL_KEY_OBJSIZE_OFFSET] = j;
|
|
jbe@0
|
1393 }
|
|
jbe@0
|
1394 /* break if object size bit is set in one key only */
|
|
jbe@0
|
1395 else if (objsize1 >= j || objsize2 >= j) break;
|
|
jbe@0
|
1396 }
|
|
jbe@0
|
1397 /* all other bits describe latitude and longitude */
|
|
jbe@0
|
1398 else {
|
|
jbe@0
|
1399 /* break if bit differs in both keys */
|
|
jbe@0
|
1400 if (PGL_KEY_LATLONBIT(dst, k)) {
|
|
jbe@0
|
1401 if (!PGL_KEY_LATLONBIT(src, k)) break;
|
|
jbe@0
|
1402 /* but set bit in destination buffer if bit is set in both keys */
|
|
jbe@0
|
1403 dstbuf[k/8] |= 0x80 >> (k%8);
|
|
jbe@0
|
1404 } else if (PGL_KEY_LATLONBIT(src, k)) break;
|
|
jbe@0
|
1405 /* increase bit counter for latitude/longitude bits */
|
|
jbe@0
|
1406 k++;
|
|
jbe@0
|
1407 }
|
|
jbe@0
|
1408 }
|
|
jbe@0
|
1409 /* set common node depth and type bit (type bit = 1) */
|
|
jbe@0
|
1410 dstbuf[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | i;
|
|
jbe@0
|
1411 /* copy contents of destination buffer to first key */
|
|
jbe@0
|
1412 memcpy(dst, dstbuf, sizeof(pgl_areakey));
|
|
jbe@0
|
1413 }
|
|
jbe@0
|
1414 /* if not, keys are point keys */
|
|
jbe@0
|
1415 else {
|
|
jbe@0
|
1416 pgl_pointkey dstbuf = { 0, }; /* destination buffer (cleared) */
|
|
jbe@0
|
1417 /* iterate through key bits */
|
|
jbe@0
|
1418 for (i=0; i<depth; i++) {
|
|
jbe@0
|
1419 /* break if bit differs in both keys */
|
|
jbe@0
|
1420 if (PGL_KEY_LATLONBIT(dst, i)) {
|
|
jbe@0
|
1421 if (!PGL_KEY_LATLONBIT(src, i)) break;
|
|
jbe@0
|
1422 /* but set bit in destination buffer if bit is set in both keys */
|
|
jbe@0
|
1423 dstbuf[i/8] |= 0x80 >> (i%8);
|
|
jbe@0
|
1424 } else if (PGL_KEY_LATLONBIT(src, i)) break;
|
|
jbe@0
|
1425 }
|
|
jbe@0
|
1426 /* set common node depth (type bit = 0) */
|
|
jbe@0
|
1427 dstbuf[PGL_KEY_NODEDEPTH_OFFSET] = i;
|
|
jbe@0
|
1428 /* copy contents of destination buffer to first key */
|
|
jbe@0
|
1429 memcpy(dst, dstbuf, sizeof(pgl_pointkey));
|
|
jbe@0
|
1430 }
|
|
jbe@0
|
1431 }
|
|
jbe@0
|
1432
|
|
jbe@0
|
1433 /* determine center(!) boundaries and radius estimation of index key */
|
|
jbe@0
|
1434 static double pgl_key_to_box(pgl_keyptr key, pgl_box *box) {
|
|
jbe@0
|
1435 int i;
|
|
jbe@0
|
1436 /* determine node depth */
|
|
jbe@0
|
1437 int depth = PGL_KEY_NODEDEPTH(key);
|
|
jbe@0
|
1438 /* center point of possible result */
|
|
jbe@0
|
1439 double lat = 0;
|
|
jbe@0
|
1440 double lon = 0;
|
|
jbe@0
|
1441 /* maximum distance of real center point from key center */
|
|
jbe@0
|
1442 double dlat = 90;
|
|
jbe@0
|
1443 double dlon = 180;
|
|
jbe@0
|
1444 /* maximum radius of contained objects */
|
|
jbe@0
|
1445 double radius = 0; /* always return zero for point index keys */
|
|
jbe@0
|
1446 /* check if key is area key */
|
|
jbe@0
|
1447 if (PGL_KEY_IS_AREAKEY(key)) {
|
|
jbe@0
|
1448 /* get logarithmic object size */
|
|
jbe@0
|
1449 int objsize = PGL_KEY_OBJSIZE(key);
|
|
jbe@0
|
1450 /* handle special cases for empty objects (universal and empty keys) */
|
|
jbe@0
|
1451 if (objsize == PGL_KEY_OBJSIZE_EMPTY) {
|
|
jbe@0
|
1452 pgl_box_set_empty(box);
|
|
jbe@0
|
1453 return 0;
|
|
jbe@0
|
1454 } else if (objsize == PGL_KEY_OBJSIZE_UNIVERSAL) {
|
|
jbe@0
|
1455 box->lat_min = -90;
|
|
jbe@0
|
1456 box->lat_max = 90;
|
|
jbe@0
|
1457 box->lon_min = -180;
|
|
jbe@0
|
1458 box->lon_max = 180;
|
|
jbe@0
|
1459 return 0; /* any value >= 0 would do */
|
|
jbe@0
|
1460 }
|
|
jbe@0
|
1461 /* calculate maximum possible radius of objects covered by the given key */
|
|
jbe@0
|
1462 if (objsize == 0) radius = INFINITY;
|
|
jbe@0
|
1463 else {
|
|
jbe@0
|
1464 radius = PGL_AREAKEY_REFOBJSIZE;
|
|
jbe@0
|
1465 while (--objsize) radius /= M_SQRT2;
|
|
jbe@0
|
1466 }
|
|
jbe@0
|
1467 /* iterate over latitude and longitude bits in key */
|
|
jbe@0
|
1468 /* (every second bit is a latitude or longitude bit) */
|
|
jbe@0
|
1469 for (i=0; i<depth/2; i++) {
|
|
jbe@0
|
1470 /* check if latitude bit */
|
|
jbe@0
|
1471 if (i%2 == 0) {
|
|
jbe@0
|
1472 /* cut latitude dimension in half */
|
|
jbe@0
|
1473 dlat /= 2;
|
|
jbe@0
|
1474 /* increase center latitude if bit is 1, otherwise decrease */
|
|
jbe@0
|
1475 if (PGL_KEY_LATLONBIT(key, i)) lat += dlat;
|
|
jbe@0
|
1476 else lat -= dlat;
|
|
jbe@0
|
1477 }
|
|
jbe@0
|
1478 /* otherwise longitude bit */
|
|
jbe@0
|
1479 else {
|
|
jbe@0
|
1480 /* cut longitude dimension in half */
|
|
jbe@0
|
1481 dlon /= 2;
|
|
jbe@0
|
1482 /* increase center longitude if bit is 1, otherwise decrease */
|
|
jbe@0
|
1483 if (PGL_KEY_LATLONBIT(key, i)) lon += dlon;
|
|
jbe@0
|
1484 else lon -= dlon;
|
|
jbe@0
|
1485 }
|
|
jbe@0
|
1486 }
|
|
jbe@0
|
1487 }
|
|
jbe@0
|
1488 /* if not, keys are point keys */
|
|
jbe@0
|
1489 else {
|
|
jbe@0
|
1490 /* iterate over all bits in key */
|
|
jbe@0
|
1491 for (i=0; i<depth; i++) {
|
|
jbe@0
|
1492 /* check if latitude bit */
|
|
jbe@0
|
1493 if (i%2 == 0) {
|
|
jbe@0
|
1494 /* cut latitude dimension in half */
|
|
jbe@0
|
1495 dlat /= 2;
|
|
jbe@0
|
1496 /* increase center latitude if bit is 1, otherwise decrease */
|
|
jbe@0
|
1497 if (PGL_KEY_LATLONBIT(key, i)) lat += dlat;
|
|
jbe@0
|
1498 else lat -= dlat;
|
|
jbe@0
|
1499 }
|
|
jbe@0
|
1500 /* otherwise longitude bit */
|
|
jbe@0
|
1501 else {
|
|
jbe@0
|
1502 /* cut longitude dimension in half */
|
|
jbe@0
|
1503 dlon /= 2;
|
|
jbe@0
|
1504 /* increase center longitude if bit is 1, otherwise decrease */
|
|
jbe@0
|
1505 if (PGL_KEY_LATLONBIT(key, i)) lon += dlon;
|
|
jbe@0
|
1506 else lon -= dlon;
|
|
jbe@0
|
1507 }
|
|
jbe@0
|
1508 }
|
|
jbe@0
|
1509 }
|
|
jbe@0
|
1510 /* calculate boundaries from center point and remaining dlat and dlon */
|
|
jbe@0
|
1511 /* (return values through pointer to box) */
|
|
jbe@0
|
1512 box->lat_min = lat - dlat;
|
|
jbe@0
|
1513 box->lat_max = lat + dlat;
|
|
jbe@0
|
1514 box->lon_min = lon - dlon;
|
|
jbe@0
|
1515 box->lon_max = lon + dlon;
|
|
jbe@0
|
1516 /* return radius (as a function return value) */
|
|
jbe@0
|
1517 return radius;
|
|
jbe@0
|
1518 }
|
|
jbe@0
|
1519
|
|
jbe@0
|
1520 /* estimator function for distance between point and index key */
|
|
jbe@16
|
1521 /* always returns a smaller value than actually correct or zero */
|
|
jbe@0
|
1522 static double pgl_estimate_key_distance(pgl_keyptr key, pgl_point *point) {
|
|
jbe@0
|
1523 pgl_box box; /* center(!) bounding box of area index key */
|
|
jbe@0
|
1524 /* calculate center(!) bounding box and maximum radius of objects covered
|
|
jbe@0
|
1525 by area index key (radius is zero for point index keys) */
|
|
jbe@0
|
1526 double distance = pgl_key_to_box(key, &box);
|
|
jbe@0
|
1527 /* calculate estimated distance between bounding box of center point of
|
|
jbe@0
|
1528 indexed object and point passed as second argument, then substract maximum
|
|
jbe@0
|
1529 radius of objects covered by index key */
|
|
jbe@16
|
1530 distance = pgl_estimate_point_box_distance(point, &box) - distance;
|
|
jbe@0
|
1531 /* truncate negative results to zero */
|
|
jbe@0
|
1532 if (distance <= 0) distance = 0;
|
|
jbe@0
|
1533 /* return result */
|
|
jbe@0
|
1534 return distance;
|
|
jbe@0
|
1535 }
|
|
jbe@0
|
1536
|
|
jbe@0
|
1537
|
|
jbe@0
|
1538 /*---------------------------------*
|
|
jbe@0
|
1539 * helper functions for text I/O *
|
|
jbe@0
|
1540 *---------------------------------*/
|
|
jbe@0
|
1541
|
|
jbe@0
|
1542 #define PGL_NUMBUFLEN 64 /* buffer size for number to string conversion */
|
|
jbe@0
|
1543
|
|
jbe@0
|
1544 /* convert floating point number to string (round-trip safe) */
|
|
jbe@0
|
1545 static void pgl_print_float(char *buf, double flt) {
|
|
jbe@0
|
1546 /* check if number is integral */
|
|
jbe@0
|
1547 if (trunc(flt) == flt) {
|
|
jbe@0
|
1548 /* for integral floats use maximum precision */
|
|
jbe@0
|
1549 snprintf(buf, PGL_NUMBUFLEN, "%.17g", flt);
|
|
jbe@0
|
1550 } else {
|
|
jbe@0
|
1551 /* otherwise check if 15, 16, or 17 digits needed (round-trip safety) */
|
|
jbe@0
|
1552 snprintf(buf, PGL_NUMBUFLEN, "%.15g", flt);
|
|
jbe@0
|
1553 if (strtod(buf, NULL) != flt) snprintf(buf, PGL_NUMBUFLEN, "%.16g", flt);
|
|
jbe@0
|
1554 if (strtod(buf, NULL) != flt) snprintf(buf, PGL_NUMBUFLEN, "%.17g", flt);
|
|
jbe@0
|
1555 }
|
|
jbe@0
|
1556 }
|
|
jbe@0
|
1557
|
|
jbe@0
|
1558 /* convert latitude floating point number (in degrees) to string */
|
|
jbe@0
|
1559 static void pgl_print_lat(char *buf, double lat) {
|
|
jbe@0
|
1560 if (signbit(lat)) {
|
|
jbe@0
|
1561 /* treat negative latitudes (including -0) as south */
|
|
jbe@0
|
1562 snprintf(buf, PGL_NUMBUFLEN, "S%015.12f", -lat);
|
|
jbe@0
|
1563 } else {
|
|
jbe@0
|
1564 /* treat positive latitudes (including +0) as north */
|
|
jbe@0
|
1565 snprintf(buf, PGL_NUMBUFLEN, "N%015.12f", lat);
|
|
jbe@0
|
1566 }
|
|
jbe@0
|
1567 }
|
|
jbe@0
|
1568
|
|
jbe@0
|
1569 /* convert longitude floating point number (in degrees) to string */
|
|
jbe@0
|
1570 static void pgl_print_lon(char *buf, double lon) {
|
|
jbe@0
|
1571 if (signbit(lon)) {
|
|
jbe@0
|
1572 /* treat negative longitudes (including -0) as west */
|
|
jbe@0
|
1573 snprintf(buf, PGL_NUMBUFLEN, "W%016.12f", -lon);
|
|
jbe@0
|
1574 } else {
|
|
jbe@0
|
1575 /* treat positive longitudes (including +0) as east */
|
|
jbe@0
|
1576 snprintf(buf, PGL_NUMBUFLEN, "E%016.12f", lon);
|
|
jbe@0
|
1577 }
|
|
jbe@0
|
1578 }
|
|
jbe@0
|
1579
|
|
jbe@0
|
1580 /* bit masks used as return value of pgl_scan() function */
|
|
jbe@0
|
1581 #define PGL_SCAN_NONE 0 /* no value has been parsed */
|
|
jbe@0
|
1582 #define PGL_SCAN_LAT (1<<0) /* latitude has been parsed */
|
|
jbe@0
|
1583 #define PGL_SCAN_LON (1<<1) /* longitude has been parsed */
|
|
jbe@0
|
1584 #define PGL_SCAN_LATLON (PGL_SCAN_LAT | PGL_SCAN_LON) /* bitwise OR of both */
|
|
jbe@0
|
1585
|
|
jbe@0
|
1586 /* parse a coordinate (can be latitude or longitude) */
|
|
jbe@0
|
1587 static int pgl_scan(char **str, double *lat, double *lon) {
|
|
jbe@0
|
1588 double val;
|
|
jbe@0
|
1589 int len;
|
|
jbe@0
|
1590 if (
|
|
jbe@0
|
1591 sscanf(*str, " N %lf %n", &val, &len) ||
|
|
jbe@0
|
1592 sscanf(*str, " n %lf %n", &val, &len)
|
|
jbe@0
|
1593 ) {
|
|
jbe@0
|
1594 *str += len; *lat = val; return PGL_SCAN_LAT;
|
|
jbe@0
|
1595 }
|
|
jbe@0
|
1596 if (
|
|
jbe@0
|
1597 sscanf(*str, " S %lf %n", &val, &len) ||
|
|
jbe@0
|
1598 sscanf(*str, " s %lf %n", &val, &len)
|
|
jbe@0
|
1599 ) {
|
|
jbe@0
|
1600 *str += len; *lat = -val; return PGL_SCAN_LAT;
|
|
jbe@0
|
1601 }
|
|
jbe@0
|
1602 if (
|
|
jbe@0
|
1603 sscanf(*str, " E %lf %n", &val, &len) ||
|
|
jbe@0
|
1604 sscanf(*str, " e %lf %n", &val, &len)
|
|
jbe@0
|
1605 ) {
|
|
jbe@0
|
1606 *str += len; *lon = val; return PGL_SCAN_LON;
|
|
jbe@0
|
1607 }
|
|
jbe@0
|
1608 if (
|
|
jbe@0
|
1609 sscanf(*str, " W %lf %n", &val, &len) ||
|
|
jbe@0
|
1610 sscanf(*str, " w %lf %n", &val, &len)
|
|
jbe@0
|
1611 ) {
|
|
jbe@0
|
1612 *str += len; *lon = -val; return PGL_SCAN_LON;
|
|
jbe@0
|
1613 }
|
|
jbe@0
|
1614 return PGL_SCAN_NONE;
|
|
jbe@0
|
1615 }
|
|
jbe@0
|
1616
|
|
jbe@0
|
1617
|
|
jbe@0
|
1618 /*-----------------*
|
|
jbe@0
|
1619 * SQL functions *
|
|
jbe@0
|
1620 *-----------------*/
|
|
jbe@0
|
1621
|
|
jbe@0
|
1622 /* Note: These function names use "epoint", "ebox", etc. notation here instead
|
|
jbe@0
|
1623 of "point", "box", etc. in order to distinguish them from any previously
|
|
jbe@0
|
1624 defined functions. */
|
|
jbe@0
|
1625
|
|
jbe@0
|
1626 /* function needed for dummy types and/or not implemented features */
|
|
jbe@0
|
1627 PG_FUNCTION_INFO_V1(pgl_notimpl);
|
|
jbe@0
|
1628 Datum pgl_notimpl(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
1629 ereport(ERROR, (errmsg("not implemented by pgLatLon")));
|
|
jbe@0
|
1630 }
|
|
jbe@0
|
1631
|
|
jbe@0
|
1632 /* set point to latitude and longitude (including checks) */
|
|
jbe@0
|
1633 static void pgl_epoint_set_latlon(pgl_point *point, double lat, double lon) {
|
|
jbe@0
|
1634 /* reject infinite or NaN values */
|
|
jbe@0
|
1635 if (!isfinite(lat) || !isfinite(lon)) {
|
|
jbe@0
|
1636 ereport(ERROR, (
|
|
jbe@0
|
1637 errcode(ERRCODE_DATA_EXCEPTION),
|
|
jbe@0
|
1638 errmsg("epoint requires finite coordinates")
|
|
jbe@0
|
1639 ));
|
|
jbe@0
|
1640 }
|
|
jbe@0
|
1641 /* check latitude bounds */
|
|
jbe@0
|
1642 if (lat < -90) {
|
|
jbe@0
|
1643 ereport(WARNING, (errmsg("latitude exceeds south pole")));
|
|
jbe@0
|
1644 lat = -90;
|
|
jbe@0
|
1645 } else if (lat > 90) {
|
|
jbe@0
|
1646 ereport(WARNING, (errmsg("latitude exceeds north pole")));
|
|
jbe@0
|
1647 lat = 90;
|
|
jbe@0
|
1648 }
|
|
jbe@70
|
1649 /* normalize longitude */
|
|
jbe@70
|
1650 lon = pgl_normalize(lon, true);
|
|
jbe@0
|
1651 /* store rounded latitude/longitude values for round-trip safety */
|
|
jbe@0
|
1652 point->lat = pgl_round(lat);
|
|
jbe@0
|
1653 point->lon = pgl_round(lon);
|
|
jbe@0
|
1654 }
|
|
jbe@0
|
1655
|
|
jbe@0
|
1656 /* create point ("epoint" in SQL) from latitude and longitude */
|
|
jbe@0
|
1657 PG_FUNCTION_INFO_V1(pgl_create_epoint);
|
|
jbe@0
|
1658 Datum pgl_create_epoint(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
1659 pgl_point *point = (pgl_point *)palloc(sizeof(pgl_point));
|
|
jbe@0
|
1660 pgl_epoint_set_latlon(point, PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1));
|
|
jbe@0
|
1661 PG_RETURN_POINTER(point);
|
|
jbe@0
|
1662 }
|
|
jbe@0
|
1663
|
|
jbe@0
|
1664 /* parse point ("epoint" in SQL) */
|
|
jbe@0
|
1665 /* format: '[NS]<float> [EW]<float>' */
|
|
jbe@0
|
1666 PG_FUNCTION_INFO_V1(pgl_epoint_in);
|
|
jbe@0
|
1667 Datum pgl_epoint_in(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
1668 char *str = PG_GETARG_CSTRING(0); /* input string */
|
|
jbe@0
|
1669 char *strptr = str; /* current position within string */
|
|
jbe@0
|
1670 int done = 0; /* bit mask storing if latitude or longitude was read */
|
|
jbe@0
|
1671 double lat, lon; /* parsed values as double precision floats */
|
|
jbe@0
|
1672 pgl_point *point; /* return value (to be palloc'ed) */
|
|
jbe@0
|
1673 /* parse two floats (each latitude or longitude) separated by white-space */
|
|
jbe@0
|
1674 done |= pgl_scan(&strptr, &lat, &lon);
|
|
jbe@0
|
1675 if (strptr != str && isspace(strptr[-1])) {
|
|
jbe@0
|
1676 done |= pgl_scan(&strptr, &lat, &lon);
|
|
jbe@0
|
1677 }
|
|
jbe@0
|
1678 /* require end of string, and latitude and longitude parsed successfully */
|
|
jbe@0
|
1679 if (strptr[0] || done != PGL_SCAN_LATLON) {
|
|
jbe@0
|
1680 ereport(ERROR, (
|
|
jbe@0
|
1681 errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
|
|
jbe@0
|
1682 errmsg("invalid input syntax for type epoint: \"%s\"", str)
|
|
jbe@0
|
1683 ));
|
|
jbe@0
|
1684 }
|
|
jbe@0
|
1685 /* allocate memory for result */
|
|
jbe@0
|
1686 point = (pgl_point *)palloc(sizeof(pgl_point));
|
|
jbe@0
|
1687 /* set latitude and longitude (and perform checks) */
|
|
jbe@0
|
1688 pgl_epoint_set_latlon(point, lat, lon);
|
|
jbe@0
|
1689 /* return result */
|
|
jbe@0
|
1690 PG_RETURN_POINTER(point);
|
|
jbe@0
|
1691 }
|
|
jbe@0
|
1692
|
|
jbe@46
|
1693 /* set sample count for numerical integration (including checks) */
|
|
jbe@46
|
1694 static void pgl_epoint_set_sample_count(pgl_point_sc *search, int32 samples) {
|
|
jbe@46
|
1695 /* require minimum of 6 samples */
|
|
jbe@46
|
1696 if (samples < 6) {
|
|
jbe@46
|
1697 ereport(ERROR, (
|
|
jbe@46
|
1698 errcode(ERRCODE_DATA_EXCEPTION),
|
|
jbe@46
|
1699 errmsg("too few sample points for numerical integration (minimum 6)")
|
|
jbe@46
|
1700 ));
|
|
jbe@46
|
1701 }
|
|
jbe@46
|
1702 /* limit sample count to avoid integer overflows on memory allocation */
|
|
jbe@46
|
1703 if (samples > PGL_CLUSTER_MAXPOINTS) {
|
|
jbe@46
|
1704 ereport(ERROR, (
|
|
jbe@46
|
1705 errcode(ERRCODE_DATA_EXCEPTION),
|
|
jbe@46
|
1706 errmsg(
|
|
jbe@46
|
1707 "too many sample points for numerical integration (maximum %i)",
|
|
jbe@46
|
1708 PGL_CLUSTER_MAXPOINTS
|
|
jbe@46
|
1709 )
|
|
jbe@46
|
1710 ));
|
|
jbe@46
|
1711 }
|
|
jbe@46
|
1712 search->samples = samples;
|
|
jbe@46
|
1713 }
|
|
jbe@46
|
1714
|
|
jbe@46
|
1715 /* create point with sample count for fair distance calculation
|
|
jbe@46
|
1716 ("epoint_with_sample_count" in SQL) from epoint and integer */
|
|
jbe@46
|
1717 PG_FUNCTION_INFO_V1(pgl_create_epoint_with_sample_count);
|
|
jbe@46
|
1718 Datum pgl_create_epoint_with_sample_count(PG_FUNCTION_ARGS) {
|
|
jbe@46
|
1719 pgl_point_sc *search = (pgl_point_sc *)palloc(sizeof(pgl_point_sc));
|
|
jbe@46
|
1720 search->point = *(pgl_point *)PG_GETARG_POINTER(0);
|
|
jbe@46
|
1721 pgl_epoint_set_sample_count(search, PG_GETARG_INT32(1));
|
|
jbe@46
|
1722 PG_RETURN_POINTER(search);
|
|
jbe@46
|
1723 }
|
|
jbe@46
|
1724
|
|
jbe@46
|
1725 /* parse point with sample count ("epoint_with_sample_count" in SQL) */
|
|
jbe@46
|
1726 /* format: '[NS]<float> [EW]<float> <integer>' */
|
|
jbe@46
|
1727 PG_FUNCTION_INFO_V1(pgl_epoint_with_sample_count_in);
|
|
jbe@46
|
1728 Datum pgl_epoint_with_sample_count_in(PG_FUNCTION_ARGS) {
|
|
jbe@46
|
1729 char *str = PG_GETARG_CSTRING(0); /* input string */
|
|
jbe@46
|
1730 char *strptr = str; /* current position within string */
|
|
jbe@46
|
1731 double lat, lon; /* parsed values for latitude and longitude */
|
|
jbe@46
|
1732 int samples; /* parsed value for sample count */
|
|
jbe@46
|
1733 int valid = 0; /* number of valid chars */
|
|
jbe@46
|
1734 int done = 0; /* stores if latitude and/or longitude was read */
|
|
jbe@46
|
1735 pgl_point_sc *search; /* return value (to be palloc'ed) */
|
|
jbe@46
|
1736 /* demand three blocks separated by whitespace */
|
|
jbe@46
|
1737 sscanf(strptr, " %*s %*s %*s %n", &valid);
|
|
jbe@46
|
1738 /* if three blocks separated by whitespace exist, parse those blocks */
|
|
jbe@46
|
1739 if (strptr[valid] == 0) {
|
|
jbe@46
|
1740 /* parse latitude and longitude */
|
|
jbe@46
|
1741 done |= pgl_scan(&strptr, &lat, &lon);
|
|
jbe@46
|
1742 done |= pgl_scan(&strptr, &lat, &lon);
|
|
jbe@46
|
1743 /* parse sample count (while incr. strptr by number of bytes parsed) */
|
|
jbe@46
|
1744 valid = 0;
|
|
jbe@46
|
1745 if (sscanf(strptr, " %d %n", &samples, &valid) == 1) strptr += valid;
|
|
jbe@46
|
1746 }
|
|
jbe@46
|
1747 /* require end of string and both latitude and longitude being parsed */
|
|
jbe@46
|
1748 if (strptr[0] || done != PGL_SCAN_LATLON) {
|
|
jbe@46
|
1749 ereport(ERROR, (
|
|
jbe@46
|
1750 errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
|
|
jbe@46
|
1751 errmsg("invalid input syntax for type ecircle: \"%s\"", str)
|
|
jbe@46
|
1752 ));
|
|
jbe@46
|
1753 }
|
|
jbe@46
|
1754 /* allocate memory for result */
|
|
jbe@46
|
1755 search = (pgl_point_sc *)palloc(sizeof(pgl_point_sc));
|
|
jbe@46
|
1756 /* set latitude, longitude, and sample count (while performing checks) */
|
|
jbe@46
|
1757 pgl_epoint_set_latlon(&search->point, lat, lon);
|
|
jbe@46
|
1758 pgl_epoint_set_sample_count(search, samples);
|
|
jbe@46
|
1759 /* return result */
|
|
jbe@46
|
1760 PG_RETURN_POINTER(search);
|
|
jbe@46
|
1761 }
|
|
jbe@46
|
1762
|
|
jbe@0
|
1763 /* create box ("ebox" in SQL) that is empty */
|
|
jbe@0
|
1764 PG_FUNCTION_INFO_V1(pgl_create_empty_ebox);
|
|
jbe@0
|
1765 Datum pgl_create_empty_ebox(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
1766 pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box));
|
|
jbe@0
|
1767 pgl_box_set_empty(box);
|
|
jbe@0
|
1768 PG_RETURN_POINTER(box);
|
|
jbe@0
|
1769 }
|
|
jbe@0
|
1770
|
|
jbe@0
|
1771 /* set box to given boundaries (including checks) */
|
|
jbe@0
|
1772 static void pgl_ebox_set_boundaries(
|
|
jbe@0
|
1773 pgl_box *box,
|
|
jbe@0
|
1774 double lat_min, double lat_max, double lon_min, double lon_max
|
|
jbe@0
|
1775 ) {
|
|
jbe@0
|
1776 /* if minimum latitude is greater than maximum latitude, return empty box */
|
|
jbe@0
|
1777 if (lat_min > lat_max) {
|
|
jbe@0
|
1778 pgl_box_set_empty(box);
|
|
jbe@0
|
1779 return;
|
|
jbe@0
|
1780 }
|
|
jbe@0
|
1781 /* otherwise reject infinite or NaN values */
|
|
jbe@0
|
1782 if (
|
|
jbe@0
|
1783 !isfinite(lat_min) || !isfinite(lat_max) ||
|
|
jbe@0
|
1784 !isfinite(lon_min) || !isfinite(lon_max)
|
|
jbe@0
|
1785 ) {
|
|
jbe@0
|
1786 ereport(ERROR, (
|
|
jbe@0
|
1787 errcode(ERRCODE_DATA_EXCEPTION),
|
|
jbe@0
|
1788 errmsg("ebox requires finite coordinates")
|
|
jbe@0
|
1789 ));
|
|
jbe@0
|
1790 }
|
|
jbe@0
|
1791 /* check latitude bounds */
|
|
jbe@0
|
1792 if (lat_max < -90) {
|
|
jbe@0
|
1793 ereport(WARNING, (errmsg("northern latitude exceeds south pole")));
|
|
jbe@0
|
1794 lat_max = -90;
|
|
jbe@0
|
1795 } else if (lat_max > 90) {
|
|
jbe@0
|
1796 ereport(WARNING, (errmsg("northern latitude exceeds north pole")));
|
|
jbe@0
|
1797 lat_max = 90;
|
|
jbe@0
|
1798 }
|
|
jbe@0
|
1799 if (lat_min < -90) {
|
|
jbe@0
|
1800 ereport(WARNING, (errmsg("southern latitude exceeds south pole")));
|
|
jbe@0
|
1801 lat_min = -90;
|
|
jbe@0
|
1802 } else if (lat_min > 90) {
|
|
jbe@0
|
1803 ereport(WARNING, (errmsg("southern latitude exceeds north pole")));
|
|
jbe@0
|
1804 lat_min = 90;
|
|
jbe@0
|
1805 }
|
|
jbe@0
|
1806 /* check if all longitudes are included */
|
|
jbe@0
|
1807 if (lon_max - lon_min >= 360) {
|
|
jbe@0
|
1808 if (lon_max - lon_min > 360) ereport(WARNING, (
|
|
jbe@0
|
1809 errmsg("longitude coverage greater than 360 degrees")
|
|
jbe@0
|
1810 ));
|
|
jbe@0
|
1811 lon_min = -180;
|
|
jbe@0
|
1812 lon_max = 180;
|
|
jbe@0
|
1813 } else {
|
|
jbe@0
|
1814 /* normalize longitude bounds */
|
|
jbe@70
|
1815 lon_min = pgl_normalize(lon_min, false);
|
|
jbe@70
|
1816 lon_max = pgl_normalize(lon_max, false);
|
|
jbe@0
|
1817 }
|
|
jbe@0
|
1818 /* store rounded latitude/longitude values for round-trip safety */
|
|
jbe@0
|
1819 box->lat_min = pgl_round(lat_min);
|
|
jbe@0
|
1820 box->lat_max = pgl_round(lat_max);
|
|
jbe@0
|
1821 box->lon_min = pgl_round(lon_min);
|
|
jbe@0
|
1822 box->lon_max = pgl_round(lon_max);
|
|
jbe@0
|
1823 /* ensure that rounding does not change orientation */
|
|
jbe@0
|
1824 if (lon_min > lon_max && box->lon_min == box->lon_max) {
|
|
jbe@0
|
1825 box->lon_min = -180;
|
|
jbe@0
|
1826 box->lon_max = 180;
|
|
jbe@0
|
1827 }
|
|
jbe@0
|
1828 }
|
|
jbe@0
|
1829
|
|
jbe@0
|
1830 /* create box ("ebox" in SQL) from min/max latitude and min/max longitude */
|
|
jbe@0
|
1831 PG_FUNCTION_INFO_V1(pgl_create_ebox);
|
|
jbe@0
|
1832 Datum pgl_create_ebox(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
1833 pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box));
|
|
jbe@0
|
1834 pgl_ebox_set_boundaries(
|
|
jbe@0
|
1835 box,
|
|
jbe@0
|
1836 PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1),
|
|
jbe@0
|
1837 PG_GETARG_FLOAT8(2), PG_GETARG_FLOAT8(3)
|
|
jbe@0
|
1838 );
|
|
jbe@0
|
1839 PG_RETURN_POINTER(box);
|
|
jbe@0
|
1840 }
|
|
jbe@0
|
1841
|
|
jbe@0
|
1842 /* create box ("ebox" in SQL) from two points ("epoint"s) */
|
|
jbe@0
|
1843 /* (can not be used to cover a longitude range of more than 120 degrees) */
|
|
jbe@0
|
1844 PG_FUNCTION_INFO_V1(pgl_create_ebox_from_epoints);
|
|
jbe@0
|
1845 Datum pgl_create_ebox_from_epoints(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
1846 pgl_point *point1 = (pgl_point *)PG_GETARG_POINTER(0);
|
|
jbe@0
|
1847 pgl_point *point2 = (pgl_point *)PG_GETARG_POINTER(1);
|
|
jbe@0
|
1848 pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box));
|
|
jbe@0
|
1849 double lat_min, lat_max, lon_min, lon_max;
|
|
jbe@0
|
1850 double dlon; /* longitude range (delta longitude) */
|
|
jbe@0
|
1851 /* order latitude and longitude boundaries */
|
|
jbe@0
|
1852 if (point2->lat < point1->lat) {
|
|
jbe@0
|
1853 lat_min = point2->lat;
|
|
jbe@0
|
1854 lat_max = point1->lat;
|
|
jbe@0
|
1855 } else {
|
|
jbe@0
|
1856 lat_min = point1->lat;
|
|
jbe@0
|
1857 lat_max = point2->lat;
|
|
jbe@0
|
1858 }
|
|
jbe@0
|
1859 if (point2->lon < point1->lon) {
|
|
jbe@0
|
1860 lon_min = point2->lon;
|
|
jbe@0
|
1861 lon_max = point1->lon;
|
|
jbe@0
|
1862 } else {
|
|
jbe@0
|
1863 lon_min = point1->lon;
|
|
jbe@0
|
1864 lon_max = point2->lon;
|
|
jbe@0
|
1865 }
|
|
jbe@0
|
1866 /* calculate longitude range (round to avoid floating point errors) */
|
|
jbe@0
|
1867 dlon = pgl_round(lon_max - lon_min);
|
|
jbe@0
|
1868 /* determine east-west direction */
|
|
jbe@0
|
1869 if (dlon >= 240) {
|
|
jbe@0
|
1870 /* assume that 180th meridian is crossed and swap min/max longitude */
|
|
jbe@0
|
1871 double swap = lon_min; lon_min = lon_max; lon_max = swap;
|
|
jbe@0
|
1872 } else if (dlon > 120) {
|
|
jbe@0
|
1873 /* unclear orientation since delta longitude > 120 */
|
|
jbe@0
|
1874 ereport(ERROR, (
|
|
jbe@0
|
1875 errcode(ERRCODE_DATA_EXCEPTION),
|
|
jbe@0
|
1876 errmsg("can not determine east/west orientation for ebox")
|
|
jbe@0
|
1877 ));
|
|
jbe@0
|
1878 }
|
|
jbe@0
|
1879 /* use boundaries to setup box (and perform checks) */
|
|
jbe@0
|
1880 pgl_ebox_set_boundaries(box, lat_min, lat_max, lon_min, lon_max);
|
|
jbe@0
|
1881 /* return result */
|
|
jbe@0
|
1882 PG_RETURN_POINTER(box);
|
|
jbe@0
|
1883 }
|
|
jbe@0
|
1884
|
|
jbe@0
|
1885 /* parse box ("ebox" in SQL) */
|
|
jbe@0
|
1886 /* format: '[NS]<float> [EW]<float> [NS]<float> [EW]<float>'
|
|
jbe@0
|
1887 or: '[NS]<float> [NS]<float> [EW]<float> [EW]<float>' */
|
|
jbe@0
|
1888 PG_FUNCTION_INFO_V1(pgl_ebox_in);
|
|
jbe@0
|
1889 Datum pgl_ebox_in(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
1890 char *str = PG_GETARG_CSTRING(0); /* input string */
|
|
jbe@0
|
1891 char *str_lower; /* lower case version of input string */
|
|
jbe@0
|
1892 char *strptr; /* current position within string */
|
|
jbe@0
|
1893 int valid; /* number of valid chars */
|
|
jbe@0
|
1894 int done; /* specifies if latitude or longitude was read */
|
|
jbe@0
|
1895 double val; /* temporary variable */
|
|
jbe@0
|
1896 int lat_count = 0; /* count of latitude values parsed */
|
|
jbe@0
|
1897 int lon_count = 0; /* count of longitufde values parsed */
|
|
jbe@0
|
1898 double lat_min, lat_max, lon_min, lon_max; /* see pgl_box struct */
|
|
jbe@0
|
1899 pgl_box *box; /* return value (to be palloc'ed) */
|
|
jbe@0
|
1900 /* lowercase input */
|
|
jbe@0
|
1901 str_lower = psprintf("%s", str);
|
|
jbe@0
|
1902 for (strptr=str_lower; *strptr; strptr++) {
|
|
jbe@0
|
1903 if (*strptr >= 'A' && *strptr <= 'Z') *strptr += 'a' - 'A';
|
|
jbe@0
|
1904 }
|
|
jbe@0
|
1905 /* reset reading position to start of (lowercase) string */
|
|
jbe@0
|
1906 strptr = str_lower;
|
|
jbe@0
|
1907 /* check if empty box */
|
|
jbe@0
|
1908 valid = 0;
|
|
jbe@0
|
1909 sscanf(strptr, " empty %n", &valid);
|
|
jbe@0
|
1910 if (valid && strptr[valid] == 0) {
|
|
jbe@0
|
1911 /* allocate and return empty box */
|
|
jbe@0
|
1912 box = (pgl_box *)palloc(sizeof(pgl_box));
|
|
jbe@0
|
1913 pgl_box_set_empty(box);
|
|
jbe@0
|
1914 PG_RETURN_POINTER(box);
|
|
jbe@0
|
1915 }
|
|
jbe@0
|
1916 /* demand four blocks separated by whitespace */
|
|
jbe@0
|
1917 valid = 0;
|
|
jbe@0
|
1918 sscanf(strptr, " %*s %*s %*s %*s %n", &valid);
|
|
jbe@0
|
1919 /* if four blocks separated by whitespace exist, parse those blocks */
|
|
jbe@0
|
1920 if (strptr[valid] == 0) while (strptr[0]) {
|
|
jbe@0
|
1921 /* parse either latitude or longitude (whichever found in input string) */
|
|
jbe@0
|
1922 done = pgl_scan(&strptr, &val, &val);
|
|
jbe@0
|
1923 /* store latitude or longitude in lat_min, lat_max, lon_min, or lon_max */
|
|
jbe@0
|
1924 if (done == PGL_SCAN_LAT) {
|
|
jbe@0
|
1925 if (!lat_count) lat_min = val; else lat_max = val;
|
|
jbe@0
|
1926 lat_count++;
|
|
jbe@0
|
1927 } else if (done == PGL_SCAN_LON) {
|
|
jbe@0
|
1928 if (!lon_count) lon_min = val; else lon_max = val;
|
|
jbe@0
|
1929 lon_count++;
|
|
jbe@0
|
1930 } else {
|
|
jbe@0
|
1931 break;
|
|
jbe@0
|
1932 }
|
|
jbe@0
|
1933 }
|
|
jbe@0
|
1934 /* require end of string, and two latitude and two longitude values */
|
|
jbe@0
|
1935 if (strptr[0] || lat_count != 2 || lon_count != 2) {
|
|
jbe@0
|
1936 ereport(ERROR, (
|
|
jbe@0
|
1937 errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
|
|
jbe@0
|
1938 errmsg("invalid input syntax for type ebox: \"%s\"", str)
|
|
jbe@0
|
1939 ));
|
|
jbe@0
|
1940 }
|
|
jbe@0
|
1941 /* free lower case string */
|
|
jbe@0
|
1942 pfree(str_lower);
|
|
jbe@0
|
1943 /* order boundaries (maximum greater than minimum) */
|
|
jbe@0
|
1944 if (lat_min > lat_max) { val = lat_min; lat_min = lat_max; lat_max = val; }
|
|
jbe@0
|
1945 if (lon_min > lon_max) { val = lon_min; lon_min = lon_max; lon_max = val; }
|
|
jbe@0
|
1946 /* allocate memory for result */
|
|
jbe@0
|
1947 box = (pgl_box *)palloc(sizeof(pgl_box));
|
|
jbe@0
|
1948 /* set boundaries (and perform checks) */
|
|
jbe@0
|
1949 pgl_ebox_set_boundaries(box, lat_min, lat_max, lon_min, lon_max);
|
|
jbe@0
|
1950 /* return result */
|
|
jbe@0
|
1951 PG_RETURN_POINTER(box);
|
|
jbe@0
|
1952 }
|
|
jbe@0
|
1953
|
|
jbe@0
|
1954 /* set circle to given latitude, longitude, and radius (including checks) */
|
|
jbe@0
|
1955 static void pgl_ecircle_set_latlon_radius(
|
|
jbe@0
|
1956 pgl_circle *circle, double lat, double lon, double radius
|
|
jbe@0
|
1957 ) {
|
|
jbe@0
|
1958 /* set center point (including checks) */
|
|
jbe@0
|
1959 pgl_epoint_set_latlon(&(circle->center), lat, lon);
|
|
jbe@0
|
1960 /* handle non-positive radius */
|
|
jbe@0
|
1961 if (isnan(radius)) {
|
|
jbe@0
|
1962 ereport(ERROR, (
|
|
jbe@0
|
1963 errcode(ERRCODE_DATA_EXCEPTION),
|
|
jbe@0
|
1964 errmsg("invalid radius for ecircle")
|
|
jbe@0
|
1965 ));
|
|
jbe@0
|
1966 }
|
|
jbe@0
|
1967 if (radius == 0) radius = 0; /* avoids -0 */
|
|
jbe@0
|
1968 else if (radius < 0) {
|
|
jbe@0
|
1969 if (isfinite(radius)) {
|
|
jbe@0
|
1970 ereport(NOTICE, (errmsg("negative radius converted to minus infinity")));
|
|
jbe@0
|
1971 }
|
|
jbe@0
|
1972 radius = -INFINITY;
|
|
jbe@0
|
1973 }
|
|
jbe@0
|
1974 /* store radius (round-trip safety is ensured by pgl_print_float) */
|
|
jbe@0
|
1975 circle->radius = radius;
|
|
jbe@0
|
1976 }
|
|
jbe@0
|
1977
|
|
jbe@0
|
1978 /* create circle ("ecircle" in SQL) from latitude, longitude, and radius */
|
|
jbe@0
|
1979 PG_FUNCTION_INFO_V1(pgl_create_ecircle);
|
|
jbe@0
|
1980 Datum pgl_create_ecircle(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
1981 pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle));
|
|
jbe@0
|
1982 pgl_ecircle_set_latlon_radius(
|
|
jbe@0
|
1983 circle, PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1), PG_GETARG_FLOAT8(2)
|
|
jbe@0
|
1984 );
|
|
jbe@0
|
1985 PG_RETURN_POINTER(circle);
|
|
jbe@0
|
1986 }
|
|
jbe@0
|
1987
|
|
jbe@0
|
1988 /* create circle ("ecircle" in SQL) from point ("epoint"), and radius */
|
|
jbe@0
|
1989 PG_FUNCTION_INFO_V1(pgl_create_ecircle_from_epoint);
|
|
jbe@0
|
1990 Datum pgl_create_ecircle_from_epoint(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
1991 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
|
|
jbe@0
|
1992 double radius = PG_GETARG_FLOAT8(1);
|
|
jbe@0
|
1993 pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle));
|
|
jbe@0
|
1994 /* set latitude, longitude, radius (and perform checks) */
|
|
jbe@0
|
1995 pgl_ecircle_set_latlon_radius(circle, point->lat, point->lon, radius);
|
|
jbe@0
|
1996 /* return result */
|
|
jbe@0
|
1997 PG_RETURN_POINTER(circle);
|
|
jbe@0
|
1998 }
|
|
jbe@0
|
1999
|
|
jbe@0
|
2000 /* parse circle ("ecircle" in SQL) */
|
|
jbe@0
|
2001 /* format: '[NS]<float> [EW]<float> <float>' */
|
|
jbe@0
|
2002 PG_FUNCTION_INFO_V1(pgl_ecircle_in);
|
|
jbe@0
|
2003 Datum pgl_ecircle_in(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2004 char *str = PG_GETARG_CSTRING(0); /* input string */
|
|
jbe@0
|
2005 char *strptr = str; /* current position within string */
|
|
jbe@46
|
2006 double lat, lon, radius; /* parsed values as double precision floats */
|
|
jbe@0
|
2007 int valid = 0; /* number of valid chars */
|
|
jbe@0
|
2008 int done = 0; /* stores if latitude and/or longitude was read */
|
|
jbe@0
|
2009 pgl_circle *circle; /* return value (to be palloc'ed) */
|
|
jbe@0
|
2010 /* demand three blocks separated by whitespace */
|
|
jbe@0
|
2011 sscanf(strptr, " %*s %*s %*s %n", &valid);
|
|
jbe@0
|
2012 /* if three blocks separated by whitespace exist, parse those blocks */
|
|
jbe@0
|
2013 if (strptr[valid] == 0) {
|
|
jbe@0
|
2014 /* parse latitude and longitude */
|
|
jbe@0
|
2015 done |= pgl_scan(&strptr, &lat, &lon);
|
|
jbe@0
|
2016 done |= pgl_scan(&strptr, &lat, &lon);
|
|
jbe@0
|
2017 /* parse radius (while incrementing strptr by number of bytes parsed) */
|
|
jbe@0
|
2018 valid = 0;
|
|
jbe@0
|
2019 if (sscanf(strptr, " %lf %n", &radius, &valid) == 1) strptr += valid;
|
|
jbe@0
|
2020 }
|
|
jbe@0
|
2021 /* require end of string and both latitude and longitude being parsed */
|
|
jbe@0
|
2022 if (strptr[0] || done != PGL_SCAN_LATLON) {
|
|
jbe@0
|
2023 ereport(ERROR, (
|
|
jbe@0
|
2024 errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
|
|
jbe@0
|
2025 errmsg("invalid input syntax for type ecircle: \"%s\"", str)
|
|
jbe@0
|
2026 ));
|
|
jbe@0
|
2027 }
|
|
jbe@0
|
2028 /* allocate memory for result */
|
|
jbe@0
|
2029 circle = (pgl_circle *)palloc(sizeof(pgl_circle));
|
|
jbe@0
|
2030 /* set latitude, longitude, radius (and perform checks) */
|
|
jbe@0
|
2031 pgl_ecircle_set_latlon_radius(circle, lat, lon, radius);
|
|
jbe@0
|
2032 /* return result */
|
|
jbe@0
|
2033 PG_RETURN_POINTER(circle);
|
|
jbe@0
|
2034 }
|
|
jbe@0
|
2035
|
|
jbe@0
|
2036 /* parse cluster ("ecluster" in SQL) */
|
|
jbe@0
|
2037 PG_FUNCTION_INFO_V1(pgl_ecluster_in);
|
|
jbe@0
|
2038 Datum pgl_ecluster_in(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2039 int i;
|
|
jbe@0
|
2040 char *str = PG_GETARG_CSTRING(0); /* input string */
|
|
jbe@0
|
2041 char *str_lower; /* lower case version of input string */
|
|
jbe@0
|
2042 char *strptr; /* pointer to current reading position of input */
|
|
jbe@0
|
2043 int npoints_total = 0; /* total number of points in cluster */
|
|
jbe@0
|
2044 int nentries = 0; /* total number of entries */
|
|
jbe@0
|
2045 pgl_newentry *entries; /* array of pgl_newentry to create pgl_cluster */
|
|
jbe@0
|
2046 int entries_buflen = 4; /* maximum number of elements in entries array */
|
|
jbe@0
|
2047 int valid; /* number of valid chars processed */
|
|
jbe@0
|
2048 double lat, lon; /* latitude and longitude of parsed point */
|
|
jbe@0
|
2049 int entrytype; /* current entry type */
|
|
jbe@0
|
2050 int npoints; /* number of points in current entry */
|
|
jbe@0
|
2051 pgl_point *points; /* array of pgl_point for pgl_newentry */
|
|
jbe@0
|
2052 int points_buflen; /* maximum number of elements in points array */
|
|
jbe@0
|
2053 int done; /* return value of pgl_scan function */
|
|
jbe@0
|
2054 pgl_cluster *cluster; /* created cluster */
|
|
jbe@0
|
2055 /* lowercase input */
|
|
jbe@0
|
2056 str_lower = psprintf("%s", str);
|
|
jbe@0
|
2057 for (strptr=str_lower; *strptr; strptr++) {
|
|
jbe@0
|
2058 if (*strptr >= 'A' && *strptr <= 'Z') *strptr += 'a' - 'A';
|
|
jbe@0
|
2059 }
|
|
jbe@0
|
2060 /* reset reading position to start of (lowercase) string */
|
|
jbe@0
|
2061 strptr = str_lower;
|
|
jbe@0
|
2062 /* allocate initial buffer for entries */
|
|
jbe@0
|
2063 entries = palloc(entries_buflen * sizeof(pgl_newentry));
|
|
jbe@0
|
2064 /* parse until end of string */
|
|
jbe@0
|
2065 while (strptr[0]) {
|
|
jbe@0
|
2066 /* require previous white-space or closing parenthesis before next token */
|
|
jbe@0
|
2067 if (strptr != str_lower && !isspace(strptr[-1]) && strptr[-1] != ')') {
|
|
jbe@0
|
2068 goto pgl_ecluster_in_error;
|
|
jbe@0
|
2069 }
|
|
jbe@0
|
2070 /* ignore token "empty" */
|
|
jbe@0
|
2071 valid = 0; sscanf(strptr, " empty %n", &valid);
|
|
jbe@0
|
2072 if (valid) { strptr += valid; continue; }
|
|
jbe@0
|
2073 /* test for "point" token */
|
|
jbe@0
|
2074 valid = 0; sscanf(strptr, " point ( %n", &valid);
|
|
jbe@0
|
2075 if (valid) {
|
|
jbe@0
|
2076 strptr += valid;
|
|
jbe@0
|
2077 entrytype = PGL_ENTRY_POINT;
|
|
jbe@0
|
2078 goto pgl_ecluster_in_type_ok;
|
|
jbe@0
|
2079 }
|
|
jbe@0
|
2080 /* test for "path" token */
|
|
jbe@0
|
2081 valid = 0; sscanf(strptr, " path ( %n", &valid);
|
|
jbe@0
|
2082 if (valid) {
|
|
jbe@0
|
2083 strptr += valid;
|
|
jbe@0
|
2084 entrytype = PGL_ENTRY_PATH;
|
|
jbe@0
|
2085 goto pgl_ecluster_in_type_ok;
|
|
jbe@0
|
2086 }
|
|
jbe@0
|
2087 /* test for "outline" token */
|
|
jbe@0
|
2088 valid = 0; sscanf(strptr, " outline ( %n", &valid);
|
|
jbe@0
|
2089 if (valid) {
|
|
jbe@0
|
2090 strptr += valid;
|
|
jbe@0
|
2091 entrytype = PGL_ENTRY_OUTLINE;
|
|
jbe@0
|
2092 goto pgl_ecluster_in_type_ok;
|
|
jbe@0
|
2093 }
|
|
jbe@0
|
2094 /* test for "polygon" token */
|
|
jbe@0
|
2095 valid = 0; sscanf(strptr, " polygon ( %n", &valid);
|
|
jbe@0
|
2096 if (valid) {
|
|
jbe@0
|
2097 strptr += valid;
|
|
jbe@0
|
2098 entrytype = PGL_ENTRY_POLYGON;
|
|
jbe@0
|
2099 goto pgl_ecluster_in_type_ok;
|
|
jbe@0
|
2100 }
|
|
jbe@0
|
2101 /* error if no valid token found */
|
|
jbe@0
|
2102 goto pgl_ecluster_in_error;
|
|
jbe@0
|
2103 pgl_ecluster_in_type_ok:
|
|
jbe@0
|
2104 /* check if pgl_newentry array needs to grow */
|
|
jbe@0
|
2105 if (nentries == entries_buflen) {
|
|
jbe@0
|
2106 pgl_newentry *newbuf;
|
|
jbe@0
|
2107 entries_buflen *= 2;
|
|
jbe@0
|
2108 newbuf = palloc(entries_buflen * sizeof(pgl_newentry));
|
|
jbe@0
|
2109 memcpy(newbuf, entries, nentries * sizeof(pgl_newentry));
|
|
jbe@0
|
2110 pfree(entries);
|
|
jbe@0
|
2111 entries = newbuf;
|
|
jbe@0
|
2112 }
|
|
jbe@0
|
2113 /* reset number of points for current entry */
|
|
jbe@0
|
2114 npoints = 0;
|
|
jbe@0
|
2115 /* allocate array for points */
|
|
jbe@0
|
2116 points_buflen = 4;
|
|
jbe@0
|
2117 points = palloc(points_buflen * sizeof(pgl_point));
|
|
jbe@0
|
2118 /* parse until closing parenthesis */
|
|
jbe@0
|
2119 while (strptr[0] != ')') {
|
|
jbe@0
|
2120 /* error on unexpected end of string */
|
|
jbe@0
|
2121 if (strptr[0] == 0) goto pgl_ecluster_in_error;
|
|
jbe@0
|
2122 /* mark neither latitude nor longitude as read */
|
|
jbe@0
|
2123 done = PGL_SCAN_NONE;
|
|
jbe@0
|
2124 /* require white-space before second, third, etc. point */
|
|
jbe@0
|
2125 if (npoints != 0 && !isspace(strptr[-1])) goto pgl_ecluster_in_error;
|
|
jbe@0
|
2126 /* scan latitude (or longitude) */
|
|
jbe@0
|
2127 done |= pgl_scan(&strptr, &lat, &lon);
|
|
jbe@0
|
2128 /* require white-space before second coordinate */
|
|
jbe@0
|
2129 if (strptr != str && !isspace(strptr[-1])) goto pgl_ecluster_in_error;
|
|
jbe@0
|
2130 /* scan longitude (or latitude) */
|
|
jbe@0
|
2131 done |= pgl_scan(&strptr, &lat, &lon);
|
|
jbe@0
|
2132 /* error unless both latitude and longitude were parsed */
|
|
jbe@0
|
2133 if (done != PGL_SCAN_LATLON) goto pgl_ecluster_in_error;
|
|
jbe@0
|
2134 /* throw error if number of points is too high */
|
|
jbe@0
|
2135 if (npoints_total == PGL_CLUSTER_MAXPOINTS) {
|
|
jbe@0
|
2136 ereport(ERROR, (
|
|
jbe@0
|
2137 errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
|
|
jbe@0
|
2138 errmsg(
|
|
jbe@0
|
2139 "too many points for ecluster entry (maximum %i)",
|
|
jbe@0
|
2140 PGL_CLUSTER_MAXPOINTS
|
|
jbe@0
|
2141 )
|
|
jbe@0
|
2142 ));
|
|
jbe@0
|
2143 }
|
|
jbe@0
|
2144 /* check if pgl_point array needs to grow */
|
|
jbe@0
|
2145 if (npoints == points_buflen) {
|
|
jbe@0
|
2146 pgl_point *newbuf;
|
|
jbe@0
|
2147 points_buflen *= 2;
|
|
jbe@0
|
2148 newbuf = palloc(points_buflen * sizeof(pgl_point));
|
|
jbe@0
|
2149 memcpy(newbuf, points, npoints * sizeof(pgl_point));
|
|
jbe@0
|
2150 pfree(points);
|
|
jbe@0
|
2151 points = newbuf;
|
|
jbe@0
|
2152 }
|
|
jbe@0
|
2153 /* append point to pgl_point array (includes checks) */
|
|
jbe@0
|
2154 pgl_epoint_set_latlon(&(points[npoints++]), lat, lon);
|
|
jbe@0
|
2155 /* increase total number of points */
|
|
jbe@0
|
2156 npoints_total++;
|
|
jbe@0
|
2157 }
|
|
jbe@0
|
2158 /* error if entry has no points */
|
|
jbe@0
|
2159 if (!npoints) goto pgl_ecluster_in_error;
|
|
jbe@0
|
2160 /* entries with one point are automatically of type "point" */
|
|
jbe@0
|
2161 if (npoints == 1) entrytype = PGL_ENTRY_POINT;
|
|
jbe@0
|
2162 /* if entries have more than one point */
|
|
jbe@0
|
2163 else {
|
|
jbe@0
|
2164 /* throw error if entry type is "point" */
|
|
jbe@0
|
2165 if (entrytype == PGL_ENTRY_POINT) {
|
|
jbe@0
|
2166 ereport(ERROR, (
|
|
jbe@0
|
2167 errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
|
|
jbe@0
|
2168 errmsg("invalid input syntax for type ecluster (point entry with more than one point)")
|
|
jbe@0
|
2169 ));
|
|
jbe@0
|
2170 }
|
|
jbe@0
|
2171 /* coerce outlines and polygons with more than 2 points to be a path */
|
|
jbe@0
|
2172 if (npoints == 2) entrytype = PGL_ENTRY_PATH;
|
|
jbe@0
|
2173 }
|
|
jbe@0
|
2174 /* append entry to pgl_newentry array */
|
|
jbe@0
|
2175 entries[nentries].entrytype = entrytype;
|
|
jbe@0
|
2176 entries[nentries].npoints = npoints;
|
|
jbe@0
|
2177 entries[nentries].points = points;
|
|
jbe@0
|
2178 nentries++;
|
|
jbe@0
|
2179 /* consume closing parenthesis */
|
|
jbe@0
|
2180 strptr++;
|
|
jbe@0
|
2181 /* consume white-space */
|
|
jbe@0
|
2182 while (isspace(strptr[0])) strptr++;
|
|
jbe@0
|
2183 }
|
|
jbe@0
|
2184 /* free lower case string */
|
|
jbe@0
|
2185 pfree(str_lower);
|
|
jbe@0
|
2186 /* create cluster from pgl_newentry array */
|
|
jbe@0
|
2187 cluster = pgl_new_cluster(nentries, entries);
|
|
jbe@0
|
2188 /* free pgl_newentry array */
|
|
jbe@0
|
2189 for (i=0; i<nentries; i++) pfree(entries[i].points);
|
|
jbe@0
|
2190 pfree(entries);
|
|
jbe@0
|
2191 /* set bounding circle of cluster and check east/west orientation */
|
|
jbe@0
|
2192 if (!pgl_finalize_cluster(cluster)) {
|
|
jbe@0
|
2193 ereport(ERROR, (
|
|
jbe@0
|
2194 errcode(ERRCODE_DATA_EXCEPTION),
|
|
jbe@0
|
2195 errmsg("can not determine east/west orientation for ecluster"),
|
|
jbe@0
|
2196 errhint("Ensure that each entry has a longitude span of less than 180 degrees.")
|
|
jbe@0
|
2197 ));
|
|
jbe@0
|
2198 }
|
|
jbe@0
|
2199 /* return cluster */
|
|
jbe@0
|
2200 PG_RETURN_POINTER(cluster);
|
|
jbe@0
|
2201 /* code to throw error */
|
|
jbe@0
|
2202 pgl_ecluster_in_error:
|
|
jbe@0
|
2203 ereport(ERROR, (
|
|
jbe@0
|
2204 errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
|
|
jbe@0
|
2205 errmsg("invalid input syntax for type ecluster: \"%s\"", str)
|
|
jbe@0
|
2206 ));
|
|
jbe@0
|
2207 }
|
|
jbe@0
|
2208
|
|
jbe@0
|
2209 /* convert point ("epoint") to string representation */
|
|
jbe@0
|
2210 PG_FUNCTION_INFO_V1(pgl_epoint_out);
|
|
jbe@0
|
2211 Datum pgl_epoint_out(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2212 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
|
|
jbe@0
|
2213 char latstr[PGL_NUMBUFLEN];
|
|
jbe@0
|
2214 char lonstr[PGL_NUMBUFLEN];
|
|
jbe@0
|
2215 pgl_print_lat(latstr, point->lat);
|
|
jbe@0
|
2216 pgl_print_lon(lonstr, point->lon);
|
|
jbe@0
|
2217 PG_RETURN_CSTRING(psprintf("%s %s", latstr, lonstr));
|
|
jbe@0
|
2218 }
|
|
jbe@0
|
2219
|
|
jbe@46
|
2220 /* convert point with sample count ("epoint_with_sample_count") to str. rep. */
|
|
jbe@46
|
2221 PG_FUNCTION_INFO_V1(pgl_epoint_with_sample_count_out);
|
|
jbe@46
|
2222 Datum pgl_epoint_with_sample_count_out(PG_FUNCTION_ARGS) {
|
|
jbe@46
|
2223 pgl_point_sc *search = (pgl_point_sc *)PG_GETARG_POINTER(0);
|
|
jbe@46
|
2224 char latstr[PGL_NUMBUFLEN];
|
|
jbe@46
|
2225 char lonstr[PGL_NUMBUFLEN];
|
|
jbe@46
|
2226 pgl_print_lat(latstr, search->point.lat);
|
|
jbe@46
|
2227 pgl_print_lon(lonstr, search->point.lon);
|
|
jbe@46
|
2228 PG_RETURN_CSTRING(
|
|
jbe@46
|
2229 psprintf("%s %s %i", latstr, lonstr, (int)search->samples)
|
|
jbe@46
|
2230 );
|
|
jbe@46
|
2231 }
|
|
jbe@46
|
2232
|
|
jbe@0
|
2233 /* convert box ("ebox") to string representation */
|
|
jbe@0
|
2234 PG_FUNCTION_INFO_V1(pgl_ebox_out);
|
|
jbe@0
|
2235 Datum pgl_ebox_out(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2236 pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
|
|
jbe@0
|
2237 double lon_min = box->lon_min;
|
|
jbe@0
|
2238 double lon_max = box->lon_max;
|
|
jbe@0
|
2239 char lat_min_str[PGL_NUMBUFLEN];
|
|
jbe@0
|
2240 char lat_max_str[PGL_NUMBUFLEN];
|
|
jbe@0
|
2241 char lon_min_str[PGL_NUMBUFLEN];
|
|
jbe@0
|
2242 char lon_max_str[PGL_NUMBUFLEN];
|
|
jbe@0
|
2243 /* return string "empty" if box is set to be empty */
|
|
jbe@0
|
2244 if (box->lat_min > box->lat_max) PG_RETURN_CSTRING("empty");
|
|
jbe@0
|
2245 /* use boundaries exceeding W180 or E180 if 180th meridian is enclosed */
|
|
jbe@0
|
2246 /* (required since pgl_box_in orders the longitude boundaries) */
|
|
jbe@0
|
2247 if (lon_min > lon_max) {
|
|
jbe@0
|
2248 if (lon_min + lon_max >= 0) lon_min -= 360;
|
|
jbe@0
|
2249 else lon_max += 360;
|
|
jbe@0
|
2250 }
|
|
jbe@0
|
2251 /* format and return result */
|
|
jbe@0
|
2252 pgl_print_lat(lat_min_str, box->lat_min);
|
|
jbe@0
|
2253 pgl_print_lat(lat_max_str, box->lat_max);
|
|
jbe@0
|
2254 pgl_print_lon(lon_min_str, lon_min);
|
|
jbe@0
|
2255 pgl_print_lon(lon_max_str, lon_max);
|
|
jbe@0
|
2256 PG_RETURN_CSTRING(psprintf(
|
|
jbe@0
|
2257 "%s %s %s %s",
|
|
jbe@0
|
2258 lat_min_str, lon_min_str, lat_max_str, lon_max_str
|
|
jbe@0
|
2259 ));
|
|
jbe@0
|
2260 }
|
|
jbe@0
|
2261
|
|
jbe@0
|
2262 /* convert circle ("ecircle") to string representation */
|
|
jbe@0
|
2263 PG_FUNCTION_INFO_V1(pgl_ecircle_out);
|
|
jbe@0
|
2264 Datum pgl_ecircle_out(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2265 pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
|
|
jbe@0
|
2266 char latstr[PGL_NUMBUFLEN];
|
|
jbe@0
|
2267 char lonstr[PGL_NUMBUFLEN];
|
|
jbe@0
|
2268 char radstr[PGL_NUMBUFLEN];
|
|
jbe@0
|
2269 pgl_print_lat(latstr, circle->center.lat);
|
|
jbe@0
|
2270 pgl_print_lon(lonstr, circle->center.lon);
|
|
jbe@0
|
2271 pgl_print_float(radstr, circle->radius);
|
|
jbe@0
|
2272 PG_RETURN_CSTRING(psprintf("%s %s %s", latstr, lonstr, radstr));
|
|
jbe@0
|
2273 }
|
|
jbe@0
|
2274
|
|
jbe@0
|
2275 /* convert cluster ("ecluster") to string representation */
|
|
jbe@0
|
2276 PG_FUNCTION_INFO_V1(pgl_ecluster_out);
|
|
jbe@0
|
2277 Datum pgl_ecluster_out(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2278 pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
|
|
jbe@0
|
2279 char latstr[PGL_NUMBUFLEN]; /* string buffer for latitude */
|
|
jbe@0
|
2280 char lonstr[PGL_NUMBUFLEN]; /* string buffer for longitude */
|
|
jbe@0
|
2281 char ***strings; /* array of array of strings */
|
|
jbe@0
|
2282 char *string; /* string of current token */
|
|
jbe@0
|
2283 char *res, *resptr; /* result and pointer to current write position */
|
|
jbe@0
|
2284 size_t reslen = 1; /* length of result (init with 1 for terminator) */
|
|
jbe@0
|
2285 int npoints; /* number of points of current entry */
|
|
jbe@0
|
2286 int i, j; /* i: entry, j: point in entry */
|
|
jbe@0
|
2287 /* handle empty clusters */
|
|
jbe@0
|
2288 if (cluster->nentries == 0) {
|
|
jbe@0
|
2289 /* free detoasted cluster (if copy) */
|
|
jbe@0
|
2290 PG_FREE_IF_COPY(cluster, 0);
|
|
jbe@0
|
2291 /* return static result */
|
|
jbe@0
|
2292 PG_RETURN_CSTRING("empty");
|
|
jbe@0
|
2293 }
|
|
jbe@0
|
2294 /* allocate array of array of strings */
|
|
jbe@0
|
2295 strings = palloc(cluster->nentries * sizeof(char **));
|
|
jbe@0
|
2296 /* iterate over all entries in cluster */
|
|
jbe@0
|
2297 for (i=0; i<cluster->nentries; i++) {
|
|
jbe@0
|
2298 /* get number of points in entry */
|
|
jbe@0
|
2299 npoints = cluster->entries[i].npoints;
|
|
jbe@0
|
2300 /* allocate array of strings (one string for each point plus two extra) */
|
|
jbe@0
|
2301 strings[i] = palloc((2 + npoints) * sizeof(char *));
|
|
jbe@0
|
2302 /* determine opening string */
|
|
jbe@0
|
2303 switch (cluster->entries[i].entrytype) {
|
|
jbe@0
|
2304 case PGL_ENTRY_POINT: string = (i==0)?"point (" :" point ("; break;
|
|
jbe@0
|
2305 case PGL_ENTRY_PATH: string = (i==0)?"path (" :" path ("; break;
|
|
jbe@0
|
2306 case PGL_ENTRY_OUTLINE: string = (i==0)?"outline (":" outline ("; break;
|
|
jbe@0
|
2307 case PGL_ENTRY_POLYGON: string = (i==0)?"polygon (":" polygon ("; break;
|
|
jbe@0
|
2308 default: string = (i==0)?"unknown" :" unknown";
|
|
jbe@0
|
2309 }
|
|
jbe@0
|
2310 /* use opening string as first string in array */
|
|
jbe@0
|
2311 strings[i][0] = string;
|
|
jbe@0
|
2312 /* update result length (for allocating result string later) */
|
|
jbe@0
|
2313 reslen += strlen(string);
|
|
jbe@0
|
2314 /* iterate over all points */
|
|
jbe@0
|
2315 for (j=0; j<npoints; j++) {
|
|
jbe@0
|
2316 /* create string representation of point */
|
|
jbe@0
|
2317 pgl_print_lat(latstr, PGL_ENTRY_POINTS(cluster, i)[j].lat);
|
|
jbe@0
|
2318 pgl_print_lon(lonstr, PGL_ENTRY_POINTS(cluster, i)[j].lon);
|
|
jbe@0
|
2319 string = psprintf((j == 0) ? "%s %s" : " %s %s", latstr, lonstr);
|
|
jbe@0
|
2320 /* copy string pointer to string array */
|
|
jbe@0
|
2321 strings[i][j+1] = string;
|
|
jbe@0
|
2322 /* update result length (for allocating result string later) */
|
|
jbe@0
|
2323 reslen += strlen(string);
|
|
jbe@0
|
2324 }
|
|
jbe@0
|
2325 /* use closing parenthesis as last string in array */
|
|
jbe@0
|
2326 strings[i][npoints+1] = ")";
|
|
jbe@0
|
2327 /* update result length (for allocating result string later) */
|
|
jbe@0
|
2328 reslen++;
|
|
jbe@0
|
2329 }
|
|
jbe@0
|
2330 /* allocate result string */
|
|
jbe@0
|
2331 res = palloc(reslen);
|
|
jbe@0
|
2332 /* set write pointer to begin of result string */
|
|
jbe@0
|
2333 resptr = res;
|
|
jbe@0
|
2334 /* copy strings into result string */
|
|
jbe@0
|
2335 for (i=0; i<cluster->nentries; i++) {
|
|
jbe@0
|
2336 npoints = cluster->entries[i].npoints;
|
|
jbe@0
|
2337 for (j=0; j<npoints+2; j++) {
|
|
jbe@0
|
2338 string = strings[i][j];
|
|
jbe@0
|
2339 strcpy(resptr, string);
|
|
jbe@0
|
2340 resptr += strlen(string);
|
|
jbe@0
|
2341 /* free strings allocated by psprintf */
|
|
jbe@0
|
2342 if (j != 0 && j != npoints+1) pfree(string);
|
|
jbe@0
|
2343 }
|
|
jbe@0
|
2344 /* free array of strings */
|
|
jbe@0
|
2345 pfree(strings[i]);
|
|
jbe@0
|
2346 }
|
|
jbe@0
|
2347 /* free array of array of strings */
|
|
jbe@0
|
2348 pfree(strings);
|
|
jbe@0
|
2349 /* free detoasted cluster (if copy) */
|
|
jbe@0
|
2350 PG_FREE_IF_COPY(cluster, 0);
|
|
jbe@0
|
2351 /* return result */
|
|
jbe@0
|
2352 PG_RETURN_CSTRING(res);
|
|
jbe@0
|
2353 }
|
|
jbe@0
|
2354
|
|
jbe@0
|
2355 /* binary input function for point ("epoint") */
|
|
jbe@0
|
2356 PG_FUNCTION_INFO_V1(pgl_epoint_recv);
|
|
jbe@0
|
2357 Datum pgl_epoint_recv(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2358 StringInfo buf = (StringInfo)PG_GETARG_POINTER(0);
|
|
jbe@0
|
2359 pgl_point *point = (pgl_point *)palloc(sizeof(pgl_point));
|
|
jbe@0
|
2360 point->lat = pq_getmsgfloat8(buf);
|
|
jbe@0
|
2361 point->lon = pq_getmsgfloat8(buf);
|
|
jbe@0
|
2362 PG_RETURN_POINTER(point);
|
|
jbe@0
|
2363 }
|
|
jbe@0
|
2364
|
|
jbe@0
|
2365 /* binary input function for box ("ebox") */
|
|
jbe@0
|
2366 PG_FUNCTION_INFO_V1(pgl_ebox_recv);
|
|
jbe@0
|
2367 Datum pgl_ebox_recv(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2368 StringInfo buf = (StringInfo)PG_GETARG_POINTER(0);
|
|
jbe@0
|
2369 pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box));
|
|
jbe@0
|
2370 box->lat_min = pq_getmsgfloat8(buf);
|
|
jbe@0
|
2371 box->lat_max = pq_getmsgfloat8(buf);
|
|
jbe@0
|
2372 box->lon_min = pq_getmsgfloat8(buf);
|
|
jbe@0
|
2373 box->lon_max = pq_getmsgfloat8(buf);
|
|
jbe@0
|
2374 PG_RETURN_POINTER(box);
|
|
jbe@0
|
2375 }
|
|
jbe@0
|
2376
|
|
jbe@0
|
2377 /* binary input function for circle ("ecircle") */
|
|
jbe@0
|
2378 PG_FUNCTION_INFO_V1(pgl_ecircle_recv);
|
|
jbe@0
|
2379 Datum pgl_ecircle_recv(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2380 StringInfo buf = (StringInfo)PG_GETARG_POINTER(0);
|
|
jbe@0
|
2381 pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle));
|
|
jbe@0
|
2382 circle->center.lat = pq_getmsgfloat8(buf);
|
|
jbe@0
|
2383 circle->center.lon = pq_getmsgfloat8(buf);
|
|
jbe@0
|
2384 circle->radius = pq_getmsgfloat8(buf);
|
|
jbe@0
|
2385 PG_RETURN_POINTER(circle);
|
|
jbe@0
|
2386 }
|
|
jbe@0
|
2387
|
|
jbe@0
|
2388 /* TODO: binary receive function for cluster */
|
|
jbe@0
|
2389
|
|
jbe@0
|
2390 /* binary output function for point ("epoint") */
|
|
jbe@0
|
2391 PG_FUNCTION_INFO_V1(pgl_epoint_send);
|
|
jbe@0
|
2392 Datum pgl_epoint_send(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2393 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
|
|
jbe@0
|
2394 StringInfoData buf;
|
|
jbe@0
|
2395 pq_begintypsend(&buf);
|
|
jbe@0
|
2396 pq_sendfloat8(&buf, point->lat);
|
|
jbe@0
|
2397 pq_sendfloat8(&buf, point->lon);
|
|
jbe@0
|
2398 PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
|
|
jbe@0
|
2399 }
|
|
jbe@0
|
2400
|
|
jbe@0
|
2401 /* binary output function for box ("ebox") */
|
|
jbe@0
|
2402 PG_FUNCTION_INFO_V1(pgl_ebox_send);
|
|
jbe@0
|
2403 Datum pgl_ebox_send(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2404 pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
|
|
jbe@0
|
2405 StringInfoData buf;
|
|
jbe@0
|
2406 pq_begintypsend(&buf);
|
|
jbe@0
|
2407 pq_sendfloat8(&buf, box->lat_min);
|
|
jbe@0
|
2408 pq_sendfloat8(&buf, box->lat_max);
|
|
jbe@0
|
2409 pq_sendfloat8(&buf, box->lon_min);
|
|
jbe@0
|
2410 pq_sendfloat8(&buf, box->lon_max);
|
|
jbe@0
|
2411 PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
|
|
jbe@0
|
2412 }
|
|
jbe@0
|
2413
|
|
jbe@0
|
2414 /* binary output function for circle ("ecircle") */
|
|
jbe@0
|
2415 PG_FUNCTION_INFO_V1(pgl_ecircle_send);
|
|
jbe@0
|
2416 Datum pgl_ecircle_send(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2417 pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
|
|
jbe@0
|
2418 StringInfoData buf;
|
|
jbe@0
|
2419 pq_begintypsend(&buf);
|
|
jbe@0
|
2420 pq_sendfloat8(&buf, circle->center.lat);
|
|
jbe@0
|
2421 pq_sendfloat8(&buf, circle->center.lon);
|
|
jbe@0
|
2422 pq_sendfloat8(&buf, circle->radius);
|
|
jbe@0
|
2423 PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
|
|
jbe@0
|
2424 }
|
|
jbe@0
|
2425
|
|
jbe@0
|
2426 /* TODO: binary send functions for cluster */
|
|
jbe@0
|
2427
|
|
jbe@0
|
2428 /* cast point ("epoint") to box ("ebox") */
|
|
jbe@0
|
2429 PG_FUNCTION_INFO_V1(pgl_epoint_to_ebox);
|
|
jbe@0
|
2430 Datum pgl_epoint_to_ebox(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2431 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
|
|
jbe@0
|
2432 pgl_box *box = palloc(sizeof(pgl_box));
|
|
jbe@0
|
2433 box->lat_min = point->lat;
|
|
jbe@0
|
2434 box->lat_max = point->lat;
|
|
jbe@0
|
2435 box->lon_min = point->lon;
|
|
jbe@0
|
2436 box->lon_max = point->lon;
|
|
jbe@0
|
2437 PG_RETURN_POINTER(box);
|
|
jbe@0
|
2438 }
|
|
jbe@0
|
2439
|
|
jbe@0
|
2440 /* cast point ("epoint") to circle ("ecircle") */
|
|
jbe@0
|
2441 PG_FUNCTION_INFO_V1(pgl_epoint_to_ecircle);
|
|
jbe@0
|
2442 Datum pgl_epoint_to_ecircle(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2443 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
|
|
jbe@0
|
2444 pgl_circle *circle = palloc(sizeof(pgl_box));
|
|
jbe@0
|
2445 circle->center = *point;
|
|
jbe@0
|
2446 circle->radius = 0;
|
|
jbe@0
|
2447 PG_RETURN_POINTER(circle);
|
|
jbe@0
|
2448 }
|
|
jbe@0
|
2449
|
|
jbe@0
|
2450 /* cast point ("epoint") to cluster ("ecluster") */
|
|
jbe@0
|
2451 PG_FUNCTION_INFO_V1(pgl_epoint_to_ecluster);
|
|
jbe@0
|
2452 Datum pgl_epoint_to_ecluster(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2453 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
|
|
jbe@0
|
2454 pgl_newentry entry;
|
|
jbe@42
|
2455 pgl_cluster *cluster;
|
|
jbe@0
|
2456 entry.entrytype = PGL_ENTRY_POINT;
|
|
jbe@0
|
2457 entry.npoints = 1;
|
|
jbe@0
|
2458 entry.points = point;
|
|
jbe@42
|
2459 cluster = pgl_new_cluster(1, &entry);
|
|
jbe@42
|
2460 pgl_finalize_cluster(cluster); /* NOTE: should not fail */
|
|
jbe@42
|
2461 PG_RETURN_POINTER(cluster);
|
|
jbe@0
|
2462 }
|
|
jbe@0
|
2463
|
|
jbe@0
|
2464 /* cast box ("ebox") to cluster ("ecluster") */
|
|
jbe@0
|
2465 #define pgl_ebox_to_ecluster_macro(i, a, b) \
|
|
jbe@0
|
2466 entries[i].entrytype = PGL_ENTRY_POLYGON; \
|
|
jbe@0
|
2467 entries[i].npoints = 4; \
|
|
jbe@0
|
2468 entries[i].points = points[i]; \
|
|
jbe@0
|
2469 points[i][0].lat = box->lat_min; \
|
|
jbe@0
|
2470 points[i][0].lon = (a); \
|
|
jbe@0
|
2471 points[i][1].lat = box->lat_min; \
|
|
jbe@0
|
2472 points[i][1].lon = (b); \
|
|
jbe@0
|
2473 points[i][2].lat = box->lat_max; \
|
|
jbe@0
|
2474 points[i][2].lon = (b); \
|
|
jbe@0
|
2475 points[i][3].lat = box->lat_max; \
|
|
jbe@0
|
2476 points[i][3].lon = (a);
|
|
jbe@0
|
2477 PG_FUNCTION_INFO_V1(pgl_ebox_to_ecluster);
|
|
jbe@0
|
2478 Datum pgl_ebox_to_ecluster(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2479 pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
|
|
jbe@0
|
2480 double lon, dlon;
|
|
jbe@0
|
2481 int nentries;
|
|
jbe@0
|
2482 pgl_newentry entries[3];
|
|
jbe@0
|
2483 pgl_point points[3][4];
|
|
jbe@42
|
2484 pgl_cluster *cluster;
|
|
jbe@0
|
2485 if (box->lat_min > box->lat_max) {
|
|
jbe@0
|
2486 nentries = 0;
|
|
jbe@0
|
2487 } else if (box->lon_min > box->lon_max) {
|
|
jbe@0
|
2488 if (box->lon_min < 0) {
|
|
jbe@0
|
2489 lon = pgl_round((box->lon_min + 180) / 2.0);
|
|
jbe@0
|
2490 nentries = 3;
|
|
jbe@0
|
2491 pgl_ebox_to_ecluster_macro(0, box->lon_min, lon);
|
|
jbe@0
|
2492 pgl_ebox_to_ecluster_macro(1, lon, 180);
|
|
jbe@0
|
2493 pgl_ebox_to_ecluster_macro(2, -180, box->lon_max);
|
|
jbe@0
|
2494 } else if (box->lon_max > 0) {
|
|
jbe@0
|
2495 lon = pgl_round((box->lon_max - 180) / 2.0);
|
|
jbe@0
|
2496 nentries = 3;
|
|
jbe@0
|
2497 pgl_ebox_to_ecluster_macro(0, box->lon_min, 180);
|
|
jbe@0
|
2498 pgl_ebox_to_ecluster_macro(1, -180, lon);
|
|
jbe@0
|
2499 pgl_ebox_to_ecluster_macro(2, lon, box->lon_max);
|
|
jbe@0
|
2500 } else {
|
|
jbe@0
|
2501 nentries = 2;
|
|
jbe@0
|
2502 pgl_ebox_to_ecluster_macro(0, box->lon_min, 180);
|
|
jbe@0
|
2503 pgl_ebox_to_ecluster_macro(1, -180, box->lon_max);
|
|
jbe@0
|
2504 }
|
|
jbe@0
|
2505 } else {
|
|
jbe@0
|
2506 dlon = pgl_round(box->lon_max - box->lon_min);
|
|
jbe@0
|
2507 if (dlon < 180) {
|
|
jbe@0
|
2508 nentries = 1;
|
|
jbe@0
|
2509 pgl_ebox_to_ecluster_macro(0, box->lon_min, box->lon_max);
|
|
jbe@0
|
2510 } else {
|
|
jbe@0
|
2511 lon = pgl_round((box->lon_min + box->lon_max) / 2.0);
|
|
jbe@0
|
2512 if (
|
|
jbe@0
|
2513 pgl_round(lon - box->lon_min) < 180 &&
|
|
jbe@0
|
2514 pgl_round(box->lon_max - lon) < 180
|
|
jbe@0
|
2515 ) {
|
|
jbe@0
|
2516 nentries = 2;
|
|
jbe@0
|
2517 pgl_ebox_to_ecluster_macro(0, box->lon_min, lon);
|
|
jbe@0
|
2518 pgl_ebox_to_ecluster_macro(1, lon, box->lon_max);
|
|
jbe@0
|
2519 } else {
|
|
jbe@0
|
2520 nentries = 3;
|
|
jbe@0
|
2521 pgl_ebox_to_ecluster_macro(0, box->lon_min, -60);
|
|
jbe@0
|
2522 pgl_ebox_to_ecluster_macro(1, -60, 60);
|
|
jbe@0
|
2523 pgl_ebox_to_ecluster_macro(2, 60, box->lon_max);
|
|
jbe@0
|
2524 }
|
|
jbe@0
|
2525 }
|
|
jbe@0
|
2526 }
|
|
jbe@42
|
2527 cluster = pgl_new_cluster(nentries, entries);
|
|
jbe@42
|
2528 pgl_finalize_cluster(cluster); /* NOTE: should not fail */
|
|
jbe@42
|
2529 PG_RETURN_POINTER(cluster);
|
|
jbe@0
|
2530 }
|
|
jbe@0
|
2531
|
|
jbe@0
|
2532 /* extract latitude from point ("epoint") */
|
|
jbe@0
|
2533 PG_FUNCTION_INFO_V1(pgl_epoint_lat);
|
|
jbe@0
|
2534 Datum pgl_epoint_lat(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2535 PG_RETURN_FLOAT8(((pgl_point *)PG_GETARG_POINTER(0))->lat);
|
|
jbe@0
|
2536 }
|
|
jbe@0
|
2537
|
|
jbe@0
|
2538 /* extract longitude from point ("epoint") */
|
|
jbe@0
|
2539 PG_FUNCTION_INFO_V1(pgl_epoint_lon);
|
|
jbe@0
|
2540 Datum pgl_epoint_lon(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2541 PG_RETURN_FLOAT8(((pgl_point *)PG_GETARG_POINTER(0))->lon);
|
|
jbe@0
|
2542 }
|
|
jbe@0
|
2543
|
|
jbe@0
|
2544 /* extract minimum latitude from box ("ebox") */
|
|
jbe@0
|
2545 PG_FUNCTION_INFO_V1(pgl_ebox_lat_min);
|
|
jbe@0
|
2546 Datum pgl_ebox_lat_min(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2547 PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lat_min);
|
|
jbe@0
|
2548 }
|
|
jbe@0
|
2549
|
|
jbe@0
|
2550 /* extract maximum latitude from box ("ebox") */
|
|
jbe@0
|
2551 PG_FUNCTION_INFO_V1(pgl_ebox_lat_max);
|
|
jbe@0
|
2552 Datum pgl_ebox_lat_max(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2553 PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lat_max);
|
|
jbe@0
|
2554 }
|
|
jbe@0
|
2555
|
|
jbe@0
|
2556 /* extract minimum longitude from box ("ebox") */
|
|
jbe@0
|
2557 PG_FUNCTION_INFO_V1(pgl_ebox_lon_min);
|
|
jbe@0
|
2558 Datum pgl_ebox_lon_min(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2559 PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lon_min);
|
|
jbe@0
|
2560 }
|
|
jbe@0
|
2561
|
|
jbe@0
|
2562 /* extract maximum longitude from box ("ebox") */
|
|
jbe@0
|
2563 PG_FUNCTION_INFO_V1(pgl_ebox_lon_max);
|
|
jbe@0
|
2564 Datum pgl_ebox_lon_max(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2565 PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lon_max);
|
|
jbe@0
|
2566 }
|
|
jbe@0
|
2567
|
|
jbe@0
|
2568 /* extract center point from circle ("ecircle") */
|
|
jbe@0
|
2569 PG_FUNCTION_INFO_V1(pgl_ecircle_center);
|
|
jbe@0
|
2570 Datum pgl_ecircle_center(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2571 PG_RETURN_POINTER(&(((pgl_circle *)PG_GETARG_POINTER(0))->center));
|
|
jbe@0
|
2572 }
|
|
jbe@0
|
2573
|
|
jbe@0
|
2574 /* extract radius from circle ("ecircle") */
|
|
jbe@0
|
2575 PG_FUNCTION_INFO_V1(pgl_ecircle_radius);
|
|
jbe@0
|
2576 Datum pgl_ecircle_radius(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2577 PG_RETURN_FLOAT8(((pgl_circle *)PG_GETARG_POINTER(0))->radius);
|
|
jbe@0
|
2578 }
|
|
jbe@0
|
2579
|
|
jbe@0
|
2580 /* check if point is inside box (overlap operator "&&") in SQL */
|
|
jbe@0
|
2581 PG_FUNCTION_INFO_V1(pgl_epoint_ebox_overlap);
|
|
jbe@0
|
2582 Datum pgl_epoint_ebox_overlap(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2583 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
|
|
jbe@0
|
2584 pgl_box *box = (pgl_box *)PG_GETARG_POINTER(1);
|
|
jbe@0
|
2585 PG_RETURN_BOOL(pgl_point_in_box(point, box));
|
|
jbe@0
|
2586 }
|
|
jbe@0
|
2587
|
|
jbe@0
|
2588 /* check if point is inside circle (overlap operator "&&") in SQL */
|
|
jbe@0
|
2589 PG_FUNCTION_INFO_V1(pgl_epoint_ecircle_overlap);
|
|
jbe@0
|
2590 Datum pgl_epoint_ecircle_overlap(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2591 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
|
|
jbe@0
|
2592 pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(1);
|
|
jbe@0
|
2593 PG_RETURN_BOOL(
|
|
jbe@0
|
2594 pgl_distance(
|
|
jbe@0
|
2595 point->lat, point->lon,
|
|
jbe@0
|
2596 circle->center.lat, circle->center.lon
|
|
jbe@0
|
2597 ) <= circle->radius
|
|
jbe@0
|
2598 );
|
|
jbe@0
|
2599 }
|
|
jbe@0
|
2600
|
|
jbe@0
|
2601 /* check if point is inside cluster (overlap operator "&&") in SQL */
|
|
jbe@0
|
2602 PG_FUNCTION_INFO_V1(pgl_epoint_ecluster_overlap);
|
|
jbe@0
|
2603 Datum pgl_epoint_ecluster_overlap(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2604 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
|
|
jbe@0
|
2605 pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
|
|
jbe@16
|
2606 bool retval;
|
|
jbe@46
|
2607 /* points outside bounding circle are always assumed to be non-overlapping */
|
|
jbe@16
|
2608 if (
|
|
jbe@16
|
2609 pgl_distance(
|
|
jbe@16
|
2610 point->lat, point->lon,
|
|
jbe@16
|
2611 cluster->bounding.center.lat, cluster->bounding.center.lon
|
|
jbe@16
|
2612 ) > cluster->bounding.radius
|
|
jbe@16
|
2613 ) retval = false;
|
|
jbe@20
|
2614 else retval = pgl_point_in_cluster(point, cluster, false);
|
|
jbe@0
|
2615 PG_FREE_IF_COPY(cluster, 1);
|
|
jbe@0
|
2616 PG_RETURN_BOOL(retval);
|
|
jbe@0
|
2617 }
|
|
jbe@0
|
2618
|
|
jbe@10
|
2619 /* check if point may be inside cluster (lossy overl. operator "&&+") in SQL */
|
|
jbe@10
|
2620 PG_FUNCTION_INFO_V1(pgl_epoint_ecluster_may_overlap);
|
|
jbe@10
|
2621 Datum pgl_epoint_ecluster_may_overlap(PG_FUNCTION_ARGS) {
|
|
jbe@10
|
2622 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
|
|
jbe@10
|
2623 pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
|
|
jbe@10
|
2624 bool retval = pgl_distance(
|
|
jbe@10
|
2625 point->lat, point->lon,
|
|
jbe@10
|
2626 cluster->bounding.center.lat, cluster->bounding.center.lon
|
|
jbe@10
|
2627 ) <= cluster->bounding.radius;
|
|
jbe@10
|
2628 PG_FREE_IF_COPY(cluster, 1);
|
|
jbe@10
|
2629 PG_RETURN_BOOL(retval);
|
|
jbe@10
|
2630 }
|
|
jbe@10
|
2631
|
|
jbe@0
|
2632 /* check if two boxes overlap (overlap operator "&&") in SQL */
|
|
jbe@0
|
2633 PG_FUNCTION_INFO_V1(pgl_ebox_overlap);
|
|
jbe@0
|
2634 Datum pgl_ebox_overlap(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2635 pgl_box *box1 = (pgl_box *)PG_GETARG_POINTER(0);
|
|
jbe@0
|
2636 pgl_box *box2 = (pgl_box *)PG_GETARG_POINTER(1);
|
|
jbe@0
|
2637 PG_RETURN_BOOL(pgl_boxes_overlap(box1, box2));
|
|
jbe@0
|
2638 }
|
|
jbe@0
|
2639
|
|
jbe@10
|
2640 /* check if box and circle may overlap (lossy overl. operator "&&+") in SQL */
|
|
jbe@10
|
2641 PG_FUNCTION_INFO_V1(pgl_ebox_ecircle_may_overlap);
|
|
jbe@10
|
2642 Datum pgl_ebox_ecircle_may_overlap(PG_FUNCTION_ARGS) {
|
|
jbe@10
|
2643 pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
|
|
jbe@10
|
2644 pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(1);
|
|
jbe@10
|
2645 PG_RETURN_BOOL(
|
|
jbe@10
|
2646 pgl_estimate_point_box_distance(&circle->center, box) <= circle->radius
|
|
jbe@10
|
2647 );
|
|
jbe@10
|
2648 }
|
|
jbe@10
|
2649
|
|
jbe@10
|
2650 /* check if box and cluster may overlap (lossy overl. operator "&&+") in SQL */
|
|
jbe@10
|
2651 PG_FUNCTION_INFO_V1(pgl_ebox_ecluster_may_overlap);
|
|
jbe@10
|
2652 Datum pgl_ebox_ecluster_may_overlap(PG_FUNCTION_ARGS) {
|
|
jbe@10
|
2653 pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
|
|
jbe@10
|
2654 pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
|
|
jbe@10
|
2655 bool retval = pgl_estimate_point_box_distance(
|
|
jbe@10
|
2656 &cluster->bounding.center,
|
|
jbe@10
|
2657 box
|
|
jbe@10
|
2658 ) <= cluster->bounding.radius;
|
|
jbe@10
|
2659 PG_FREE_IF_COPY(cluster, 1);
|
|
jbe@10
|
2660 PG_RETURN_BOOL(retval);
|
|
jbe@10
|
2661 }
|
|
jbe@10
|
2662
|
|
jbe@0
|
2663 /* check if two circles overlap (overlap operator "&&") in SQL */
|
|
jbe@0
|
2664 PG_FUNCTION_INFO_V1(pgl_ecircle_overlap);
|
|
jbe@0
|
2665 Datum pgl_ecircle_overlap(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2666 pgl_circle *circle1 = (pgl_circle *)PG_GETARG_POINTER(0);
|
|
jbe@0
|
2667 pgl_circle *circle2 = (pgl_circle *)PG_GETARG_POINTER(1);
|
|
jbe@0
|
2668 PG_RETURN_BOOL(
|
|
jbe@0
|
2669 pgl_distance(
|
|
jbe@0
|
2670 circle1->center.lat, circle1->center.lon,
|
|
jbe@0
|
2671 circle2->center.lat, circle2->center.lon
|
|
jbe@0
|
2672 ) <= circle1->radius + circle2->radius
|
|
jbe@0
|
2673 );
|
|
jbe@0
|
2674 }
|
|
jbe@0
|
2675
|
|
jbe@0
|
2676 /* check if circle and cluster overlap (overlap operator "&&") in SQL */
|
|
jbe@0
|
2677 PG_FUNCTION_INFO_V1(pgl_ecircle_ecluster_overlap);
|
|
jbe@0
|
2678 Datum pgl_ecircle_ecluster_overlap(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2679 pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
|
|
jbe@0
|
2680 pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
|
|
jbe@46
|
2681 bool retval;
|
|
jbe@46
|
2682 /* points outside bounding circle (with margin for flattening) are always
|
|
jbe@46
|
2683 assumed to be non-overlapping */
|
|
jbe@46
|
2684 if (
|
|
jbe@46
|
2685 (1.0-PGL_SPHEROID_F) * /* margin for flattening and approximation */
|
|
jbe@46
|
2686 pgl_distance(
|
|
jbe@46
|
2687 circle->center.lat, circle->center.lon,
|
|
jbe@46
|
2688 cluster->bounding.center.lat, cluster->bounding.center.lon
|
|
jbe@46
|
2689 ) > circle->radius + cluster->bounding.radius
|
|
jbe@46
|
2690 ) retval = false;
|
|
jbe@46
|
2691 else retval = pgl_point_cluster_distance(&(circle->center), cluster) <= circle->radius;
|
|
jbe@0
|
2692 PG_FREE_IF_COPY(cluster, 1);
|
|
jbe@0
|
2693 PG_RETURN_BOOL(retval);
|
|
jbe@0
|
2694 }
|
|
jbe@0
|
2695
|
|
jbe@17
|
2696 /* check if circle and cluster may overlap (l. ov. operator "&&+") in SQL */
|
|
jbe@10
|
2697 PG_FUNCTION_INFO_V1(pgl_ecircle_ecluster_may_overlap);
|
|
jbe@10
|
2698 Datum pgl_ecircle_ecluster_may_overlap(PG_FUNCTION_ARGS) {
|
|
jbe@10
|
2699 pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
|
|
jbe@10
|
2700 pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
|
|
jbe@46
|
2701 bool retval = (
|
|
jbe@46
|
2702 (1.0-PGL_SPHEROID_F) * /* margin for flattening and approximation */
|
|
jbe@46
|
2703 pgl_distance(
|
|
jbe@46
|
2704 circle->center.lat, circle->center.lon,
|
|
jbe@46
|
2705 cluster->bounding.center.lat, cluster->bounding.center.lon
|
|
jbe@46
|
2706 )
|
|
jbe@10
|
2707 ) <= circle->radius + cluster->bounding.radius;
|
|
jbe@10
|
2708 PG_FREE_IF_COPY(cluster, 1);
|
|
jbe@10
|
2709 PG_RETURN_BOOL(retval);
|
|
jbe@10
|
2710 }
|
|
jbe@10
|
2711
|
|
jbe@16
|
2712 /* check if two clusters overlap (overlap operator "&&") in SQL */
|
|
jbe@16
|
2713 PG_FUNCTION_INFO_V1(pgl_ecluster_overlap);
|
|
jbe@16
|
2714 Datum pgl_ecluster_overlap(PG_FUNCTION_ARGS) {
|
|
jbe@16
|
2715 pgl_cluster *cluster1 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
|
|
jbe@16
|
2716 pgl_cluster *cluster2 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
|
|
jbe@16
|
2717 bool retval;
|
|
jbe@46
|
2718 /* clusters with non-touching bounding circles (with margin for flattening)
|
|
jbe@46
|
2719 are always assumed to be non-overlapping */
|
|
jbe@16
|
2720 if (
|
|
jbe@46
|
2721 (1.0-PGL_SPHEROID_F) * /* margin for flattening and approximation */
|
|
jbe@16
|
2722 pgl_distance(
|
|
jbe@16
|
2723 cluster1->bounding.center.lat, cluster1->bounding.center.lon,
|
|
jbe@16
|
2724 cluster2->bounding.center.lat, cluster2->bounding.center.lon
|
|
jbe@16
|
2725 ) > cluster1->bounding.radius + cluster2->bounding.radius
|
|
jbe@16
|
2726 ) retval = false;
|
|
jbe@16
|
2727 else retval = pgl_clusters_overlap(cluster1, cluster2);
|
|
jbe@16
|
2728 PG_FREE_IF_COPY(cluster1, 0);
|
|
jbe@16
|
2729 PG_FREE_IF_COPY(cluster2, 1);
|
|
jbe@16
|
2730 PG_RETURN_BOOL(retval);
|
|
jbe@16
|
2731 }
|
|
jbe@16
|
2732
|
|
jbe@10
|
2733 /* check if two clusters may overlap (lossy overlap operator "&&+") in SQL */
|
|
jbe@10
|
2734 PG_FUNCTION_INFO_V1(pgl_ecluster_may_overlap);
|
|
jbe@10
|
2735 Datum pgl_ecluster_may_overlap(PG_FUNCTION_ARGS) {
|
|
jbe@10
|
2736 pgl_cluster *cluster1 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
|
|
jbe@10
|
2737 pgl_cluster *cluster2 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
|
|
jbe@46
|
2738 bool retval = (
|
|
jbe@46
|
2739 (1.0-PGL_SPHEROID_F) * /* margin for flattening and approximation */
|
|
jbe@46
|
2740 pgl_distance(
|
|
jbe@46
|
2741 cluster1->bounding.center.lat, cluster1->bounding.center.lon,
|
|
jbe@46
|
2742 cluster2->bounding.center.lat, cluster2->bounding.center.lon
|
|
jbe@46
|
2743 )
|
|
jbe@10
|
2744 ) <= cluster1->bounding.radius + cluster2->bounding.radius;
|
|
jbe@10
|
2745 PG_FREE_IF_COPY(cluster1, 0);
|
|
jbe@10
|
2746 PG_FREE_IF_COPY(cluster2, 1);
|
|
jbe@10
|
2747 PG_RETURN_BOOL(retval);
|
|
jbe@10
|
2748 }
|
|
jbe@10
|
2749
|
|
jbe@16
|
2750 /* check if second cluster is in first cluster (cont. operator "@>) in SQL */
|
|
jbe@16
|
2751 PG_FUNCTION_INFO_V1(pgl_ecluster_contains);
|
|
jbe@16
|
2752 Datum pgl_ecluster_contains(PG_FUNCTION_ARGS) {
|
|
jbe@16
|
2753 pgl_cluster *outer = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
|
|
jbe@16
|
2754 pgl_cluster *inner = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
|
|
jbe@16
|
2755 bool retval;
|
|
jbe@16
|
2756 /* clusters with non-touching bounding circles are always assumed to be
|
|
jbe@46
|
2757 non-overlapping */
|
|
jbe@16
|
2758 if (
|
|
jbe@46
|
2759 (1.0-PGL_SPHEROID_F) * /* margin for flattening and approximation */
|
|
jbe@16
|
2760 pgl_distance(
|
|
jbe@16
|
2761 outer->bounding.center.lat, outer->bounding.center.lon,
|
|
jbe@16
|
2762 inner->bounding.center.lat, inner->bounding.center.lon
|
|
jbe@16
|
2763 ) > outer->bounding.radius + inner->bounding.radius
|
|
jbe@16
|
2764 ) retval = false;
|
|
jbe@16
|
2765 else retval = pgl_cluster_in_cluster(outer, inner);
|
|
jbe@16
|
2766 PG_FREE_IF_COPY(outer, 0);
|
|
jbe@16
|
2767 PG_FREE_IF_COPY(inner, 1);
|
|
jbe@16
|
2768 PG_RETURN_BOOL(retval);
|
|
jbe@16
|
2769 }
|
|
jbe@16
|
2770
|
|
jbe@0
|
2771 /* calculate distance between two points ("<->" operator) in SQL */
|
|
jbe@0
|
2772 PG_FUNCTION_INFO_V1(pgl_epoint_distance);
|
|
jbe@0
|
2773 Datum pgl_epoint_distance(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2774 pgl_point *point1 = (pgl_point *)PG_GETARG_POINTER(0);
|
|
jbe@0
|
2775 pgl_point *point2 = (pgl_point *)PG_GETARG_POINTER(1);
|
|
jbe@0
|
2776 PG_RETURN_FLOAT8(pgl_distance(
|
|
jbe@0
|
2777 point1->lat, point1->lon, point2->lat, point2->lon
|
|
jbe@0
|
2778 ));
|
|
jbe@0
|
2779 }
|
|
jbe@0
|
2780
|
|
jbe@0
|
2781 /* calculate point to circle distance ("<->" operator) in SQL */
|
|
jbe@0
|
2782 PG_FUNCTION_INFO_V1(pgl_epoint_ecircle_distance);
|
|
jbe@0
|
2783 Datum pgl_epoint_ecircle_distance(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2784 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
|
|
jbe@0
|
2785 pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(1);
|
|
jbe@0
|
2786 double distance = pgl_distance(
|
|
jbe@0
|
2787 point->lat, point->lon, circle->center.lat, circle->center.lon
|
|
jbe@0
|
2788 ) - circle->radius;
|
|
jbe@0
|
2789 PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance);
|
|
jbe@0
|
2790 }
|
|
jbe@0
|
2791
|
|
jbe@0
|
2792 /* calculate point to cluster distance ("<->" operator) in SQL */
|
|
jbe@0
|
2793 PG_FUNCTION_INFO_V1(pgl_epoint_ecluster_distance);
|
|
jbe@0
|
2794 Datum pgl_epoint_ecluster_distance(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2795 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
|
|
jbe@0
|
2796 pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
|
|
jbe@0
|
2797 double distance = pgl_point_cluster_distance(point, cluster);
|
|
jbe@0
|
2798 PG_FREE_IF_COPY(cluster, 1);
|
|
jbe@0
|
2799 PG_RETURN_FLOAT8(distance);
|
|
jbe@0
|
2800 }
|
|
jbe@0
|
2801
|
|
jbe@0
|
2802 /* calculate distance between two circles ("<->" operator) in SQL */
|
|
jbe@0
|
2803 PG_FUNCTION_INFO_V1(pgl_ecircle_distance);
|
|
jbe@0
|
2804 Datum pgl_ecircle_distance(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2805 pgl_circle *circle1 = (pgl_circle *)PG_GETARG_POINTER(0);
|
|
jbe@0
|
2806 pgl_circle *circle2 = (pgl_circle *)PG_GETARG_POINTER(1);
|
|
jbe@0
|
2807 double distance = pgl_distance(
|
|
jbe@0
|
2808 circle1->center.lat, circle1->center.lon,
|
|
jbe@0
|
2809 circle2->center.lat, circle2->center.lon
|
|
jbe@0
|
2810 ) - (circle1->radius + circle2->radius);
|
|
jbe@0
|
2811 PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance);
|
|
jbe@0
|
2812 }
|
|
jbe@0
|
2813
|
|
jbe@0
|
2814 /* calculate circle to cluster distance ("<->" operator) in SQL */
|
|
jbe@0
|
2815 PG_FUNCTION_INFO_V1(pgl_ecircle_ecluster_distance);
|
|
jbe@0
|
2816 Datum pgl_ecircle_ecluster_distance(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2817 pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
|
|
jbe@0
|
2818 pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
|
|
jbe@0
|
2819 double distance = (
|
|
jbe@0
|
2820 pgl_point_cluster_distance(&(circle->center), cluster) - circle->radius
|
|
jbe@0
|
2821 );
|
|
jbe@0
|
2822 PG_FREE_IF_COPY(cluster, 1);
|
|
jbe@0
|
2823 PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance);
|
|
jbe@0
|
2824 }
|
|
jbe@0
|
2825
|
|
jbe@16
|
2826 /* calculate distance between two clusters ("<->" operator) in SQL */
|
|
jbe@16
|
2827 PG_FUNCTION_INFO_V1(pgl_ecluster_distance);
|
|
jbe@16
|
2828 Datum pgl_ecluster_distance(PG_FUNCTION_ARGS) {
|
|
jbe@16
|
2829 pgl_cluster *cluster1 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
|
|
jbe@16
|
2830 pgl_cluster *cluster2 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
|
|
jbe@16
|
2831 double retval = pgl_cluster_distance(cluster1, cluster2);
|
|
jbe@16
|
2832 PG_FREE_IF_COPY(cluster1, 0);
|
|
jbe@16
|
2833 PG_FREE_IF_COPY(cluster2, 1);
|
|
jbe@16
|
2834 PG_RETURN_FLOAT8(retval);
|
|
jbe@16
|
2835 }
|
|
jbe@16
|
2836
|
|
jbe@46
|
2837 /* calculate fair distance (see README) between cluster and point with
|
|
jbe@46
|
2838 precision denoted by sample count ("<=> operator) in SQL */
|
|
jbe@46
|
2839 PG_FUNCTION_INFO_V1(pgl_ecluster_epoint_sc_fair_distance);
|
|
jbe@46
|
2840 Datum pgl_ecluster_epoint_sc_fair_distance(PG_FUNCTION_ARGS) {
|
|
jbe@42
|
2841 pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
|
|
jbe@46
|
2842 pgl_point_sc *search = (pgl_point_sc *)PG_GETARG_POINTER(1);
|
|
jbe@46
|
2843 double retval = pgl_fair_distance(
|
|
jbe@46
|
2844 &search->point, cluster, search->samples
|
|
jbe@46
|
2845 );
|
|
jbe@42
|
2846 PG_FREE_IF_COPY(cluster, 0);
|
|
jbe@42
|
2847 PG_RETURN_FLOAT8(retval);
|
|
jbe@42
|
2848 }
|
|
jbe@42
|
2849
|
|
jbe@0
|
2850
|
|
jbe@0
|
2851 /*-----------------------------------------------------------*
|
|
jbe@0
|
2852 * B-tree comparison operators and index support functions *
|
|
jbe@0
|
2853 *-----------------------------------------------------------*/
|
|
jbe@0
|
2854
|
|
jbe@0
|
2855 /* macro for a B-tree operator (without detoasting) */
|
|
jbe@0
|
2856 #define PGL_BTREE_OPER(func, type, cmpfunc, oper) \
|
|
jbe@0
|
2857 PG_FUNCTION_INFO_V1(func); \
|
|
jbe@0
|
2858 Datum func(PG_FUNCTION_ARGS) { \
|
|
jbe@0
|
2859 type *a = (type *)PG_GETARG_POINTER(0); \
|
|
jbe@0
|
2860 type *b = (type *)PG_GETARG_POINTER(1); \
|
|
jbe@0
|
2861 PG_RETURN_BOOL(cmpfunc(a, b) oper 0); \
|
|
jbe@0
|
2862 }
|
|
jbe@0
|
2863
|
|
jbe@0
|
2864 /* macro for a B-tree comparison function (without detoasting) */
|
|
jbe@0
|
2865 #define PGL_BTREE_CMP(func, type, cmpfunc) \
|
|
jbe@0
|
2866 PG_FUNCTION_INFO_V1(func); \
|
|
jbe@0
|
2867 Datum func(PG_FUNCTION_ARGS) { \
|
|
jbe@0
|
2868 type *a = (type *)PG_GETARG_POINTER(0); \
|
|
jbe@0
|
2869 type *b = (type *)PG_GETARG_POINTER(1); \
|
|
jbe@0
|
2870 PG_RETURN_INT32(cmpfunc(a, b)); \
|
|
jbe@0
|
2871 }
|
|
jbe@0
|
2872
|
|
jbe@0
|
2873 /* macro for a B-tree operator (with detoasting) */
|
|
jbe@0
|
2874 #define PGL_BTREE_OPER_DETOAST(func, type, cmpfunc, oper) \
|
|
jbe@0
|
2875 PG_FUNCTION_INFO_V1(func); \
|
|
jbe@0
|
2876 Datum func(PG_FUNCTION_ARGS) { \
|
|
jbe@0
|
2877 bool res; \
|
|
jbe@0
|
2878 type *a = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); \
|
|
jbe@0
|
2879 type *b = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); \
|
|
jbe@0
|
2880 res = cmpfunc(a, b) oper 0; \
|
|
jbe@0
|
2881 PG_FREE_IF_COPY(a, 0); \
|
|
jbe@0
|
2882 PG_FREE_IF_COPY(b, 1); \
|
|
jbe@0
|
2883 PG_RETURN_BOOL(res); \
|
|
jbe@0
|
2884 }
|
|
jbe@0
|
2885
|
|
jbe@0
|
2886 /* macro for a B-tree comparison function (with detoasting) */
|
|
jbe@0
|
2887 #define PGL_BTREE_CMP_DETOAST(func, type, cmpfunc) \
|
|
jbe@0
|
2888 PG_FUNCTION_INFO_V1(func); \
|
|
jbe@0
|
2889 Datum func(PG_FUNCTION_ARGS) { \
|
|
jbe@0
|
2890 int32_t res; \
|
|
jbe@0
|
2891 type *a = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); \
|
|
jbe@0
|
2892 type *b = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); \
|
|
jbe@0
|
2893 res = cmpfunc(a, b); \
|
|
jbe@0
|
2894 PG_FREE_IF_COPY(a, 0); \
|
|
jbe@0
|
2895 PG_FREE_IF_COPY(b, 1); \
|
|
jbe@0
|
2896 PG_RETURN_INT32(res); \
|
|
jbe@0
|
2897 }
|
|
jbe@0
|
2898
|
|
jbe@0
|
2899 /* B-tree operators and comparison function for point */
|
|
jbe@0
|
2900 PGL_BTREE_OPER(pgl_btree_epoint_lt, pgl_point, pgl_point_cmp, <)
|
|
jbe@0
|
2901 PGL_BTREE_OPER(pgl_btree_epoint_le, pgl_point, pgl_point_cmp, <=)
|
|
jbe@0
|
2902 PGL_BTREE_OPER(pgl_btree_epoint_eq, pgl_point, pgl_point_cmp, ==)
|
|
jbe@0
|
2903 PGL_BTREE_OPER(pgl_btree_epoint_ne, pgl_point, pgl_point_cmp, !=)
|
|
jbe@0
|
2904 PGL_BTREE_OPER(pgl_btree_epoint_ge, pgl_point, pgl_point_cmp, >=)
|
|
jbe@0
|
2905 PGL_BTREE_OPER(pgl_btree_epoint_gt, pgl_point, pgl_point_cmp, >)
|
|
jbe@0
|
2906 PGL_BTREE_CMP(pgl_btree_epoint_cmp, pgl_point, pgl_point_cmp)
|
|
jbe@0
|
2907
|
|
jbe@0
|
2908 /* B-tree operators and comparison function for box */
|
|
jbe@0
|
2909 PGL_BTREE_OPER(pgl_btree_ebox_lt, pgl_box, pgl_box_cmp, <)
|
|
jbe@0
|
2910 PGL_BTREE_OPER(pgl_btree_ebox_le, pgl_box, pgl_box_cmp, <=)
|
|
jbe@0
|
2911 PGL_BTREE_OPER(pgl_btree_ebox_eq, pgl_box, pgl_box_cmp, ==)
|
|
jbe@0
|
2912 PGL_BTREE_OPER(pgl_btree_ebox_ne, pgl_box, pgl_box_cmp, !=)
|
|
jbe@0
|
2913 PGL_BTREE_OPER(pgl_btree_ebox_ge, pgl_box, pgl_box_cmp, >=)
|
|
jbe@0
|
2914 PGL_BTREE_OPER(pgl_btree_ebox_gt, pgl_box, pgl_box_cmp, >)
|
|
jbe@0
|
2915 PGL_BTREE_CMP(pgl_btree_ebox_cmp, pgl_box, pgl_box_cmp)
|
|
jbe@0
|
2916
|
|
jbe@0
|
2917 /* B-tree operators and comparison function for circle */
|
|
jbe@0
|
2918 PGL_BTREE_OPER(pgl_btree_ecircle_lt, pgl_circle, pgl_circle_cmp, <)
|
|
jbe@0
|
2919 PGL_BTREE_OPER(pgl_btree_ecircle_le, pgl_circle, pgl_circle_cmp, <=)
|
|
jbe@0
|
2920 PGL_BTREE_OPER(pgl_btree_ecircle_eq, pgl_circle, pgl_circle_cmp, ==)
|
|
jbe@0
|
2921 PGL_BTREE_OPER(pgl_btree_ecircle_ne, pgl_circle, pgl_circle_cmp, !=)
|
|
jbe@0
|
2922 PGL_BTREE_OPER(pgl_btree_ecircle_ge, pgl_circle, pgl_circle_cmp, >=)
|
|
jbe@0
|
2923 PGL_BTREE_OPER(pgl_btree_ecircle_gt, pgl_circle, pgl_circle_cmp, >)
|
|
jbe@0
|
2924 PGL_BTREE_CMP(pgl_btree_ecircle_cmp, pgl_circle, pgl_circle_cmp)
|
|
jbe@0
|
2925
|
|
jbe@0
|
2926
|
|
jbe@0
|
2927 /*--------------------------------*
|
|
jbe@0
|
2928 * GiST index support functions *
|
|
jbe@0
|
2929 *--------------------------------*/
|
|
jbe@0
|
2930
|
|
jbe@0
|
2931 /* GiST "consistent" support function */
|
|
jbe@0
|
2932 PG_FUNCTION_INFO_V1(pgl_gist_consistent);
|
|
jbe@0
|
2933 Datum pgl_gist_consistent(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
2934 GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
|
|
jbe@0
|
2935 pgl_keyptr key = (pgl_keyptr)DatumGetPointer(entry->key);
|
|
jbe@0
|
2936 StrategyNumber strategy = (StrategyNumber)PG_GETARG_UINT16(2);
|
|
jbe@0
|
2937 bool *recheck = (bool *)PG_GETARG_POINTER(4);
|
|
jbe@0
|
2938 /* demand recheck because index and query methods are lossy */
|
|
jbe@0
|
2939 *recheck = true;
|
|
jbe@10
|
2940 /* strategy number aliases for different operators using the same strategy */
|
|
jbe@10
|
2941 strategy %= 100;
|
|
jbe@0
|
2942 /* strategy number 11: equality of two points */
|
|
jbe@0
|
2943 if (strategy == 11) {
|
|
jbe@0
|
2944 /* query datum is another point */
|
|
jbe@0
|
2945 pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1);
|
|
jbe@0
|
2946 /* convert other point to key */
|
|
jbe@0
|
2947 pgl_pointkey querykey;
|
|
jbe@0
|
2948 pgl_point_to_key(query, querykey);
|
|
jbe@0
|
2949 /* return true if both keys overlap */
|
|
jbe@0
|
2950 PG_RETURN_BOOL(pgl_keys_overlap(key, querykey));
|
|
jbe@0
|
2951 }
|
|
jbe@0
|
2952 /* strategy number 13: equality of two circles */
|
|
jbe@0
|
2953 if (strategy == 13) {
|
|
jbe@0
|
2954 /* query datum is another circle */
|
|
jbe@0
|
2955 pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1);
|
|
jbe@0
|
2956 /* convert other circle to key */
|
|
jbe@0
|
2957 pgl_areakey querykey;
|
|
jbe@0
|
2958 pgl_circle_to_key(query, querykey);
|
|
jbe@0
|
2959 /* return true if both keys overlap */
|
|
jbe@0
|
2960 PG_RETURN_BOOL(pgl_keys_overlap(key, querykey));
|
|
jbe@0
|
2961 }
|
|
jbe@0
|
2962 /* for all remaining strategies, keys on empty objects produce no match */
|
|
jbe@0
|
2963 /* (check necessary because query radius may be infinite) */
|
|
jbe@0
|
2964 if (PGL_KEY_IS_EMPTY(key)) PG_RETURN_BOOL(false);
|
|
jbe@0
|
2965 /* strategy number 21: overlapping with point */
|
|
jbe@0
|
2966 if (strategy == 21) {
|
|
jbe@0
|
2967 /* query datum is a point */
|
|
jbe@0
|
2968 pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1);
|
|
jbe@0
|
2969 /* return true if estimated distance (allowed to be smaller than real
|
|
jbe@0
|
2970 distance) between index key and point is zero */
|
|
jbe@0
|
2971 PG_RETURN_BOOL(pgl_estimate_key_distance(key, query) == 0);
|
|
jbe@0
|
2972 }
|
|
jbe@0
|
2973 /* strategy number 22: (point) overlapping with box */
|
|
jbe@0
|
2974 if (strategy == 22) {
|
|
jbe@0
|
2975 /* query datum is a box */
|
|
jbe@0
|
2976 pgl_box *query = (pgl_box *)PG_GETARG_POINTER(1);
|
|
jbe@0
|
2977 /* determine bounding box of indexed key */
|
|
jbe@0
|
2978 pgl_box keybox;
|
|
jbe@0
|
2979 pgl_key_to_box(key, &keybox);
|
|
jbe@0
|
2980 /* return true if query box overlaps with bounding box of indexed key */
|
|
jbe@0
|
2981 PG_RETURN_BOOL(pgl_boxes_overlap(query, &keybox));
|
|
jbe@0
|
2982 }
|
|
jbe@0
|
2983 /* strategy number 23: overlapping with circle */
|
|
jbe@0
|
2984 if (strategy == 23) {
|
|
jbe@0
|
2985 /* query datum is a circle */
|
|
jbe@0
|
2986 pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1);
|
|
jbe@0
|
2987 /* return true if estimated distance (allowed to be smaller than real
|
|
jbe@0
|
2988 distance) between index key and circle center is smaller than radius */
|
|
jbe@0
|
2989 PG_RETURN_BOOL(
|
|
jbe@46
|
2990 (1.0-PGL_SPHEROID_F) * /* safety margin for lossy operator */
|
|
jbe@46
|
2991 pgl_estimate_key_distance(key, &(query->center))
|
|
jbe@46
|
2992 <= query->radius
|
|
jbe@0
|
2993 );
|
|
jbe@0
|
2994 }
|
|
jbe@0
|
2995 /* strategy number 24: overlapping with cluster */
|
|
jbe@0
|
2996 if (strategy == 24) {
|
|
jbe@0
|
2997 bool retval; /* return value */
|
|
jbe@0
|
2998 /* query datum is a cluster */
|
|
jbe@0
|
2999 pgl_cluster *query = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
|
|
jbe@0
|
3000 /* return true if estimated distance (allowed to be smaller than real
|
|
jbe@0
|
3001 distance) between index key and circle center is smaller than radius */
|
|
jbe@0
|
3002 retval = (
|
|
jbe@46
|
3003 (1.0-PGL_SPHEROID_F) * /* safety margin for lossy operator */
|
|
jbe@46
|
3004 pgl_estimate_key_distance(key, &(query->bounding.center))
|
|
jbe@46
|
3005 <= query->bounding.radius
|
|
jbe@0
|
3006 );
|
|
jbe@0
|
3007 PG_FREE_IF_COPY(query, 1); /* free detoasted cluster (if copy) */
|
|
jbe@0
|
3008 PG_RETURN_BOOL(retval);
|
|
jbe@0
|
3009 }
|
|
jbe@0
|
3010 /* throw error for any unknown strategy number */
|
|
jbe@0
|
3011 elog(ERROR, "unrecognized strategy number: %d", strategy);
|
|
jbe@0
|
3012 }
|
|
jbe@0
|
3013
|
|
jbe@0
|
3014 /* GiST "union" support function */
|
|
jbe@0
|
3015 PG_FUNCTION_INFO_V1(pgl_gist_union);
|
|
jbe@0
|
3016 Datum pgl_gist_union(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
3017 GistEntryVector *entryvec = (GistEntryVector *)PG_GETARG_POINTER(0);
|
|
jbe@0
|
3018 pgl_keyptr out; /* return value (to be palloc'ed) */
|
|
jbe@0
|
3019 int i;
|
|
jbe@0
|
3020 /* determine key size */
|
|
jbe@0
|
3021 size_t keysize = PGL_KEY_IS_AREAKEY(
|
|
jbe@0
|
3022 (pgl_keyptr)DatumGetPointer(entryvec->vector[0].key)
|
|
jbe@0
|
3023 ) ? sizeof (pgl_areakey) : sizeof(pgl_pointkey);
|
|
jbe@0
|
3024 /* begin with first key as result */
|
|
jbe@0
|
3025 out = palloc(keysize);
|
|
jbe@0
|
3026 memcpy(out, (pgl_keyptr)DatumGetPointer(entryvec->vector[0].key), keysize);
|
|
jbe@0
|
3027 /* unite current result with second, third, etc. key */
|
|
jbe@0
|
3028 for (i=1; i<entryvec->n; i++) {
|
|
jbe@0
|
3029 pgl_unite_keys(out, (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key));
|
|
jbe@0
|
3030 }
|
|
jbe@0
|
3031 /* return result */
|
|
jbe@0
|
3032 PG_RETURN_POINTER(out);
|
|
jbe@0
|
3033 }
|
|
jbe@0
|
3034
|
|
jbe@0
|
3035 /* GiST "compress" support function for indicis on points */
|
|
jbe@0
|
3036 PG_FUNCTION_INFO_V1(pgl_gist_compress_epoint);
|
|
jbe@0
|
3037 Datum pgl_gist_compress_epoint(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
3038 GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
|
|
jbe@0
|
3039 GISTENTRY *retval; /* return value (to be palloc'ed unless set to entry) */
|
|
jbe@0
|
3040 /* only transform new leaves */
|
|
jbe@0
|
3041 if (entry->leafkey) {
|
|
jbe@0
|
3042 /* get point to be transformed */
|
|
jbe@0
|
3043 pgl_point *point = (pgl_point *)DatumGetPointer(entry->key);
|
|
jbe@0
|
3044 /* allocate memory for key */
|
|
jbe@0
|
3045 pgl_keyptr key = palloc(sizeof(pgl_pointkey));
|
|
jbe@0
|
3046 /* transform point to key */
|
|
jbe@0
|
3047 pgl_point_to_key(point, key);
|
|
jbe@0
|
3048 /* create new GISTENTRY structure as return value */
|
|
jbe@0
|
3049 retval = palloc(sizeof(GISTENTRY));
|
|
jbe@0
|
3050 gistentryinit(
|
|
jbe@0
|
3051 *retval, PointerGetDatum(key),
|
|
jbe@60
|
3052 entry->rel, entry->page, entry->offset, false
|
|
jbe@0
|
3053 );
|
|
jbe@0
|
3054 } else {
|
|
jbe@0
|
3055 /* inner nodes have already been transformed */
|
|
jbe@0
|
3056 retval = entry;
|
|
jbe@0
|
3057 }
|
|
jbe@0
|
3058 /* return pointer to old or new GISTENTRY structure */
|
|
jbe@0
|
3059 PG_RETURN_POINTER(retval);
|
|
jbe@0
|
3060 }
|
|
jbe@0
|
3061
|
|
jbe@0
|
3062 /* GiST "compress" support function for indicis on circles */
|
|
jbe@0
|
3063 PG_FUNCTION_INFO_V1(pgl_gist_compress_ecircle);
|
|
jbe@0
|
3064 Datum pgl_gist_compress_ecircle(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
3065 GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
|
|
jbe@0
|
3066 GISTENTRY *retval; /* return value (to be palloc'ed unless set to entry) */
|
|
jbe@0
|
3067 /* only transform new leaves */
|
|
jbe@0
|
3068 if (entry->leafkey) {
|
|
jbe@0
|
3069 /* get circle to be transformed */
|
|
jbe@0
|
3070 pgl_circle *circle = (pgl_circle *)DatumGetPointer(entry->key);
|
|
jbe@0
|
3071 /* allocate memory for key */
|
|
jbe@0
|
3072 pgl_keyptr key = palloc(sizeof(pgl_areakey));
|
|
jbe@0
|
3073 /* transform circle to key */
|
|
jbe@0
|
3074 pgl_circle_to_key(circle, key);
|
|
jbe@0
|
3075 /* create new GISTENTRY structure as return value */
|
|
jbe@0
|
3076 retval = palloc(sizeof(GISTENTRY));
|
|
jbe@0
|
3077 gistentryinit(
|
|
jbe@0
|
3078 *retval, PointerGetDatum(key),
|
|
jbe@60
|
3079 entry->rel, entry->page, entry->offset, false
|
|
jbe@0
|
3080 );
|
|
jbe@0
|
3081 } else {
|
|
jbe@0
|
3082 /* inner nodes have already been transformed */
|
|
jbe@0
|
3083 retval = entry;
|
|
jbe@0
|
3084 }
|
|
jbe@0
|
3085 /* return pointer to old or new GISTENTRY structure */
|
|
jbe@0
|
3086 PG_RETURN_POINTER(retval);
|
|
jbe@0
|
3087 }
|
|
jbe@0
|
3088
|
|
jbe@0
|
3089 /* GiST "compress" support function for indices on clusters */
|
|
jbe@0
|
3090 PG_FUNCTION_INFO_V1(pgl_gist_compress_ecluster);
|
|
jbe@0
|
3091 Datum pgl_gist_compress_ecluster(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
3092 GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
|
|
jbe@0
|
3093 GISTENTRY *retval; /* return value (to be palloc'ed unless set to entry) */
|
|
jbe@0
|
3094 /* only transform new leaves */
|
|
jbe@0
|
3095 if (entry->leafkey) {
|
|
jbe@0
|
3096 /* get cluster to be transformed (detoasting necessary!) */
|
|
jbe@0
|
3097 pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(entry->key);
|
|
jbe@0
|
3098 /* allocate memory for key */
|
|
jbe@0
|
3099 pgl_keyptr key = palloc(sizeof(pgl_areakey));
|
|
jbe@0
|
3100 /* transform cluster to key */
|
|
jbe@0
|
3101 pgl_circle_to_key(&(cluster->bounding), key);
|
|
jbe@0
|
3102 /* create new GISTENTRY structure as return value */
|
|
jbe@0
|
3103 retval = palloc(sizeof(GISTENTRY));
|
|
jbe@0
|
3104 gistentryinit(
|
|
jbe@0
|
3105 *retval, PointerGetDatum(key),
|
|
jbe@60
|
3106 entry->rel, entry->page, entry->offset, false
|
|
jbe@0
|
3107 );
|
|
jbe@0
|
3108 /* free detoasted datum */
|
|
jbe@0
|
3109 if ((void *)cluster != (void *)DatumGetPointer(entry->key)) pfree(cluster);
|
|
jbe@0
|
3110 } else {
|
|
jbe@0
|
3111 /* inner nodes have already been transformed */
|
|
jbe@0
|
3112 retval = entry;
|
|
jbe@0
|
3113 }
|
|
jbe@0
|
3114 /* return pointer to old or new GISTENTRY structure */
|
|
jbe@0
|
3115 PG_RETURN_POINTER(retval);
|
|
jbe@0
|
3116 }
|
|
jbe@0
|
3117
|
|
jbe@0
|
3118 /* GiST "decompress" support function for indices */
|
|
jbe@0
|
3119 PG_FUNCTION_INFO_V1(pgl_gist_decompress);
|
|
jbe@0
|
3120 Datum pgl_gist_decompress(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
3121 /* return passed pointer without transformation */
|
|
jbe@0
|
3122 PG_RETURN_POINTER(PG_GETARG_POINTER(0));
|
|
jbe@0
|
3123 }
|
|
jbe@0
|
3124
|
|
jbe@0
|
3125 /* GiST "penalty" support function */
|
|
jbe@0
|
3126 PG_FUNCTION_INFO_V1(pgl_gist_penalty);
|
|
jbe@0
|
3127 Datum pgl_gist_penalty(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
3128 GISTENTRY *origentry = (GISTENTRY *)PG_GETARG_POINTER(0);
|
|
jbe@0
|
3129 GISTENTRY *newentry = (GISTENTRY *)PG_GETARG_POINTER(1);
|
|
jbe@0
|
3130 float *penalty = (float *)PG_GETARG_POINTER(2);
|
|
jbe@0
|
3131 /* get original key and key to insert */
|
|
jbe@0
|
3132 pgl_keyptr orig = (pgl_keyptr)DatumGetPointer(origentry->key);
|
|
jbe@0
|
3133 pgl_keyptr new = (pgl_keyptr)DatumGetPointer(newentry->key);
|
|
jbe@0
|
3134 /* copy original key */
|
|
jbe@0
|
3135 union { pgl_pointkey pointkey; pgl_areakey areakey; } union_key;
|
|
jbe@0
|
3136 if (PGL_KEY_IS_AREAKEY(orig)) {
|
|
jbe@0
|
3137 memcpy(union_key.areakey, orig, sizeof(union_key.areakey));
|
|
jbe@0
|
3138 } else {
|
|
jbe@0
|
3139 memcpy(union_key.pointkey, orig, sizeof(union_key.pointkey));
|
|
jbe@0
|
3140 }
|
|
jbe@0
|
3141 /* calculate union of both keys */
|
|
jbe@0
|
3142 pgl_unite_keys((pgl_keyptr)&union_key, new);
|
|
jbe@0
|
3143 /* penalty equal to reduction of key length (logarithm of added area) */
|
|
jbe@0
|
3144 /* (return value by setting referenced value and returning pointer) */
|
|
jbe@0
|
3145 *penalty = (
|
|
jbe@0
|
3146 PGL_KEY_NODEDEPTH(orig) - PGL_KEY_NODEDEPTH((pgl_keyptr)&union_key)
|
|
jbe@0
|
3147 );
|
|
jbe@0
|
3148 PG_RETURN_POINTER(penalty);
|
|
jbe@0
|
3149 }
|
|
jbe@0
|
3150
|
|
jbe@0
|
3151 /* GiST "picksplit" support function */
|
|
jbe@0
|
3152 PG_FUNCTION_INFO_V1(pgl_gist_picksplit);
|
|
jbe@0
|
3153 Datum pgl_gist_picksplit(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
3154 GistEntryVector *entryvec = (GistEntryVector *)PG_GETARG_POINTER(0);
|
|
jbe@0
|
3155 GIST_SPLITVEC *v = (GIST_SPLITVEC *)PG_GETARG_POINTER(1);
|
|
jbe@57
|
3156 OffsetNumber i; /* between FirstOffsetNumber and entryvec->n (exclusive) */
|
|
jbe@0
|
3157 union {
|
|
jbe@0
|
3158 pgl_pointkey pointkey;
|
|
jbe@0
|
3159 pgl_areakey areakey;
|
|
jbe@0
|
3160 } union_all; /* union of all keys (to be calculated from scratch)
|
|
jbe@0
|
3161 (later cut in half) */
|
|
jbe@0
|
3162 int is_areakey = PGL_KEY_IS_AREAKEY(
|
|
jbe@0
|
3163 (pgl_keyptr)DatumGetPointer(entryvec->vector[FirstOffsetNumber].key)
|
|
jbe@0
|
3164 );
|
|
jbe@0
|
3165 int keysize = is_areakey ? sizeof(pgl_areakey) : sizeof(pgl_pointkey);
|
|
jbe@0
|
3166 pgl_keyptr unionL = palloc(keysize); /* union of keys that go left */
|
|
jbe@0
|
3167 pgl_keyptr unionR = palloc(keysize); /* union of keys that go right */
|
|
jbe@0
|
3168 pgl_keyptr key; /* current key to be processed */
|
|
jbe@0
|
3169 /* allocate memory for array of left and right keys, set counts to zero */
|
|
jbe@0
|
3170 v->spl_left = (OffsetNumber *)palloc(entryvec->n * sizeof(OffsetNumber));
|
|
jbe@0
|
3171 v->spl_nleft = 0;
|
|
jbe@0
|
3172 v->spl_right = (OffsetNumber *)palloc(entryvec->n * sizeof(OffsetNumber));
|
|
jbe@0
|
3173 v->spl_nright = 0;
|
|
jbe@0
|
3174 /* calculate union of all keys from scratch */
|
|
jbe@0
|
3175 memcpy(
|
|
jbe@0
|
3176 (pgl_keyptr)&union_all,
|
|
jbe@0
|
3177 (pgl_keyptr)DatumGetPointer(entryvec->vector[FirstOffsetNumber].key),
|
|
jbe@0
|
3178 keysize
|
|
jbe@0
|
3179 );
|
|
jbe@0
|
3180 for (i=FirstOffsetNumber+1; i<entryvec->n; i=OffsetNumberNext(i)) {
|
|
jbe@0
|
3181 pgl_unite_keys(
|
|
jbe@0
|
3182 (pgl_keyptr)&union_all,
|
|
jbe@0
|
3183 (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key)
|
|
jbe@0
|
3184 );
|
|
jbe@0
|
3185 }
|
|
jbe@0
|
3186 /* check if trivial split is necessary due to exhausted key length */
|
|
jbe@0
|
3187 /* (Note: keys for empty objects must have node depth set to maximum) */
|
|
jbe@0
|
3188 if (PGL_KEY_NODEDEPTH((pgl_keyptr)&union_all) == (
|
|
jbe@0
|
3189 is_areakey ? PGL_AREAKEY_MAXDEPTH : PGL_POINTKEY_MAXDEPTH
|
|
jbe@0
|
3190 )) {
|
|
jbe@0
|
3191 /* half of all keys go left */
|
|
jbe@0
|
3192 for (
|
|
jbe@0
|
3193 i=FirstOffsetNumber;
|
|
jbe@0
|
3194 i<FirstOffsetNumber+(entryvec->n - FirstOffsetNumber)/2;
|
|
jbe@0
|
3195 i=OffsetNumberNext(i)
|
|
jbe@0
|
3196 ) {
|
|
jbe@0
|
3197 /* pointer to current key */
|
|
jbe@0
|
3198 key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key);
|
|
jbe@0
|
3199 /* update unionL */
|
|
jbe@0
|
3200 /* check if key is first key that goes left */
|
|
jbe@0
|
3201 if (!v->spl_nleft) {
|
|
jbe@0
|
3202 /* first key that goes left is just copied to unionL */
|
|
jbe@0
|
3203 memcpy(unionL, key, keysize);
|
|
jbe@0
|
3204 } else {
|
|
jbe@0
|
3205 /* unite current value and next key */
|
|
jbe@0
|
3206 pgl_unite_keys(unionL, key);
|
|
jbe@0
|
3207 }
|
|
jbe@0
|
3208 /* append offset number to list of keys that go left */
|
|
jbe@0
|
3209 v->spl_left[v->spl_nleft++] = i;
|
|
jbe@0
|
3210 }
|
|
jbe@0
|
3211 /* other half goes right */
|
|
jbe@0
|
3212 for (
|
|
jbe@0
|
3213 i=FirstOffsetNumber+(entryvec->n - FirstOffsetNumber)/2;
|
|
jbe@0
|
3214 i<entryvec->n;
|
|
jbe@0
|
3215 i=OffsetNumberNext(i)
|
|
jbe@0
|
3216 ) {
|
|
jbe@0
|
3217 /* pointer to current key */
|
|
jbe@0
|
3218 key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key);
|
|
jbe@0
|
3219 /* update unionR */
|
|
jbe@0
|
3220 /* check if key is first key that goes right */
|
|
jbe@0
|
3221 if (!v->spl_nright) {
|
|
jbe@0
|
3222 /* first key that goes right is just copied to unionR */
|
|
jbe@0
|
3223 memcpy(unionR, key, keysize);
|
|
jbe@0
|
3224 } else {
|
|
jbe@0
|
3225 /* unite current value and next key */
|
|
jbe@0
|
3226 pgl_unite_keys(unionR, key);
|
|
jbe@0
|
3227 }
|
|
jbe@0
|
3228 /* append offset number to list of keys that go right */
|
|
jbe@0
|
3229 v->spl_right[v->spl_nright++] = i;
|
|
jbe@0
|
3230 }
|
|
jbe@0
|
3231 }
|
|
jbe@0
|
3232 /* otherwise, a non-trivial split is possible */
|
|
jbe@0
|
3233 else {
|
|
jbe@0
|
3234 /* cut covered area in half */
|
|
jbe@0
|
3235 /* (union_all then refers to area of keys that go left) */
|
|
jbe@0
|
3236 /* check if union of all keys covers empty and non-empty objects */
|
|
jbe@0
|
3237 if (PGL_KEY_IS_UNIVERSAL((pgl_keyptr)&union_all)) {
|
|
jbe@0
|
3238 /* if yes, split into empty and non-empty objects */
|
|
jbe@0
|
3239 pgl_key_set_empty((pgl_keyptr)&union_all);
|
|
jbe@0
|
3240 } else {
|
|
jbe@0
|
3241 /* otherwise split by next bit */
|
|
jbe@0
|
3242 ((pgl_keyptr)&union_all)[PGL_KEY_NODEDEPTH_OFFSET]++;
|
|
jbe@0
|
3243 /* NOTE: type bit conserved */
|
|
jbe@0
|
3244 }
|
|
jbe@0
|
3245 /* determine for each key if it goes left or right */
|
|
jbe@0
|
3246 for (i=FirstOffsetNumber; i<entryvec->n; i=OffsetNumberNext(i)) {
|
|
jbe@0
|
3247 /* pointer to current key */
|
|
jbe@0
|
3248 key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key);
|
|
jbe@0
|
3249 /* keys within one half of the area go left */
|
|
jbe@0
|
3250 if (pgl_keys_overlap((pgl_keyptr)&union_all, key)) {
|
|
jbe@0
|
3251 /* update unionL */
|
|
jbe@0
|
3252 /* check if key is first key that goes left */
|
|
jbe@0
|
3253 if (!v->spl_nleft) {
|
|
jbe@0
|
3254 /* first key that goes left is just copied to unionL */
|
|
jbe@0
|
3255 memcpy(unionL, key, keysize);
|
|
jbe@0
|
3256 } else {
|
|
jbe@0
|
3257 /* unite current value of unionL and processed key */
|
|
jbe@0
|
3258 pgl_unite_keys(unionL, key);
|
|
jbe@0
|
3259 }
|
|
jbe@0
|
3260 /* append offset number to list of keys that go left */
|
|
jbe@0
|
3261 v->spl_left[v->spl_nleft++] = i;
|
|
jbe@0
|
3262 }
|
|
jbe@0
|
3263 /* the other keys go right */
|
|
jbe@0
|
3264 else {
|
|
jbe@0
|
3265 /* update unionR */
|
|
jbe@0
|
3266 /* check if key is first key that goes right */
|
|
jbe@0
|
3267 if (!v->spl_nright) {
|
|
jbe@0
|
3268 /* first key that goes right is just copied to unionR */
|
|
jbe@0
|
3269 memcpy(unionR, key, keysize);
|
|
jbe@0
|
3270 } else {
|
|
jbe@0
|
3271 /* unite current value of unionR and processed key */
|
|
jbe@0
|
3272 pgl_unite_keys(unionR, key);
|
|
jbe@0
|
3273 }
|
|
jbe@0
|
3274 /* append offset number to list of keys that go right */
|
|
jbe@0
|
3275 v->spl_right[v->spl_nright++] = i;
|
|
jbe@0
|
3276 }
|
|
jbe@0
|
3277 }
|
|
jbe@0
|
3278 }
|
|
jbe@0
|
3279 /* store unions in return value */
|
|
jbe@0
|
3280 v->spl_ldatum = PointerGetDatum(unionL);
|
|
jbe@0
|
3281 v->spl_rdatum = PointerGetDatum(unionR);
|
|
jbe@0
|
3282 /* return all results */
|
|
jbe@0
|
3283 PG_RETURN_POINTER(v);
|
|
jbe@0
|
3284 }
|
|
jbe@0
|
3285
|
|
jbe@0
|
3286 /* GiST "same"/"equal" support function */
|
|
jbe@0
|
3287 PG_FUNCTION_INFO_V1(pgl_gist_same);
|
|
jbe@0
|
3288 Datum pgl_gist_same(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
3289 pgl_keyptr key1 = (pgl_keyptr)PG_GETARG_POINTER(0);
|
|
jbe@0
|
3290 pgl_keyptr key2 = (pgl_keyptr)PG_GETARG_POINTER(1);
|
|
jbe@0
|
3291 bool *boolptr = (bool *)PG_GETARG_POINTER(2);
|
|
jbe@0
|
3292 /* two keys are equal if they are binary equal */
|
|
jbe@0
|
3293 /* (return result by setting referenced boolean and returning pointer) */
|
|
jbe@0
|
3294 *boolptr = !memcmp(
|
|
jbe@0
|
3295 key1,
|
|
jbe@0
|
3296 key2,
|
|
jbe@0
|
3297 PGL_KEY_IS_AREAKEY(key1) ? sizeof(pgl_areakey) : sizeof(pgl_pointkey)
|
|
jbe@0
|
3298 );
|
|
jbe@0
|
3299 PG_RETURN_POINTER(boolptr);
|
|
jbe@0
|
3300 }
|
|
jbe@0
|
3301
|
|
jbe@0
|
3302 /* GiST "distance" support function */
|
|
jbe@0
|
3303 PG_FUNCTION_INFO_V1(pgl_gist_distance);
|
|
jbe@0
|
3304 Datum pgl_gist_distance(PG_FUNCTION_ARGS) {
|
|
jbe@0
|
3305 GISTENTRY *entry = (GISTENTRY *)PG_GETARG_POINTER(0);
|
|
jbe@0
|
3306 pgl_keyptr key = (pgl_keyptr)DatumGetPointer(entry->key);
|
|
jbe@0
|
3307 StrategyNumber strategy = (StrategyNumber)PG_GETARG_UINT16(2);
|
|
jbe@0
|
3308 bool *recheck = (bool *)PG_GETARG_POINTER(4);
|
|
jbe@0
|
3309 double distance; /* return value */
|
|
jbe@0
|
3310 /* demand recheck because distance is just an estimation */
|
|
jbe@0
|
3311 /* (real distance may be bigger) */
|
|
jbe@0
|
3312 *recheck = true;
|
|
jbe@10
|
3313 /* strategy number aliases for different operators using the same strategy */
|
|
jbe@10
|
3314 strategy %= 100;
|
|
jbe@0
|
3315 /* strategy number 31: distance to point */
|
|
jbe@0
|
3316 if (strategy == 31) {
|
|
jbe@0
|
3317 /* query datum is a point */
|
|
jbe@0
|
3318 pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1);
|
|
jbe@0
|
3319 /* use pgl_estimate_pointkey_distance() function to compute result */
|
|
jbe@0
|
3320 distance = pgl_estimate_key_distance(key, query);
|
|
jbe@0
|
3321 /* avoid infinity (reserved!) */
|
|
jbe@0
|
3322 if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE;
|
|
jbe@0
|
3323 /* return result */
|
|
jbe@0
|
3324 PG_RETURN_FLOAT8(distance);
|
|
jbe@0
|
3325 }
|
|
jbe@0
|
3326 /* strategy number 33: distance to circle */
|
|
jbe@0
|
3327 if (strategy == 33) {
|
|
jbe@0
|
3328 /* query datum is a circle */
|
|
jbe@0
|
3329 pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1);
|
|
jbe@0
|
3330 /* estimate distance to circle center and substract circle radius */
|
|
jbe@0
|
3331 distance = (
|
|
jbe@0
|
3332 pgl_estimate_key_distance(key, &(query->center)) - query->radius
|
|
jbe@0
|
3333 );
|
|
jbe@0
|
3334 /* convert non-positive values to zero and avoid infinity (reserved!) */
|
|
jbe@0
|
3335 if (distance <= 0) distance = 0;
|
|
jbe@0
|
3336 else if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE;
|
|
jbe@0
|
3337 /* return result */
|
|
jbe@0
|
3338 PG_RETURN_FLOAT8(distance);
|
|
jbe@0
|
3339 }
|
|
jbe@0
|
3340 /* strategy number 34: distance to cluster */
|
|
jbe@0
|
3341 if (strategy == 34) {
|
|
jbe@0
|
3342 /* query datum is a cluster */
|
|
jbe@0
|
3343 pgl_cluster *query = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
|
|
jbe@0
|
3344 /* estimate distance to bounding center and substract bounding radius */
|
|
jbe@0
|
3345 distance = (
|
|
jbe@0
|
3346 pgl_estimate_key_distance(key, &(query->bounding.center)) -
|
|
jbe@0
|
3347 query->bounding.radius
|
|
jbe@0
|
3348 );
|
|
jbe@0
|
3349 /* convert non-positive values to zero and avoid infinity (reserved!) */
|
|
jbe@0
|
3350 if (distance <= 0) distance = 0;
|
|
jbe@0
|
3351 else if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE;
|
|
jbe@0
|
3352 /* free detoasted cluster (if copy) */
|
|
jbe@0
|
3353 PG_FREE_IF_COPY(query, 1);
|
|
jbe@0
|
3354 /* return result */
|
|
jbe@0
|
3355 PG_RETURN_FLOAT8(distance);
|
|
jbe@0
|
3356 }
|
|
jbe@0
|
3357 /* throw error for any unknown strategy number */
|
|
jbe@0
|
3358 elog(ERROR, "unrecognized strategy number: %d", strategy);
|
|
jbe@0
|
3359 }
|
|
jbe@0
|
3360
|