granicus.if.org Git - postgresql/blob - src/backend/utils/adt/geo_ops.c

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