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 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 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 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 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 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 * 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 2118 lseg_invsl(LSEG *lseg) 2119 { 2120 return point_invsl(&lseg->p[0], &lseg->p[1]); 2121 } 2122 2123 2124 Datum 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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