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