Line data Source code
1 : /*-------------------------------------------------------------------------
2 : *
3 : * geo_ops.c
4 : * 2D geometric operations
5 : *
6 : * Portions Copyright (c) 1996-2017, PostgreSQL Global Development Group
7 : * Portions Copyright (c) 1994, Regents of the University of California
8 : *
9 : *
10 : * IDENTIFICATION
11 : * src/backend/utils/adt/geo_ops.c
12 : *
13 : *-------------------------------------------------------------------------
14 : */
15 : #include "postgres.h"
16 :
17 : #include <math.h>
18 : #include <limits.h>
19 : #include <float.h>
20 : #include <ctype.h>
21 :
22 : #include "libpq/pqformat.h"
23 : #include "miscadmin.h"
24 : #include "utils/builtins.h"
25 : #include "utils/geo_decls.h"
26 :
27 : #ifndef M_PI
28 : #define M_PI 3.14159265358979323846
29 : #endif
30 :
31 :
32 : /*
33 : * Internal routines
34 : */
35 :
36 : enum path_delim
37 : {
38 : PATH_NONE, PATH_OPEN, PATH_CLOSED
39 : };
40 :
41 : static int point_inside(Point *p, int npts, Point *plist);
42 : static int lseg_crossing(double x, double y, double px, double py);
43 : static BOX *box_construct(double x1, double x2, double y1, double y2);
44 : static BOX *box_copy(BOX *box);
45 : static BOX *box_fill(BOX *result, double x1, double x2, double y1, double y2);
46 : static bool box_ov(BOX *box1, BOX *box2);
47 : static double box_ht(BOX *box);
48 : static double box_wd(BOX *box);
49 : static double circle_ar(CIRCLE *circle);
50 : static CIRCLE *circle_copy(CIRCLE *circle);
51 : static LINE *line_construct_pm(Point *pt, double m);
52 : static void line_construct_pts(LINE *line, Point *pt1, Point *pt2);
53 : static bool lseg_intersect_internal(LSEG *l1, LSEG *l2);
54 : static double lseg_dt(LSEG *l1, LSEG *l2);
55 : static bool on_ps_internal(Point *pt, LSEG *lseg);
56 : static void make_bound_box(POLYGON *poly);
57 : static bool plist_same(int npts, Point *p1, Point *p2);
58 : static Point *point_construct(double x, double y);
59 : static Point *point_copy(Point *pt);
60 : static double single_decode(char *num, char **endptr_p,
61 : const char *type_name, const char *orig_string);
62 : static void single_encode(float8 x, StringInfo str);
63 : static void pair_decode(char *str, double *x, double *y, char **endptr_p,
64 : const char *type_name, const char *orig_string);
65 : static void pair_encode(float8 x, float8 y, StringInfo str);
66 : static int pair_count(char *s, char delim);
67 : static void path_decode(char *str, bool opentype, int npts, Point *p,
68 : bool *isopen, char **endptr_p,
69 : const char *type_name, const char *orig_string);
70 : static char *path_encode(enum path_delim path_delim, int npts, Point *pt);
71 : static void statlseg_construct(LSEG *lseg, Point *pt1, Point *pt2);
72 : static double box_ar(BOX *box);
73 : static void box_cn(Point *center, BOX *box);
74 : static Point *interpt_sl(LSEG *lseg, LINE *line);
75 : static bool has_interpt_sl(LSEG *lseg, LINE *line);
76 : static double dist_pl_internal(Point *pt, LINE *line);
77 : static double dist_ps_internal(Point *pt, LSEG *lseg);
78 : static Point *line_interpt_internal(LINE *l1, LINE *l2);
79 : static bool lseg_inside_poly(Point *a, Point *b, POLYGON *poly, int start);
80 : static Point *lseg_interpt_internal(LSEG *l1, LSEG *l2);
81 : static double dist_ppoly_internal(Point *pt, POLYGON *poly);
82 :
83 :
84 : /*
85 : * Delimiters for input and output strings.
86 : * LDELIM, RDELIM, and DELIM are left, right, and separator delimiters, respectively.
87 : * LDELIM_EP, RDELIM_EP are left and right delimiters for paths with endpoints.
88 : */
89 :
90 : #define LDELIM '('
91 : #define RDELIM ')'
92 : #define DELIM ','
93 : #define LDELIM_EP '['
94 : #define RDELIM_EP ']'
95 : #define LDELIM_C '<'
96 : #define RDELIM_C '>'
97 :
98 :
99 : /*
100 : * Geometric data types are composed of points.
101 : * This code tries to support a common format throughout the data types,
102 : * to allow for more predictable usage and data type conversion.
103 : * The fundamental unit is the point. Other units are line segments,
104 : * open paths, boxes, closed paths, and polygons (which should be considered
105 : * non-intersecting closed paths).
106 : *
107 : * Data representation is as follows:
108 : * point: (x,y)
109 : * line segment: [(x1,y1),(x2,y2)]
110 : * box: (x1,y1),(x2,y2)
111 : * open path: [(x1,y1),...,(xn,yn)]
112 : * closed path: ((x1,y1),...,(xn,yn))
113 : * polygon: ((x1,y1),...,(xn,yn))
114 : *
115 : * For boxes, the points are opposite corners with the first point at the top right.
116 : * For closed paths and polygons, the points should be reordered to allow
117 : * fast and correct equality comparisons.
118 : *
119 : * XXX perhaps points in complex shapes should be reordered internally
120 : * to allow faster internal operations, but should keep track of input order
121 : * and restore that order for text output - tgl 97/01/16
122 : */
123 :
124 : static double
125 56 : single_decode(char *num, char **endptr_p,
126 : const char *type_name, const char *orig_string)
127 : {
128 56 : return float8in_internal(num, endptr_p, type_name, orig_string);
129 : } /* single_decode() */
130 :
131 : static void
132 74 : single_encode(float8 x, StringInfo str)
133 : {
134 74 : char *xstr = float8out_internal(x);
135 :
136 74 : appendStringInfoString(str, xstr);
137 74 : pfree(xstr);
138 74 : } /* single_encode() */
139 :
140 : static void
141 18595 : pair_decode(char *str, double *x, double *y, char **endptr_p,
142 : const char *type_name, const char *orig_string)
143 : {
144 : bool has_delim;
145 :
146 37191 : while (isspace((unsigned char) *str))
147 1 : str++;
148 18595 : if ((has_delim = (*str == LDELIM)))
149 12214 : str++;
150 :
151 18595 : *x = float8in_internal(str, &str, type_name, orig_string);
152 :
153 18589 : if (*str++ != DELIM)
154 4 : ereport(ERROR,
155 : (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
156 : errmsg("invalid input syntax for type %s: \"%s\"",
157 : type_name, orig_string)));
158 :
159 18585 : *y = float8in_internal(str, &str, type_name, orig_string);
160 :
161 18584 : if (has_delim)
162 : {
163 12209 : if (*str++ != RDELIM)
164 1 : ereport(ERROR,
165 : (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
166 : errmsg("invalid input syntax for type %s: \"%s\"",
167 : type_name, orig_string)));
168 24416 : while (isspace((unsigned char) *str))
169 0 : str++;
170 : }
171 :
172 : /* report stopping point if wanted, else complain if not end of string */
173 18583 : if (endptr_p)
174 18396 : *endptr_p = str;
175 187 : else if (*str != '\0')
176 0 : ereport(ERROR,
177 : (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
178 : errmsg("invalid input syntax for type %s: \"%s\"",
179 : type_name, orig_string)));
180 18583 : }
181 :
182 : static void
183 3129 : pair_encode(float8 x, float8 y, StringInfo str)
184 : {
185 3129 : char *xstr = float8out_internal(x);
186 3129 : char *ystr = float8out_internal(y);
187 :
188 3129 : appendStringInfo(str, "%s,%s", xstr, ystr);
189 3129 : pfree(xstr);
190 3129 : pfree(ystr);
191 3129 : }
192 :
193 : static void
194 8560 : path_decode(char *str, bool opentype, int npts, Point *p,
195 : bool *isopen, char **endptr_p,
196 : const char *type_name, const char *orig_string)
197 : {
198 8560 : int depth = 0;
199 : char *cp;
200 : int i;
201 :
202 17120 : while (isspace((unsigned char) *str))
203 0 : str++;
204 8560 : if ((*isopen = (*str == LDELIM_EP)))
205 : {
206 : /* no open delimiter allowed? */
207 5114 : if (!opentype)
208 0 : ereport(ERROR,
209 : (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
210 : errmsg("invalid input syntax for type %s: \"%s\"",
211 : type_name, orig_string)));
212 5114 : depth++;
213 5114 : str++;
214 : }
215 3446 : else if (*str == LDELIM)
216 : {
217 3439 : cp = (str + 1);
218 6878 : while (isspace((unsigned char) *cp))
219 0 : cp++;
220 3439 : if (*cp == LDELIM)
221 : {
222 181 : depth++;
223 181 : str = cp;
224 : }
225 3258 : else if (strrchr(str, LDELIM) == str)
226 : {
227 3171 : depth++;
228 3171 : str = cp;
229 : }
230 : }
231 :
232 26906 : for (i = 0; i < npts; i++)
233 : {
234 18353 : pair_decode(str, &(p->x), &(p->y), &str, type_name, orig_string);
235 18346 : if (*str == DELIM)
236 9787 : str++;
237 18346 : p++;
238 : }
239 :
240 17106 : while (isspace((unsigned char) *str))
241 0 : str++;
242 25560 : while (depth > 0)
243 : {
244 16914 : if ((*str == RDELIM)
245 5112 : || ((*str == RDELIM_EP) && (*isopen) && (depth == 1)))
246 : {
247 8454 : depth--;
248 8454 : str++;
249 16908 : while (isspace((unsigned char) *str))
250 0 : str++;
251 : }
252 : else
253 6 : ereport(ERROR,
254 : (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
255 : errmsg("invalid input syntax for type %s: \"%s\"",
256 : type_name, orig_string)));
257 : }
258 :
259 : /* report stopping point if wanted, else complain if not end of string */
260 8547 : if (endptr_p)
261 5140 : *endptr_p = str;
262 3407 : else if (*str != '\0')
263 0 : ereport(ERROR,
264 : (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
265 : errmsg("invalid input syntax for type %s: \"%s\"",
266 : type_name, orig_string)));
267 8547 : } /* path_decode() */
268 :
269 : static char *
270 1479 : path_encode(enum path_delim path_delim, int npts, Point *pt)
271 : {
272 : StringInfoData str;
273 : int i;
274 :
275 1479 : initStringInfo(&str);
276 :
277 1479 : switch (path_delim)
278 : {
279 : case PATH_CLOSED:
280 166 : appendStringInfoChar(&str, LDELIM);
281 166 : break;
282 : case PATH_OPEN:
283 397 : appendStringInfoChar(&str, LDELIM_EP);
284 397 : break;
285 : case PATH_NONE:
286 916 : break;
287 : }
288 :
289 4534 : for (i = 0; i < npts; i++)
290 : {
291 3055 : if (i > 0)
292 1576 : appendStringInfoChar(&str, DELIM);
293 3055 : appendStringInfoChar(&str, LDELIM);
294 3055 : pair_encode(pt->x, pt->y, &str);
295 3055 : appendStringInfoChar(&str, RDELIM);
296 3055 : pt++;
297 : }
298 :
299 1479 : switch (path_delim)
300 : {
301 : case PATH_CLOSED:
302 166 : appendStringInfoChar(&str, RDELIM);
303 166 : break;
304 : case PATH_OPEN:
305 397 : appendStringInfoChar(&str, RDELIM_EP);
306 397 : break;
307 : case PATH_NONE:
308 916 : break;
309 : }
310 :
311 1479 : return str.data;
312 : } /* path_encode() */
313 :
314 : /*-------------------------------------------------------------
315 : * pair_count - count the number of points
316 : * allow the following notation:
317 : * '((1,2),(3,4))'
318 : * '(1,3,2,4)'
319 : * require an odd number of delim characters in the string
320 : *-------------------------------------------------------------*/
321 : static int
322 5217 : pair_count(char *s, char delim)
323 : {
324 5217 : int ndelim = 0;
325 :
326 28553 : while ((s = strchr(s, delim)) != NULL)
327 : {
328 18119 : ndelim++;
329 18119 : s++;
330 : }
331 5217 : return (ndelim % 2) ? ((ndelim + 1) / 2) : -1;
332 : }
333 :
334 :
335 : /***********************************************************************
336 : **
337 : ** Routines for two-dimensional boxes.
338 : **
339 : ***********************************************************************/
340 :
341 : /*----------------------------------------------------------
342 : * Formatting and conversion routines.
343 : *---------------------------------------------------------*/
344 :
345 : /* box_in - convert a string to internal form.
346 : *
347 : * External format: (two corners of box)
348 : * "(f8, f8), (f8, f8)"
349 : * also supports the older style "(f8, f8, f8, f8)"
350 : */
351 : Datum
352 3271 : box_in(PG_FUNCTION_ARGS)
353 : {
354 3271 : char *str = PG_GETARG_CSTRING(0);
355 3271 : BOX *box = (BOX *) palloc(sizeof(BOX));
356 : bool isopen;
357 : double x,
358 : y;
359 :
360 3271 : path_decode(str, false, 2, &(box->high), &isopen, NULL, "box", str);
361 :
362 : /* reorder corners if necessary... */
363 3269 : if (box->high.x < box->low.x)
364 : {
365 137 : x = box->high.x;
366 137 : box->high.x = box->low.x;
367 137 : box->low.x = x;
368 : }
369 3269 : if (box->high.y < box->low.y)
370 : {
371 140 : y = box->high.y;
372 140 : box->high.y = box->low.y;
373 140 : box->low.y = y;
374 : }
375 :
376 3269 : PG_RETURN_BOX_P(box);
377 : }
378 :
379 : /* box_out - convert a box to external form.
380 : */
381 : Datum
382 362 : box_out(PG_FUNCTION_ARGS)
383 : {
384 362 : BOX *box = PG_GETARG_BOX_P(0);
385 :
386 362 : PG_RETURN_CSTRING(path_encode(PATH_NONE, 2, &(box->high)));
387 : }
388 :
389 : /*
390 : * box_recv - converts external binary format to box
391 : */
392 : Datum
393 0 : box_recv(PG_FUNCTION_ARGS)
394 : {
395 0 : StringInfo buf = (StringInfo) PG_GETARG_POINTER(0);
396 : BOX *box;
397 : double x,
398 : y;
399 :
400 0 : box = (BOX *) palloc(sizeof(BOX));
401 :
402 0 : box->high.x = pq_getmsgfloat8(buf);
403 0 : box->high.y = pq_getmsgfloat8(buf);
404 0 : box->low.x = pq_getmsgfloat8(buf);
405 0 : box->low.y = pq_getmsgfloat8(buf);
406 :
407 : /* reorder corners if necessary... */
408 0 : if (box->high.x < box->low.x)
409 : {
410 0 : x = box->high.x;
411 0 : box->high.x = box->low.x;
412 0 : box->low.x = x;
413 : }
414 0 : if (box->high.y < box->low.y)
415 : {
416 0 : y = box->high.y;
417 0 : box->high.y = box->low.y;
418 0 : box->low.y = y;
419 : }
420 :
421 0 : PG_RETURN_BOX_P(box);
422 : }
423 :
424 : /*
425 : * box_send - converts box to binary format
426 : */
427 : Datum
428 0 : box_send(PG_FUNCTION_ARGS)
429 : {
430 0 : BOX *box = PG_GETARG_BOX_P(0);
431 : StringInfoData buf;
432 :
433 0 : pq_begintypsend(&buf);
434 0 : pq_sendfloat8(&buf, box->high.x);
435 0 : pq_sendfloat8(&buf, box->high.y);
436 0 : pq_sendfloat8(&buf, box->low.x);
437 0 : pq_sendfloat8(&buf, box->low.y);
438 0 : PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
439 : }
440 :
441 :
442 : /* box_construct - fill in a new box.
443 : */
444 : static BOX *
445 20270 : box_construct(double x1, double x2, double y1, double y2)
446 : {
447 20270 : BOX *result = (BOX *) palloc(sizeof(BOX));
448 :
449 20270 : return box_fill(result, x1, x2, y1, y2);
450 : }
451 :
452 :
453 : /* box_fill - fill in a given box struct
454 : */
455 : static BOX *
456 23464 : box_fill(BOX *result, double x1, double x2, double y1, double y2)
457 : {
458 23464 : if (x1 > x2)
459 : {
460 3120 : result->high.x = x1;
461 3120 : result->low.x = x2;
462 : }
463 : else
464 : {
465 20344 : result->high.x = x2;
466 20344 : result->low.x = x1;
467 : }
468 23464 : if (y1 > y2)
469 : {
470 3124 : result->high.y = y1;
471 3124 : result->low.y = y2;
472 : }
473 : else
474 : {
475 20340 : result->high.y = y2;
476 20340 : result->low.y = y1;
477 : }
478 :
479 23464 : return result;
480 : }
481 :
482 :
483 : /* box_copy - copy a box
484 : */
485 : static BOX *
486 0 : box_copy(BOX *box)
487 : {
488 0 : BOX *result = (BOX *) palloc(sizeof(BOX));
489 :
490 0 : memcpy((char *) result, (char *) box, sizeof(BOX));
491 :
492 0 : return result;
493 : }
494 :
495 :
496 : /*----------------------------------------------------------
497 : * Relational operators for BOXes.
498 : * <, >, <=, >=, and == are based on box area.
499 : *---------------------------------------------------------*/
500 :
501 : /* box_same - are two boxes identical?
502 : */
503 : Datum
504 1194 : box_same(PG_FUNCTION_ARGS)
505 : {
506 1194 : BOX *box1 = PG_GETARG_BOX_P(0);
507 1194 : BOX *box2 = PG_GETARG_BOX_P(1);
508 :
509 1194 : PG_RETURN_BOOL(FPeq(box1->high.x, box2->high.x) &&
510 : FPeq(box1->low.x, box2->low.x) &&
511 : FPeq(box1->high.y, box2->high.y) &&
512 : FPeq(box1->low.y, box2->low.y));
513 : }
514 :
515 : /* box_overlap - does box1 overlap box2?
516 : */
517 : Datum
518 6636 : box_overlap(PG_FUNCTION_ARGS)
519 : {
520 6636 : BOX *box1 = PG_GETARG_BOX_P(0);
521 6636 : BOX *box2 = PG_GETARG_BOX_P(1);
522 :
523 6636 : PG_RETURN_BOOL(box_ov(box1, box2));
524 : }
525 :
526 : static bool
527 239903 : box_ov(BOX *box1, BOX *box2)
528 : {
529 601293 : return (FPle(box1->low.x, box2->high.x) &&
530 138393 : FPle(box2->low.x, box1->high.x) &&
531 269626 : FPle(box1->low.y, box2->high.y) &&
532 12817 : FPle(box2->low.y, box1->high.y));
533 : }
534 :
535 : /* box_left - is box1 strictly left of box2?
536 : */
537 : Datum
538 1750 : box_left(PG_FUNCTION_ARGS)
539 : {
540 1750 : BOX *box1 = PG_GETARG_BOX_P(0);
541 1750 : BOX *box2 = PG_GETARG_BOX_P(1);
542 :
543 1750 : PG_RETURN_BOOL(FPlt(box1->high.x, box2->low.x));
544 : }
545 :
546 : /* box_overleft - is the right edge of box1 at or left of
547 : * the right edge of box2?
548 : *
549 : * This is "less than or equal" for the end of a time range,
550 : * when time ranges are stored as rectangles.
551 : */
552 : Datum
553 5546 : box_overleft(PG_FUNCTION_ARGS)
554 : {
555 5546 : BOX *box1 = PG_GETARG_BOX_P(0);
556 5546 : BOX *box2 = PG_GETARG_BOX_P(1);
557 :
558 5546 : PG_RETURN_BOOL(FPle(box1->high.x, box2->high.x));
559 : }
560 :
561 : /* box_right - is box1 strictly right of box2?
562 : */
563 : Datum
564 15849 : box_right(PG_FUNCTION_ARGS)
565 : {
566 15849 : BOX *box1 = PG_GETARG_BOX_P(0);
567 15849 : BOX *box2 = PG_GETARG_BOX_P(1);
568 :
569 15849 : PG_RETURN_BOOL(FPgt(box1->low.x, box2->high.x));
570 : }
571 :
572 : /* box_overright - is the left edge of box1 at or right of
573 : * the left edge of box2?
574 : *
575 : * This is "greater than or equal" for time ranges, when time ranges
576 : * are stored as rectangles.
577 : */
578 : Datum
579 10693 : box_overright(PG_FUNCTION_ARGS)
580 : {
581 10693 : BOX *box1 = PG_GETARG_BOX_P(0);
582 10693 : BOX *box2 = PG_GETARG_BOX_P(1);
583 :
584 10693 : PG_RETURN_BOOL(FPge(box1->low.x, box2->low.x));
585 : }
586 :
587 : /* box_below - is box1 strictly below box2?
588 : */
589 : Datum
590 2467 : box_below(PG_FUNCTION_ARGS)
591 : {
592 2467 : BOX *box1 = PG_GETARG_BOX_P(0);
593 2467 : BOX *box2 = PG_GETARG_BOX_P(1);
594 :
595 2467 : PG_RETURN_BOOL(FPlt(box1->high.y, box2->low.y));
596 : }
597 :
598 : /* box_overbelow - is the upper edge of box1 at or below
599 : * the upper edge of box2?
600 : */
601 : Datum
602 6456 : box_overbelow(PG_FUNCTION_ARGS)
603 : {
604 6456 : BOX *box1 = PG_GETARG_BOX_P(0);
605 6456 : BOX *box2 = PG_GETARG_BOX_P(1);
606 :
607 6456 : PG_RETURN_BOOL(FPle(box1->high.y, box2->high.y));
608 : }
609 :
610 : /* box_above - is box1 strictly above box2?
611 : */
612 : Datum
613 5390 : box_above(PG_FUNCTION_ARGS)
614 : {
615 5390 : BOX *box1 = PG_GETARG_BOX_P(0);
616 5390 : BOX *box2 = PG_GETARG_BOX_P(1);
617 :
618 5390 : PG_RETURN_BOOL(FPgt(box1->low.y, box2->high.y));
619 : }
620 :
621 : /* box_overabove - is the lower edge of box1 at or above
622 : * the lower edge of box2?
623 : */
624 : Datum
625 9506 : box_overabove(PG_FUNCTION_ARGS)
626 : {
627 9506 : BOX *box1 = PG_GETARG_BOX_P(0);
628 9506 : BOX *box2 = PG_GETARG_BOX_P(1);
629 :
630 9506 : PG_RETURN_BOOL(FPge(box1->low.y, box2->low.y));
631 : }
632 :
633 : /* box_contained - is box1 contained by box2?
634 : */
635 : Datum
636 5747 : box_contained(PG_FUNCTION_ARGS)
637 : {
638 5747 : BOX *box1 = PG_GETARG_BOX_P(0);
639 5747 : BOX *box2 = PG_GETARG_BOX_P(1);
640 :
641 5747 : PG_RETURN_BOOL(FPle(box1->high.x, box2->high.x) &&
642 : FPge(box1->low.x, box2->low.x) &&
643 : FPle(box1->high.y, box2->high.y) &&
644 : FPge(box1->low.y, box2->low.y));
645 : }
646 :
647 : /* box_contain - does box1 contain box2?
648 : */
649 : Datum
650 1645 : box_contain(PG_FUNCTION_ARGS)
651 : {
652 1645 : BOX *box1 = PG_GETARG_BOX_P(0);
653 1645 : BOX *box2 = PG_GETARG_BOX_P(1);
654 :
655 1645 : PG_RETURN_BOOL(FPge(box1->high.x, box2->high.x) &&
656 : FPle(box1->low.x, box2->low.x) &&
657 : FPge(box1->high.y, box2->high.y) &&
658 : FPle(box1->low.y, box2->low.y));
659 : }
660 :
661 :
662 : /* box_positionop -
663 : * is box1 entirely {above,below} box2?
664 : *
665 : * box_below_eq and box_above_eq are obsolete versions that (probably
666 : * erroneously) accept the equal-boundaries case. Since these are not
667 : * in sync with the box_left and box_right code, they are deprecated and
668 : * not supported in the PG 8.1 rtree operator class extension.
669 : */
670 : Datum
671 0 : box_below_eq(PG_FUNCTION_ARGS)
672 : {
673 0 : BOX *box1 = PG_GETARG_BOX_P(0);
674 0 : BOX *box2 = PG_GETARG_BOX_P(1);
675 :
676 0 : PG_RETURN_BOOL(FPle(box1->high.y, box2->low.y));
677 : }
678 :
679 : Datum
680 0 : box_above_eq(PG_FUNCTION_ARGS)
681 : {
682 0 : BOX *box1 = PG_GETARG_BOX_P(0);
683 0 : BOX *box2 = PG_GETARG_BOX_P(1);
684 :
685 0 : PG_RETURN_BOOL(FPge(box1->low.y, box2->high.y));
686 : }
687 :
688 :
689 : /* box_relop - is area(box1) relop area(box2), within
690 : * our accuracy constraint?
691 : */
692 : Datum
693 4 : box_lt(PG_FUNCTION_ARGS)
694 : {
695 4 : BOX *box1 = PG_GETARG_BOX_P(0);
696 4 : BOX *box2 = PG_GETARG_BOX_P(1);
697 :
698 4 : PG_RETURN_BOOL(FPlt(box_ar(box1), box_ar(box2)));
699 : }
700 :
701 : Datum
702 4 : box_gt(PG_FUNCTION_ARGS)
703 : {
704 4 : BOX *box1 = PG_GETARG_BOX_P(0);
705 4 : BOX *box2 = PG_GETARG_BOX_P(1);
706 :
707 4 : PG_RETURN_BOOL(FPgt(box_ar(box1), box_ar(box2)));
708 : }
709 :
710 : Datum
711 4 : box_eq(PG_FUNCTION_ARGS)
712 : {
713 4 : BOX *box1 = PG_GETARG_BOX_P(0);
714 4 : BOX *box2 = PG_GETARG_BOX_P(1);
715 :
716 4 : PG_RETURN_BOOL(FPeq(box_ar(box1), box_ar(box2)));
717 : }
718 :
719 : Datum
720 4 : box_le(PG_FUNCTION_ARGS)
721 : {
722 4 : BOX *box1 = PG_GETARG_BOX_P(0);
723 4 : BOX *box2 = PG_GETARG_BOX_P(1);
724 :
725 4 : PG_RETURN_BOOL(FPle(box_ar(box1), box_ar(box2)));
726 : }
727 :
728 : Datum
729 4 : box_ge(PG_FUNCTION_ARGS)
730 : {
731 4 : BOX *box1 = PG_GETARG_BOX_P(0);
732 4 : BOX *box2 = PG_GETARG_BOX_P(1);
733 :
734 4 : PG_RETURN_BOOL(FPge(box_ar(box1), box_ar(box2)));
735 : }
736 :
737 :
738 : /*----------------------------------------------------------
739 : * "Arithmetic" operators on boxes.
740 : *---------------------------------------------------------*/
741 :
742 : /* box_area - returns the area of the box.
743 : */
744 : Datum
745 4 : box_area(PG_FUNCTION_ARGS)
746 : {
747 4 : BOX *box = PG_GETARG_BOX_P(0);
748 :
749 4 : PG_RETURN_FLOAT8(box_ar(box));
750 : }
751 :
752 :
753 : /* box_width - returns the width of the box
754 : * (horizontal magnitude).
755 : */
756 : Datum
757 4 : box_width(PG_FUNCTION_ARGS)
758 : {
759 4 : BOX *box = PG_GETARG_BOX_P(0);
760 :
761 4 : PG_RETURN_FLOAT8(box->high.x - box->low.x);
762 : }
763 :
764 :
765 : /* box_height - returns the height of the box
766 : * (vertical magnitude).
767 : */
768 : Datum
769 4 : box_height(PG_FUNCTION_ARGS)
770 : {
771 4 : BOX *box = PG_GETARG_BOX_P(0);
772 :
773 4 : PG_RETURN_FLOAT8(box->high.y - box->low.y);
774 : }
775 :
776 :
777 : /* box_distance - returns the distance between the
778 : * center points of two boxes.
779 : */
780 : Datum
781 0 : box_distance(PG_FUNCTION_ARGS)
782 : {
783 0 : BOX *box1 = PG_GETARG_BOX_P(0);
784 0 : BOX *box2 = PG_GETARG_BOX_P(1);
785 : Point a,
786 : b;
787 :
788 0 : box_cn(&a, box1);
789 0 : box_cn(&b, box2);
790 :
791 0 : PG_RETURN_FLOAT8(HYPOT(a.x - b.x, a.y - b.y));
792 : }
793 :
794 :
795 : /* box_center - returns the center point of the box.
796 : */
797 : Datum
798 12 : box_center(PG_FUNCTION_ARGS)
799 : {
800 12 : BOX *box = PG_GETARG_BOX_P(0);
801 12 : Point *result = (Point *) palloc(sizeof(Point));
802 :
803 12 : box_cn(result, box);
804 :
805 12 : PG_RETURN_POINT_P(result);
806 : }
807 :
808 :
809 : /* box_ar - returns the area of the box.
810 : */
811 : static double
812 44 : box_ar(BOX *box)
813 : {
814 44 : return box_wd(box) * box_ht(box);
815 : }
816 :
817 :
818 : /* box_cn - stores the centerpoint of the box into *center.
819 : */
820 : static void
821 12 : box_cn(Point *center, BOX *box)
822 : {
823 12 : center->x = (box->high.x + box->low.x) / 2.0;
824 12 : center->y = (box->high.y + box->low.y) / 2.0;
825 12 : }
826 :
827 :
828 : /* box_wd - returns the width (length) of the box
829 : * (horizontal magnitude).
830 : */
831 : static double
832 44 : box_wd(BOX *box)
833 : {
834 44 : return box->high.x - box->low.x;
835 : }
836 :
837 :
838 : /* box_ht - returns the height of the box
839 : * (vertical magnitude).
840 : */
841 : static double
842 44 : box_ht(BOX *box)
843 : {
844 44 : return box->high.y - box->low.y;
845 : }
846 :
847 :
848 : /*----------------------------------------------------------
849 : * Funky operations.
850 : *---------------------------------------------------------*/
851 :
852 : /* box_intersect -
853 : * returns the overlapping portion of two boxes,
854 : * or NULL if they do not intersect.
855 : */
856 : Datum
857 0 : box_intersect(PG_FUNCTION_ARGS)
858 : {
859 0 : BOX *box1 = PG_GETARG_BOX_P(0);
860 0 : BOX *box2 = PG_GETARG_BOX_P(1);
861 : BOX *result;
862 :
863 0 : if (!box_ov(box1, box2))
864 0 : PG_RETURN_NULL();
865 :
866 0 : result = (BOX *) palloc(sizeof(BOX));
867 :
868 0 : result->high.x = Min(box1->high.x, box2->high.x);
869 0 : result->low.x = Max(box1->low.x, box2->low.x);
870 0 : result->high.y = Min(box1->high.y, box2->high.y);
871 0 : result->low.y = Max(box1->low.y, box2->low.y);
872 :
873 0 : PG_RETURN_BOX_P(result);
874 : }
875 :
876 :
877 : /* box_diagonal -
878 : * returns a line segment which happens to be the
879 : * positive-slope diagonal of "box".
880 : */
881 : Datum
882 0 : box_diagonal(PG_FUNCTION_ARGS)
883 : {
884 0 : BOX *box = PG_GETARG_BOX_P(0);
885 0 : LSEG *result = (LSEG *) palloc(sizeof(LSEG));
886 :
887 0 : statlseg_construct(result, &box->high, &box->low);
888 :
889 0 : PG_RETURN_LSEG_P(result);
890 : }
891 :
892 : /***********************************************************************
893 : **
894 : ** Routines for 2D lines.
895 : **
896 : ***********************************************************************/
897 :
898 : static bool
899 2 : line_decode(char *s, const char *str, LINE *line)
900 : {
901 : /* s was already advanced over leading '{' */
902 2 : line->A = single_decode(s, &s, "line", str);
903 2 : if (*s++ != DELIM)
904 0 : return false;
905 2 : line->B = single_decode(s, &s, "line", str);
906 2 : if (*s++ != DELIM)
907 0 : return false;
908 2 : line->C = single_decode(s, &s, "line", str);
909 2 : if (*s++ != '}')
910 0 : return false;
911 4 : while (isspace((unsigned char) *s))
912 0 : s++;
913 2 : if (*s != '\0')
914 0 : return false;
915 2 : return true;
916 : }
917 :
918 : Datum
919 56 : line_in(PG_FUNCTION_ARGS)
920 : {
921 56 : char *str = PG_GETARG_CSTRING(0);
922 56 : LINE *line = (LINE *) palloc(sizeof(LINE));
923 : LSEG lseg;
924 : bool isopen;
925 : char *s;
926 :
927 56 : s = str;
928 112 : while (isspace((unsigned char) *s))
929 0 : s++;
930 56 : if (*s == '{')
931 : {
932 2 : if (!line_decode(s + 1, str, line))
933 0 : ereport(ERROR,
934 : (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
935 : errmsg("invalid input syntax for type %s: \"%s\"",
936 : "line", str)));
937 2 : if (FPzero(line->A) && FPzero(line->B))
938 1 : ereport(ERROR,
939 : (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
940 : errmsg("invalid line specification: A and B cannot both be zero")));
941 : }
942 : else
943 : {
944 54 : path_decode(s, true, 2, &(lseg.p[0]), &isopen, NULL, "line", str);
945 50 : if (FPeq(lseg.p[0].x, lseg.p[1].x) && FPeq(lseg.p[0].y, lseg.p[1].y))
946 1 : ereport(ERROR,
947 : (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
948 : errmsg("invalid line specification: must be two distinct points")));
949 49 : line_construct_pts(line, &lseg.p[0], &lseg.p[1]);
950 : }
951 :
952 50 : PG_RETURN_LINE_P(line);
953 : }
954 :
955 :
956 : Datum
957 19 : line_out(PG_FUNCTION_ARGS)
958 : {
959 19 : LINE *line = PG_GETARG_LINE_P(0);
960 19 : char *astr = float8out_internal(line->A);
961 19 : char *bstr = float8out_internal(line->B);
962 19 : char *cstr = float8out_internal(line->C);
963 :
964 19 : PG_RETURN_CSTRING(psprintf("{%s,%s,%s}", astr, bstr, cstr));
965 : }
966 :
967 : /*
968 : * line_recv - converts external binary format to line
969 : */
970 : Datum
971 0 : line_recv(PG_FUNCTION_ARGS)
972 : {
973 0 : StringInfo buf = (StringInfo) PG_GETARG_POINTER(0);
974 : LINE *line;
975 :
976 0 : line = (LINE *) palloc(sizeof(LINE));
977 :
978 0 : line->A = pq_getmsgfloat8(buf);
979 0 : line->B = pq_getmsgfloat8(buf);
980 0 : line->C = pq_getmsgfloat8(buf);
981 :
982 0 : PG_RETURN_LINE_P(line);
983 : }
984 :
985 : /*
986 : * line_send - converts line to binary format
987 : */
988 : Datum
989 0 : line_send(PG_FUNCTION_ARGS)
990 : {
991 0 : LINE *line = PG_GETARG_LINE_P(0);
992 : StringInfoData buf;
993 :
994 0 : pq_begintypsend(&buf);
995 0 : pq_sendfloat8(&buf, line->A);
996 0 : pq_sendfloat8(&buf, line->B);
997 0 : pq_sendfloat8(&buf, line->C);
998 0 : PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
999 : }
1000 :
1001 :
1002 : /*----------------------------------------------------------
1003 : * Conversion routines from one line formula to internal.
1004 : * Internal form: Ax+By+C=0
1005 : *---------------------------------------------------------*/
1006 :
1007 : /* line_construct_pm()
1008 : * point-slope
1009 : */
1010 : static LINE *
1011 12564 : line_construct_pm(Point *pt, double m)
1012 : {
1013 12564 : LINE *result = (LINE *) palloc(sizeof(LINE));
1014 :
1015 12564 : if (m == DBL_MAX)
1016 : {
1017 : /* vertical - use "x = C" */
1018 6168 : result->A = -1;
1019 6168 : result->B = 0;
1020 6168 : result->C = pt->x;
1021 : }
1022 : else
1023 : {
1024 : /* use "mx - y + yinter = 0" */
1025 6396 : result->A = m;
1026 6396 : result->B = -1.0;
1027 6396 : result->C = pt->y - m * pt->x;
1028 : }
1029 :
1030 12564 : return result;
1031 : }
1032 :
1033 : /*
1034 : * Fill already-allocated LINE struct from two points on the line
1035 : */
1036 : static void
1037 85958 : line_construct_pts(LINE *line, Point *pt1, Point *pt2)
1038 : {
1039 85958 : if (FPeq(pt1->x, pt2->x))
1040 : { /* vertical */
1041 : /* use "x = C" */
1042 17119 : line->A = -1;
1043 17119 : line->B = 0;
1044 17119 : line->C = pt1->x;
1045 : #ifdef GEODEBUG
1046 : printf("line_construct_pts- line is vertical\n");
1047 : #endif
1048 : }
1049 68839 : else if (FPeq(pt1->y, pt2->y))
1050 : { /* horizontal */
1051 : /* use "y = C" */
1052 18025 : line->A = 0;
1053 18025 : line->B = -1;
1054 18025 : line->C = pt1->y;
1055 : #ifdef GEODEBUG
1056 : printf("line_construct_pts- line is horizontal\n");
1057 : #endif
1058 : }
1059 : else
1060 : {
1061 : /* use "mx - y + yinter = 0" */
1062 50814 : line->A = (pt2->y - pt1->y) / (pt2->x - pt1->x);
1063 50814 : line->B = -1.0;
1064 50814 : line->C = pt1->y - line->A * pt1->x;
1065 : /* on some platforms, the preceding expression tends to produce -0 */
1066 50814 : if (line->C == 0.0)
1067 76 : line->C = 0.0;
1068 : #ifdef GEODEBUG
1069 : printf("line_construct_pts- line is neither vertical nor horizontal (diffs x=%.*g, y=%.*g\n",
1070 : DBL_DIG, (pt2->x - pt1->x), DBL_DIG, (pt2->y - pt1->y));
1071 : #endif
1072 : }
1073 85958 : }
1074 :
1075 : /* line_construct_pp()
1076 : * two points
1077 : */
1078 : Datum
1079 1 : line_construct_pp(PG_FUNCTION_ARGS)
1080 : {
1081 1 : Point *pt1 = PG_GETARG_POINT_P(0);
1082 1 : Point *pt2 = PG_GETARG_POINT_P(1);
1083 1 : LINE *result = (LINE *) palloc(sizeof(LINE));
1084 :
1085 1 : line_construct_pts(result, pt1, pt2);
1086 1 : PG_RETURN_LINE_P(result);
1087 : }
1088 :
1089 :
1090 : /*----------------------------------------------------------
1091 : * Relative position routines.
1092 : *---------------------------------------------------------*/
1093 :
1094 : Datum
1095 2 : line_intersect(PG_FUNCTION_ARGS)
1096 : {
1097 2 : LINE *l1 = PG_GETARG_LINE_P(0);
1098 2 : LINE *l2 = PG_GETARG_LINE_P(1);
1099 :
1100 2 : PG_RETURN_BOOL(!DatumGetBool(DirectFunctionCall2(line_parallel,
1101 : LinePGetDatum(l1),
1102 : LinePGetDatum(l2))));
1103 : }
1104 :
1105 : Datum
1106 49234 : line_parallel(PG_FUNCTION_ARGS)
1107 : {
1108 49234 : LINE *l1 = PG_GETARG_LINE_P(0);
1109 49234 : LINE *l2 = PG_GETARG_LINE_P(1);
1110 :
1111 49234 : if (FPzero(l1->B))
1112 16033 : PG_RETURN_BOOL(FPzero(l2->B));
1113 :
1114 33201 : PG_RETURN_BOOL(FPeq(l2->A, l1->A * (l2->B / l1->B)));
1115 : }
1116 :
1117 : Datum
1118 2 : line_perp(PG_FUNCTION_ARGS)
1119 : {
1120 2 : LINE *l1 = PG_GETARG_LINE_P(0);
1121 2 : LINE *l2 = PG_GETARG_LINE_P(1);
1122 :
1123 2 : if (FPzero(l1->A))
1124 1 : PG_RETURN_BOOL(FPzero(l2->B));
1125 1 : else if (FPzero(l1->B))
1126 0 : PG_RETURN_BOOL(FPzero(l2->A));
1127 :
1128 1 : PG_RETURN_BOOL(FPeq(((l1->A * l2->B) / (l1->B * l2->A)), -1.0));
1129 : }
1130 :
1131 : Datum
1132 2 : line_vertical(PG_FUNCTION_ARGS)
1133 : {
1134 2 : LINE *line = PG_GETARG_LINE_P(0);
1135 :
1136 2 : PG_RETURN_BOOL(FPzero(line->B));
1137 : }
1138 :
1139 : Datum
1140 2 : line_horizontal(PG_FUNCTION_ARGS)
1141 : {
1142 2 : LINE *line = PG_GETARG_LINE_P(0);
1143 :
1144 2 : PG_RETURN_BOOL(FPzero(line->A));
1145 : }
1146 :
1147 : Datum
1148 2 : line_eq(PG_FUNCTION_ARGS)
1149 : {
1150 2 : LINE *l1 = PG_GETARG_LINE_P(0);
1151 2 : LINE *l2 = PG_GETARG_LINE_P(1);
1152 : double k;
1153 :
1154 2 : if (!FPzero(l2->A))
1155 1 : k = l1->A / l2->A;
1156 1 : else if (!FPzero(l2->B))
1157 1 : k = l1->B / l2->B;
1158 0 : else if (!FPzero(l2->C))
1159 0 : k = l1->C / l2->C;
1160 : else
1161 0 : k = 1.0;
1162 :
1163 2 : PG_RETURN_BOOL(FPeq(l1->A, k * l2->A) &&
1164 : FPeq(l1->B, k * l2->B) &&
1165 : FPeq(l1->C, k * l2->C));
1166 : }
1167 :
1168 :
1169 : /*----------------------------------------------------------
1170 : * Line arithmetic routines.
1171 : *---------------------------------------------------------*/
1172 :
1173 : /* line_distance()
1174 : * Distance between two lines.
1175 : */
1176 : Datum
1177 8 : line_distance(PG_FUNCTION_ARGS)
1178 : {
1179 8 : LINE *l1 = PG_GETARG_LINE_P(0);
1180 8 : LINE *l2 = PG_GETARG_LINE_P(1);
1181 : float8 result;
1182 : Point *tmp;
1183 :
1184 8 : if (!DatumGetBool(DirectFunctionCall2(line_parallel,
1185 : LinePGetDatum(l1),
1186 : LinePGetDatum(l2))))
1187 4 : PG_RETURN_FLOAT8(0.0);
1188 4 : if (FPzero(l1->B)) /* vertical? */
1189 0 : PG_RETURN_FLOAT8(fabs(l1->C - l2->C));
1190 4 : tmp = point_construct(0.0, l1->C);
1191 4 : result = dist_pl_internal(tmp, l2);
1192 4 : PG_RETURN_FLOAT8(result);
1193 : }
1194 :
1195 : /* line_interpt()
1196 : * Point where two lines l1, l2 intersect (if any)
1197 : */
1198 : Datum
1199 2 : line_interpt(PG_FUNCTION_ARGS)
1200 : {
1201 2 : LINE *l1 = PG_GETARG_LINE_P(0);
1202 2 : LINE *l2 = PG_GETARG_LINE_P(1);
1203 : Point *result;
1204 :
1205 2 : result = line_interpt_internal(l1, l2);
1206 :
1207 2 : if (result == NULL)
1208 1 : PG_RETURN_NULL();
1209 1 : PG_RETURN_POINT_P(result);
1210 : }
1211 :
1212 : /*
1213 : * Internal version of line_interpt
1214 : *
1215 : * returns a NULL pointer if no intersection point
1216 : */
1217 : static Point *
1218 49222 : line_interpt_internal(LINE *l1, LINE *l2)
1219 : {
1220 : Point *result;
1221 : double x,
1222 : y;
1223 :
1224 : /*
1225 : * NOTE: if the lines are identical then we will find they are parallel
1226 : * and report "no intersection". This is a little weird, but since
1227 : * there's no *unique* intersection, maybe it's appropriate behavior.
1228 : */
1229 49222 : if (DatumGetBool(DirectFunctionCall2(line_parallel,
1230 : LinePGetDatum(l1),
1231 : LinePGetDatum(l2))))
1232 1637 : return NULL;
1233 :
1234 47585 : if (FPzero(l1->B)) /* l1 vertical? */
1235 : {
1236 15642 : x = l1->C;
1237 15642 : y = (l2->A * x + l2->C);
1238 : }
1239 31943 : else if (FPzero(l2->B)) /* l2 vertical? */
1240 : {
1241 6860 : x = l2->C;
1242 6860 : y = (l1->A * x + l1->C);
1243 : }
1244 : else
1245 : {
1246 25083 : x = (l1->C - l2->C) / (l2->A - l1->A);
1247 25083 : y = (l1->A * x + l1->C);
1248 : }
1249 47585 : result = point_construct(x, y);
1250 :
1251 : #ifdef GEODEBUG
1252 : printf("line_interpt- lines are A=%.*g, B=%.*g, C=%.*g, A=%.*g, B=%.*g, C=%.*g\n",
1253 : DBL_DIG, l1->A, DBL_DIG, l1->B, DBL_DIG, l1->C, DBL_DIG, l2->A, DBL_DIG, l2->B, DBL_DIG, l2->C);
1254 : printf("line_interpt- lines intersect at (%.*g,%.*g)\n", DBL_DIG, x, DBL_DIG, y);
1255 : #endif
1256 :
1257 47585 : return result;
1258 : }
1259 :
1260 :
1261 : /***********************************************************************
1262 : **
1263 : ** Routines for 2D paths (sequences of line segments, also
1264 : ** called `polylines').
1265 : **
1266 : ** This is not a general package for geometric paths,
1267 : ** which of course include polygons; the emphasis here
1268 : ** is on (for example) usefulness in wire layout.
1269 : **
1270 : ***********************************************************************/
1271 :
1272 : /*----------------------------------------------------------
1273 : * String to path / path to string conversion.
1274 : * External format:
1275 : * "((xcoord, ycoord),... )"
1276 : * "[(xcoord, ycoord),... ]"
1277 : * "(xcoord, ycoord),... "
1278 : * "[xcoord, ycoord,... ]"
1279 : * Also support older format:
1280 : * "(closed, npts, xcoord, ycoord,... )"
1281 : *---------------------------------------------------------*/
1282 :
1283 : Datum
1284 0 : path_area(PG_FUNCTION_ARGS)
1285 : {
1286 0 : PATH *path = PG_GETARG_PATH_P(0);
1287 0 : double area = 0.0;
1288 : int i,
1289 : j;
1290 :
1291 0 : if (!path->closed)
1292 0 : PG_RETURN_NULL();
1293 :
1294 0 : for (i = 0; i < path->npts; i++)
1295 : {
1296 0 : j = (i + 1) % path->npts;
1297 0 : area += path->p[i].x * path->p[j].y;
1298 0 : area -= path->p[i].y * path->p[j].x;
1299 : }
1300 :
1301 0 : area *= 0.5;
1302 0 : PG_RETURN_FLOAT8(area < 0.0 ? -area : area);
1303 : }
1304 :
1305 :
1306 : Datum
1307 5142 : path_in(PG_FUNCTION_ARGS)
1308 : {
1309 5142 : char *str = PG_GETARG_CSTRING(0);
1310 : PATH *path;
1311 : bool isopen;
1312 : char *s;
1313 : int npts;
1314 : int size;
1315 : int base_size;
1316 5142 : int depth = 0;
1317 :
1318 5142 : if ((npts = pair_count(str, ',')) <= 0)
1319 0 : ereport(ERROR,
1320 : (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
1321 : errmsg("invalid input syntax for type %s: \"%s\"",
1322 : "path", str)));
1323 :
1324 5142 : s = str;
1325 10284 : while (isspace((unsigned char) *s))
1326 0 : s++;
1327 :
1328 : /* skip single leading paren */
1329 5142 : if ((*s == LDELIM) && (strrchr(s, LDELIM) == s))
1330 : {
1331 3 : s++;
1332 3 : depth++;
1333 : }
1334 :
1335 5142 : base_size = sizeof(path->p[0]) * npts;
1336 5142 : size = offsetof(PATH, p) + base_size;
1337 :
1338 : /* Check for integer overflow */
1339 5142 : if (base_size / npts != sizeof(path->p[0]) || size <= base_size)
1340 0 : ereport(ERROR,
1341 : (errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED),
1342 : errmsg("too many points requested")));
1343 :
1344 5142 : path = (PATH *) palloc(size);
1345 :
1346 5142 : SET_VARSIZE(path, size);
1347 5142 : path->npts = npts;
1348 :
1349 5142 : path_decode(s, true, npts, &(path->p[0]), &isopen, &s, "path", str);
1350 :
1351 5140 : if (depth >= 1)
1352 : {
1353 3 : if (*s++ != RDELIM)
1354 0 : ereport(ERROR,
1355 : (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
1356 : errmsg("invalid input syntax for type %s: \"%s\"",
1357 : "path", str)));
1358 6 : while (isspace((unsigned char) *s))
1359 0 : s++;
1360 : }
1361 5140 : if (*s != '\0')
1362 0 : ereport(ERROR,
1363 : (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
1364 : errmsg("invalid input syntax for type %s: \"%s\"",
1365 : "path", str)));
1366 :
1367 5140 : path->closed = (!isopen);
1368 : /* prevent instability in unused pad bytes */
1369 5140 : path->dummy = 0;
1370 :
1371 5140 : PG_RETURN_PATH_P(path);
1372 : }
1373 :
1374 :
1375 : Datum
1376 399 : path_out(PG_FUNCTION_ARGS)
1377 : {
1378 399 : PATH *path = PG_GETARG_PATH_P(0);
1379 :
1380 399 : PG_RETURN_CSTRING(path_encode(path->closed ? PATH_CLOSED : PATH_OPEN, path->npts, path->p));
1381 : }
1382 :
1383 : /*
1384 : * path_recv - converts external binary format to path
1385 : *
1386 : * External representation is closed flag (a boolean byte), int32 number
1387 : * of points, and the points.
1388 : */
1389 : Datum
1390 0 : path_recv(PG_FUNCTION_ARGS)
1391 : {
1392 0 : StringInfo buf = (StringInfo) PG_GETARG_POINTER(0);
1393 : PATH *path;
1394 : int closed;
1395 : int32 npts;
1396 : int32 i;
1397 : int size;
1398 :
1399 0 : closed = pq_getmsgbyte(buf);
1400 0 : npts = pq_getmsgint(buf, sizeof(int32));
1401 0 : if (npts <= 0 || npts >= (int32) ((INT_MAX - offsetof(PATH, p)) / sizeof(Point)))
1402 0 : ereport(ERROR,
1403 : (errcode(ERRCODE_INVALID_BINARY_REPRESENTATION),
1404 : errmsg("invalid number of points in external \"path\" value")));
1405 :
1406 0 : size = offsetof(PATH, p) + sizeof(path->p[0]) * npts;
1407 0 : path = (PATH *) palloc(size);
1408 :
1409 0 : SET_VARSIZE(path, size);
1410 0 : path->npts = npts;
1411 0 : path->closed = (closed ? 1 : 0);
1412 : /* prevent instability in unused pad bytes */
1413 0 : path->dummy = 0;
1414 :
1415 0 : for (i = 0; i < npts; i++)
1416 : {
1417 0 : path->p[i].x = pq_getmsgfloat8(buf);
1418 0 : path->p[i].y = pq_getmsgfloat8(buf);
1419 : }
1420 :
1421 0 : PG_RETURN_PATH_P(path);
1422 : }
1423 :
1424 : /*
1425 : * path_send - converts path to binary format
1426 : */
1427 : Datum
1428 0 : path_send(PG_FUNCTION_ARGS)
1429 : {
1430 0 : PATH *path = PG_GETARG_PATH_P(0);
1431 : StringInfoData buf;
1432 : int32 i;
1433 :
1434 0 : pq_begintypsend(&buf);
1435 0 : pq_sendbyte(&buf, path->closed ? 1 : 0);
1436 0 : pq_sendint(&buf, path->npts, sizeof(int32));
1437 0 : for (i = 0; i < path->npts; i++)
1438 : {
1439 0 : pq_sendfloat8(&buf, path->p[i].x);
1440 0 : pq_sendfloat8(&buf, path->p[i].y);
1441 : }
1442 0 : PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
1443 : }
1444 :
1445 :
1446 : /*----------------------------------------------------------
1447 : * Relational operators.
1448 : * These are based on the path cardinality,
1449 : * as stupid as that sounds.
1450 : *
1451 : * Better relops and access methods coming soon.
1452 : *---------------------------------------------------------*/
1453 :
1454 : Datum
1455 0 : path_n_lt(PG_FUNCTION_ARGS)
1456 : {
1457 0 : PATH *p1 = PG_GETARG_PATH_P(0);
1458 0 : PATH *p2 = PG_GETARG_PATH_P(1);
1459 :
1460 0 : PG_RETURN_BOOL(p1->npts < p2->npts);
1461 : }
1462 :
1463 : Datum
1464 0 : path_n_gt(PG_FUNCTION_ARGS)
1465 : {
1466 0 : PATH *p1 = PG_GETARG_PATH_P(0);
1467 0 : PATH *p2 = PG_GETARG_PATH_P(1);
1468 :
1469 0 : PG_RETURN_BOOL(p1->npts > p2->npts);
1470 : }
1471 :
1472 : Datum
1473 0 : path_n_eq(PG_FUNCTION_ARGS)
1474 : {
1475 0 : PATH *p1 = PG_GETARG_PATH_P(0);
1476 0 : PATH *p2 = PG_GETARG_PATH_P(1);
1477 :
1478 0 : PG_RETURN_BOOL(p1->npts == p2->npts);
1479 : }
1480 :
1481 : Datum
1482 0 : path_n_le(PG_FUNCTION_ARGS)
1483 : {
1484 0 : PATH *p1 = PG_GETARG_PATH_P(0);
1485 0 : PATH *p2 = PG_GETARG_PATH_P(1);
1486 :
1487 0 : PG_RETURN_BOOL(p1->npts <= p2->npts);
1488 : }
1489 :
1490 : Datum
1491 0 : path_n_ge(PG_FUNCTION_ARGS)
1492 : {
1493 0 : PATH *p1 = PG_GETARG_PATH_P(0);
1494 0 : PATH *p2 = PG_GETARG_PATH_P(1);
1495 :
1496 0 : PG_RETURN_BOOL(p1->npts >= p2->npts);
1497 : }
1498 :
1499 : /*----------------------------------------------------------
1500 : * Conversion operators.
1501 : *---------------------------------------------------------*/
1502 :
1503 : Datum
1504 16 : path_isclosed(PG_FUNCTION_ARGS)
1505 : {
1506 16 : PATH *path = PG_GETARG_PATH_P(0);
1507 :
1508 16 : PG_RETURN_BOOL(path->closed);
1509 : }
1510 :
1511 : Datum
1512 16 : path_isopen(PG_FUNCTION_ARGS)
1513 : {
1514 16 : PATH *path = PG_GETARG_PATH_P(0);
1515 :
1516 16 : PG_RETURN_BOOL(!path->closed);
1517 : }
1518 :
1519 : Datum
1520 904 : path_npoints(PG_FUNCTION_ARGS)
1521 : {
1522 904 : PATH *path = PG_GETARG_PATH_P(0);
1523 :
1524 904 : PG_RETURN_INT32(path->npts);
1525 : }
1526 :
1527 :
1528 : Datum
1529 12 : path_close(PG_FUNCTION_ARGS)
1530 : {
1531 12 : PATH *path = PG_GETARG_PATH_P_COPY(0);
1532 :
1533 12 : path->closed = TRUE;
1534 :
1535 12 : PG_RETURN_PATH_P(path);
1536 : }
1537 :
1538 : Datum
1539 8 : path_open(PG_FUNCTION_ARGS)
1540 : {
1541 8 : PATH *path = PG_GETARG_PATH_P_COPY(0);
1542 :
1543 8 : path->closed = FALSE;
1544 :
1545 8 : PG_RETURN_PATH_P(path);
1546 : }
1547 :
1548 :
1549 : /* path_inter -
1550 : * Does p1 intersect p2 at any point?
1551 : * Use bounding boxes for a quick (O(n)) check, then do a
1552 : * O(n^2) iterative edge check.
1553 : */
1554 : Datum
1555 230153 : path_inter(PG_FUNCTION_ARGS)
1556 : {
1557 230153 : PATH *p1 = PG_GETARG_PATH_P(0);
1558 230153 : PATH *p2 = PG_GETARG_PATH_P(1);
1559 : BOX b1,
1560 : b2;
1561 : int i,
1562 : j;
1563 : LSEG seg1,
1564 : seg2;
1565 :
1566 230153 : if (p1->npts <= 0 || p2->npts <= 0)
1567 0 : PG_RETURN_BOOL(false);
1568 :
1569 230153 : b1.high.x = b1.low.x = p1->p[0].x;
1570 230153 : b1.high.y = b1.low.y = p1->p[0].y;
1571 875350 : for (i = 1; i < p1->npts; i++)
1572 : {
1573 645197 : b1.high.x = Max(p1->p[i].x, b1.high.x);
1574 645197 : b1.high.y = Max(p1->p[i].y, b1.high.y);
1575 645197 : b1.low.x = Min(p1->p[i].x, b1.low.x);
1576 645197 : b1.low.y = Min(p1->p[i].y, b1.low.y);
1577 : }
1578 230153 : b2.high.x = b2.low.x = p2->p[0].x;
1579 230153 : b2.high.y = b2.low.y = p2->p[0].y;
1580 658554 : for (i = 1; i < p2->npts; i++)
1581 : {
1582 428401 : b2.high.x = Max(p2->p[i].x, b2.high.x);
1583 428401 : b2.high.y = Max(p2->p[i].y, b2.high.y);
1584 428401 : b2.low.x = Min(p2->p[i].x, b2.low.x);
1585 428401 : b2.low.y = Min(p2->p[i].y, b2.low.y);
1586 : }
1587 230153 : if (!box_ov(&b1, &b2))
1588 224570 : PG_RETURN_BOOL(false);
1589 :
1590 : /* pairwise check lseg intersections */
1591 27502 : for (i = 0; i < p1->npts; i++)
1592 : {
1593 : int iprev;
1594 :
1595 23148 : if (i > 0)
1596 17565 : iprev = i - 1;
1597 : else
1598 : {
1599 5583 : if (!p1->closed)
1600 1331 : continue;
1601 4252 : iprev = p1->npts - 1; /* include the closure segment */
1602 : }
1603 :
1604 72611 : for (j = 0; j < p2->npts; j++)
1605 : {
1606 : int jprev;
1607 :
1608 52023 : if (j > 0)
1609 30206 : jprev = j - 1;
1610 : else
1611 : {
1612 21817 : if (!p2->closed)
1613 20623 : continue;
1614 1194 : jprev = p2->npts - 1; /* include the closure segment */
1615 : }
1616 :
1617 31400 : statlseg_construct(&seg1, &p1->p[iprev], &p1->p[i]);
1618 31400 : statlseg_construct(&seg2, &p2->p[jprev], &p2->p[j]);
1619 31400 : if (lseg_intersect_internal(&seg1, &seg2))
1620 1229 : PG_RETURN_BOOL(true);
1621 : }
1622 : }
1623 :
1624 : /* if we dropped through, no two segs intersected */
1625 4354 : PG_RETURN_BOOL(false);
1626 : }
1627 :
1628 : /* path_distance()
1629 : * This essentially does a cartesian product of the lsegs in the
1630 : * two paths, and finds the min distance between any two lsegs
1631 : */
1632 : Datum
1633 0 : path_distance(PG_FUNCTION_ARGS)
1634 : {
1635 0 : PATH *p1 = PG_GETARG_PATH_P(0);
1636 0 : PATH *p2 = PG_GETARG_PATH_P(1);
1637 0 : float8 min = 0.0; /* initialize to keep compiler quiet */
1638 0 : bool have_min = false;
1639 : float8 tmp;
1640 : int i,
1641 : j;
1642 : LSEG seg1,
1643 : seg2;
1644 :
1645 0 : for (i = 0; i < p1->npts; i++)
1646 : {
1647 : int iprev;
1648 :
1649 0 : if (i > 0)
1650 0 : iprev = i - 1;
1651 : else
1652 : {
1653 0 : if (!p1->closed)
1654 0 : continue;
1655 0 : iprev = p1->npts - 1; /* include the closure segment */
1656 : }
1657 :
1658 0 : for (j = 0; j < p2->npts; j++)
1659 : {
1660 : int jprev;
1661 :
1662 0 : if (j > 0)
1663 0 : jprev = j - 1;
1664 : else
1665 : {
1666 0 : if (!p2->closed)
1667 0 : continue;
1668 0 : jprev = p2->npts - 1; /* include the closure segment */
1669 : }
1670 :
1671 0 : statlseg_construct(&seg1, &p1->p[iprev], &p1->p[i]);
1672 0 : statlseg_construct(&seg2, &p2->p[jprev], &p2->p[j]);
1673 :
1674 0 : tmp = DatumGetFloat8(DirectFunctionCall2(lseg_distance,
1675 : LsegPGetDatum(&seg1),
1676 : LsegPGetDatum(&seg2)));
1677 0 : if (!have_min || tmp < min)
1678 : {
1679 0 : min = tmp;
1680 0 : have_min = true;
1681 : }
1682 : }
1683 : }
1684 :
1685 0 : if (!have_min)
1686 0 : PG_RETURN_NULL();
1687 :
1688 0 : PG_RETURN_FLOAT8(min);
1689 : }
1690 :
1691 :
1692 : /*----------------------------------------------------------
1693 : * "Arithmetic" operations.
1694 : *---------------------------------------------------------*/
1695 :
1696 : Datum
1697 0 : path_length(PG_FUNCTION_ARGS)
1698 : {
1699 0 : PATH *path = PG_GETARG_PATH_P(0);
1700 0 : float8 result = 0.0;
1701 : int i;
1702 :
1703 0 : for (i = 0; i < path->npts; i++)
1704 : {
1705 : int iprev;
1706 :
1707 0 : if (i > 0)
1708 0 : iprev = i - 1;
1709 : else
1710 : {
1711 0 : if (!path->closed)
1712 0 : continue;
1713 0 : iprev = path->npts - 1; /* include the closure segment */
1714 : }
1715 :
1716 0 : result += point_dt(&path->p[iprev], &path->p[i]);
1717 : }
1718 :
1719 0 : PG_RETURN_FLOAT8(result);
1720 : }
1721 :
1722 : /***********************************************************************
1723 : **
1724 : ** Routines for 2D points.
1725 : **
1726 : ***********************************************************************/
1727 :
1728 : /*----------------------------------------------------------
1729 : * String to point, point to string conversion.
1730 : * External format:
1731 : * "(x,y)"
1732 : * "x,y"
1733 : *---------------------------------------------------------*/
1734 :
1735 : Datum
1736 190 : point_in(PG_FUNCTION_ARGS)
1737 : {
1738 190 : char *str = PG_GETARG_CSTRING(0);
1739 190 : Point *point = (Point *) palloc(sizeof(Point));
1740 :
1741 190 : pair_decode(str, &point->x, &point->y, NULL, "point", str);
1742 187 : PG_RETURN_POINT_P(point);
1743 : }
1744 :
1745 : Datum
1746 554 : point_out(PG_FUNCTION_ARGS)
1747 : {
1748 554 : Point *pt = PG_GETARG_POINT_P(0);
1749 :
1750 554 : PG_RETURN_CSTRING(path_encode(PATH_NONE, 1, pt));
1751 : }
1752 :
1753 : /*
1754 : * point_recv - converts external binary format to point
1755 : */
1756 : Datum
1757 3 : point_recv(PG_FUNCTION_ARGS)
1758 : {
1759 3 : StringInfo buf = (StringInfo) PG_GETARG_POINTER(0);
1760 : Point *point;
1761 :
1762 3 : point = (Point *) palloc(sizeof(Point));
1763 3 : point->x = pq_getmsgfloat8(buf);
1764 3 : point->y = pq_getmsgfloat8(buf);
1765 3 : PG_RETURN_POINT_P(point);
1766 : }
1767 :
1768 : /*
1769 : * point_send - converts point to binary format
1770 : */
1771 : Datum
1772 3 : point_send(PG_FUNCTION_ARGS)
1773 : {
1774 3 : Point *pt = PG_GETARG_POINT_P(0);
1775 : StringInfoData buf;
1776 :
1777 3 : pq_begintypsend(&buf);
1778 3 : pq_sendfloat8(&buf, pt->x);
1779 3 : pq_sendfloat8(&buf, pt->y);
1780 3 : PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
1781 : }
1782 :
1783 :
1784 : static Point *
1785 157964 : point_construct(double x, double y)
1786 : {
1787 157964 : Point *result = (Point *) palloc(sizeof(Point));
1788 :
1789 157964 : result->x = x;
1790 157964 : result->y = y;
1791 157964 : return result;
1792 : }
1793 :
1794 :
1795 : static Point *
1796 19 : point_copy(Point *pt)
1797 : {
1798 : Point *result;
1799 :
1800 19 : if (!PointerIsValid(pt))
1801 0 : return NULL;
1802 :
1803 19 : result = (Point *) palloc(sizeof(Point));
1804 :
1805 19 : result->x = pt->x;
1806 19 : result->y = pt->y;
1807 19 : return result;
1808 : }
1809 :
1810 :
1811 : /*----------------------------------------------------------
1812 : * Relational operators for Points.
1813 : * Since we do have a sense of coordinates being
1814 : * "equal" to a given accuracy (point_vert, point_horiz),
1815 : * the other ops must preserve that sense. This means
1816 : * that results may, strictly speaking, be a lie (unless
1817 : * EPSILON = 0.0).
1818 : *---------------------------------------------------------*/
1819 :
1820 : Datum
1821 88111 : point_left(PG_FUNCTION_ARGS)
1822 : {
1823 88111 : Point *pt1 = PG_GETARG_POINT_P(0);
1824 88111 : Point *pt2 = PG_GETARG_POINT_P(1);
1825 :
1826 88111 : PG_RETURN_BOOL(FPlt(pt1->x, pt2->x));
1827 : }
1828 :
1829 : Datum
1830 1273977 : point_right(PG_FUNCTION_ARGS)
1831 : {
1832 1273977 : Point *pt1 = PG_GETARG_POINT_P(0);
1833 1273977 : Point *pt2 = PG_GETARG_POINT_P(1);
1834 :
1835 1273977 : PG_RETURN_BOOL(FPgt(pt1->x, pt2->x));
1836 : }
1837 :
1838 : Datum
1839 1308987 : point_above(PG_FUNCTION_ARGS)
1840 : {
1841 1308987 : Point *pt1 = PG_GETARG_POINT_P(0);
1842 1308987 : Point *pt2 = PG_GETARG_POINT_P(1);
1843 :
1844 1308987 : PG_RETURN_BOOL(FPgt(pt1->y, pt2->y));
1845 : }
1846 :
1847 : Datum
1848 146923 : point_below(PG_FUNCTION_ARGS)
1849 : {
1850 146923 : Point *pt1 = PG_GETARG_POINT_P(0);
1851 146923 : Point *pt2 = PG_GETARG_POINT_P(1);
1852 :
1853 146923 : PG_RETURN_BOOL(FPlt(pt1->y, pt2->y));
1854 : }
1855 :
1856 : Datum
1857 52452 : point_vert(PG_FUNCTION_ARGS)
1858 : {
1859 52452 : Point *pt1 = PG_GETARG_POINT_P(0);
1860 52452 : Point *pt2 = PG_GETARG_POINT_P(1);
1861 :
1862 52452 : PG_RETURN_BOOL(FPeq(pt1->x, pt2->x));
1863 : }
1864 :
1865 : Datum
1866 55215 : point_horiz(PG_FUNCTION_ARGS)
1867 : {
1868 55215 : Point *pt1 = PG_GETARG_POINT_P(0);
1869 55215 : Point *pt2 = PG_GETARG_POINT_P(1);
1870 :
1871 55215 : PG_RETURN_BOOL(FPeq(pt1->y, pt2->y));
1872 : }
1873 :
1874 : Datum
1875 13210 : point_eq(PG_FUNCTION_ARGS)
1876 : {
1877 13210 : Point *pt1 = PG_GETARG_POINT_P(0);
1878 13210 : Point *pt2 = PG_GETARG_POINT_P(1);
1879 :
1880 13210 : PG_RETURN_BOOL(FPeq(pt1->x, pt2->x) && FPeq(pt1->y, pt2->y));
1881 : }
1882 :
1883 : Datum
1884 17 : point_ne(PG_FUNCTION_ARGS)
1885 : {
1886 17 : Point *pt1 = PG_GETARG_POINT_P(0);
1887 17 : Point *pt2 = PG_GETARG_POINT_P(1);
1888 :
1889 17 : PG_RETURN_BOOL(FPne(pt1->x, pt2->x) || FPne(pt1->y, pt2->y));
1890 : }
1891 :
1892 : /*----------------------------------------------------------
1893 : * "Arithmetic" operators on points.
1894 : *---------------------------------------------------------*/
1895 :
1896 : Datum
1897 1436 : point_distance(PG_FUNCTION_ARGS)
1898 : {
1899 1436 : Point *pt1 = PG_GETARG_POINT_P(0);
1900 1436 : Point *pt2 = PG_GETARG_POINT_P(1);
1901 :
1902 1436 : PG_RETURN_FLOAT8(HYPOT(pt1->x - pt2->x, pt1->y - pt2->y));
1903 : }
1904 :
1905 : double
1906 217435 : point_dt(Point *pt1, Point *pt2)
1907 : {
1908 : #ifdef GEODEBUG
1909 : printf("point_dt- segment (%f,%f),(%f,%f) length is %f\n",
1910 : pt1->x, pt1->y, pt2->x, pt2->y, HYPOT(pt1->x - pt2->x, pt1->y - pt2->y));
1911 : #endif
1912 217435 : return HYPOT(pt1->x - pt2->x, pt1->y - pt2->y);
1913 : }
1914 :
1915 : Datum
1916 0 : point_slope(PG_FUNCTION_ARGS)
1917 : {
1918 0 : Point *pt1 = PG_GETARG_POINT_P(0);
1919 0 : Point *pt2 = PG_GETARG_POINT_P(1);
1920 :
1921 0 : PG_RETURN_FLOAT8(point_sl(pt1, pt2));
1922 : }
1923 :
1924 :
1925 : double
1926 30 : point_sl(Point *pt1, Point *pt2)
1927 : {
1928 150 : return (FPeq(pt1->x, pt2->x)
1929 : ? (double) DBL_MAX
1930 120 : : (pt1->y - pt2->y) / (pt1->x - pt2->x));
1931 : }
1932 :
1933 :
1934 : /***********************************************************************
1935 : **
1936 : ** Routines for 2D line segments.
1937 : **
1938 : ***********************************************************************/
1939 :
1940 : /*----------------------------------------------------------
1941 : * String to lseg, lseg to string conversion.
1942 : * External forms: "[(x1, y1), (x2, y2)]"
1943 : * "(x1, y1), (x2, y2)"
1944 : * "x1, y1, x2, y2"
1945 : * closed form ok "((x1, y1), (x2, y2))"
1946 : * (old form) "(x1, y1, x2, y2)"
1947 : *---------------------------------------------------------*/
1948 :
1949 : Datum
1950 22 : lseg_in(PG_FUNCTION_ARGS)
1951 : {
1952 22 : char *str = PG_GETARG_CSTRING(0);
1953 22 : LSEG *lseg = (LSEG *) palloc(sizeof(LSEG));
1954 : bool isopen;
1955 :
1956 22 : path_decode(str, true, 2, &(lseg->p[0]), &isopen, NULL, "lseg", str);
1957 18 : PG_RETURN_LSEG_P(lseg);
1958 : }
1959 :
1960 :
1961 : Datum
1962 41 : lseg_out(PG_FUNCTION_ARGS)
1963 : {
1964 41 : LSEG *ls = PG_GETARG_LSEG_P(0);
1965 :
1966 41 : PG_RETURN_CSTRING(path_encode(PATH_OPEN, 2, (Point *) &(ls->p[0])));
1967 : }
1968 :
1969 : /*
1970 : * lseg_recv - converts external binary format to lseg
1971 : */
1972 : Datum
1973 0 : lseg_recv(PG_FUNCTION_ARGS)
1974 : {
1975 0 : StringInfo buf = (StringInfo) PG_GETARG_POINTER(0);
1976 : LSEG *lseg;
1977 :
1978 0 : lseg = (LSEG *) palloc(sizeof(LSEG));
1979 :
1980 0 : lseg->p[0].x = pq_getmsgfloat8(buf);
1981 0 : lseg->p[0].y = pq_getmsgfloat8(buf);
1982 0 : lseg->p[1].x = pq_getmsgfloat8(buf);
1983 0 : lseg->p[1].y = pq_getmsgfloat8(buf);
1984 :
1985 0 : PG_RETURN_LSEG_P(lseg);
1986 : }
1987 :
1988 : /*
1989 : * lseg_send - converts lseg to binary format
1990 : */
1991 : Datum
1992 0 : lseg_send(PG_FUNCTION_ARGS)
1993 : {
1994 0 : LSEG *ls = PG_GETARG_LSEG_P(0);
1995 : StringInfoData buf;
1996 :
1997 0 : pq_begintypsend(&buf);
1998 0 : pq_sendfloat8(&buf, ls->p[0].x);
1999 0 : pq_sendfloat8(&buf, ls->p[0].y);
2000 0 : pq_sendfloat8(&buf, ls->p[1].x);
2001 0 : pq_sendfloat8(&buf, ls->p[1].y);
2002 0 : PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
2003 : }
2004 :
2005 :
2006 : /* lseg_construct -
2007 : * form a LSEG from two Points.
2008 : */
2009 : Datum
2010 0 : lseg_construct(PG_FUNCTION_ARGS)
2011 : {
2012 0 : Point *pt1 = PG_GETARG_POINT_P(0);
2013 0 : Point *pt2 = PG_GETARG_POINT_P(1);
2014 0 : LSEG *result = (LSEG *) palloc(sizeof(LSEG));
2015 :
2016 0 : result->p[0].x = pt1->x;
2017 0 : result->p[0].y = pt1->y;
2018 0 : result->p[1].x = pt2->x;
2019 0 : result->p[1].y = pt2->y;
2020 :
2021 0 : PG_RETURN_LSEG_P(result);
2022 : }
2023 :
2024 : /* like lseg_construct, but assume space already allocated */
2025 : static void
2026 62805 : statlseg_construct(LSEG *lseg, Point *pt1, Point *pt2)
2027 : {
2028 62805 : lseg->p[0].x = pt1->x;
2029 62805 : lseg->p[0].y = pt1->y;
2030 62805 : lseg->p[1].x = pt2->x;
2031 62805 : lseg->p[1].y = pt2->y;
2032 62805 : }
2033 :
2034 : Datum
2035 0 : lseg_length(PG_FUNCTION_ARGS)
2036 : {
2037 0 : LSEG *lseg = PG_GETARG_LSEG_P(0);
2038 :
2039 0 : PG_RETURN_FLOAT8(point_dt(&lseg->p[0], &lseg->p[1]));
2040 : }
2041 :
2042 : /*----------------------------------------------------------
2043 : * Relative position routines.
2044 : *---------------------------------------------------------*/
2045 :
2046 : /*
2047 : ** find intersection of the two lines, and see if it falls on
2048 : ** both segments.
2049 : */
2050 : Datum
2051 4228 : lseg_intersect(PG_FUNCTION_ARGS)
2052 : {
2053 4228 : LSEG *l1 = PG_GETARG_LSEG_P(0);
2054 4228 : LSEG *l2 = PG_GETARG_LSEG_P(1);
2055 :
2056 4228 : PG_RETURN_BOOL(lseg_intersect_internal(l1, l2));
2057 : }
2058 :
2059 : static bool
2060 35719 : lseg_intersect_internal(LSEG *l1, LSEG *l2)
2061 : {
2062 : LINE ln;
2063 : Point *interpt;
2064 : bool retval;
2065 :
2066 35719 : line_construct_pts(&ln, &l2->p[0], &l2->p[1]);
2067 35719 : interpt = interpt_sl(l1, &ln);
2068 :
2069 35719 : if (interpt != NULL && on_ps_internal(interpt, l2))
2070 2132 : retval = true; /* interpt on l1 and l2 */
2071 : else
2072 33587 : retval = false;
2073 35719 : return retval;
2074 : }
2075 :
2076 : Datum
2077 0 : lseg_parallel(PG_FUNCTION_ARGS)
2078 : {
2079 0 : LSEG *l1 = PG_GETARG_LSEG_P(0);
2080 0 : LSEG *l2 = PG_GETARG_LSEG_P(1);
2081 :
2082 0 : PG_RETURN_BOOL(FPeq(point_sl(&l1->p[0], &l1->p[1]),
2083 : point_sl(&l2->p[0], &l2->p[1])));
2084 : }
2085 :
2086 : /* lseg_perp()
2087 : * Determine if two line segments are perpendicular.
2088 : *
2089 : * This code did not get the correct answer for
2090 : * '((0,0),(0,1))'::lseg ?-| '((0,0),(1,0))'::lseg
2091 : * So, modified it to check explicitly for slope of vertical line
2092 : * returned by point_sl() and the results seem better.
2093 : * - thomas 1998-01-31
2094 : */
2095 : Datum
2096 0 : lseg_perp(PG_FUNCTION_ARGS)
2097 : {
2098 0 : LSEG *l1 = PG_GETARG_LSEG_P(0);
2099 0 : LSEG *l2 = PG_GETARG_LSEG_P(1);
2100 : double m1,
2101 : m2;
2102 :
2103 0 : m1 = point_sl(&(l1->p[0]), &(l1->p[1]));
2104 0 : m2 = point_sl(&(l2->p[0]), &(l2->p[1]));
2105 :
2106 : #ifdef GEODEBUG
2107 : printf("lseg_perp- slopes are %g and %g\n", m1, m2);
2108 : #endif
2109 0 : if (FPzero(m1))
2110 0 : PG_RETURN_BOOL(FPeq(m2, DBL_MAX));
2111 0 : else if (FPzero(m2))
2112 0 : PG_RETURN_BOOL(FPeq(m1, DBL_MAX));
2113 :
2114 0 : PG_RETURN_BOOL(FPeq(m1 / m2, -1.0));
2115 : }
2116 :
2117 : Datum
2118 0 : lseg_vertical(PG_FUNCTION_ARGS)
2119 : {
2120 0 : LSEG *lseg = PG_GETARG_LSEG_P(0);
2121 :
2122 0 : PG_RETURN_BOOL(FPeq(lseg->p[0].x, lseg->p[1].x));
2123 : }
2124 :
2125 : Datum
2126 0 : lseg_horizontal(PG_FUNCTION_ARGS)
2127 : {
2128 0 : LSEG *lseg = PG_GETARG_LSEG_P(0);
2129 :
2130 0 : PG_RETURN_BOOL(FPeq(lseg->p[0].y, lseg->p[1].y));
2131 : }
2132 :
2133 :
2134 : Datum
2135 0 : lseg_eq(PG_FUNCTION_ARGS)
2136 : {
2137 0 : LSEG *l1 = PG_GETARG_LSEG_P(0);
2138 0 : LSEG *l2 = PG_GETARG_LSEG_P(1);
2139 :
2140 0 : PG_RETURN_BOOL(FPeq(l1->p[0].x, l2->p[0].x) &&
2141 : FPeq(l1->p[0].y, l2->p[0].y) &&
2142 : FPeq(l1->p[1].x, l2->p[1].x) &&
2143 : FPeq(l1->p[1].y, l2->p[1].y));
2144 : }
2145 :
2146 : Datum
2147 0 : lseg_ne(PG_FUNCTION_ARGS)
2148 : {
2149 0 : LSEG *l1 = PG_GETARG_LSEG_P(0);
2150 0 : LSEG *l2 = PG_GETARG_LSEG_P(1);
2151 :
2152 0 : PG_RETURN_BOOL(!FPeq(l1->p[0].x, l2->p[0].x) ||
2153 : !FPeq(l1->p[0].y, l2->p[0].y) ||
2154 : !FPeq(l1->p[1].x, l2->p[1].x) ||
2155 : !FPeq(l1->p[1].y, l2->p[1].y));
2156 : }
2157 :
2158 : Datum
2159 0 : lseg_lt(PG_FUNCTION_ARGS)
2160 : {
2161 0 : LSEG *l1 = PG_GETARG_LSEG_P(0);
2162 0 : LSEG *l2 = PG_GETARG_LSEG_P(1);
2163 :
2164 0 : PG_RETURN_BOOL(FPlt(point_dt(&l1->p[0], &l1->p[1]),
2165 : point_dt(&l2->p[0], &l2->p[1])));
2166 : }
2167 :
2168 : Datum
2169 5 : lseg_le(PG_FUNCTION_ARGS)
2170 : {
2171 5 : LSEG *l1 = PG_GETARG_LSEG_P(0);
2172 5 : LSEG *l2 = PG_GETARG_LSEG_P(1);
2173 :
2174 5 : PG_RETURN_BOOL(FPle(point_dt(&l1->p[0], &l1->p[1]),
2175 : point_dt(&l2->p[0], &l2->p[1])));
2176 : }
2177 :
2178 : Datum
2179 0 : lseg_gt(PG_FUNCTION_ARGS)
2180 : {
2181 0 : LSEG *l1 = PG_GETARG_LSEG_P(0);
2182 0 : LSEG *l2 = PG_GETARG_LSEG_P(1);
2183 :
2184 0 : PG_RETURN_BOOL(FPgt(point_dt(&l1->p[0], &l1->p[1]),
2185 : point_dt(&l2->p[0], &l2->p[1])));
2186 : }
2187 :
2188 : Datum
2189 0 : lseg_ge(PG_FUNCTION_ARGS)
2190 : {
2191 0 : LSEG *l1 = PG_GETARG_LSEG_P(0);
2192 0 : LSEG *l2 = PG_GETARG_LSEG_P(1);
2193 :
2194 0 : PG_RETURN_BOOL(FPge(point_dt(&l1->p[0], &l1->p[1]),
2195 : point_dt(&l2->p[0], &l2->p[1])));
2196 : }
2197 :
2198 :
2199 : /*----------------------------------------------------------
2200 : * Line arithmetic routines.
2201 : *---------------------------------------------------------*/
2202 :
2203 : /* lseg_distance -
2204 : * If two segments don't intersect, then the closest
2205 : * point will be from one of the endpoints to the other
2206 : * segment.
2207 : */
2208 : Datum
2209 5 : lseg_distance(PG_FUNCTION_ARGS)
2210 : {
2211 5 : LSEG *l1 = PG_GETARG_LSEG_P(0);
2212 5 : LSEG *l2 = PG_GETARG_LSEG_P(1);
2213 :
2214 5 : PG_RETURN_FLOAT8(lseg_dt(l1, l2));
2215 : }
2216 :
2217 : /* lseg_dt()
2218 : * Distance between two line segments.
2219 : * Must check both sets of endpoints to ensure minimum distance is found.
2220 : * - thomas 1998-02-01
2221 : */
2222 : static double
2223 5 : lseg_dt(LSEG *l1, LSEG *l2)
2224 : {
2225 : double result,
2226 : d;
2227 :
2228 5 : if (lseg_intersect_internal(l1, l2))
2229 0 : return 0.0;
2230 :
2231 5 : d = dist_ps_internal(&l1->p[0], l2);
2232 5 : result = d;
2233 5 : d = dist_ps_internal(&l1->p[1], l2);
2234 5 : result = Min(result, d);
2235 5 : d = dist_ps_internal(&l2->p[0], l1);
2236 5 : result = Min(result, d);
2237 5 : d = dist_ps_internal(&l2->p[1], l1);
2238 5 : result = Min(result, d);
2239 :
2240 5 : return result;
2241 : }
2242 :
2243 :
2244 : Datum
2245 0 : lseg_center(PG_FUNCTION_ARGS)
2246 : {
2247 0 : LSEG *lseg = PG_GETARG_LSEG_P(0);
2248 : Point *result;
2249 :
2250 0 : result = (Point *) palloc(sizeof(Point));
2251 :
2252 0 : result->x = (lseg->p[0].x + lseg->p[1].x) / 2.0;
2253 0 : result->y = (lseg->p[0].y + lseg->p[1].y) / 2.0;
2254 :
2255 0 : PG_RETURN_POINT_P(result);
2256 : }
2257 :
2258 : static Point *
2259 970 : lseg_interpt_internal(LSEG *l1, LSEG *l2)
2260 : {
2261 : Point *result;
2262 : LINE tmp1,
2263 : tmp2;
2264 :
2265 : /*
2266 : * Find the intersection of the appropriate lines, if any.
2267 : */
2268 970 : line_construct_pts(&tmp1, &l1->p[0], &l1->p[1]);
2269 970 : line_construct_pts(&tmp2, &l2->p[0], &l2->p[1]);
2270 970 : result = line_interpt_internal(&tmp1, &tmp2);
2271 970 : if (!PointerIsValid(result))
2272 25 : return NULL;
2273 :
2274 : /*
2275 : * If the line intersection point isn't within l1 (or equivalently l2),
2276 : * there is no valid segment intersection point at all.
2277 : */
2278 1847 : if (!on_ps_internal(result, l1) ||
2279 902 : !on_ps_internal(result, l2))
2280 : {
2281 45 : pfree(result);
2282 45 : return NULL;
2283 : }
2284 :
2285 : /*
2286 : * If there is an intersection, then check explicitly for matching
2287 : * endpoints since there may be rounding effects with annoying lsb
2288 : * residue. - tgl 1997-07-09
2289 : */
2290 1693 : if ((FPeq(l1->p[0].x, l2->p[0].x) && FPeq(l1->p[0].y, l2->p[0].y)) ||
2291 902 : (FPeq(l1->p[0].x, l2->p[1].x) && FPeq(l1->p[0].y, l2->p[1].y)))
2292 : {
2293 201 : result->x = l1->p[0].x;
2294 201 : result->y = l1->p[0].y;
2295 : }
2296 1054 : else if ((FPeq(l1->p[1].x, l2->p[0].x) && FPeq(l1->p[1].y, l2->p[0].y)) ||
2297 645 : (FPeq(l1->p[1].x, l2->p[1].x) && FPeq(l1->p[1].y, l2->p[1].y)))
2298 : {
2299 627 : result->x = l1->p[1].x;
2300 627 : result->y = l1->p[1].y;
2301 : }
2302 :
2303 900 : return result;
2304 : }
2305 :
2306 : /* lseg_interpt -
2307 : * Find the intersection point of two segments (if any).
2308 : */
2309 : Datum
2310 894 : lseg_interpt(PG_FUNCTION_ARGS)
2311 : {
2312 894 : LSEG *l1 = PG_GETARG_LSEG_P(0);
2313 894 : LSEG *l2 = PG_GETARG_LSEG_P(1);
2314 : Point *result;
2315 :
2316 894 : result = lseg_interpt_internal(l1, l2);
2317 894 : if (!PointerIsValid(result))
2318 0 : PG_RETURN_NULL();
2319 :
2320 894 : PG_RETURN_POINT_P(result);
2321 : }
2322 :
2323 : /***********************************************************************
2324 : **
2325 : ** Routines for position comparisons of differently-typed
2326 : ** 2D objects.
2327 : **
2328 : ***********************************************************************/
2329 :
2330 : /*---------------------------------------------------------------------
2331 : * dist_
2332 : * Minimum distance from one object to another.
2333 : *-------------------------------------------------------------------*/
2334 :
2335 : /*
2336 : * Distance from a point to a line
2337 : */
2338 : Datum
2339 8 : dist_pl(PG_FUNCTION_ARGS)
2340 : {
2341 8 : Point *pt = PG_GETARG_POINT_P(0);
2342 8 : LINE *line = PG_GETARG_LINE_P(1);
2343 :
2344 8 : PG_RETURN_FLOAT8(dist_pl_internal(pt, line));
2345 : }
2346 :
2347 : static double
2348 30 : dist_pl_internal(Point *pt, LINE *line)
2349 : {
2350 60 : return fabs((line->A * pt->x + line->B * pt->y + line->C) /
2351 30 : HYPOT(line->A, line->B));
2352 : }
2353 :
2354 : /*
2355 : * Distance from a point to a lseg
2356 : */
2357 : Datum
2358 0 : dist_ps(PG_FUNCTION_ARGS)
2359 : {
2360 0 : Point *pt = PG_GETARG_POINT_P(0);
2361 0 : LSEG *lseg = PG_GETARG_LSEG_P(1);
2362 :
2363 0 : PG_RETURN_FLOAT8(dist_ps_internal(pt, lseg));
2364 : }
2365 :
2366 : static double
2367 12502 : dist_ps_internal(Point *pt, LSEG *lseg)
2368 : {
2369 : double m; /* slope of perp. */
2370 : LINE *ln;
2371 : double result,
2372 : tmpdist;
2373 : Point *ip;
2374 :
2375 : /*
2376 : * Construct a line perpendicular to the input segment and through the
2377 : * input point
2378 : */
2379 12502 : if (lseg->p[1].x == lseg->p[0].x)
2380 6306 : m = 0;
2381 6196 : else if (lseg->p[1].y == lseg->p[0].y)
2382 6168 : m = (double) DBL_MAX; /* slope is infinite */
2383 : else
2384 28 : m = (lseg->p[0].x - lseg->p[1].x) / (lseg->p[1].y - lseg->p[0].y);
2385 12502 : ln = line_construct_pm(pt, m);
2386 :
2387 : #ifdef GEODEBUG
2388 : printf("dist_ps- line is A=%g B=%g C=%g from (point) slope (%f,%f) %g\n",
2389 : ln->A, ln->B, ln->C, pt->x, pt->y, m);
2390 : #endif
2391 :
2392 : /*
2393 : * Calculate distance to the line segment or to the nearest endpoint of
2394 : * the segment.
2395 : */
2396 :
2397 : /* intersection is on the line segment? */
2398 12502 : if ((ip = interpt_sl(lseg, ln)) != NULL)
2399 : {
2400 : /* yes, so use distance to the intersection point */
2401 25 : result = point_dt(pt, ip);
2402 : #ifdef GEODEBUG
2403 : printf("dist_ps- distance is %f to intersection point is (%f,%f)\n",
2404 : result, ip->x, ip->y);
2405 : #endif
2406 : }
2407 : else
2408 : {
2409 : /* no, so use distance to the nearer endpoint */
2410 12477 : result = point_dt(pt, &lseg->p[0]);
2411 12477 : tmpdist = point_dt(pt, &lseg->p[1]);
2412 12477 : if (tmpdist < result)
2413 3088 : result = tmpdist;
2414 : }
2415 :
2416 12502 : return result;
2417 : }
2418 :
2419 : /*
2420 : * Distance from a point to a path
2421 : */
2422 : Datum
2423 0 : dist_ppath(PG_FUNCTION_ARGS)
2424 : {
2425 0 : Point *pt = PG_GETARG_POINT_P(0);
2426 0 : PATH *path = PG_GETARG_PATH_P(1);
2427 0 : float8 result = 0.0; /* keep compiler quiet */
2428 0 : bool have_min = false;
2429 : float8 tmp;
2430 : int i;
2431 : LSEG lseg;
2432 :
2433 0 : switch (path->npts)
2434 : {
2435 : case 0:
2436 : /* no points in path? then result is undefined... */
2437 0 : PG_RETURN_NULL();
2438 : case 1:
2439 : /* one point in path? then get distance between two points... */
2440 0 : result = point_dt(pt, &path->p[0]);
2441 0 : break;
2442 : default:
2443 : /* make sure the path makes sense... */
2444 0 : Assert(path->npts > 1);
2445 :
2446 : /*
2447 : * the distance from a point to a path is the smallest distance
2448 : * from the point to any of its constituent segments.
2449 : */
2450 0 : for (i = 0; i < path->npts; i++)
2451 : {
2452 : int iprev;
2453 :
2454 0 : if (i > 0)
2455 0 : iprev = i - 1;
2456 : else
2457 : {
2458 0 : if (!path->closed)
2459 0 : continue;
2460 0 : iprev = path->npts - 1; /* include the closure segment */
2461 : }
2462 :
2463 0 : statlseg_construct(&lseg, &path->p[iprev], &path->p[i]);
2464 0 : tmp = dist_ps_internal(pt, &lseg);
2465 0 : if (!have_min || tmp < result)
2466 : {
2467 0 : result = tmp;
2468 0 : have_min = true;
2469 : }
2470 : }
2471 0 : break;
2472 : }
2473 0 : PG_RETURN_FLOAT8(result);
2474 : }
2475 :
2476 : /*
2477 : * Distance from a point to a box
2478 : */
2479 : Datum
2480 0 : dist_pb(PG_FUNCTION_ARGS)
2481 : {
2482 0 : Point *pt = PG_GETARG_POINT_P(0);
2483 0 : BOX *box = PG_GETARG_BOX_P(1);
2484 : float8 result;
2485 : Point *near;
2486 :
2487 0 : near = DatumGetPointP(DirectFunctionCall2(close_pb,
2488 : PointPGetDatum(pt),
2489 : BoxPGetDatum(box)));
2490 0 : result = point_dt(near, pt);
2491 :
2492 0 : PG_RETURN_FLOAT8(result);
2493 : }
2494 :
2495 : /*
2496 : * Distance from a lseg to a line
2497 : */
2498 : Datum
2499 8 : dist_sl(PG_FUNCTION_ARGS)
2500 : {
2501 8 : LSEG *lseg = PG_GETARG_LSEG_P(0);
2502 8 : LINE *line = PG_GETARG_LINE_P(1);
2503 : float8 result,
2504 : d2;
2505 :
2506 8 : if (has_interpt_sl(lseg, line))
2507 0 : result = 0.0;
2508 : else
2509 : {
2510 8 : result = dist_pl_internal(&lseg->p[0], line);
2511 8 : d2 = dist_pl_internal(&lseg->p[1], line);
2512 : /* XXX shouldn't we take the min not max? */
2513 8 : if (d2 > result)
2514 1 : result = d2;
2515 : }
2516 :
2517 8 : PG_RETURN_FLOAT8(result);
2518 : }
2519 :
2520 : /*
2521 : * Distance from a lseg to a box
2522 : */
2523 : Datum
2524 0 : dist_sb(PG_FUNCTION_ARGS)
2525 : {
2526 0 : LSEG *lseg = PG_GETARG_LSEG_P(0);
2527 0 : BOX *box = PG_GETARG_BOX_P(1);
2528 : Point *tmp;
2529 : Datum result;
2530 :
2531 0 : tmp = DatumGetPointP(DirectFunctionCall2(close_sb,
2532 : LsegPGetDatum(lseg),
2533 : BoxPGetDatum(box)));
2534 0 : result = DirectFunctionCall2(dist_pb,
2535 : PointPGetDatum(tmp),
2536 : BoxPGetDatum(box));
2537 :
2538 0 : PG_RETURN_DATUM(result);
2539 : }
2540 :
2541 : /*
2542 : * Distance from a line to a box
2543 : */
2544 : Datum
2545 0 : dist_lb(PG_FUNCTION_ARGS)
2546 : {
2547 : #ifdef NOT_USED
2548 : LINE *line = PG_GETARG_LINE_P(0);
2549 : BOX *box = PG_GETARG_BOX_P(1);
2550 : #endif
2551 :
2552 : /* need to think about this one for a while */
2553 0 : ereport(ERROR,
2554 : (errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
2555 : errmsg("function \"dist_lb\" not implemented")));
2556 :
2557 : PG_RETURN_NULL();
2558 : }
2559 :
2560 : /*
2561 : * Distance from a circle to a polygon
2562 : */
2563 : Datum
2564 0 : dist_cpoly(PG_FUNCTION_ARGS)
2565 : {
2566 0 : CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
2567 0 : POLYGON *poly = PG_GETARG_POLYGON_P(1);
2568 : float8 result;
2569 :
2570 : /* calculate distance to center, and subtract radius */
2571 0 : result = dist_ppoly_internal(&circle->center, poly);
2572 :
2573 0 : result -= circle->radius;
2574 0 : if (result < 0)
2575 0 : result = 0;
2576 :
2577 0 : PG_RETURN_FLOAT8(result);
2578 : }
2579 :
2580 : /*
2581 : * Distance from a point to a polygon
2582 : */
2583 : Datum
2584 5 : dist_ppoly(PG_FUNCTION_ARGS)
2585 : {
2586 5 : Point *point = PG_GETARG_POINT_P(0);
2587 5 : POLYGON *poly = PG_GETARG_POLYGON_P(1);
2588 : float8 result;
2589 :
2590 5 : result = dist_ppoly_internal(point, poly);
2591 :
2592 5 : PG_RETURN_FLOAT8(result);
2593 : }
2594 :
2595 : Datum
2596 3122 : dist_polyp(PG_FUNCTION_ARGS)
2597 : {
2598 3122 : POLYGON *poly = PG_GETARG_POLYGON_P(0);
2599 3122 : Point *point = PG_GETARG_POINT_P(1);
2600 : float8 result;
2601 :
2602 3122 : result = dist_ppoly_internal(point, poly);
2603 :
2604 3122 : PG_RETURN_FLOAT8(result);
2605 : }
2606 :
2607 : static double
2608 3127 : dist_ppoly_internal(Point *pt, POLYGON *poly)
2609 : {
2610 : float8 result;
2611 : float8 d;
2612 : int i;
2613 : LSEG seg;
2614 :
2615 3127 : if (point_inside(pt, poly->npts, poly->p) != 0)
2616 : {
2617 : #ifdef GEODEBUG
2618 : printf("dist_ppoly_internal- point inside of polygon\n");
2619 : #endif
2620 3 : return 0.0;
2621 : }
2622 :
2623 : /* initialize distance with segment between first and last points */
2624 3124 : seg.p[0].x = poly->p[0].x;
2625 3124 : seg.p[0].y = poly->p[0].y;
2626 3124 : seg.p[1].x = poly->p[poly->npts - 1].x;
2627 3124 : seg.p[1].y = poly->p[poly->npts - 1].y;
2628 3124 : result = dist_ps_internal(pt, &seg);
2629 : #ifdef GEODEBUG
2630 : printf("dist_ppoly_internal- segment 0/n distance is %f\n", result);
2631 : #endif
2632 :
2633 : /* check distances for other segments */
2634 12482 : for (i = 0; (i < poly->npts - 1); i++)
2635 : {
2636 9358 : seg.p[0].x = poly->p[i].x;
2637 9358 : seg.p[0].y = poly->p[i].y;
2638 9358 : seg.p[1].x = poly->p[i + 1].x;
2639 9358 : seg.p[1].y = poly->p[i + 1].y;
2640 9358 : d = dist_ps_internal(pt, &seg);
2641 : #ifdef GEODEBUG
2642 : printf("dist_ppoly_internal- segment %d distance is %f\n", (i + 1), d);
2643 : #endif
2644 9358 : if (d < result)
2645 1 : result = d;
2646 : }
2647 :
2648 3124 : return result;
2649 : }
2650 :
2651 :
2652 : /*---------------------------------------------------------------------
2653 : * interpt_
2654 : * Intersection point of objects.
2655 : * We choose to ignore the "point" of intersection between
2656 : * lines and boxes, since there are typically two.
2657 : *-------------------------------------------------------------------*/
2658 :
2659 : /* Get intersection point of lseg and line; returns NULL if no intersection */
2660 : static Point *
2661 48249 : interpt_sl(LSEG *lseg, LINE *line)
2662 : {
2663 : LINE tmp;
2664 : Point *p;
2665 :
2666 48249 : line_construct_pts(&tmp, &lseg->p[0], &lseg->p[1]);
2667 48249 : p = line_interpt_internal(&tmp, line);
2668 : #ifdef GEODEBUG
2669 : printf("interpt_sl- segment is (%.*g %.*g) (%.*g %.*g)\n",
2670 : DBL_DIG, lseg->p[0].x, DBL_DIG, lseg->p[0].y, DBL_DIG, lseg->p[1].x, DBL_DIG, lseg->p[1].y);
2671 : printf("interpt_sl- segment becomes line A=%.*g B=%.*g C=%.*g\n",
2672 : DBL_DIG, tmp.A, DBL_DIG, tmp.B, DBL_DIG, tmp.C);
2673 : #endif
2674 48249 : if (PointerIsValid(p))
2675 : {
2676 : #ifdef GEODEBUG
2677 : printf("interpt_sl- intersection point is (%.*g %.*g)\n", DBL_DIG, p->x, DBL_DIG, p->y);
2678 : #endif
2679 46638 : if (on_ps_internal(p, lseg))
2680 : {
2681 : #ifdef GEODEBUG
2682 : printf("interpt_sl- intersection point is on segment\n");
2683 : #endif
2684 : }
2685 : else
2686 34326 : p = NULL;
2687 : }
2688 :
2689 48249 : return p;
2690 : }
2691 :
2692 : /* variant: just indicate if intersection point exists */
2693 : static bool
2694 15 : has_interpt_sl(LSEG *lseg, LINE *line)
2695 : {
2696 : Point *tmp;
2697 :
2698 15 : tmp = interpt_sl(lseg, line);
2699 15 : if (tmp)
2700 2 : return true;
2701 13 : return false;
2702 : }
2703 :
2704 : /*---------------------------------------------------------------------
2705 : * close_
2706 : * Point of closest proximity between objects.
2707 : *-------------------------------------------------------------------*/
2708 :
2709 : /* close_pl -
2710 : * The intersection point of a perpendicular of the line
2711 : * through the point.
2712 : */
2713 : Datum
2714 1 : close_pl(PG_FUNCTION_ARGS)
2715 : {
2716 1 : Point *pt = PG_GETARG_POINT_P(0);
2717 1 : LINE *line = PG_GETARG_LINE_P(1);
2718 : Point *result;
2719 : LINE *tmp;
2720 : double invm;
2721 :
2722 1 : result = (Point *) palloc(sizeof(Point));
2723 :
2724 1 : if (FPzero(line->B)) /* vertical? */
2725 : {
2726 0 : result->x = line->C;
2727 0 : result->y = pt->y;
2728 0 : PG_RETURN_POINT_P(result);
2729 : }
2730 1 : if (FPzero(line->A)) /* horizontal? */
2731 : {
2732 0 : result->x = pt->x;
2733 0 : result->y = line->C;
2734 0 : PG_RETURN_POINT_P(result);
2735 : }
2736 : /* drop a perpendicular and find the intersection point */
2737 :
2738 : /* invert and flip the sign on the slope to get a perpendicular */
2739 1 : invm = line->B / line->A;
2740 1 : tmp = line_construct_pm(pt, invm);
2741 1 : result = line_interpt_internal(tmp, line);
2742 1 : Assert(result != NULL);
2743 1 : PG_RETURN_POINT_P(result);
2744 : }
2745 :
2746 :
2747 : /* close_ps()
2748 : * Closest point on line segment to specified point.
2749 : * Take the closest endpoint if the point is left, right,
2750 : * above, or below the segment, otherwise find the intersection
2751 : * point of the segment and its perpendicular through the point.
2752 : *
2753 : * Some tricky code here, relying on boolean expressions
2754 : * evaluating to only zero or one to use as an array index.
2755 : * bug fixes by gthaker@atl.lmco.com; May 1, 1998
2756 : */
2757 : Datum
2758 30 : close_ps(PG_FUNCTION_ARGS)
2759 : {
2760 30 : Point *pt = PG_GETARG_POINT_P(0);
2761 30 : LSEG *lseg = PG_GETARG_LSEG_P(1);
2762 30 : Point *result = NULL;
2763 : LINE *tmp;
2764 : double invm;
2765 : int xh,
2766 : yh;
2767 :
2768 : #ifdef GEODEBUG
2769 : printf("close_sp:pt->x %f pt->y %f\nlseg(0).x %f lseg(0).y %f lseg(1).x %f lseg(1).y %f\n",
2770 : pt->x, pt->y, lseg->p[0].x, lseg->p[0].y,
2771 : lseg->p[1].x, lseg->p[1].y);
2772 : #endif
2773 :
2774 : /* xh (or yh) is the index of upper x( or y) end point of lseg */
2775 : /* !xh (or !yh) is the index of lower x( or y) end point of lseg */
2776 30 : xh = lseg->p[0].x < lseg->p[1].x;
2777 30 : yh = lseg->p[0].y < lseg->p[1].y;
2778 :
2779 30 : if (FPeq(lseg->p[0].x, lseg->p[1].x)) /* vertical? */
2780 : {
2781 : #ifdef GEODEBUG
2782 : printf("close_ps- segment is vertical\n");
2783 : #endif
2784 : /* first check if point is below or above the entire lseg. */
2785 0 : if (pt->y < lseg->p[!yh].y)
2786 0 : result = point_copy(&lseg->p[!yh]); /* below the lseg */
2787 0 : else if (pt->y > lseg->p[yh].y)
2788 0 : result = point_copy(&lseg->p[yh]); /* above the lseg */
2789 0 : if (result != NULL)
2790 0 : PG_RETURN_POINT_P(result);
2791 :
2792 : /* point lines along (to left or right) of the vertical lseg. */
2793 :
2794 0 : result = (Point *) palloc(sizeof(Point));
2795 0 : result->x = lseg->p[0].x;
2796 0 : result->y = pt->y;
2797 0 : PG_RETURN_POINT_P(result);
2798 : }
2799 30 : else if (FPeq(lseg->p[0].y, lseg->p[1].y)) /* horizontal? */
2800 : {
2801 : #ifdef GEODEBUG
2802 : printf("close_ps- segment is horizontal\n");
2803 : #endif
2804 : /* first check if point is left or right of the entire lseg. */
2805 0 : if (pt->x < lseg->p[!xh].x)
2806 0 : result = point_copy(&lseg->p[!xh]); /* left of the lseg */
2807 0 : else if (pt->x > lseg->p[xh].x)
2808 0 : result = point_copy(&lseg->p[xh]); /* right of the lseg */
2809 0 : if (result != NULL)
2810 0 : PG_RETURN_POINT_P(result);
2811 :
2812 : /* point lines along (at top or below) the horiz. lseg. */
2813 0 : result = (Point *) palloc(sizeof(Point));
2814 0 : result->x = pt->x;
2815 0 : result->y = lseg->p[0].y;
2816 0 : PG_RETURN_POINT_P(result);
2817 : }
2818 :
2819 : /*
2820 : * vert. and horiz. cases are down, now check if the closest point is one
2821 : * of the end points or someplace on the lseg.
2822 : */
2823 :
2824 30 : invm = -1.0 / point_sl(&(lseg->p[0]), &(lseg->p[1]));
2825 30 : tmp = line_construct_pm(&lseg->p[!yh], invm); /* lower edge of the
2826 : * "band" */
2827 30 : if (pt->y < (tmp->A * pt->x + tmp->C))
2828 : { /* we are below the lower edge */
2829 11 : result = point_copy(&lseg->p[!yh]); /* below the lseg, take lower end
2830 : * pt */
2831 : #ifdef GEODEBUG
2832 : printf("close_ps below: tmp A %f B %f C %f\n",
2833 : tmp->A, tmp->B, tmp->C);
2834 : #endif
2835 11 : PG_RETURN_POINT_P(result);
2836 : }
2837 19 : tmp = line_construct_pm(&lseg->p[yh], invm); /* upper edge of the
2838 : * "band" */
2839 19 : if (pt->y > (tmp->A * pt->x + tmp->C))
2840 : { /* we are below the lower edge */
2841 7 : result = point_copy(&lseg->p[yh]); /* above the lseg, take higher end
2842 : * pt */
2843 : #ifdef GEODEBUG
2844 : printf("close_ps above: tmp A %f B %f C %f\n",
2845 : tmp->A, tmp->B, tmp->C);
2846 : #endif
2847 7 : PG_RETURN_POINT_P(result);
2848 : }
2849 :
2850 : /*
2851 : * at this point the "normal" from point will hit lseg. The closest point
2852 : * will be somewhere on the lseg
2853 : */
2854 12 : tmp = line_construct_pm(pt, invm);
2855 : #ifdef GEODEBUG
2856 : printf("close_ps- tmp A %f B %f C %f\n",
2857 : tmp->A, tmp->B, tmp->C);
2858 : #endif
2859 12 : result = interpt_sl(lseg, tmp);
2860 :
2861 : /*
2862 : * ordinarily we should always find an intersection point, but that could
2863 : * fail in the presence of NaN coordinates, and perhaps even from simple
2864 : * roundoff issues. Return a SQL NULL if so.
2865 : */
2866 12 : if (result == NULL)
2867 0 : PG_RETURN_NULL();
2868 :
2869 : #ifdef GEODEBUG
2870 : printf("close_ps- result.x %f result.y %f\n", result->x, result->y);
2871 : #endif
2872 12 : PG_RETURN_POINT_P(result);
2873 : }
2874 :
2875 :
2876 : /* close_lseg()
2877 : * Closest point to l1 on l2.
2878 : */
2879 : Datum
2880 0 : close_lseg(PG_FUNCTION_ARGS)
2881 : {
2882 0 : LSEG *l1 = PG_GETARG_LSEG_P(0);
2883 0 : LSEG *l2 = PG_GETARG_LSEG_P(1);
2884 0 : Point *result = NULL;
2885 : Point point;
2886 : double dist;
2887 : double d;
2888 :
2889 0 : d = dist_ps_internal(&l1->p[0], l2);
2890 0 : dist = d;
2891 0 : memcpy(&point, &l1->p[0], sizeof(Point));
2892 :
2893 0 : if ((d = dist_ps_internal(&l1->p[1], l2)) < dist)
2894 : {
2895 0 : dist = d;
2896 0 : memcpy(&point, &l1->p[1], sizeof(Point));
2897 : }
2898 :
2899 0 : if (dist_ps_internal(&l2->p[0], l1) < dist)
2900 : {
2901 0 : result = DatumGetPointP(DirectFunctionCall2(close_ps,
2902 : PointPGetDatum(&l2->p[0]),
2903 : LsegPGetDatum(l1)));
2904 0 : memcpy(&point, result, sizeof(Point));
2905 0 : result = DatumGetPointP(DirectFunctionCall2(close_ps,
2906 : PointPGetDatum(&point),
2907 : LsegPGetDatum(l2)));
2908 : }
2909 :
2910 0 : if (dist_ps_internal(&l2->p[1], l1) < dist)
2911 : {
2912 0 : result = DatumGetPointP(DirectFunctionCall2(close_ps,
2913 : PointPGetDatum(&l2->p[1]),
2914 : LsegPGetDatum(l1)));
2915 0 : memcpy(&point, result, sizeof(Point));
2916 0 : result = DatumGetPointP(DirectFunctionCall2(close_ps,
2917 : PointPGetDatum(&point),
2918 : LsegPGetDatum(l2)));
2919 : }
2920 :
2921 0 : if (result == NULL)
2922 0 : result = point_copy(&point);
2923 :
2924 0 : PG_RETURN_POINT_P(result);
2925 : }
2926 :
2927 : /* close_pb()
2928 : * Closest point on or in box to specified point.
2929 : */
2930 : Datum
2931 0 : close_pb(PG_FUNCTION_ARGS)
2932 : {
2933 0 : Point *pt = PG_GETARG_POINT_P(0);
2934 0 : BOX *box = PG_GETARG_BOX_P(1);
2935 : LSEG lseg,
2936 : seg;
2937 : Point point;
2938 : double dist,
2939 : d;
2940 :
2941 0 : if (DatumGetBool(DirectFunctionCall2(on_pb,
2942 : PointPGetDatum(pt),
2943 : BoxPGetDatum(box))))
2944 0 : PG_RETURN_POINT_P(pt);
2945 :
2946 : /* pairwise check lseg distances */
2947 0 : point.x = box->low.x;
2948 0 : point.y = box->high.y;
2949 0 : statlseg_construct(&lseg, &box->low, &point);
2950 0 : dist = dist_ps_internal(pt, &lseg);
2951 :
2952 0 : statlseg_construct(&seg, &box->high, &point);
2953 0 : if ((d = dist_ps_internal(pt, &seg)) < dist)
2954 : {
2955 0 : dist = d;
2956 0 : memcpy(&lseg, &seg, sizeof(lseg));
2957 : }
2958 :
2959 0 : point.x = box->high.x;
2960 0 : point.y = box->low.y;
2961 0 : statlseg_construct(&seg, &box->low, &point);
2962 0 : if ((d = dist_ps_internal(pt, &seg)) < dist)
2963 : {
2964 0 : dist = d;
2965 0 : memcpy(&lseg, &seg, sizeof(lseg));
2966 : }
2967 :
2968 0 : statlseg_construct(&seg, &box->high, &point);
2969 0 : if ((d = dist_ps_internal(pt, &seg)) < dist)
2970 : {
2971 0 : dist = d;
2972 0 : memcpy(&lseg, &seg, sizeof(lseg));
2973 : }
2974 :
2975 0 : PG_RETURN_DATUM(DirectFunctionCall2(close_ps,
2976 : PointPGetDatum(pt),
2977 : LsegPGetDatum(&lseg)));
2978 : }
2979 :
2980 : /* close_sl()
2981 : * Closest point on line to line segment.
2982 : *
2983 : * XXX THIS CODE IS WRONG
2984 : * The code is actually calculating the point on the line segment
2985 : * which is backwards from the routine naming convention.
2986 : * Copied code to new routine close_ls() but haven't fixed this one yet.
2987 : * - thomas 1998-01-31
2988 : */
2989 : Datum
2990 0 : close_sl(PG_FUNCTION_ARGS)
2991 : {
2992 : #ifdef NOT_USED
2993 : LSEG *lseg = PG_GETARG_LSEG_P(0);
2994 : LINE *line = PG_GETARG_LINE_P(1);
2995 : Point *result;
2996 : float8 d1,
2997 : d2;
2998 :
2999 : result = interpt_sl(lseg, line);
3000 : if (result)
3001 : PG_RETURN_POINT_P(result);
3002 :
3003 : d1 = dist_pl_internal(&lseg->p[0], line);
3004 : d2 = dist_pl_internal(&lseg->p[1], line);
3005 : if (d1 < d2)
3006 : result = point_copy(&lseg->p[0]);
3007 : else
3008 : result = point_copy(&lseg->p[1]);
3009 :
3010 : PG_RETURN_POINT_P(result);
3011 : #endif
3012 :
3013 0 : ereport(ERROR,
3014 : (errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
3015 : errmsg("function \"close_sl\" not implemented")));
3016 :
3017 : PG_RETURN_NULL();
3018 : }
3019 :
3020 : /* close_ls()
3021 : * Closest point on line segment to line.
3022 : */
3023 : Datum
3024 1 : close_ls(PG_FUNCTION_ARGS)
3025 : {
3026 1 : LINE *line = PG_GETARG_LINE_P(0);
3027 1 : LSEG *lseg = PG_GETARG_LSEG_P(1);
3028 : Point *result;
3029 : float8 d1,
3030 : d2;
3031 :
3032 1 : result = interpt_sl(lseg, line);
3033 1 : if (result)
3034 0 : PG_RETURN_POINT_P(result);
3035 :
3036 1 : d1 = dist_pl_internal(&lseg->p[0], line);
3037 1 : d2 = dist_pl_internal(&lseg->p[1], line);
3038 1 : if (d1 < d2)
3039 1 : result = point_copy(&lseg->p[0]);
3040 : else
3041 0 : result = point_copy(&lseg->p[1]);
3042 :
3043 1 : PG_RETURN_POINT_P(result);
3044 : }
3045 :
3046 : /* close_sb()
3047 : * Closest point on or in box to line segment.
3048 : */
3049 : Datum
3050 0 : close_sb(PG_FUNCTION_ARGS)
3051 : {
3052 0 : LSEG *lseg = PG_GETARG_LSEG_P(0);
3053 0 : BOX *box = PG_GETARG_BOX_P(1);
3054 : Point point;
3055 : LSEG bseg,
3056 : seg;
3057 : double dist,
3058 : d;
3059 :
3060 : /* segment intersects box? then just return closest point to center */
3061 0 : if (DatumGetBool(DirectFunctionCall2(inter_sb,
3062 : LsegPGetDatum(lseg),
3063 : BoxPGetDatum(box))))
3064 : {
3065 0 : box_cn(&point, box);
3066 0 : PG_RETURN_DATUM(DirectFunctionCall2(close_ps,
3067 : PointPGetDatum(&point),
3068 : LsegPGetDatum(lseg)));
3069 : }
3070 :
3071 : /* pairwise check lseg distances */
3072 0 : point.x = box->low.x;
3073 0 : point.y = box->high.y;
3074 0 : statlseg_construct(&bseg, &box->low, &point);
3075 0 : dist = lseg_dt(lseg, &bseg);
3076 :
3077 0 : statlseg_construct(&seg, &box->high, &point);
3078 0 : if ((d = lseg_dt(lseg, &seg)) < dist)
3079 : {
3080 0 : dist = d;
3081 0 : memcpy(&bseg, &seg, sizeof(bseg));
3082 : }
3083 :
3084 0 : point.x = box->high.x;
3085 0 : point.y = box->low.y;
3086 0 : statlseg_construct(&seg, &box->low, &point);
3087 0 : if ((d = lseg_dt(lseg, &seg)) < dist)
3088 : {
3089 0 : dist = d;
3090 0 : memcpy(&bseg, &seg, sizeof(bseg));
3091 : }
3092 :
3093 0 : statlseg_construct(&seg, &box->high, &point);
3094 0 : if ((d = lseg_dt(lseg, &seg)) < dist)
3095 : {
3096 0 : dist = d;
3097 0 : memcpy(&bseg, &seg, sizeof(bseg));
3098 : }
3099 :
3100 : /* OK, we now have the closest line segment on the box boundary */
3101 0 : PG_RETURN_DATUM(DirectFunctionCall2(close_lseg,
3102 : LsegPGetDatum(lseg),
3103 : LsegPGetDatum(&bseg)));
3104 : }
3105 :
3106 : Datum
3107 0 : close_lb(PG_FUNCTION_ARGS)
3108 : {
3109 : #ifdef NOT_USED
3110 : LINE *line = PG_GETARG_LINE_P(0);
3111 : BOX *box = PG_GETARG_BOX_P(1);
3112 : #endif
3113 :
3114 : /* think about this one for a while */
3115 0 : ereport(ERROR,
3116 : (errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
3117 : errmsg("function \"close_lb\" not implemented")));
3118 :
3119 : PG_RETURN_NULL();
3120 : }
3121 :
3122 : /*---------------------------------------------------------------------
3123 : * on_
3124 : * Whether one object lies completely within another.
3125 : *-------------------------------------------------------------------*/
3126 :
3127 : /* on_pl -
3128 : * Does the point satisfy the equation?
3129 : */
3130 : Datum
3131 10 : on_pl(PG_FUNCTION_ARGS)
3132 : {
3133 10 : Point *pt = PG_GETARG_POINT_P(0);
3134 10 : LINE *line = PG_GETARG_LINE_P(1);
3135 :
3136 10 : PG_RETURN_BOOL(FPzero(line->A * pt->x + line->B * pt->y + line->C));
3137 : }
3138 :
3139 :
3140 : /* on_ps -
3141 : * Determine colinearity by detecting a triangle inequality.
3142 : * This algorithm seems to behave nicely even with lsb residues - tgl 1997-07-09
3143 : */
3144 : Datum
3145 0 : on_ps(PG_FUNCTION_ARGS)
3146 : {
3147 0 : Point *pt = PG_GETARG_POINT_P(0);
3148 0 : LSEG *lseg = PG_GETARG_LSEG_P(1);
3149 :
3150 0 : PG_RETURN_BOOL(on_ps_internal(pt, lseg));
3151 : }
3152 :
3153 : static bool
3154 60986 : on_ps_internal(Point *pt, LSEG *lseg)
3155 : {
3156 60986 : return FPeq(point_dt(pt, &lseg->p[0]) + point_dt(pt, &lseg->p[1]),
3157 : point_dt(&lseg->p[0], &lseg->p[1]));
3158 : }
3159 :
3160 : Datum
3161 12026 : on_pb(PG_FUNCTION_ARGS)
3162 : {
3163 12026 : Point *pt = PG_GETARG_POINT_P(0);
3164 12026 : BOX *box = PG_GETARG_BOX_P(1);
3165 :
3166 12026 : PG_RETURN_BOOL(pt->x <= box->high.x && pt->x >= box->low.x &&
3167 : pt->y <= box->high.y && pt->y >= box->low.y);
3168 : }
3169 :
3170 : Datum
3171 20873 : box_contain_pt(PG_FUNCTION_ARGS)
3172 : {
3173 20873 : BOX *box = PG_GETARG_BOX_P(0);
3174 20873 : Point *pt = PG_GETARG_POINT_P(1);
3175 :
3176 20873 : PG_RETURN_BOOL(pt->x <= box->high.x && pt->x >= box->low.x &&
3177 : pt->y <= box->high.y && pt->y >= box->low.y);
3178 : }
3179 :
3180 : /* on_ppath -
3181 : * Whether a point lies within (on) a polyline.
3182 : * If open, we have to (groan) check each segment.
3183 : * (uses same algorithm as for point intersecting segment - tgl 1997-07-09)
3184 : * If closed, we use the old O(n) ray method for point-in-polygon.
3185 : * The ray is horizontal, from pt out to the right.
3186 : * Each segment that crosses the ray counts as an
3187 : * intersection; note that an endpoint or edge may touch
3188 : * but not cross.
3189 : * (we can do p-in-p in lg(n), but it takes preprocessing)
3190 : */
3191 : Datum
3192 6 : on_ppath(PG_FUNCTION_ARGS)
3193 : {
3194 6 : Point *pt = PG_GETARG_POINT_P(0);
3195 6 : PATH *path = PG_GETARG_PATH_P(1);
3196 : int i,
3197 : n;
3198 : double a,
3199 : b;
3200 :
3201 : /*-- OPEN --*/
3202 6 : if (!path->closed)
3203 : {
3204 6 : n = path->npts - 1;
3205 6 : a = point_dt(pt, &path->p[0]);
3206 14 : for (i = 0; i < n; i++)
3207 : {
3208 10 : b = point_dt(pt, &path->p[i + 1]);
3209 10 : if (FPeq(a + b,
3210 : point_dt(&path->p[i], &path->p[i + 1])))
3211 2 : PG_RETURN_BOOL(true);
3212 8 : a = b;
3213 : }
3214 4 : PG_RETURN_BOOL(false);
3215 : }
3216 :
3217 : /*-- CLOSED --*/
3218 0 : PG_RETURN_BOOL(point_inside(pt, path->npts, path->p) != 0);
3219 : }
3220 :
3221 : Datum
3222 4 : on_sl(PG_FUNCTION_ARGS)
3223 : {
3224 4 : LSEG *lseg = PG_GETARG_LSEG_P(0);
3225 4 : LINE *line = PG_GETARG_LINE_P(1);
3226 :
3227 4 : PG_RETURN_BOOL(DatumGetBool(DirectFunctionCall2(on_pl,
3228 : PointPGetDatum(&lseg->p[0]),
3229 : LinePGetDatum(line))) &&
3230 : DatumGetBool(DirectFunctionCall2(on_pl,
3231 : PointPGetDatum(&lseg->p[1]),
3232 : LinePGetDatum(line))));
3233 : }
3234 :
3235 : Datum
3236 0 : on_sb(PG_FUNCTION_ARGS)
3237 : {
3238 0 : LSEG *lseg = PG_GETARG_LSEG_P(0);
3239 0 : BOX *box = PG_GETARG_BOX_P(1);
3240 :
3241 0 : PG_RETURN_BOOL(DatumGetBool(DirectFunctionCall2(on_pb,
3242 : PointPGetDatum(&lseg->p[0]),
3243 : BoxPGetDatum(box))) &&
3244 : DatumGetBool(DirectFunctionCall2(on_pb,
3245 : PointPGetDatum(&lseg->p[1]),
3246 : BoxPGetDatum(box))));
3247 : }
3248 :
3249 : /*---------------------------------------------------------------------
3250 : * inter_
3251 : * Whether one object intersects another.
3252 : *-------------------------------------------------------------------*/
3253 :
3254 : Datum
3255 2 : inter_sl(PG_FUNCTION_ARGS)
3256 : {
3257 2 : LSEG *lseg = PG_GETARG_LSEG_P(0);
3258 2 : LINE *line = PG_GETARG_LINE_P(1);
3259 :
3260 2 : PG_RETURN_BOOL(has_interpt_sl(lseg, line));
3261 : }
3262 :
3263 : /* inter_sb()
3264 : * Do line segment and box intersect?
3265 : *
3266 : * Segment completely inside box counts as intersection.
3267 : * If you want only segments crossing box boundaries,
3268 : * try converting box to path first.
3269 : *
3270 : * Optimize for non-intersection by checking for box intersection first.
3271 : * - thomas 1998-01-30
3272 : */
3273 : Datum
3274 0 : inter_sb(PG_FUNCTION_ARGS)
3275 : {
3276 0 : LSEG *lseg = PG_GETARG_LSEG_P(0);
3277 0 : BOX *box = PG_GETARG_BOX_P(1);
3278 : BOX lbox;
3279 : LSEG bseg;
3280 : Point point;
3281 :
3282 0 : lbox.low.x = Min(lseg->p[0].x, lseg->p[1].x);
3283 0 : lbox.low.y = Min(lseg->p[0].y, lseg->p[1].y);
3284 0 : lbox.high.x = Max(lseg->p[0].x, lseg->p[1].x);
3285 0 : lbox.high.y = Max(lseg->p[0].y, lseg->p[1].y);
3286 :
3287 : /* nothing close to overlap? then not going to intersect */
3288 0 : if (!box_ov(&lbox, box))
3289 0 : PG_RETURN_BOOL(false);
3290 :
3291 : /* an endpoint of segment is inside box? then clearly intersects */
3292 0 : if (DatumGetBool(DirectFunctionCall2(on_pb,
3293 : PointPGetDatum(&lseg->p[0]),
3294 0 : BoxPGetDatum(box))) ||
3295 0 : DatumGetBool(DirectFunctionCall2(on_pb,
3296 : PointPGetDatum(&lseg->p[1]),
3297 : BoxPGetDatum(box))))
3298 0 : PG_RETURN_BOOL(true);
3299 :
3300 : /* pairwise check lseg intersections */
3301 0 : point.x = box->low.x;
3302 0 : point.y = box->high.y;
3303 0 : statlseg_construct(&bseg, &box->low, &point);
3304 0 : if (lseg_intersect_internal(&bseg, lseg))
3305 0 : PG_RETURN_BOOL(true);
3306 :
3307 0 : statlseg_construct(&bseg, &box->high, &point);
3308 0 : if (lseg_intersect_internal(&bseg, lseg))
3309 0 : PG_RETURN_BOOL(true);
3310 :
3311 0 : point.x = box->high.x;
3312 0 : point.y = box->low.y;
3313 0 : statlseg_construct(&bseg, &box->low, &point);
3314 0 : if (lseg_intersect_internal(&bseg, lseg))
3315 0 : PG_RETURN_BOOL(true);
3316 :
3317 0 : statlseg_construct(&bseg, &box->high, &point);
3318 0 : if (lseg_intersect_internal(&bseg, lseg))
3319 0 : PG_RETURN_BOOL(true);
3320 :
3321 : /* if we dropped through, no two segs intersected */
3322 0 : PG_RETURN_BOOL(false);
3323 : }
3324 :
3325 : /* inter_lb()
3326 : * Do line and box intersect?
3327 : */
3328 : Datum
3329 2 : inter_lb(PG_FUNCTION_ARGS)
3330 : {
3331 2 : LINE *line = PG_GETARG_LINE_P(0);
3332 2 : BOX *box = PG_GETARG_BOX_P(1);
3333 : LSEG bseg;
3334 : Point p1,
3335 : p2;
3336 :
3337 : /* pairwise check lseg intersections */
3338 2 : p1.x = box->low.x;
3339 2 : p1.y = box->low.y;
3340 2 : p2.x = box->low.x;
3341 2 : p2.y = box->high.y;
3342 2 : statlseg_construct(&bseg, &p1, &p2);
3343 2 : if (has_interpt_sl(&bseg, line))
3344 1 : PG_RETURN_BOOL(true);
3345 1 : p1.x = box->high.x;
3346 1 : p1.y = box->high.y;
3347 1 : statlseg_construct(&bseg, &p1, &p2);
3348 1 : if (has_interpt_sl(&bseg, line))
3349 0 : PG_RETURN_BOOL(true);
3350 1 : p2.x = box->high.x;
3351 1 : p2.y = box->low.y;
3352 1 : statlseg_construct(&bseg, &p1, &p2);
3353 1 : if (has_interpt_sl(&bseg, line))
3354 0 : PG_RETURN_BOOL(true);
3355 1 : p1.x = box->low.x;
3356 1 : p1.y = box->low.y;
3357 1 : statlseg_construct(&bseg, &p1, &p2);
3358 1 : if (has_interpt_sl(&bseg, line))
3359 0 : PG_RETURN_BOOL(true);
3360 :
3361 : /* if we dropped through, no intersection */
3362 1 : PG_RETURN_BOOL(false);
3363 : }
3364 :
3365 : /*------------------------------------------------------------------
3366 : * The following routines define a data type and operator class for
3367 : * POLYGONS .... Part of which (the polygon's bounding box) is built on
3368 : * top of the BOX data type.
3369 : *
3370 : * make_bound_box - create the bounding box for the input polygon
3371 : *------------------------------------------------------------------*/
3372 :
3373 : /*---------------------------------------------------------------------
3374 : * Make the smallest bounding box for the given polygon.
3375 : *---------------------------------------------------------------------*/
3376 : static void
3377 90 : make_bound_box(POLYGON *poly)
3378 : {
3379 : int i;
3380 : double x1,
3381 : y1,
3382 : x2,
3383 : y2;
3384 :
3385 90 : if (poly->npts > 0)
3386 : {
3387 90 : x2 = x1 = poly->p[0].x;
3388 90 : y2 = y1 = poly->p[0].y;
3389 363 : for (i = 1; i < poly->npts; i++)
3390 : {
3391 273 : if (poly->p[i].x < x1)
3392 34 : x1 = poly->p[i].x;
3393 273 : if (poly->p[i].x > x2)
3394 108 : x2 = poly->p[i].x;
3395 273 : if (poly->p[i].y < y1)
3396 68 : y1 = poly->p[i].y;
3397 273 : if (poly->p[i].y > y2)
3398 97 : y2 = poly->p[i].y;
3399 : }
3400 :
3401 90 : box_fill(&(poly->boundbox), x1, x2, y1, y2);
3402 : }
3403 : else
3404 0 : ereport(ERROR,
3405 : (errcode(ERRCODE_INVALID_PARAMETER_VALUE),
3406 : errmsg("cannot create bounding box for empty polygon")));
3407 90 : }
3408 :
3409 : /*------------------------------------------------------------------
3410 : * poly_in - read in the polygon from a string specification
3411 : *
3412 : * External format:
3413 : * "((x0,y0),...,(xn,yn))"
3414 : * "x0,y0,...,xn,yn"
3415 : * also supports the older style "(x1,...,xn,y1,...yn)"
3416 : *------------------------------------------------------------------*/
3417 : Datum
3418 75 : poly_in(PG_FUNCTION_ARGS)
3419 : {
3420 75 : char *str = PG_GETARG_CSTRING(0);
3421 : POLYGON *poly;
3422 : int npts;
3423 : int size;
3424 : int base_size;
3425 : bool isopen;
3426 :
3427 75 : if ((npts = pair_count(str, ',')) <= 0)
3428 4 : ereport(ERROR,
3429 : (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
3430 : errmsg("invalid input syntax for type %s: \"%s\"",
3431 : "polygon", str)));
3432 :
3433 71 : base_size = sizeof(poly->p[0]) * npts;
3434 71 : size = offsetof(POLYGON, p) + base_size;
3435 :
3436 : /* Check for integer overflow */
3437 71 : if (base_size / npts != sizeof(poly->p[0]) || size <= base_size)
3438 0 : ereport(ERROR,
3439 : (errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED),
3440 : errmsg("too many points requested")));
3441 :
3442 71 : poly = (POLYGON *) palloc0(size); /* zero any holes */
3443 :
3444 71 : SET_VARSIZE(poly, size);
3445 71 : poly->npts = npts;
3446 :
3447 71 : path_decode(str, false, npts, &(poly->p[0]), &isopen, NULL, "polygon", str);
3448 :
3449 70 : make_bound_box(poly);
3450 :
3451 70 : PG_RETURN_POLYGON_P(poly);
3452 : }
3453 :
3454 : /*---------------------------------------------------------------
3455 : * poly_out - convert internal POLYGON representation to the
3456 : * character string format "((f8,f8),...,(f8,f8))"
3457 : *---------------------------------------------------------------*/
3458 : Datum
3459 123 : poly_out(PG_FUNCTION_ARGS)
3460 : {
3461 123 : POLYGON *poly = PG_GETARG_POLYGON_P(0);
3462 :
3463 123 : PG_RETURN_CSTRING(path_encode(PATH_CLOSED, poly->npts, poly->p));
3464 : }
3465 :
3466 : /*
3467 : * poly_recv - converts external binary format to polygon
3468 : *
3469 : * External representation is int32 number of points, and the points.
3470 : * We recompute the bounding box on read, instead of trusting it to
3471 : * be valid. (Checking it would take just as long, so may as well
3472 : * omit it from external representation.)
3473 : */
3474 : Datum
3475 0 : poly_recv(PG_FUNCTION_ARGS)
3476 : {
3477 0 : StringInfo buf = (StringInfo) PG_GETARG_POINTER(0);
3478 : POLYGON *poly;
3479 : int32 npts;
3480 : int32 i;
3481 : int size;
3482 :
3483 0 : npts = pq_getmsgint(buf, sizeof(int32));
3484 0 : if (npts <= 0 || npts >= (int32) ((INT_MAX - offsetof(POLYGON, p)) / sizeof(Point)))
3485 0 : ereport(ERROR,
3486 : (errcode(ERRCODE_INVALID_BINARY_REPRESENTATION),
3487 : errmsg("invalid number of points in external \"polygon\" value")));
3488 :
3489 0 : size = offsetof(POLYGON, p) + sizeof(poly->p[0]) * npts;
3490 0 : poly = (POLYGON *) palloc0(size); /* zero any holes */
3491 :
3492 0 : SET_VARSIZE(poly, size);
3493 0 : poly->npts = npts;
3494 :
3495 0 : for (i = 0; i < npts; i++)
3496 : {
3497 0 : poly->p[i].x = pq_getmsgfloat8(buf);
3498 0 : poly->p[i].y = pq_getmsgfloat8(buf);
3499 : }
3500 :
3501 0 : make_bound_box(poly);
3502 :
3503 0 : PG_RETURN_POLYGON_P(poly);
3504 : }
3505 :
3506 : /*
3507 : * poly_send - converts polygon to binary format
3508 : */
3509 : Datum
3510 0 : poly_send(PG_FUNCTION_ARGS)
3511 : {
3512 0 : POLYGON *poly = PG_GETARG_POLYGON_P(0);
3513 : StringInfoData buf;
3514 : int32 i;
3515 :
3516 0 : pq_begintypsend(&buf);
3517 0 : pq_sendint(&buf, poly->npts, sizeof(int32));
3518 0 : for (i = 0; i < poly->npts; i++)
3519 : {
3520 0 : pq_sendfloat8(&buf, poly->p[i].x);
3521 0 : pq_sendfloat8(&buf, poly->p[i].y);
3522 : }
3523 0 : PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
3524 : }
3525 :
3526 :
3527 : /*-------------------------------------------------------
3528 : * Is polygon A strictly left of polygon B? i.e. is
3529 : * the right most point of A left of the left most point
3530 : * of B?
3531 : *-------------------------------------------------------*/
3532 : Datum
3533 6 : poly_left(PG_FUNCTION_ARGS)
3534 : {
3535 6 : POLYGON *polya = PG_GETARG_POLYGON_P(0);
3536 6 : POLYGON *polyb = PG_GETARG_POLYGON_P(1);
3537 : bool result;
3538 :
3539 6 : result = polya->boundbox.high.x < polyb->boundbox.low.x;
3540 :
3541 : /*
3542 : * Avoid leaking memory for toasted inputs ... needed for rtree indexes
3543 : */
3544 6 : PG_FREE_IF_COPY(polya, 0);
3545 6 : PG_FREE_IF_COPY(polyb, 1);
3546 :
3547 6 : PG_RETURN_BOOL(result);
3548 : }
3549 :
3550 : /*-------------------------------------------------------
3551 : * Is polygon A overlapping or left of polygon B? i.e. is
3552 : * the right most point of A at or left of the right most point
3553 : * of B?
3554 : *-------------------------------------------------------*/
3555 : Datum
3556 4 : poly_overleft(PG_FUNCTION_ARGS)
3557 : {
3558 4 : POLYGON *polya = PG_GETARG_POLYGON_P(0);
3559 4 : POLYGON *polyb = PG_GETARG_POLYGON_P(1);
3560 : bool result;
3561 :
3562 4 : result = polya->boundbox.high.x <= polyb->boundbox.high.x;
3563 :
3564 : /*
3565 : * Avoid leaking memory for toasted inputs ... needed for rtree indexes
3566 : */
3567 4 : PG_FREE_IF_COPY(polya, 0);
3568 4 : PG_FREE_IF_COPY(polyb, 1);
3569 :
3570 4 : PG_RETURN_BOOL(result);
3571 : }
3572 :
3573 : /*-------------------------------------------------------
3574 : * Is polygon A strictly right of polygon B? i.e. is
3575 : * the left most point of A right of the right most point
3576 : * of B?
3577 : *-------------------------------------------------------*/
3578 : Datum
3579 5 : poly_right(PG_FUNCTION_ARGS)
3580 : {
3581 5 : POLYGON *polya = PG_GETARG_POLYGON_P(0);
3582 5 : POLYGON *polyb = PG_GETARG_POLYGON_P(1);
3583 : bool result;
3584 :
3585 5 : result = polya->boundbox.low.x > polyb->boundbox.high.x;
3586 :
3587 : /*
3588 : * Avoid leaking memory for toasted inputs ... needed for rtree indexes
3589 : */
3590 5 : PG_FREE_IF_COPY(polya, 0);
3591 5 : PG_FREE_IF_COPY(polyb, 1);
3592 :
3593 5 : PG_RETURN_BOOL(result);
3594 : }
3595 :
3596 : /*-------------------------------------------------------
3597 : * Is polygon A overlapping or right of polygon B? i.e. is
3598 : * the left most point of A at or right of the left most point
3599 : * of B?
3600 : *-------------------------------------------------------*/
3601 : Datum
3602 5 : poly_overright(PG_FUNCTION_ARGS)
3603 : {
3604 5 : POLYGON *polya = PG_GETARG_POLYGON_P(0);
3605 5 : POLYGON *polyb = PG_GETARG_POLYGON_P(1);
3606 : bool result;
3607 :
3608 5 : result = polya->boundbox.low.x >= polyb->boundbox.low.x;
3609 :
3610 : /*
3611 : * Avoid leaking memory for toasted inputs ... needed for rtree indexes
3612 : */
3613 5 : PG_FREE_IF_COPY(polya, 0);
3614 5 : PG_FREE_IF_COPY(polyb, 1);
3615 :
3616 5 : PG_RETURN_BOOL(result);
3617 : }
3618 :
3619 : /*-------------------------------------------------------
3620 : * Is polygon A strictly below polygon B? i.e. is
3621 : * the upper most point of A below the lower most point
3622 : * of B?
3623 : *-------------------------------------------------------*/
3624 : Datum
3625 0 : poly_below(PG_FUNCTION_ARGS)
3626 : {
3627 0 : POLYGON *polya = PG_GETARG_POLYGON_P(0);
3628 0 : POLYGON *polyb = PG_GETARG_POLYGON_P(1);
3629 : bool result;
3630 :
3631 0 : result = polya->boundbox.high.y < polyb->boundbox.low.y;
3632 :
3633 : /*
3634 : * Avoid leaking memory for toasted inputs ... needed for rtree indexes
3635 : */
3636 0 : PG_FREE_IF_COPY(polya, 0);
3637 0 : PG_FREE_IF_COPY(polyb, 1);
3638 :
3639 0 : PG_RETURN_BOOL(result);
3640 : }
3641 :
3642 : /*-------------------------------------------------------
3643 : * Is polygon A overlapping or below polygon B? i.e. is
3644 : * the upper most point of A at or below the upper most point
3645 : * of B?
3646 : *-------------------------------------------------------*/
3647 : Datum
3648 0 : poly_overbelow(PG_FUNCTION_ARGS)
3649 : {
3650 0 : POLYGON *polya = PG_GETARG_POLYGON_P(0);
3651 0 : POLYGON *polyb = PG_GETARG_POLYGON_P(1);
3652 : bool result;
3653 :
3654 0 : result = polya->boundbox.high.y <= polyb->boundbox.high.y;
3655 :
3656 : /*
3657 : * Avoid leaking memory for toasted inputs ... needed for rtree indexes
3658 : */
3659 0 : PG_FREE_IF_COPY(polya, 0);
3660 0 : PG_FREE_IF_COPY(polyb, 1);
3661 :
3662 0 : PG_RETURN_BOOL(result);
3663 : }
3664 :
3665 : /*-------------------------------------------------------
3666 : * Is polygon A strictly above polygon B? i.e. is
3667 : * the lower most point of A above the upper most point
3668 : * of B?
3669 : *-------------------------------------------------------*/
3670 : Datum
3671 0 : poly_above(PG_FUNCTION_ARGS)
3672 : {
3673 0 : POLYGON *polya = PG_GETARG_POLYGON_P(0);
3674 0 : POLYGON *polyb = PG_GETARG_POLYGON_P(1);
3675 : bool result;
3676 :
3677 0 : result = polya->boundbox.low.y > polyb->boundbox.high.y;
3678 :
3679 : /*
3680 : * Avoid leaking memory for toasted inputs ... needed for rtree indexes
3681 : */
3682 0 : PG_FREE_IF_COPY(polya, 0);
3683 0 : PG_FREE_IF_COPY(polyb, 1);
3684 :
3685 0 : PG_RETURN_BOOL(result);
3686 : }
3687 :
3688 : /*-------------------------------------------------------
3689 : * Is polygon A overlapping or above polygon B? i.e. is
3690 : * the lower most point of A at or above the lower most point
3691 : * of B?
3692 : *-------------------------------------------------------*/
3693 : Datum
3694 0 : poly_overabove(PG_FUNCTION_ARGS)
3695 : {
3696 0 : POLYGON *polya = PG_GETARG_POLYGON_P(0);
3697 0 : POLYGON *polyb = PG_GETARG_POLYGON_P(1);
3698 : bool result;
3699 :
3700 0 : result = polya->boundbox.low.y >= polyb->boundbox.low.y;
3701 :
3702 : /*
3703 : * Avoid leaking memory for toasted inputs ... needed for rtree indexes
3704 : */
3705 0 : PG_FREE_IF_COPY(polya, 0);
3706 0 : PG_FREE_IF_COPY(polyb, 1);
3707 :
3708 0 : PG_RETURN_BOOL(result);
3709 : }
3710 :
3711 :
3712 : /*-------------------------------------------------------
3713 : * Is polygon A the same as polygon B? i.e. are all the
3714 : * points the same?
3715 : * Check all points for matches in both forward and reverse
3716 : * direction since polygons are non-directional and are
3717 : * closed shapes.
3718 : *-------------------------------------------------------*/
3719 : Datum
3720 5 : poly_same(PG_FUNCTION_ARGS)
3721 : {
3722 5 : POLYGON *polya = PG_GETARG_POLYGON_P(0);
3723 5 : POLYGON *polyb = PG_GETARG_POLYGON_P(1);
3724 : bool result;
3725 :
3726 5 : if (polya->npts != polyb->npts)
3727 2 : result = false;
3728 : else
3729 3 : result = plist_same(polya->npts, polya->p, polyb->p);
3730 :
3731 : /*
3732 : * Avoid leaking memory for toasted inputs ... needed for rtree indexes
3733 : */
3734 5 : PG_FREE_IF_COPY(polya, 0);
3735 5 : PG_FREE_IF_COPY(polyb, 1);
3736 :
3737 5 : PG_RETURN_BOOL(result);
3738 : }
3739 :
3740 : /*-----------------------------------------------------------------
3741 : * Determine if polygon A overlaps polygon B
3742 : *-----------------------------------------------------------------*/
3743 : Datum
3744 3114 : poly_overlap(PG_FUNCTION_ARGS)
3745 : {
3746 3114 : POLYGON *polya = PG_GETARG_POLYGON_P(0);
3747 3114 : POLYGON *polyb = PG_GETARG_POLYGON_P(1);
3748 : bool result;
3749 :
3750 : /* Quick check by bounding box */
3751 9342 : result = (polya->npts > 0 && polyb->npts > 0 &&
3752 6228 : box_ov(&polya->boundbox, &polyb->boundbox)) ? true : false;
3753 :
3754 : /*
3755 : * Brute-force algorithm - try to find intersected edges, if so then
3756 : * polygons are overlapped else check is one polygon inside other or not
3757 : * by testing single point of them.
3758 : */
3759 3114 : if (result)
3760 : {
3761 : int ia,
3762 : ib;
3763 : LSEG sa,
3764 : sb;
3765 :
3766 : /* Init first of polya's edge with last point */
3767 14 : sa.p[0] = polya->p[polya->npts - 1];
3768 14 : result = false;
3769 :
3770 48 : for (ia = 0; ia < polya->npts && result == false; ia++)
3771 : {
3772 : /* Second point of polya's edge is a current one */
3773 34 : sa.p[1] = polya->p[ia];
3774 :
3775 : /* Init first of polyb's edge with last point */
3776 34 : sb.p[0] = polyb->p[polyb->npts - 1];
3777 :
3778 120 : for (ib = 0; ib < polyb->npts && result == false; ib++)
3779 : {
3780 86 : sb.p[1] = polyb->p[ib];
3781 86 : result = lseg_intersect_internal(&sa, &sb);
3782 86 : sb.p[0] = sb.p[1];
3783 : }
3784 :
3785 : /*
3786 : * move current endpoint to the first point of next edge
3787 : */
3788 34 : sa.p[0] = sa.p[1];
3789 : }
3790 :
3791 14 : if (result == false)
3792 : {
3793 10 : result = (point_inside(polya->p, polyb->npts, polyb->p)
3794 10 : ||
3795 5 : point_inside(polyb->p, polya->npts, polya->p));
3796 : }
3797 : }
3798 :
3799 : /*
3800 : * Avoid leaking memory for toasted inputs ... needed for rtree indexes
3801 : */
3802 3114 : PG_FREE_IF_COPY(polya, 0);
3803 3114 : PG_FREE_IF_COPY(polyb, 1);
3804 :
3805 3114 : PG_RETURN_BOOL(result);
3806 : }
3807 :
3808 : /*
3809 : * Tests special kind of segment for in/out of polygon.
3810 : * Special kind means:
3811 : * - point a should be on segment s
3812 : * - segment (a,b) should not be contained by s
3813 : * Returns true if:
3814 : * - segment (a,b) is collinear to s and (a,b) is in polygon
3815 : * - segment (a,b) s not collinear to s. Note: that doesn't
3816 : * mean that segment is in polygon!
3817 : */
3818 :
3819 : static bool
3820 16 : touched_lseg_inside_poly(Point *a, Point *b, LSEG *s, POLYGON *poly, int start)
3821 : {
3822 : /* point a is on s, b is not */
3823 : LSEG t;
3824 :
3825 16 : t.p[0] = *a;
3826 16 : t.p[1] = *b;
3827 :
3828 : #define POINTEQ(pt1, pt2) (FPeq((pt1)->x, (pt2)->x) && FPeq((pt1)->y, (pt2)->y))
3829 16 : if (POINTEQ(a, s->p))
3830 : {
3831 8 : if (on_ps_internal(s->p + 1, &t))
3832 0 : return lseg_inside_poly(b, s->p + 1, poly, start);
3833 : }
3834 12 : else if (POINTEQ(a, s->p + 1))
3835 : {
3836 8 : if (on_ps_internal(s->p, &t))
3837 0 : return lseg_inside_poly(b, s->p, poly, start);
3838 : }
3839 8 : else if (on_ps_internal(s->p, &t))
3840 : {
3841 0 : return lseg_inside_poly(b, s->p, poly, start);
3842 : }
3843 8 : else if (on_ps_internal(s->p + 1, &t))
3844 : {
3845 0 : return lseg_inside_poly(b, s->p + 1, poly, start);
3846 : }
3847 :
3848 16 : return true; /* may be not true, but that will check later */
3849 : }
3850 :
3851 : /*
3852 : * Returns true if segment (a,b) is in polygon, option
3853 : * start is used for optimization - function checks
3854 : * polygon's edges started from start
3855 : */
3856 : static bool
3857 37 : lseg_inside_poly(Point *a, Point *b, POLYGON *poly, int start)
3858 : {
3859 : LSEG s,
3860 : t;
3861 : int i;
3862 37 : bool res = true,
3863 37 : intersection = false;
3864 :
3865 37 : t.p[0] = *a;
3866 37 : t.p[1] = *b;
3867 37 : s.p[0] = poly->p[(start == 0) ? (poly->npts - 1) : (start - 1)];
3868 :
3869 129 : for (i = start; i < poly->npts && res; i++)
3870 : {
3871 : Point *interpt;
3872 :
3873 102 : CHECK_FOR_INTERRUPTS();
3874 :
3875 102 : s.p[1] = poly->p[i];
3876 :
3877 102 : if (on_ps_internal(t.p, &s))
3878 : {
3879 18 : if (on_ps_internal(t.p + 1, &s))
3880 10 : return true; /* t is contained by s */
3881 :
3882 : /* Y-cross */
3883 8 : res = touched_lseg_inside_poly(t.p, t.p + 1, &s, poly, i + 1);
3884 : }
3885 84 : else if (on_ps_internal(t.p + 1, &s))
3886 : {
3887 : /* Y-cross */
3888 8 : res = touched_lseg_inside_poly(t.p + 1, t.p, &s, poly, i + 1);
3889 : }
3890 76 : else if ((interpt = lseg_interpt_internal(&t, &s)) != NULL)
3891 : {
3892 : /*
3893 : * segments are X-crossing, go to check each subsegment
3894 : */
3895 :
3896 6 : intersection = true;
3897 6 : res = lseg_inside_poly(t.p, interpt, poly, i + 1);
3898 6 : if (res)
3899 5 : res = lseg_inside_poly(t.p + 1, interpt, poly, i + 1);
3900 6 : pfree(interpt);
3901 : }
3902 :
3903 92 : s.p[0] = s.p[1];
3904 : }
3905 :
3906 27 : if (res && !intersection)
3907 : {
3908 : Point p;
3909 :
3910 : /*
3911 : * if X-intersection wasn't found then check central point of tested
3912 : * segment. In opposite case we already check all subsegments
3913 : */
3914 22 : p.x = (t.p[0].x + t.p[1].x) / 2.0;
3915 22 : p.y = (t.p[0].y + t.p[1].y) / 2.0;
3916 :
3917 22 : res = point_inside(&p, poly->npts, poly->p);
3918 : }
3919 :
3920 27 : return res;
3921 : }
3922 :
3923 : /*-----------------------------------------------------------------
3924 : * Determine if polygon A contains polygon B.
3925 : *-----------------------------------------------------------------*/
3926 : Datum
3927 20 : poly_contain(PG_FUNCTION_ARGS)
3928 : {
3929 20 : POLYGON *polya = PG_GETARG_POLYGON_P(0);
3930 20 : POLYGON *polyb = PG_GETARG_POLYGON_P(1);
3931 : bool result;
3932 :
3933 : /*
3934 : * Quick check to see if bounding box is contained.
3935 : */
3936 40 : if (polya->npts > 0 && polyb->npts > 0 &&
3937 20 : DatumGetBool(DirectFunctionCall2(box_contain,
3938 : BoxPGetDatum(&polya->boundbox),
3939 : BoxPGetDatum(&polyb->boundbox))))
3940 10 : {
3941 : int i;
3942 : LSEG s;
3943 :
3944 10 : s.p[0] = polyb->p[polyb->npts - 1];
3945 10 : result = true;
3946 :
3947 36 : for (i = 0; i < polyb->npts && result; i++)
3948 : {
3949 26 : s.p[1] = polyb->p[i];
3950 26 : result = lseg_inside_poly(s.p, s.p + 1, polya, 0);
3951 26 : s.p[0] = s.p[1];
3952 : }
3953 : }
3954 : else
3955 : {
3956 10 : result = false;
3957 : }
3958 :
3959 : /*
3960 : * Avoid leaking memory for toasted inputs ... needed for rtree indexes
3961 : */
3962 20 : PG_FREE_IF_COPY(polya, 0);
3963 20 : PG_FREE_IF_COPY(polyb, 1);
3964 :
3965 20 : PG_RETURN_BOOL(result);
3966 : }
3967 :
3968 :
3969 : /*-----------------------------------------------------------------
3970 : * Determine if polygon A is contained by polygon B
3971 : *-----------------------------------------------------------------*/
3972 : Datum
3973 5 : poly_contained(PG_FUNCTION_ARGS)
3974 : {
3975 5 : Datum polya = PG_GETARG_DATUM(0);
3976 5 : Datum polyb = PG_GETARG_DATUM(1);
3977 :
3978 : /* Just switch the arguments and pass it off to poly_contain */
3979 5 : PG_RETURN_DATUM(DirectFunctionCall2(poly_contain, polyb, polya));
3980 : }
3981 :
3982 :
3983 : Datum
3984 27 : poly_contain_pt(PG_FUNCTION_ARGS)
3985 : {
3986 27 : POLYGON *poly = PG_GETARG_POLYGON_P(0);
3987 27 : Point *p = PG_GETARG_POINT_P(1);
3988 :
3989 27 : PG_RETURN_BOOL(point_inside(p, poly->npts, poly->p) != 0);
3990 : }
3991 :
3992 : Datum
3993 30 : pt_contained_poly(PG_FUNCTION_ARGS)
3994 : {
3995 30 : Point *p = PG_GETARG_POINT_P(0);
3996 30 : POLYGON *poly = PG_GETARG_POLYGON_P(1);
3997 :
3998 30 : PG_RETURN_BOOL(point_inside(p, poly->npts, poly->p) != 0);
3999 : }
4000 :
4001 :
4002 : Datum
4003 0 : poly_distance(PG_FUNCTION_ARGS)
4004 : {
4005 : #ifdef NOT_USED
4006 : POLYGON *polya = PG_GETARG_POLYGON_P(0);
4007 : POLYGON *polyb = PG_GETARG_POLYGON_P(1);
4008 : #endif
4009 :
4010 0 : ereport(ERROR,
4011 : (errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
4012 : errmsg("function \"poly_distance\" not implemented")));
4013 :
4014 : PG_RETURN_NULL();
4015 : }
4016 :
4017 :
4018 : /***********************************************************************
4019 : **
4020 : ** Routines for 2D points.
4021 : **
4022 : ***********************************************************************/
4023 :
4024 : Datum
4025 110375 : construct_point(PG_FUNCTION_ARGS)
4026 : {
4027 110375 : float8 x = PG_GETARG_FLOAT8(0);
4028 110375 : float8 y = PG_GETARG_FLOAT8(1);
4029 :
4030 110375 : PG_RETURN_POINT_P(point_construct(x, y));
4031 : }
4032 :
4033 : Datum
4034 0 : point_add(PG_FUNCTION_ARGS)
4035 : {
4036 0 : Point *p1 = PG_GETARG_POINT_P(0);
4037 0 : Point *p2 = PG_GETARG_POINT_P(1);
4038 : Point *result;
4039 :
4040 0 : result = (Point *) palloc(sizeof(Point));
4041 :
4042 0 : result->x = (p1->x + p2->x);
4043 0 : result->y = (p1->y + p2->y);
4044 :
4045 0 : PG_RETURN_POINT_P(result);
4046 : }
4047 :
4048 : Datum
4049 0 : point_sub(PG_FUNCTION_ARGS)
4050 : {
4051 0 : Point *p1 = PG_GETARG_POINT_P(0);
4052 0 : Point *p2 = PG_GETARG_POINT_P(1);
4053 : Point *result;
4054 :
4055 0 : result = (Point *) palloc(sizeof(Point));
4056 :
4057 0 : result->x = (p1->x - p2->x);
4058 0 : result->y = (p1->y - p2->y);
4059 :
4060 0 : PG_RETURN_POINT_P(result);
4061 : }
4062 :
4063 : Datum
4064 66 : point_mul(PG_FUNCTION_ARGS)
4065 : {
4066 66 : Point *p1 = PG_GETARG_POINT_P(0);
4067 66 : Point *p2 = PG_GETARG_POINT_P(1);
4068 : Point *result;
4069 :
4070 66 : result = (Point *) palloc(sizeof(Point));
4071 :
4072 66 : result->x = (p1->x * p2->x) - (p1->y * p2->y);
4073 66 : result->y = (p1->x * p2->y) + (p1->y * p2->x);
4074 :
4075 66 : PG_RETURN_POINT_P(result);
4076 : }
4077 :
4078 : Datum
4079 40 : point_div(PG_FUNCTION_ARGS)
4080 : {
4081 40 : Point *p1 = PG_GETARG_POINT_P(0);
4082 40 : Point *p2 = PG_GETARG_POINT_P(1);
4083 : Point *result;
4084 : double div;
4085 :
4086 40 : result = (Point *) palloc(sizeof(Point));
4087 :
4088 40 : div = (p2->x * p2->x) + (p2->y * p2->y);
4089 :
4090 40 : if (div == 0.0)
4091 0 : ereport(ERROR,
4092 : (errcode(ERRCODE_DIVISION_BY_ZERO),
4093 : errmsg("division by zero")));
4094 :
4095 40 : result->x = ((p1->x * p2->x) + (p1->y * p2->y)) / div;
4096 40 : result->y = ((p2->x * p1->y) - (p2->y * p1->x)) / div;
4097 :
4098 40 : PG_RETURN_POINT_P(result);
4099 : }
4100 :
4101 :
4102 : /***********************************************************************
4103 : **
4104 : ** Routines for 2D boxes.
4105 : **
4106 : ***********************************************************************/
4107 :
4108 : Datum
4109 20178 : points_box(PG_FUNCTION_ARGS)
4110 : {
4111 20178 : Point *p1 = PG_GETARG_POINT_P(0);
4112 20178 : Point *p2 = PG_GETARG_POINT_P(1);
4113 :
4114 20178 : PG_RETURN_BOX_P(box_construct(p1->x, p2->x, p1->y, p2->y));
4115 : }
4116 :
4117 : Datum
4118 24 : box_add(PG_FUNCTION_ARGS)
4119 : {
4120 24 : BOX *box = PG_GETARG_BOX_P(0);
4121 24 : Point *p = PG_GETARG_POINT_P(1);
4122 :
4123 24 : PG_RETURN_BOX_P(box_construct((box->high.x + p->x),
4124 : (box->low.x + p->x),
4125 : (box->high.y + p->y),
4126 : (box->low.y + p->y)));
4127 : }
4128 :
4129 : Datum
4130 24 : box_sub(PG_FUNCTION_ARGS)
4131 : {
4132 24 : BOX *box = PG_GETARG_BOX_P(0);
4133 24 : Point *p = PG_GETARG_POINT_P(1);
4134 :
4135 24 : PG_RETURN_BOX_P(box_construct((box->high.x - p->x),
4136 : (box->low.x - p->x),
4137 : (box->high.y - p->y),
4138 : (box->low.y - p->y)));
4139 : }
4140 :
4141 : Datum
4142 24 : box_mul(PG_FUNCTION_ARGS)
4143 : {
4144 24 : BOX *box = PG_GETARG_BOX_P(0);
4145 24 : Point *p = PG_GETARG_POINT_P(1);
4146 : BOX *result;
4147 : Point *high,
4148 : *low;
4149 :
4150 24 : high = DatumGetPointP(DirectFunctionCall2(point_mul,
4151 : PointPGetDatum(&box->high),
4152 : PointPGetDatum(p)));
4153 24 : low = DatumGetPointP(DirectFunctionCall2(point_mul,
4154 : PointPGetDatum(&box->low),
4155 : PointPGetDatum(p)));
4156 :
4157 24 : result = box_construct(high->x, low->x, high->y, low->y);
4158 :
4159 24 : PG_RETURN_BOX_P(result);
4160 : }
4161 :
4162 : Datum
4163 20 : box_div(PG_FUNCTION_ARGS)
4164 : {
4165 20 : BOX *box = PG_GETARG_BOX_P(0);
4166 20 : Point *p = PG_GETARG_POINT_P(1);
4167 : BOX *result;
4168 : Point *high,
4169 : *low;
4170 :
4171 20 : high = DatumGetPointP(DirectFunctionCall2(point_div,
4172 : PointPGetDatum(&box->high),
4173 : PointPGetDatum(p)));
4174 20 : low = DatumGetPointP(DirectFunctionCall2(point_div,
4175 : PointPGetDatum(&box->low),
4176 : PointPGetDatum(p)));
4177 :
4178 20 : result = box_construct(high->x, low->x, high->y, low->y);
4179 :
4180 20 : PG_RETURN_BOX_P(result);
4181 : }
4182 :
4183 : /*
4184 : * Convert point to empty box
4185 : */
4186 : Datum
4187 6 : point_box(PG_FUNCTION_ARGS)
4188 : {
4189 6 : Point *pt = PG_GETARG_POINT_P(0);
4190 : BOX *box;
4191 :
4192 6 : box = (BOX *) palloc(sizeof(BOX));
4193 :
4194 6 : box->high.x = pt->x;
4195 6 : box->low.x = pt->x;
4196 6 : box->high.y = pt->y;
4197 6 : box->low.y = pt->y;
4198 :
4199 6 : PG_RETURN_BOX_P(box);
4200 : }
4201 :
4202 : /*
4203 : * Smallest bounding box that includes both of the given boxes
4204 : */
4205 : Datum
4206 16 : boxes_bound_box(PG_FUNCTION_ARGS)
4207 : {
4208 16 : BOX *box1 = PG_GETARG_BOX_P(0),
4209 16 : *box2 = PG_GETARG_BOX_P(1),
4210 : *container;
4211 :
4212 16 : container = (BOX *) palloc(sizeof(BOX));
4213 :
4214 16 : container->high.x = Max(box1->high.x, box2->high.x);
4215 16 : container->low.x = Min(box1->low.x, box2->low.x);
4216 16 : container->high.y = Max(box1->high.y, box2->high.y);
4217 16 : container->low.y = Min(box1->low.y, box2->low.y);
4218 :
4219 16 : PG_RETURN_BOX_P(container);
4220 : }
4221 :
4222 :
4223 : /***********************************************************************
4224 : **
4225 : ** Routines for 2D paths.
4226 : **
4227 : ***********************************************************************/
4228 :
4229 : /* path_add()
4230 : * Concatenate two paths (only if they are both open).
4231 : */
4232 : Datum
4233 0 : path_add(PG_FUNCTION_ARGS)
4234 : {
4235 0 : PATH *p1 = PG_GETARG_PATH_P(0);
4236 0 : PATH *p2 = PG_GETARG_PATH_P(1);
4237 : PATH *result;
4238 : int size,
4239 : base_size;
4240 : int i;
4241 :
4242 0 : if (p1->closed || p2->closed)
4243 0 : PG_RETURN_NULL();
4244 :
4245 0 : base_size = sizeof(p1->p[0]) * (p1->npts + p2->npts);
4246 0 : size = offsetof(PATH, p) + base_size;
4247 :
4248 : /* Check for integer overflow */
4249 0 : if (base_size / sizeof(p1->p[0]) != (p1->npts + p2->npts) ||
4250 : size <= base_size)
4251 0 : ereport(ERROR,
4252 : (errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED),
4253 : errmsg("too many points requested")));
4254 :
4255 0 : result = (PATH *) palloc(size);
4256 :
4257 0 : SET_VARSIZE(result, size);
4258 0 : result->npts = (p1->npts + p2->npts);
4259 0 : result->closed = p1->closed;
4260 : /* prevent instability in unused pad bytes */
4261 0 : result->dummy = 0;
4262 :
4263 0 : for (i = 0; i < p1->npts; i++)
4264 : {
4265 0 : result->p[i].x = p1->p[i].x;
4266 0 : result->p[i].y = p1->p[i].y;
4267 : }
4268 0 : for (i = 0; i < p2->npts; i++)
4269 : {
4270 0 : result->p[i + p1->npts].x = p2->p[i].x;
4271 0 : result->p[i + p1->npts].y = p2->p[i].y;
4272 : }
4273 :
4274 0 : PG_RETURN_PATH_P(result);
4275 : }
4276 :
4277 : /* path_add_pt()
4278 : * Translation operators.
4279 : */
4280 : Datum
4281 8 : path_add_pt(PG_FUNCTION_ARGS)
4282 : {
4283 8 : PATH *path = PG_GETARG_PATH_P_COPY(0);
4284 8 : Point *point = PG_GETARG_POINT_P(1);
4285 : int i;
4286 :
4287 26 : for (i = 0; i < path->npts; i++)
4288 : {
4289 18 : path->p[i].x += point->x;
4290 18 : path->p[i].y += point->y;
4291 : }
4292 :
4293 8 : PG_RETURN_PATH_P(path);
4294 : }
4295 :
4296 : Datum
4297 0 : path_sub_pt(PG_FUNCTION_ARGS)
4298 : {
4299 0 : PATH *path = PG_GETARG_PATH_P_COPY(0);
4300 0 : Point *point = PG_GETARG_POINT_P(1);
4301 : int i;
4302 :
4303 0 : for (i = 0; i < path->npts; i++)
4304 : {
4305 0 : path->p[i].x -= point->x;
4306 0 : path->p[i].y -= point->y;
4307 : }
4308 :
4309 0 : PG_RETURN_PATH_P(path);
4310 : }
4311 :
4312 : /* path_mul_pt()
4313 : * Rotation and scaling operators.
4314 : */
4315 : Datum
4316 8 : path_mul_pt(PG_FUNCTION_ARGS)
4317 : {
4318 8 : PATH *path = PG_GETARG_PATH_P_COPY(0);
4319 8 : Point *point = PG_GETARG_POINT_P(1);
4320 : Point *p;
4321 : int i;
4322 :
4323 26 : for (i = 0; i < path->npts; i++)
4324 : {
4325 18 : p = DatumGetPointP(DirectFunctionCall2(point_mul,
4326 : PointPGetDatum(&path->p[i]),
4327 : PointPGetDatum(point)));
4328 18 : path->p[i].x = p->x;
4329 18 : path->p[i].y = p->y;
4330 : }
4331 :
4332 8 : PG_RETURN_PATH_P(path);
4333 : }
4334 :
4335 : Datum
4336 0 : path_div_pt(PG_FUNCTION_ARGS)
4337 : {
4338 0 : PATH *path = PG_GETARG_PATH_P_COPY(0);
4339 0 : Point *point = PG_GETARG_POINT_P(1);
4340 : Point *p;
4341 : int i;
4342 :
4343 0 : for (i = 0; i < path->npts; i++)
4344 : {
4345 0 : p = DatumGetPointP(DirectFunctionCall2(point_div,
4346 : PointPGetDatum(&path->p[i]),
4347 : PointPGetDatum(point)));
4348 0 : path->p[i].x = p->x;
4349 0 : path->p[i].y = p->y;
4350 : }
4351 :
4352 0 : PG_RETURN_PATH_P(path);
4353 : }
4354 :
4355 :
4356 : Datum
4357 0 : path_center(PG_FUNCTION_ARGS)
4358 : {
4359 : #ifdef NOT_USED
4360 : PATH *path = PG_GETARG_PATH_P(0);
4361 : #endif
4362 :
4363 0 : ereport(ERROR,
4364 : (errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
4365 : errmsg("function \"path_center\" not implemented")));
4366 :
4367 : PG_RETURN_NULL();
4368 : }
4369 :
4370 : Datum
4371 8 : path_poly(PG_FUNCTION_ARGS)
4372 : {
4373 8 : PATH *path = PG_GETARG_PATH_P(0);
4374 : POLYGON *poly;
4375 : int size;
4376 : int i;
4377 :
4378 : /* This is not very consistent --- other similar cases return NULL ... */
4379 8 : if (!path->closed)
4380 0 : ereport(ERROR,
4381 : (errcode(ERRCODE_INVALID_PARAMETER_VALUE),
4382 : errmsg("open path cannot be converted to polygon")));
4383 :
4384 : /*
4385 : * Never overflows: the old size fit in MaxAllocSize, and the new size is
4386 : * just a small constant larger.
4387 : */
4388 8 : size = offsetof(POLYGON, p) + sizeof(poly->p[0]) * path->npts;
4389 8 : poly = (POLYGON *) palloc(size);
4390 :
4391 8 : SET_VARSIZE(poly, size);
4392 8 : poly->npts = path->npts;
4393 :
4394 26 : for (i = 0; i < path->npts; i++)
4395 : {
4396 18 : poly->p[i].x = path->p[i].x;
4397 18 : poly->p[i].y = path->p[i].y;
4398 : }
4399 :
4400 8 : make_bound_box(poly);
4401 :
4402 8 : PG_RETURN_POLYGON_P(poly);
4403 : }
4404 :
4405 :
4406 : /***********************************************************************
4407 : **
4408 : ** Routines for 2D polygons.
4409 : **
4410 : ***********************************************************************/
4411 :
4412 : Datum
4413 12 : poly_npoints(PG_FUNCTION_ARGS)
4414 : {
4415 12 : POLYGON *poly = PG_GETARG_POLYGON_P(0);
4416 :
4417 12 : PG_RETURN_INT32(poly->npts);
4418 : }
4419 :
4420 :
4421 : Datum
4422 4 : poly_center(PG_FUNCTION_ARGS)
4423 : {
4424 4 : POLYGON *poly = PG_GETARG_POLYGON_P(0);
4425 : Datum result;
4426 : CIRCLE *circle;
4427 :
4428 4 : circle = DatumGetCircleP(DirectFunctionCall1(poly_circle,
4429 : PolygonPGetDatum(poly)));
4430 4 : result = DirectFunctionCall1(circle_center,
4431 : CirclePGetDatum(circle));
4432 :
4433 4 : PG_RETURN_DATUM(result);
4434 : }
4435 :
4436 :
4437 : Datum
4438 0 : poly_box(PG_FUNCTION_ARGS)
4439 : {
4440 0 : POLYGON *poly = PG_GETARG_POLYGON_P(0);
4441 : BOX *box;
4442 :
4443 0 : if (poly->npts < 1)
4444 0 : PG_RETURN_NULL();
4445 :
4446 0 : box = box_copy(&poly->boundbox);
4447 :
4448 0 : PG_RETURN_BOX_P(box);
4449 : }
4450 :
4451 :
4452 : /* box_poly()
4453 : * Convert a box to a polygon.
4454 : */
4455 : Datum
4456 3104 : box_poly(PG_FUNCTION_ARGS)
4457 : {
4458 3104 : BOX *box = PG_GETARG_BOX_P(0);
4459 : POLYGON *poly;
4460 : int size;
4461 :
4462 : /* map four corners of the box to a polygon */
4463 3104 : size = offsetof(POLYGON, p) + sizeof(poly->p[0]) * 4;
4464 3104 : poly = (POLYGON *) palloc(size);
4465 :
4466 3104 : SET_VARSIZE(poly, size);
4467 3104 : poly->npts = 4;
4468 :
4469 3104 : poly->p[0].x = box->low.x;
4470 3104 : poly->p[0].y = box->low.y;
4471 3104 : poly->p[1].x = box->low.x;
4472 3104 : poly->p[1].y = box->high.y;
4473 3104 : poly->p[2].x = box->high.x;
4474 3104 : poly->p[2].y = box->high.y;
4475 3104 : poly->p[3].x = box->high.x;
4476 3104 : poly->p[3].y = box->low.y;
4477 :
4478 3104 : box_fill(&poly->boundbox, box->high.x, box->low.x,
4479 : box->high.y, box->low.y);
4480 :
4481 3104 : PG_RETURN_POLYGON_P(poly);
4482 : }
4483 :
4484 :
4485 : Datum
4486 4 : poly_path(PG_FUNCTION_ARGS)
4487 : {
4488 4 : POLYGON *poly = PG_GETARG_POLYGON_P(0);
4489 : PATH *path;
4490 : int size;
4491 : int i;
4492 :
4493 : /*
4494 : * Never overflows: the old size fit in MaxAllocSize, and the new size is
4495 : * smaller by a small constant.
4496 : */
4497 4 : size = offsetof(PATH, p) + sizeof(path->p[0]) * poly->npts;
4498 4 : path = (PATH *) palloc(size);
4499 :
4500 4 : SET_VARSIZE(path, size);
4501 4 : path->npts = poly->npts;
4502 4 : path->closed = TRUE;
4503 : /* prevent instability in unused pad bytes */
4504 4 : path->dummy = 0;
4505 :
4506 13 : for (i = 0; i < poly->npts; i++)
4507 : {
4508 9 : path->p[i].x = poly->p[i].x;
4509 9 : path->p[i].y = poly->p[i].y;
4510 : }
4511 :
4512 4 : PG_RETURN_PATH_P(path);
4513 : }
4514 :
4515 :
4516 : /***********************************************************************
4517 : **
4518 : ** Routines for circles.
4519 : **
4520 : ***********************************************************************/
4521 :
4522 : /*----------------------------------------------------------
4523 : * Formatting and conversion routines.
4524 : *---------------------------------------------------------*/
4525 :
4526 : /* circle_in - convert a string to internal form.
4527 : *
4528 : * External format: (center and radius of circle)
4529 : * "((f8,f8)<f8>)"
4530 : * also supports quick entry style "(f8,f8,f8)"
4531 : */
4532 : Datum
4533 52 : circle_in(PG_FUNCTION_ARGS)
4534 : {
4535 52 : char *str = PG_GETARG_CSTRING(0);
4536 52 : CIRCLE *circle = (CIRCLE *) palloc(sizeof(CIRCLE));
4537 : char *s,
4538 : *cp;
4539 52 : int depth = 0;
4540 :
4541 52 : s = str;
4542 104 : while (isspace((unsigned char) *s))
4543 0 : s++;
4544 52 : if ((*s == LDELIM_C) || (*s == LDELIM))
4545 : {
4546 50 : depth++;
4547 50 : cp = (s + 1);
4548 100 : while (isspace((unsigned char) *cp))
4549 0 : cp++;
4550 50 : if (*cp == LDELIM)
4551 49 : s = cp;
4552 : }
4553 :
4554 52 : pair_decode(s, &circle->center.x, &circle->center.y, &s, "circle", str);
4555 :
4556 50 : if (*s == DELIM)
4557 50 : s++;
4558 :
4559 50 : circle->radius = single_decode(s, &s, "circle", str);
4560 50 : if (circle->radius < 0)
4561 1 : ereport(ERROR,
4562 : (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
4563 : errmsg("invalid input syntax for type %s: \"%s\"",
4564 : "circle", str)));
4565 :
4566 146 : while (depth > 0)
4567 : {
4568 96 : if ((*s == RDELIM)
4569 47 : || ((*s == RDELIM_C) && (depth == 1)))
4570 : {
4571 48 : depth--;
4572 48 : s++;
4573 96 : while (isspace((unsigned char) *s))
4574 0 : s++;
4575 : }
4576 : else
4577 0 : ereport(ERROR,
4578 : (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
4579 : errmsg("invalid input syntax for type %s: \"%s\"",
4580 : "circle", str)));
4581 : }
4582 :
4583 49 : if (*s != '\0')
4584 0 : ereport(ERROR,
4585 : (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
4586 : errmsg("invalid input syntax for type %s: \"%s\"",
4587 : "circle", str)));
4588 :
4589 49 : PG_RETURN_CIRCLE_P(circle);
4590 : }
4591 :
4592 : /* circle_out - convert a circle to external form.
4593 : */
4594 : Datum
4595 74 : circle_out(PG_FUNCTION_ARGS)
4596 : {
4597 74 : CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
4598 : StringInfoData str;
4599 :
4600 74 : initStringInfo(&str);
4601 :
4602 74 : appendStringInfoChar(&str, LDELIM_C);
4603 74 : appendStringInfoChar(&str, LDELIM);
4604 74 : pair_encode(circle->center.x, circle->center.y, &str);
4605 74 : appendStringInfoChar(&str, RDELIM);
4606 74 : appendStringInfoChar(&str, DELIM);
4607 74 : single_encode(circle->radius, &str);
4608 74 : appendStringInfoChar(&str, RDELIM_C);
4609 :
4610 74 : PG_RETURN_CSTRING(str.data);
4611 : }
4612 :
4613 : /*
4614 : * circle_recv - converts external binary format to circle
4615 : */
4616 : Datum
4617 0 : circle_recv(PG_FUNCTION_ARGS)
4618 : {
4619 0 : StringInfo buf = (StringInfo) PG_GETARG_POINTER(0);
4620 : CIRCLE *circle;
4621 :
4622 0 : circle = (CIRCLE *) palloc(sizeof(CIRCLE));
4623 :
4624 0 : circle->center.x = pq_getmsgfloat8(buf);
4625 0 : circle->center.y = pq_getmsgfloat8(buf);
4626 0 : circle->radius = pq_getmsgfloat8(buf);
4627 :
4628 0 : if (circle->radius < 0)
4629 0 : ereport(ERROR,
4630 : (errcode(ERRCODE_INVALID_BINARY_REPRESENTATION),
4631 : errmsg("invalid radius in external \"circle\" value")));
4632 :
4633 0 : PG_RETURN_CIRCLE_P(circle);
4634 : }
4635 :
4636 : /*
4637 : * circle_send - converts circle to binary format
4638 : */
4639 : Datum
4640 0 : circle_send(PG_FUNCTION_ARGS)
4641 : {
4642 0 : CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
4643 : StringInfoData buf;
4644 :
4645 0 : pq_begintypsend(&buf);
4646 0 : pq_sendfloat8(&buf, circle->center.x);
4647 0 : pq_sendfloat8(&buf, circle->center.y);
4648 0 : pq_sendfloat8(&buf, circle->radius);
4649 0 : PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
4650 : }
4651 :
4652 :
4653 : /*----------------------------------------------------------
4654 : * Relational operators for CIRCLEs.
4655 : * <, >, <=, >=, and == are based on circle area.
4656 : *---------------------------------------------------------*/
4657 :
4658 : /* circles identical?
4659 : */
4660 : Datum
4661 0 : circle_same(PG_FUNCTION_ARGS)
4662 : {
4663 0 : CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4664 0 : CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4665 :
4666 0 : PG_RETURN_BOOL(FPeq(circle1->radius, circle2->radius) &&
4667 : FPeq(circle1->center.x, circle2->center.x) &&
4668 : FPeq(circle1->center.y, circle2->center.y));
4669 : }
4670 :
4671 : /* circle_overlap - does circle1 overlap circle2?
4672 : */
4673 : Datum
4674 3134 : circle_overlap(PG_FUNCTION_ARGS)
4675 : {
4676 3134 : CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4677 3134 : CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4678 :
4679 3134 : PG_RETURN_BOOL(FPle(point_dt(&circle1->center, &circle2->center),
4680 : circle1->radius + circle2->radius));
4681 : }
4682 :
4683 : /* circle_overleft - is the right edge of circle1 at or left of
4684 : * the right edge of circle2?
4685 : */
4686 : Datum
4687 0 : circle_overleft(PG_FUNCTION_ARGS)
4688 : {
4689 0 : CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4690 0 : CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4691 :
4692 0 : PG_RETURN_BOOL(FPle((circle1->center.x + circle1->radius),
4693 : (circle2->center.x + circle2->radius)));
4694 : }
4695 :
4696 : /* circle_left - is circle1 strictly left of circle2?
4697 : */
4698 : Datum
4699 0 : circle_left(PG_FUNCTION_ARGS)
4700 : {
4701 0 : CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4702 0 : CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4703 :
4704 0 : PG_RETURN_BOOL(FPlt((circle1->center.x + circle1->radius),
4705 : (circle2->center.x - circle2->radius)));
4706 : }
4707 :
4708 : /* circle_right - is circle1 strictly right of circle2?
4709 : */
4710 : Datum
4711 0 : circle_right(PG_FUNCTION_ARGS)
4712 : {
4713 0 : CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4714 0 : CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4715 :
4716 0 : PG_RETURN_BOOL(FPgt((circle1->center.x - circle1->radius),
4717 : (circle2->center.x + circle2->radius)));
4718 : }
4719 :
4720 : /* circle_overright - is the left edge of circle1 at or right of
4721 : * the left edge of circle2?
4722 : */
4723 : Datum
4724 0 : circle_overright(PG_FUNCTION_ARGS)
4725 : {
4726 0 : CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4727 0 : CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4728 :
4729 0 : PG_RETURN_BOOL(FPge((circle1->center.x - circle1->radius),
4730 : (circle2->center.x - circle2->radius)));
4731 : }
4732 :
4733 : /* circle_contained - is circle1 contained by circle2?
4734 : */
4735 : Datum
4736 0 : circle_contained(PG_FUNCTION_ARGS)
4737 : {
4738 0 : CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4739 0 : CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4740 :
4741 0 : PG_RETURN_BOOL(FPle((point_dt(&circle1->center, &circle2->center) + circle1->radius), circle2->radius));
4742 : }
4743 :
4744 : /* circle_contain - does circle1 contain circle2?
4745 : */
4746 : Datum
4747 0 : circle_contain(PG_FUNCTION_ARGS)
4748 : {
4749 0 : CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4750 0 : CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4751 :
4752 0 : PG_RETURN_BOOL(FPle((point_dt(&circle1->center, &circle2->center) + circle2->radius), circle1->radius));
4753 : }
4754 :
4755 :
4756 : /* circle_below - is circle1 strictly below circle2?
4757 : */
4758 : Datum
4759 0 : circle_below(PG_FUNCTION_ARGS)
4760 : {
4761 0 : CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4762 0 : CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4763 :
4764 0 : PG_RETURN_BOOL(FPlt((circle1->center.y + circle1->radius),
4765 : (circle2->center.y - circle2->radius)));
4766 : }
4767 :
4768 : /* circle_above - is circle1 strictly above circle2?
4769 : */
4770 : Datum
4771 0 : circle_above(PG_FUNCTION_ARGS)
4772 : {
4773 0 : CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4774 0 : CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4775 :
4776 0 : PG_RETURN_BOOL(FPgt((circle1->center.y - circle1->radius),
4777 : (circle2->center.y + circle2->radius)));
4778 : }
4779 :
4780 : /* circle_overbelow - is the upper edge of circle1 at or below
4781 : * the upper edge of circle2?
4782 : */
4783 : Datum
4784 0 : circle_overbelow(PG_FUNCTION_ARGS)
4785 : {
4786 0 : CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4787 0 : CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4788 :
4789 0 : PG_RETURN_BOOL(FPle((circle1->center.y + circle1->radius),
4790 : (circle2->center.y + circle2->radius)));
4791 : }
4792 :
4793 : /* circle_overabove - is the lower edge of circle1 at or above
4794 : * the lower edge of circle2?
4795 : */
4796 : Datum
4797 0 : circle_overabove(PG_FUNCTION_ARGS)
4798 : {
4799 0 : CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4800 0 : CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4801 :
4802 0 : PG_RETURN_BOOL(FPge((circle1->center.y - circle1->radius),
4803 : (circle2->center.y - circle2->radius)));
4804 : }
4805 :
4806 :
4807 : /* circle_relop - is area(circle1) relop area(circle2), within
4808 : * our accuracy constraint?
4809 : */
4810 : Datum
4811 0 : circle_eq(PG_FUNCTION_ARGS)
4812 : {
4813 0 : CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4814 0 : CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4815 :
4816 0 : PG_RETURN_BOOL(FPeq(circle_ar(circle1), circle_ar(circle2)));
4817 : }
4818 :
4819 : Datum
4820 0 : circle_ne(PG_FUNCTION_ARGS)
4821 : {
4822 0 : CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4823 0 : CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4824 :
4825 0 : PG_RETURN_BOOL(FPne(circle_ar(circle1), circle_ar(circle2)));
4826 : }
4827 :
4828 : Datum
4829 36 : circle_lt(PG_FUNCTION_ARGS)
4830 : {
4831 36 : CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4832 36 : CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4833 :
4834 36 : PG_RETURN_BOOL(FPlt(circle_ar(circle1), circle_ar(circle2)));
4835 : }
4836 :
4837 : Datum
4838 0 : circle_gt(PG_FUNCTION_ARGS)
4839 : {
4840 0 : CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4841 0 : CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4842 :
4843 0 : PG_RETURN_BOOL(FPgt(circle_ar(circle1), circle_ar(circle2)));
4844 : }
4845 :
4846 : Datum
4847 0 : circle_le(PG_FUNCTION_ARGS)
4848 : {
4849 0 : CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4850 0 : CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4851 :
4852 0 : PG_RETURN_BOOL(FPle(circle_ar(circle1), circle_ar(circle2)));
4853 : }
4854 :
4855 : Datum
4856 0 : circle_ge(PG_FUNCTION_ARGS)
4857 : {
4858 0 : CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4859 0 : CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4860 :
4861 0 : PG_RETURN_BOOL(FPge(circle_ar(circle1), circle_ar(circle2)));
4862 : }
4863 :
4864 :
4865 : /*----------------------------------------------------------
4866 : * "Arithmetic" operators on circles.
4867 : *---------------------------------------------------------*/
4868 :
4869 : static CIRCLE *
4870 0 : circle_copy(CIRCLE *circle)
4871 : {
4872 : CIRCLE *result;
4873 :
4874 0 : if (!PointerIsValid(circle))
4875 0 : return NULL;
4876 :
4877 0 : result = (CIRCLE *) palloc(sizeof(CIRCLE));
4878 0 : memcpy((char *) result, (char *) circle, sizeof(CIRCLE));
4879 0 : return result;
4880 : }
4881 :
4882 :
4883 : /* circle_add_pt()
4884 : * Translation operator.
4885 : */
4886 : Datum
4887 0 : circle_add_pt(PG_FUNCTION_ARGS)
4888 : {
4889 0 : CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
4890 0 : Point *point = PG_GETARG_POINT_P(1);
4891 : CIRCLE *result;
4892 :
4893 0 : result = circle_copy(circle);
4894 :
4895 0 : result->center.x += point->x;
4896 0 : result->center.y += point->y;
4897 :
4898 0 : PG_RETURN_CIRCLE_P(result);
4899 : }
4900 :
4901 : Datum
4902 0 : circle_sub_pt(PG_FUNCTION_ARGS)
4903 : {
4904 0 : CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
4905 0 : Point *point = PG_GETARG_POINT_P(1);
4906 : CIRCLE *result;
4907 :
4908 0 : result = circle_copy(circle);
4909 :
4910 0 : result->center.x -= point->x;
4911 0 : result->center.y -= point->y;
4912 :
4913 0 : PG_RETURN_CIRCLE_P(result);
4914 : }
4915 :
4916 :
4917 : /* circle_mul_pt()
4918 : * Rotation and scaling operators.
4919 : */
4920 : Datum
4921 0 : circle_mul_pt(PG_FUNCTION_ARGS)
4922 : {
4923 0 : CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
4924 0 : Point *point = PG_GETARG_POINT_P(1);
4925 : CIRCLE *result;
4926 : Point *p;
4927 :
4928 0 : result = circle_copy(circle);
4929 :
4930 0 : p = DatumGetPointP(DirectFunctionCall2(point_mul,
4931 : PointPGetDatum(&circle->center),
4932 : PointPGetDatum(point)));
4933 0 : result->center.x = p->x;
4934 0 : result->center.y = p->y;
4935 0 : result->radius *= HYPOT(point->x, point->y);
4936 :
4937 0 : PG_RETURN_CIRCLE_P(result);
4938 : }
4939 :
4940 : Datum
4941 0 : circle_div_pt(PG_FUNCTION_ARGS)
4942 : {
4943 0 : CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
4944 0 : Point *point = PG_GETARG_POINT_P(1);
4945 : CIRCLE *result;
4946 : Point *p;
4947 :
4948 0 : result = circle_copy(circle);
4949 :
4950 0 : p = DatumGetPointP(DirectFunctionCall2(point_div,
4951 : PointPGetDatum(&circle->center),
4952 : PointPGetDatum(point)));
4953 0 : result->center.x = p->x;
4954 0 : result->center.y = p->y;
4955 0 : result->radius /= HYPOT(point->x, point->y);
4956 :
4957 0 : PG_RETURN_CIRCLE_P(result);
4958 : }
4959 :
4960 :
4961 : /* circle_area - returns the area of the circle.
4962 : */
4963 : Datum
4964 39 : circle_area(PG_FUNCTION_ARGS)
4965 : {
4966 39 : CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
4967 :
4968 39 : PG_RETURN_FLOAT8(circle_ar(circle));
4969 : }
4970 :
4971 :
4972 : /* circle_diameter - returns the diameter of the circle.
4973 : */
4974 : Datum
4975 12 : circle_diameter(PG_FUNCTION_ARGS)
4976 : {
4977 12 : CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
4978 :
4979 12 : PG_RETURN_FLOAT8(2 * circle->radius);
4980 : }
4981 :
4982 :
4983 : /* circle_radius - returns the radius of the circle.
4984 : */
4985 : Datum
4986 3122 : circle_radius(PG_FUNCTION_ARGS)
4987 : {
4988 3122 : CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
4989 :
4990 3122 : PG_RETURN_FLOAT8(circle->radius);
4991 : }
4992 :
4993 :
4994 : /* circle_distance - returns the distance between
4995 : * two circles.
4996 : */
4997 : Datum
4998 19 : circle_distance(PG_FUNCTION_ARGS)
4999 : {
5000 19 : CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
5001 19 : CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
5002 : float8 result;
5003 :
5004 57 : result = point_dt(&circle1->center, &circle2->center)
5005 38 : - (circle1->radius + circle2->radius);
5006 19 : if (result < 0)
5007 9 : result = 0;
5008 19 : PG_RETURN_FLOAT8(result);
5009 : }
5010 :
5011 :
5012 : Datum
5013 3 : circle_contain_pt(PG_FUNCTION_ARGS)
5014 : {
5015 3 : CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
5016 3 : Point *point = PG_GETARG_POINT_P(1);
5017 : double d;
5018 :
5019 3 : d = point_dt(&circle->center, point);
5020 3 : PG_RETURN_BOOL(d <= circle->radius);
5021 : }
5022 :
5023 :
5024 : Datum
5025 6 : pt_contained_circle(PG_FUNCTION_ARGS)
5026 : {
5027 6 : Point *point = PG_GETARG_POINT_P(0);
5028 6 : CIRCLE *circle = PG_GETARG_CIRCLE_P(1);
5029 : double d;
5030 :
5031 6 : d = point_dt(&circle->center, point);
5032 6 : PG_RETURN_BOOL(d <= circle->radius);
5033 : }
5034 :
5035 :
5036 : /* dist_pc - returns the distance between
5037 : * a point and a circle.
5038 : */
5039 : Datum
5040 57 : dist_pc(PG_FUNCTION_ARGS)
5041 : {
5042 57 : Point *point = PG_GETARG_POINT_P(0);
5043 57 : CIRCLE *circle = PG_GETARG_CIRCLE_P(1);
5044 : float8 result;
5045 :
5046 57 : result = point_dt(point, &circle->center) - circle->radius;
5047 57 : if (result < 0)
5048 15 : result = 0;
5049 57 : PG_RETURN_FLOAT8(result);
5050 : }
5051 :
5052 : /*
5053 : * Distance from a circle to a point
5054 : */
5055 : Datum
5056 3121 : dist_cpoint(PG_FUNCTION_ARGS)
5057 : {
5058 3121 : CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
5059 3121 : Point *point = PG_GETARG_POINT_P(1);
5060 : float8 result;
5061 :
5062 3121 : result = point_dt(point, &circle->center) - circle->radius;
5063 3121 : if (result < 0)
5064 0 : result = 0;
5065 3121 : PG_RETURN_FLOAT8(result);
5066 : }
5067 :
5068 : /* circle_center - returns the center point of the circle.
5069 : */
5070 : Datum
5071 3149 : circle_center(PG_FUNCTION_ARGS)
5072 : {
5073 3149 : CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
5074 : Point *result;
5075 :
5076 3149 : result = (Point *) palloc(sizeof(Point));
5077 3149 : result->x = circle->center.x;
5078 3149 : result->y = circle->center.y;
5079 :
5080 3149 : PG_RETURN_POINT_P(result);
5081 : }
5082 :
5083 :
5084 : /* circle_ar - returns the area of the circle.
5085 : */
5086 : static double
5087 111 : circle_ar(CIRCLE *circle)
5088 : {
5089 111 : return M_PI * (circle->radius * circle->radius);
5090 : }
5091 :
5092 :
5093 : /*----------------------------------------------------------
5094 : * Conversion operators.
5095 : *---------------------------------------------------------*/
5096 :
5097 : Datum
5098 10010 : cr_circle(PG_FUNCTION_ARGS)
5099 : {
5100 10010 : Point *center = PG_GETARG_POINT_P(0);
5101 10010 : float8 radius = PG_GETARG_FLOAT8(1);
5102 : CIRCLE *result;
5103 :
5104 10010 : result = (CIRCLE *) palloc(sizeof(CIRCLE));
5105 :
5106 10010 : result->center.x = center->x;
5107 10010 : result->center.y = center->y;
5108 10010 : result->radius = radius;
5109 :
5110 10010 : PG_RETURN_CIRCLE_P(result);
5111 : }
5112 :
5113 : Datum
5114 6 : circle_box(PG_FUNCTION_ARGS)
5115 : {
5116 6 : CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
5117 : BOX *box;
5118 : double delta;
5119 :
5120 6 : box = (BOX *) palloc(sizeof(BOX));
5121 :
5122 6 : delta = circle->radius / sqrt(2.0);
5123 :
5124 6 : box->high.x = circle->center.x + delta;
5125 6 : box->low.x = circle->center.x - delta;
5126 6 : box->high.y = circle->center.y + delta;
5127 6 : box->low.y = circle->center.y - delta;
5128 :
5129 6 : PG_RETURN_BOX_P(box);
5130 : }
5131 :
5132 : /* box_circle()
5133 : * Convert a box to a circle.
5134 : */
5135 : Datum
5136 3104 : box_circle(PG_FUNCTION_ARGS)
5137 : {
5138 3104 : BOX *box = PG_GETARG_BOX_P(0);
5139 : CIRCLE *circle;
5140 :
5141 3104 : circle = (CIRCLE *) palloc(sizeof(CIRCLE));
5142 :
5143 3104 : circle->center.x = (box->high.x + box->low.x) / 2;
5144 3104 : circle->center.y = (box->high.y + box->low.y) / 2;
5145 :
5146 3104 : circle->radius = point_dt(&circle->center, &box->high);
5147 :
5148 3104 : PG_RETURN_CIRCLE_P(circle);
5149 : }
5150 :
5151 :
5152 : Datum
5153 12 : circle_poly(PG_FUNCTION_ARGS)
5154 : {
5155 12 : int32 npts = PG_GETARG_INT32(0);
5156 12 : CIRCLE *circle = PG_GETARG_CIRCLE_P(1);
5157 : POLYGON *poly;
5158 : int base_size,
5159 : size;
5160 : int i;
5161 : double angle;
5162 : double anglestep;
5163 :
5164 12 : if (FPzero(circle->radius))
5165 0 : ereport(ERROR,
5166 : (errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
5167 : errmsg("cannot convert circle with radius zero to polygon")));
5168 :
5169 12 : if (npts < 2)
5170 0 : ereport(ERROR,
5171 : (errcode(ERRCODE_INVALID_PARAMETER_VALUE),
5172 : errmsg("must request at least 2 points")));
5173 :
5174 12 : base_size = sizeof(poly->p[0]) * npts;
5175 12 : size = offsetof(POLYGON, p) + base_size;
5176 :
5177 : /* Check for integer overflow */
5178 12 : if (base_size / npts != sizeof(poly->p[0]) || size <= base_size)
5179 0 : ereport(ERROR,
5180 : (errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED),
5181 : errmsg("too many points requested")));
5182 :
5183 12 : poly = (POLYGON *) palloc0(size); /* zero any holes */
5184 12 : SET_VARSIZE(poly, size);
5185 12 : poly->npts = npts;
5186 :
5187 12 : anglestep = (2.0 * M_PI) / npts;
5188 :
5189 132 : for (i = 0; i < npts; i++)
5190 : {
5191 120 : angle = i * anglestep;
5192 120 : poly->p[i].x = circle->center.x - (circle->radius * cos(angle));
5193 120 : poly->p[i].y = circle->center.y + (circle->radius * sin(angle));
5194 : }
5195 :
5196 12 : make_bound_box(poly);
5197 :
5198 12 : PG_RETURN_POLYGON_P(poly);
5199 : }
5200 :
5201 : /* poly_circle - convert polygon to circle
5202 : *
5203 : * XXX This algorithm should use weighted means of line segments
5204 : * rather than straight average values of points - tgl 97/01/21.
5205 : */
5206 : Datum
5207 6 : poly_circle(PG_FUNCTION_ARGS)
5208 : {
5209 6 : POLYGON *poly = PG_GETARG_POLYGON_P(0);
5210 : CIRCLE *circle;
5211 : int i;
5212 :
5213 6 : if (poly->npts < 1)
5214 0 : ereport(ERROR,
5215 : (errcode(ERRCODE_INVALID_PARAMETER_VALUE),
5216 : errmsg("cannot convert empty polygon to circle")));
5217 :
5218 6 : circle = (CIRCLE *) palloc(sizeof(CIRCLE));
5219 :
5220 6 : circle->center.x = 0;
5221 6 : circle->center.y = 0;
5222 6 : circle->radius = 0;
5223 :
5224 24 : for (i = 0; i < poly->npts; i++)
5225 : {
5226 18 : circle->center.x += poly->p[i].x;
5227 18 : circle->center.y += poly->p[i].y;
5228 : }
5229 6 : circle->center.x /= poly->npts;
5230 6 : circle->center.y /= poly->npts;
5231 :
5232 24 : for (i = 0; i < poly->npts; i++)
5233 18 : circle->radius += point_dt(&poly->p[i], &circle->center);
5234 6 : circle->radius /= poly->npts;
5235 :
5236 6 : PG_RETURN_CIRCLE_P(circle);
5237 : }
5238 :
5239 :
5240 : /***********************************************************************
5241 : **
5242 : ** Private routines for multiple types.
5243 : **
5244 : ***********************************************************************/
5245 :
5246 : /*
5247 : * Test to see if the point is inside the polygon, returns 1/0, or 2 if
5248 : * the point is on the polygon.
5249 : * Code adapted but not copied from integer-based routines in WN: A
5250 : * Server for the HTTP
5251 : * version 1.15.1, file wn/image.c
5252 : * http://hopf.math.northwestern.edu/index.html
5253 : * Description of algorithm: http://www.linuxjournal.com/article/2197
5254 : * http://www.linuxjournal.com/article/2029
5255 : */
5256 :
5257 : #define POINT_ON_POLYGON INT_MAX
5258 :
5259 : static int
5260 3216 : point_inside(Point *p, int npts, Point *plist)
5261 : {
5262 : double x0,
5263 : y0;
5264 : double prev_x,
5265 : prev_y;
5266 3216 : int i = 0;
5267 : double x,
5268 : y;
5269 : int cross,
5270 3216 : total_cross = 0;
5271 :
5272 3216 : if (npts <= 0)
5273 0 : return 0;
5274 :
5275 : /* compute first polygon point relative to single point */
5276 3216 : x0 = plist[0].x - p->x;
5277 3216 : y0 = plist[0].y - p->y;
5278 :
5279 3216 : prev_x = x0;
5280 3216 : prev_y = y0;
5281 : /* loop over polygon points and aggregate total_cross */
5282 12777 : for (i = 1; i < npts; i++)
5283 : {
5284 : /* compute next polygon point relative to single point */
5285 9566 : x = plist[i].x - p->x;
5286 9566 : y = plist[i].y - p->y;
5287 :
5288 : /* compute previous to current point crossing */
5289 9566 : if ((cross = lseg_crossing(x, y, prev_x, prev_y)) == POINT_ON_POLYGON)
5290 5 : return 2;
5291 9561 : total_cross += cross;
5292 :
5293 9561 : prev_x = x;
5294 9561 : prev_y = y;
5295 : }
5296 :
5297 : /* now do the first point */
5298 3211 : if ((cross = lseg_crossing(x0, y0, prev_x, prev_y)) == POINT_ON_POLYGON)
5299 3 : return 2;
5300 3208 : total_cross += cross;
5301 :
5302 3208 : if (total_cross != 0)
5303 23 : return 1;
5304 3185 : return 0;
5305 : }
5306 :
5307 :
5308 : /* lseg_crossing()
5309 : * Returns +/-2 if line segment crosses the positive X-axis in a +/- direction.
5310 : * Returns +/-1 if one point is on the positive X-axis.
5311 : * Returns 0 if both points are on the positive X-axis, or there is no crossing.
5312 : * Returns POINT_ON_POLYGON if the segment contains (0,0).
5313 : * Wow, that is one confusing API, but it is used above, and when summed,
5314 : * can tell is if a point is in a polygon.
5315 : */
5316 :
5317 : static int
5318 12777 : lseg_crossing(double x, double y, double prev_x, double prev_y)
5319 : {
5320 : double z;
5321 : int y_sign;
5322 :
5323 12777 : if (FPzero(y))
5324 : { /* y == 0, on X axis */
5325 35 : if (FPzero(x)) /* (x,y) is (0,0)? */
5326 7 : return POINT_ON_POLYGON;
5327 28 : else if (FPgt(x, 0))
5328 : { /* x > 0 */
5329 27 : if (FPzero(prev_y)) /* y and prev_y are zero */
5330 : /* prev_x > 0? */
5331 7 : return FPgt(prev_x, 0) ? 0 : POINT_ON_POLYGON;
5332 20 : return FPlt(prev_y, 0) ? 1 : -1;
5333 : }
5334 : else
5335 : { /* x < 0, x not on positive X axis */
5336 1 : if (FPzero(prev_y))
5337 : /* prev_x < 0? */
5338 0 : return FPlt(prev_x, 0) ? 0 : POINT_ON_POLYGON;
5339 1 : return 0;
5340 : }
5341 : }
5342 : else
5343 : { /* y != 0 */
5344 : /* compute y crossing direction from previous point */
5345 12742 : y_sign = FPgt(y, 0) ? 1 : -1;
5346 :
5347 12742 : if (FPzero(prev_y))
5348 : /* previous point was on X axis, so new point is either off or on */
5349 24 : return FPlt(prev_x, 0) ? 0 : y_sign;
5350 12718 : else if (FPgt(y_sign * prev_y, 0))
5351 : /* both above or below X axis */
5352 12657 : return 0; /* same sign */
5353 : else
5354 : { /* y and prev_y cross X-axis */
5355 61 : if (FPge(x, 0) && FPgt(prev_x, 0))
5356 : /* both non-negative so cross positive X-axis */
5357 17 : return 2 * y_sign;
5358 44 : if (FPlt(x, 0) && FPle(prev_x, 0))
5359 : /* both non-positive so do not cross positive X-axis */
5360 22 : return 0;
5361 :
5362 : /* x and y cross axises, see URL above point_inside() */
5363 22 : z = (x - prev_x) * y - (y - prev_y) * x;
5364 22 : if (FPzero(z))
5365 1 : return POINT_ON_POLYGON;
5366 21 : return FPgt((y_sign * z), 0) ? 0 : 2 * y_sign;
5367 : }
5368 : }
5369 : }
5370 :
5371 :
5372 : static bool
5373 3 : plist_same(int npts, Point *p1, Point *p2)
5374 : {
5375 : int i,
5376 : ii,
5377 : j;
5378 :
5379 : /* find match for first point */
5380 9 : for (i = 0; i < npts; i++)
5381 : {
5382 7 : if ((FPeq(p2[i].x, p1[0].x))
5383 1 : && (FPeq(p2[i].y, p1[0].y)))
5384 : {
5385 :
5386 : /* match found? then look forward through remaining points */
5387 3 : for (ii = 1, j = i + 1; ii < npts; ii++, j++)
5388 : {
5389 2 : if (j >= npts)
5390 0 : j = 0;
5391 2 : if ((!FPeq(p2[j].x, p1[ii].x))
5392 2 : || (!FPeq(p2[j].y, p1[ii].y)))
5393 : {
5394 : #ifdef GEODEBUG
5395 : printf("plist_same- %d failed forward match with %d\n", j, ii);
5396 : #endif
5397 : break;
5398 : }
5399 : }
5400 : #ifdef GEODEBUG
5401 : printf("plist_same- ii = %d/%d after forward match\n", ii, npts);
5402 : #endif
5403 1 : if (ii == npts)
5404 1 : return TRUE;
5405 :
5406 : /* match not found forwards? then look backwards */
5407 0 : for (ii = 1, j = i - 1; ii < npts; ii++, j--)
5408 : {
5409 0 : if (j < 0)
5410 0 : j = (npts - 1);
5411 0 : if ((!FPeq(p2[j].x, p1[ii].x))
5412 0 : || (!FPeq(p2[j].y, p1[ii].y)))
5413 : {
5414 : #ifdef GEODEBUG
5415 : printf("plist_same- %d failed reverse match with %d\n", j, ii);
5416 : #endif
5417 : break;
5418 : }
5419 : }
5420 : #ifdef GEODEBUG
5421 : printf("plist_same- ii = %d/%d after reverse match\n", ii, npts);
5422 : #endif
5423 0 : if (ii == npts)
5424 0 : return TRUE;
5425 : }
5426 : }
5427 :
5428 2 : return FALSE;
5429 : }
5430 :
5431 :
5432 : /*-------------------------------------------------------------------------
5433 : * Determine the hypotenuse.
5434 : *
5435 : * If required, x and y are swapped to make x the larger number. The
5436 : * traditional formula of x^2+y^2 is rearranged to factor x outside the
5437 : * sqrt. This allows computation of the hypotenuse for significantly
5438 : * larger values, and with a higher precision than when using the naive
5439 : * formula. In particular, this cannot overflow unless the final result
5440 : * would be out-of-range.
5441 : *
5442 : * sqrt( x^2 + y^2 ) = sqrt( x^2( 1 + y^2/x^2) )
5443 : * = x * sqrt( 1 + y^2/x^2 )
5444 : * = x * sqrt( 1 + y/x * y/x )
5445 : *
5446 : * It is expected that this routine will eventually be replaced with the
5447 : * C99 hypot() function.
5448 : *
5449 : * This implementation conforms to IEEE Std 1003.1 and GLIBC, in that the
5450 : * case of hypot(inf,nan) results in INF, and not NAN.
5451 : *-----------------------------------------------------------------------
5452 : */
5453 : double
5454 218901 : pg_hypot(double x, double y)
5455 : {
5456 : double yx;
5457 :
5458 : /* Handle INF and NaN properly */
5459 218901 : if (isinf(x) || isinf(y))
5460 0 : return get_float8_infinity();
5461 :
5462 218901 : if (isnan(x) || isnan(y))
5463 0 : return get_float8_nan();
5464 :
5465 : /* Else, drop any minus signs */
5466 218901 : x = fabs(x);
5467 218901 : y = fabs(y);
5468 :
5469 : /* Swap x and y if needed to make x the larger one */
5470 218901 : if (x < y)
5471 : {
5472 146520 : double temp = x;
5473 :
5474 146520 : x = y;
5475 146520 : y = temp;
5476 : }
5477 :
5478 : /*
5479 : * If y is zero, the hypotenuse is x. This test saves a few cycles in
5480 : * such cases, but more importantly it also protects against
5481 : * divide-by-zero errors, since now x >= y.
5482 : */
5483 218901 : if (y == 0.0)
5484 101346 : return x;
5485 :
5486 : /* Determine the hypotenuse */
5487 117555 : yx = y / x;
5488 117555 : return x * sqrt(1.0 + (yx * yx));
5489 : }
|
