LCOV - code coverage report
Current view: top level - src/backend/utils/adt - geo_ops.c (source / functions) Hit Total Coverage
Test: PostgreSQL 15beta1 Lines: 1841 1990 92.5 %
Date: 2022-05-18 03:10:05 Functions: 254 267 95.1 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.14