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

Generated by: LCOV version 1.13