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