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

Generated by: LCOV version 1.13