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