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