LCOV - code coverage report
Current view: top level - src/backend/utils/adt - geo_ops.c (source / functions) Hit Total Coverage
Test: PostgreSQL 13devel Lines: 1813 1957 92.6 %
Date: 2019-09-22 08:06:49 Functions: 258 271 95.2 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /*-------------------------------------------------------------------------
       2             :  *
       3             :  * geo_ops.c
       4             :  *    2D geometric operations
       5             :  *
       6             :  * This module implements the geometric functions and operators.  The
       7             :  * geometric types are (from simple to more complicated):
       8             :  *
       9             :  * - point
      10             :  * - line
      11             :  * - line segment
      12             :  * - box
      13             :  * - circle
      14             :  * - polygon
      15             :  *
      16             :  * Portions Copyright (c) 1996-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       74552 : 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      149152 :     while (isspace((unsigned char) *str))
     211          48 :         str++;
     212       74552 :     if ((has_delim = (*str == LDELIM)))
     213       48664 :         str++;
     214             : 
     215       74552 :     *x = float8in_internal(str, &str, type_name, orig_string);
     216             : 
     217       74528 :     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       74512 :     *y = float8in_internal(str, &str, type_name, orig_string);
     224             : 
     225       74504 :     if (has_delim)
     226             :     {
     227       48640 :         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       97312 :         while (isspace((unsigned char) *str))
     233          40 :             str++;
     234             :     }
     235             : 
     236             :     /* report stopping point if wanted, else complain if not end of string */
     237       74500 :     if (endptr_p)
     238       73402 :         *endptr_p = str;
     239        1098 :     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       74496 : }
     245             : 
     246             : static void
     247      903044 : pair_encode(float8 x, float8 y, StringInfo str)
     248             : {
     249      903044 :     char       *xstr = float8out_internal(x);
     250      903044 :     char       *ystr = float8out_internal(y);
     251             : 
     252      903044 :     appendStringInfo(str, "%s,%s", xstr, ystr);
     253      903044 :     pfree(xstr);
     254      903044 :     pfree(ystr);
     255      903044 : }
     256             : 
     257             : static void
     258       34146 : path_decode(char *str, bool opentype, int npts, Point *p,
     259             :             bool *isopen, char **endptr_p,
     260             :             const char *type_name, const char *orig_string)
     261             : {
     262       34146 :     int         depth = 0;
     263             :     char       *cp;
     264             :     int         i;
     265             : 
     266       68296 :     while (isspace((unsigned char) *str))
     267           4 :         str++;
     268       34146 :     if ((*isopen = (*str == LDELIM_EP)))
     269             :     {
     270             :         /* no open delimiter allowed? */
     271       20254 :         if (!opentype)
     272           4 :             ereport(ERROR,
     273             :                     (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
     274             :                      errmsg("invalid input syntax for type %s: \"%s\"",
     275             :                             type_name, orig_string)));
     276       20250 :         depth++;
     277       20250 :         str++;
     278             :     }
     279       13892 :     else if (*str == LDELIM)
     280             :     {
     281       13796 :         cp = (str + 1);
     282       27608 :         while (isspace((unsigned char) *cp))
     283          16 :             cp++;
     284       13796 :         if (*cp == LDELIM)
     285             :         {
     286         804 :             depth++;
     287         804 :             str = cp;
     288             :         }
     289       12992 :         else if (strrchr(str, LDELIM) == str)
     290             :         {
     291       12696 :             depth++;
     292       12696 :             str = cp;
     293             :         }
     294             :     }
     295             : 
     296      107300 :     for (i = 0; i < npts; i++)
     297             :     {
     298       73186 :         pair_decode(str, &(p->x), &(p->y), &str, type_name, orig_string);
     299       73158 :         if (*str == DELIM)
     300       39020 :             str++;
     301       73158 :         p++;
     302             :     }
     303             : 
     304      101926 :     while (depth > 0)
     305             :     {
     306       67424 :         if (*str == RDELIM || (*str == RDELIM_EP && *isopen && depth == 1))
     307             :         {
     308       33698 :             depth--;
     309       33698 :             str++;
     310       67408 :             while (isspace((unsigned char) *str))
     311          12 :                 str++;
     312             :         }
     313             :         else
     314          28 :             ereport(ERROR,
     315             :                     (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
     316             :                      errmsg("invalid input syntax for type %s: \"%s\"",
     317             :                             type_name, orig_string)));
     318             :     }
     319             : 
     320             :     /* report stopping point if wanted, else complain if not end of string */
     321       34086 :     if (endptr_p)
     322       20582 :         *endptr_p = str;
     323       13504 :     else if (*str != '\0')
     324           4 :         ereport(ERROR,
     325             :                 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
     326             :                  errmsg("invalid input syntax for type %s: \"%s\"",
     327             :                         type_name, orig_string)));
     328       34082 : }                               /* path_decode() */
     329             : 
     330             : static char *
     331      303860 : path_encode(enum path_delim path_delim, int npts, Point *pt)
     332             : {
     333             :     StringInfoData str;
     334             :     int         i;
     335             : 
     336      303860 :     initStringInfo(&str);
     337             : 
     338      303860 :     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      220198 :             break;
     348             :     }
     349             : 
     350     1201052 :     for (i = 0; i < npts; i++)
     351             :     {
     352      897192 :         if (i > 0)
     353      593332 :             appendStringInfoChar(&str, DELIM);
     354      897192 :         appendStringInfoChar(&str, LDELIM);
     355      897192 :         pair_encode(pt->x, pt->y, &str);
     356      897192 :         appendStringInfoChar(&str, RDELIM);
     357      897192 :         pt++;
     358             :     }
     359             : 
     360      303860 :     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      220198 :             break;
     370             :     }
     371             : 
     372      303860 :     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       13212 : box_in(PG_FUNCTION_ARGS)
     414             : {
     415       13212 :     char       *str = PG_GETARG_CSTRING(0);
     416       13212 :     BOX        *box = (BOX *) palloc(sizeof(BOX));
     417             :     bool        isopen;
     418             :     float8      x,
     419             :                 y;
     420             : 
     421       13212 :     path_decode(str, false, 2, &(box->high), &isopen, NULL, "box", str);
     422             : 
     423             :     /* reorder corners if necessary... */
     424       13192 :     if (float8_lt(box->high.x, box->low.x))
     425             :     {
     426         600 :         x = box->high.x;
     427         600 :         box->high.x = box->low.x;
     428         600 :         box->low.x = x;
     429             :     }
     430       13192 :     if (float8_lt(box->high.y, box->low.y))
     431             :     {
     432         608 :         y = box->high.y;
     433         608 :         box->high.y = box->low.y;
     434         608 :         box->low.y = y;
     435             :     }
     436             : 
     437       13192 :     PG_RETURN_BOX_P(box);
     438             : }
     439             : 
     440             : /*      box_out -       convert a box to external form.
     441             :  */
     442             : Datum
     443       89188 : box_out(PG_FUNCTION_ARGS)
     444             : {
     445       89188 :     BOX        *box = PG_GETARG_BOX_P(0);
     446             : 
     447       89188 :     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      173508 : box_construct(BOX *result, Point *pt1, Point *pt2)
     507             : {
     508      173508 :     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      161136 :         result->high.x = pt2->x;
     516      161136 :         result->low.x = pt1->x;
     517             :     }
     518      173508 :     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      161156 :         result->high.y = pt2->y;
     526      161156 :         result->low.y = pt1->y;
     527             :     }
     528      173508 : }
     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       43180 : box_overlap(PG_FUNCTION_ARGS)
     552             : {
     553       43180 :     BOX        *box1 = PG_GETARG_BOX_P(0);
     554       43180 :     BOX        *box2 = PG_GETARG_BOX_P(1);
     555             : 
     556       43180 :     PG_RETURN_BOOL(box_ov(box1, box2));
     557             : }
     558             : 
     559             : static bool
     560      983948 : box_ov(BOX *box1, BOX *box2)
     561             : {
     562     1486380 :     return (FPle(box1->low.x, box2->high.x) &&
     563      584268 :             FPle(box2->low.x, box1->high.x) &&
     564     1131324 :             FPle(box1->low.y, box2->high.y) &&
     565       65540 :             FPle(box2->low.y, box1->high.y));
     566             : }
     567             : 
     568             : /*      box_left        -       is box1 strictly left of box2?
     569             :  */
     570             : Datum
     571       33348 : box_left(PG_FUNCTION_ARGS)
     572             : {
     573       33348 :     BOX        *box1 = PG_GETARG_BOX_P(0);
     574       33348 :     BOX        *box2 = PG_GETARG_BOX_P(1);
     575             : 
     576       33348 :     PG_RETURN_BOOL(FPlt(box1->high.x, box2->low.x));
     577             : }
     578             : 
     579             : /*      box_overleft    -       is the right edge of box1 at or left of
     580             :  *                              the right edge of box2?
     581             :  *
     582             :  *      This is "less than or equal" for the end of a time range,
     583             :  *      when time ranges are stored as rectangles.
     584             :  */
     585             : Datum
     586       66472 : box_overleft(PG_FUNCTION_ARGS)
     587             : {
     588       66472 :     BOX        *box1 = PG_GETARG_BOX_P(0);
     589       66472 :     BOX        *box2 = PG_GETARG_BOX_P(1);
     590             : 
     591       66472 :     PG_RETURN_BOOL(FPle(box1->high.x, box2->high.x));
     592             : }
     593             : 
     594             : /*      box_right       -       is box1 strictly right of box2?
     595             :  */
     596             : Datum
     597       87828 : box_right(PG_FUNCTION_ARGS)
     598             : {
     599       87828 :     BOX        *box1 = PG_GETARG_BOX_P(0);
     600       87828 :     BOX        *box2 = PG_GETARG_BOX_P(1);
     601             : 
     602       87828 :     PG_RETURN_BOOL(FPgt(box1->low.x, box2->high.x));
     603             : }
     604             : 
     605             : /*      box_overright   -       is the left edge of box1 at or right of
     606             :  *                              the left edge of box2?
     607             :  *
     608             :  *      This is "greater than or equal" for time ranges, when time ranges
     609             :  *      are stored as rectangles.
     610             :  */
     611             : Datum
     612       74932 : box_overright(PG_FUNCTION_ARGS)
     613             : {
     614       74932 :     BOX        *box1 = PG_GETARG_BOX_P(0);
     615       74932 :     BOX        *box2 = PG_GETARG_BOX_P(1);
     616             : 
     617       74932 :     PG_RETURN_BOOL(FPge(box1->low.x, box2->low.x));
     618             : }
     619             : 
     620             : /*      box_below       -       is box1 strictly below box2?
     621             :  */
     622             : Datum
     623       18452 : box_below(PG_FUNCTION_ARGS)
     624             : {
     625       18452 :     BOX        *box1 = PG_GETARG_BOX_P(0);
     626       18452 :     BOX        *box2 = PG_GETARG_BOX_P(1);
     627             : 
     628       18452 :     PG_RETURN_BOOL(FPlt(box1->high.y, box2->low.y));
     629             : }
     630             : 
     631             : /*      box_overbelow   -       is the upper edge of box1 at or below
     632             :  *                              the upper edge of box2?
     633             :  */
     634             : Datum
     635       53624 : box_overbelow(PG_FUNCTION_ARGS)
     636             : {
     637       53624 :     BOX        *box1 = PG_GETARG_BOX_P(0);
     638       53624 :     BOX        *box2 = PG_GETARG_BOX_P(1);
     639             : 
     640       53624 :     PG_RETURN_BOOL(FPle(box1->high.y, box2->high.y));
     641             : }
     642             : 
     643             : /*      box_above       -       is box1 strictly above box2?
     644             :  */
     645             : Datum
     646       38480 : box_above(PG_FUNCTION_ARGS)
     647             : {
     648       38480 :     BOX        *box1 = PG_GETARG_BOX_P(0);
     649       38480 :     BOX        *box2 = PG_GETARG_BOX_P(1);
     650             : 
     651       38480 :     PG_RETURN_BOOL(FPgt(box1->low.y, box2->high.y));
     652             : }
     653             : 
     654             : /*      box_overabove   -       is the lower edge of box1 at or above
     655             :  *                              the lower edge of box2?
     656             :  */
     657             : Datum
     658       74100 : box_overabove(PG_FUNCTION_ARGS)
     659             : {
     660       74100 :     BOX        *box1 = PG_GETARG_BOX_P(0);
     661       74100 :     BOX        *box2 = PG_GETARG_BOX_P(1);
     662             : 
     663       74100 :     PG_RETURN_BOOL(FPge(box1->low.y, box2->low.y));
     664             : }
     665             : 
     666             : /*      box_contained   -       is box1 contained by box2?
     667             :  */
     668             : Datum
     669      105528 : box_contained(PG_FUNCTION_ARGS)
     670             : {
     671      105528 :     BOX        *box1 = PG_GETARG_BOX_P(0);
     672      105528 :     BOX        *box2 = PG_GETARG_BOX_P(1);
     673             : 
     674      105528 :     PG_RETURN_BOOL(box_contain_box(box2, box1));
     675             : }
     676             : 
     677             : /*      box_contain     -       does box1 contain box2?
     678             :  */
     679             : Datum
     680        7804 : box_contain(PG_FUNCTION_ARGS)
     681             : {
     682        7804 :     BOX        *box1 = PG_GETARG_BOX_P(0);
     683        7804 :     BOX        *box2 = PG_GETARG_BOX_P(1);
     684             : 
     685        7804 :     PG_RETURN_BOOL(box_contain_box(box1, box2));
     686             : }
     687             : 
     688             : /*
     689             :  * Check whether the second box is in the first box or on its border
     690             :  */
     691             : static bool
     692      169936 : box_contain_box(BOX *contains_box, BOX *contained_box)
     693             : {
     694      277608 :     return FPge(contains_box->high.x, contained_box->high.x) &&
     695      184860 :         FPle(contains_box->low.x, contained_box->low.x) &&
     696      312124 :         FPge(contains_box->high.y, contained_box->high.y) &&
     697       65000 :         FPle(contains_box->low.y, contained_box->low.y);
     698             : }
     699             : 
     700             : 
     701             : /*      box_positionop  -
     702             :  *              is box1 entirely {above,below} box2?
     703             :  *
     704             :  * box_below_eq and box_above_eq are obsolete versions that (probably
     705             :  * erroneously) accept the equal-boundaries case.  Since these are not
     706             :  * in sync with the box_left and box_right code, they are deprecated and
     707             :  * not supported in the PG 8.1 rtree operator class extension.
     708             :  */
     709             : Datum
     710         100 : box_below_eq(PG_FUNCTION_ARGS)
     711             : {
     712         100 :     BOX        *box1 = PG_GETARG_BOX_P(0);
     713         100 :     BOX        *box2 = PG_GETARG_BOX_P(1);
     714             : 
     715         100 :     PG_RETURN_BOOL(FPle(box1->high.y, box2->low.y));
     716             : }
     717             : 
     718             : Datum
     719         100 : box_above_eq(PG_FUNCTION_ARGS)
     720             : {
     721         100 :     BOX        *box1 = PG_GETARG_BOX_P(0);
     722         100 :     BOX        *box2 = PG_GETARG_BOX_P(1);
     723             : 
     724         100 :     PG_RETURN_BOOL(FPge(box1->low.y, box2->high.y));
     725             : }
     726             : 
     727             : 
     728             : /*      box_relop       -       is area(box1) relop area(box2), within
     729             :  *                              our accuracy constraint?
     730             :  */
     731             : Datum
     732          20 : box_lt(PG_FUNCTION_ARGS)
     733             : {
     734          20 :     BOX        *box1 = PG_GETARG_BOX_P(0);
     735          20 :     BOX        *box2 = PG_GETARG_BOX_P(1);
     736             : 
     737          20 :     PG_RETURN_BOOL(FPlt(box_ar(box1), box_ar(box2)));
     738             : }
     739             : 
     740             : Datum
     741          20 : box_gt(PG_FUNCTION_ARGS)
     742             : {
     743          20 :     BOX        *box1 = PG_GETARG_BOX_P(0);
     744          20 :     BOX        *box2 = PG_GETARG_BOX_P(1);
     745             : 
     746          20 :     PG_RETURN_BOOL(FPgt(box_ar(box1), box_ar(box2)));
     747             : }
     748             : 
     749             : Datum
     750          20 : box_eq(PG_FUNCTION_ARGS)
     751             : {
     752          20 :     BOX        *box1 = PG_GETARG_BOX_P(0);
     753          20 :     BOX        *box2 = PG_GETARG_BOX_P(1);
     754             : 
     755          20 :     PG_RETURN_BOOL(FPeq(box_ar(box1), box_ar(box2)));
     756             : }
     757             : 
     758             : Datum
     759          20 : box_le(PG_FUNCTION_ARGS)
     760             : {
     761          20 :     BOX        *box1 = PG_GETARG_BOX_P(0);
     762          20 :     BOX        *box2 = PG_GETARG_BOX_P(1);
     763             : 
     764          20 :     PG_RETURN_BOOL(FPle(box_ar(box1), box_ar(box2)));
     765             : }
     766             : 
     767             : Datum
     768          20 : box_ge(PG_FUNCTION_ARGS)
     769             : {
     770          20 :     BOX        *box1 = PG_GETARG_BOX_P(0);
     771          20 :     BOX        *box2 = PG_GETARG_BOX_P(1);
     772             : 
     773          20 :     PG_RETURN_BOOL(FPge(box_ar(box1), box_ar(box2)));
     774             : }
     775             : 
     776             : 
     777             : /*----------------------------------------------------------
     778             :  *  "Arithmetic" operators on boxes.
     779             :  *---------------------------------------------------------*/
     780             : 
     781             : /*      box_area        -       returns the area of the box.
     782             :  */
     783             : Datum
     784          20 : box_area(PG_FUNCTION_ARGS)
     785             : {
     786          20 :     BOX        *box = PG_GETARG_BOX_P(0);
     787             : 
     788          20 :     PG_RETURN_FLOAT8(box_ar(box));
     789             : }
     790             : 
     791             : 
     792             : /*      box_width       -       returns the width of the box
     793             :  *                                (horizontal magnitude).
     794             :  */
     795             : Datum
     796          20 : box_width(PG_FUNCTION_ARGS)
     797             : {
     798          20 :     BOX        *box = PG_GETARG_BOX_P(0);
     799             : 
     800          20 :     PG_RETURN_FLOAT8(box_wd(box));
     801             : }
     802             : 
     803             : 
     804             : /*      box_height      -       returns the height of the box
     805             :  *                                (vertical magnitude).
     806             :  */
     807             : Datum
     808          20 : box_height(PG_FUNCTION_ARGS)
     809             : {
     810          20 :     BOX        *box = PG_GETARG_BOX_P(0);
     811             : 
     812          20 :     PG_RETURN_FLOAT8(box_ht(box));
     813             : }
     814             : 
     815             : 
     816             : /*      box_distance    -       returns the distance between the
     817             :  *                                center points of two boxes.
     818             :  */
     819             : Datum
     820         100 : box_distance(PG_FUNCTION_ARGS)
     821             : {
     822         100 :     BOX        *box1 = PG_GETARG_BOX_P(0);
     823         100 :     BOX        *box2 = PG_GETARG_BOX_P(1);
     824             :     Point       a,
     825             :                 b;
     826             : 
     827         100 :     box_cn(&a, box1);
     828         100 :     box_cn(&b, box2);
     829             : 
     830         100 :     PG_RETURN_FLOAT8(point_dt(&a, &b));
     831             : }
     832             : 
     833             : 
     834             : /*      box_center      -       returns the center point of the box.
     835             :  */
     836             : Datum
     837          60 : box_center(PG_FUNCTION_ARGS)
     838             : {
     839          60 :     BOX        *box = PG_GETARG_BOX_P(0);
     840          60 :     Point      *result = (Point *) palloc(sizeof(Point));
     841             : 
     842          60 :     box_cn(result, box);
     843             : 
     844          60 :     PG_RETURN_POINT_P(result);
     845             : }
     846             : 
     847             : 
     848             : /*      box_ar  -       returns the area of the box.
     849             :  */
     850             : static float8
     851         220 : box_ar(BOX *box)
     852             : {
     853         220 :     return float8_mul(box_wd(box), box_ht(box));
     854             : }
     855             : 
     856             : 
     857             : /*      box_cn  -       stores the centerpoint of the box into *center.
     858             :  */
     859             : static void
     860         316 : box_cn(Point *center, BOX *box)
     861             : {
     862         316 :     center->x = float8_div(float8_pl(box->high.x, box->low.x), 2.0);
     863         316 :     center->y = float8_div(float8_pl(box->high.y, box->low.y), 2.0);
     864         316 : }
     865             : 
     866             : 
     867             : /*      box_wd  -       returns the width (length) of the box
     868             :  *                                (horizontal magnitude).
     869             :  */
     870             : static float8
     871         240 : box_wd(BOX *box)
     872             : {
     873         240 :     return float8_mi(box->high.x, box->low.x);
     874             : }
     875             : 
     876             : 
     877             : /*      box_ht  -       returns the height of the box
     878             :  *                                (vertical magnitude).
     879             :  */
     880             : static float8
     881         240 : box_ht(BOX *box)
     882             : {
     883         240 :     return float8_mi(box->high.y, box->low.y);
     884             : }
     885             : 
     886             : 
     887             : /*----------------------------------------------------------
     888             :  *  Funky operations.
     889             :  *---------------------------------------------------------*/
     890             : 
     891             : /*      box_intersect   -
     892             :  *              returns the overlapping portion of two boxes,
     893             :  *                or NULL if they do not intersect.
     894             :  */
     895             : Datum
     896         100 : box_intersect(PG_FUNCTION_ARGS)
     897             : {
     898         100 :     BOX        *box1 = PG_GETARG_BOX_P(0);
     899         100 :     BOX        *box2 = PG_GETARG_BOX_P(1);
     900             :     BOX        *result;
     901             : 
     902         100 :     if (!box_ov(box1, box2))
     903          56 :         PG_RETURN_NULL();
     904             : 
     905          44 :     result = (BOX *) palloc(sizeof(BOX));
     906             : 
     907          44 :     result->high.x = float8_min(box1->high.x, box2->high.x);
     908          44 :     result->low.x = float8_max(box1->low.x, box2->low.x);
     909          44 :     result->high.y = float8_min(box1->high.y, box2->high.y);
     910          44 :     result->low.y = float8_max(box1->low.y, box2->low.y);
     911             : 
     912          44 :     PG_RETURN_BOX_P(result);
     913             : }
     914             : 
     915             : 
     916             : /*      box_diagonal    -
     917             :  *              returns a line segment which happens to be the
     918             :  *                positive-slope diagonal of "box".
     919             :  */
     920             : Datum
     921          20 : box_diagonal(PG_FUNCTION_ARGS)
     922             : {
     923          20 :     BOX        *box = PG_GETARG_BOX_P(0);
     924          20 :     LSEG       *result = (LSEG *) palloc(sizeof(LSEG));
     925             : 
     926          20 :     statlseg_construct(result, &box->high, &box->low);
     927             : 
     928          20 :     PG_RETURN_LSEG_P(result);
     929             : }
     930             : 
     931             : /***********************************************************************
     932             :  **
     933             :  **     Routines for 2D lines.
     934             :  **
     935             :  ***********************************************************************/
     936             : 
     937             : static bool
     938          60 : line_decode(char *s, const char *str, LINE *line)
     939             : {
     940             :     /* s was already advanced over leading '{' */
     941          60 :     line->A = single_decode(s, &s, "line", str);
     942          56 :     if (*s++ != DELIM)
     943           4 :         return false;
     944          52 :     line->B = single_decode(s, &s, "line", str);
     945          52 :     if (*s++ != DELIM)
     946           4 :         return false;
     947          48 :     line->C = single_decode(s, &s, "line", str);
     948          48 :     if (*s++ != RDELIM_L)
     949           4 :         return false;
     950          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     4347908 : line_construct(LINE *result, Point *pt, float8 m)
    1057             : {
    1058     4347908 :     if (m == DBL_MAX)
    1059             :     {
    1060             :         /* vertical - use "x = C" */
    1061      991196 :         result->A = -1.0;
    1062      991196 :         result->B = 0.0;
    1063      991196 :         result->C = pt->x;
    1064             :     }
    1065             :     else
    1066             :     {
    1067             :         /* use "mx - y + yinter = 0" */
    1068     3356712 :         result->A = m;
    1069     3356712 :         result->B = -1.0;
    1070     3356712 :         result->C = float8_mi(pt->y, float8_mul(m, pt->x));
    1071             :         /* on some platforms, the preceding expression tends to produce -0 */
    1072     3356712 :         if (result->C == 0.0)
    1073       30352 :             result->C = 0.0;
    1074             :     }
    1075     4347908 : }
    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     1189640 : line_invsl(LINE *line)
    1210             : {
    1211     1189640 :     if (FPzero(line->A))
    1212      474228 :         return DBL_MAX;
    1213      715412 :     if (FPzero(line->B))
    1214      463192 :         return 0.0;
    1215      252220 :     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     2771312 : line_interpt_line(Point *result, LINE *l1, LINE *l2)
    1277             : {
    1278             :     float8      x,
    1279             :                 y;
    1280             : 
    1281     2771312 :     if (!FPzero(l1->B))
    1282             :     {
    1283     2018256 :         if (FPeq(l2->A, float8_mul(l1->A, float8_div(l2->B, l1->B))))
    1284        6124 :             return false;
    1285             : 
    1286     2012132 :         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     2012132 :         y = float8_div(-float8_pl(float8_mul(l1->A, x), l1->C), l1->B);
    1291             :     }
    1292      753056 :     else if (!FPzero(l2->B))
    1293             :     {
    1294      751176 :         if (FPeq(l1->A, float8_mul(l2->A, float8_div(l1->B, l2->B))))
    1295           0 :             return false;
    1296             : 
    1297      751176 :         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      751176 :         y = float8_div(-float8_pl(float8_mul(l2->A, x), l2->C), l2->B);
    1302             :     }
    1303             :     else
    1304        1880 :         return false;
    1305             : 
    1306             :     /* On some platforms, the preceding expressions tend to produce -0. */
    1307     2763308 :     if (x == 0.0)
    1308       76608 :         x = 0.0;
    1309     2763308 :     if (y == 0.0)
    1310       79292 :         y = 0.0;
    1311             : 
    1312     2763308 :     if (result != NULL)
    1313     2762300 :         point_construct(result, x, y);
    1314             : 
    1315     2763308 :     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        1114 : point_in(PG_FUNCTION_ARGS)
    1791             : {
    1792        1114 :     char       *str = PG_GETARG_CSTRING(0);
    1793        1114 :     Point      *point = (Point *) palloc(sizeof(Point));
    1794             : 
    1795        1114 :     pair_decode(str, &point->x, &point->y, NULL, "point", str);
    1796        1094 :     PG_RETURN_POINT_P(point);
    1797             : }
    1798             : 
    1799             : Datum
    1800      131010 : point_out(PG_FUNCTION_ARGS)
    1801             : {
    1802      131010 :     Point      *pt = PG_GETARG_POINT_P(0);
    1803             : 
    1804      131010 :     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     3491118 : point_construct(Point *result, float8 x, float8 y)
    1843             : {
    1844     3491118 :     result->x = x;
    1845     3491118 :     result->y = y;
    1846     3491118 : }
    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      434532 : point_left(PG_FUNCTION_ARGS)
    1860             : {
    1861      434532 :     Point      *pt1 = PG_GETARG_POINT_P(0);
    1862      434532 :     Point      *pt2 = PG_GETARG_POINT_P(1);
    1863             : 
    1864      434532 :     PG_RETURN_BOOL(FPlt(pt1->x, pt2->x));
    1865             : }
    1866             : 
    1867             : Datum
    1868    11012732 : point_right(PG_FUNCTION_ARGS)
    1869             : {
    1870    11012732 :     Point      *pt1 = PG_GETARG_POINT_P(0);
    1871    11012732 :     Point      *pt2 = PG_GETARG_POINT_P(1);
    1872             : 
    1873    11012732 :     PG_RETURN_BOOL(FPgt(pt1->x, pt2->x));
    1874             : }
    1875             : 
    1876             : Datum
    1877    11158264 : point_above(PG_FUNCTION_ARGS)
    1878             : {
    1879    11158264 :     Point      *pt1 = PG_GETARG_POINT_P(0);
    1880    11158264 :     Point      *pt2 = PG_GETARG_POINT_P(1);
    1881             : 
    1882    11158264 :     PG_RETURN_BOOL(FPgt(pt1->y, pt2->y));
    1883             : }
    1884             : 
    1885             : Datum
    1886      746884 : point_below(PG_FUNCTION_ARGS)
    1887             : {
    1888      746884 :     Point      *pt1 = PG_GETARG_POINT_P(0);
    1889      746884 :     Point      *pt2 = PG_GETARG_POINT_P(1);
    1890             : 
    1891      746884 :     PG_RETURN_BOOL(FPlt(pt1->y, pt2->y));
    1892             : }
    1893             : 
    1894             : Datum
    1895      291808 : point_vert(PG_FUNCTION_ARGS)
    1896             : {
    1897      291808 :     Point      *pt1 = PG_GETARG_POINT_P(0);
    1898      291808 :     Point      *pt2 = PG_GETARG_POINT_P(1);
    1899             : 
    1900      291808 :     PG_RETURN_BOOL(FPeq(pt1->x, pt2->x));
    1901             : }
    1902             : 
    1903             : Datum
    1904      301340 : point_horiz(PG_FUNCTION_ARGS)
    1905             : {
    1906      301340 :     Point      *pt1 = PG_GETARG_POINT_P(0);
    1907      301340 :     Point      *pt2 = PG_GETARG_POINT_P(1);
    1908             : 
    1909      301340 :     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      199856 : point_eq_point(Point *pt1, Point *pt2)
    1939             : {
    1940      399728 :     return ((FPeq(pt1->x, pt2->x) && FPeq(pt1->y, pt2->y)) ||
    1941      177060 :             (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      470880 : point_distance(PG_FUNCTION_ARGS)
    1951             : {
    1952      470880 :     Point      *pt1 = PG_GETARG_POINT_P(0);
    1953      470880 :     Point      *pt2 = PG_GETARG_POINT_P(1);
    1954             : 
    1955      470880 :     PG_RETURN_FLOAT8(point_dt(pt1, pt2));
    1956             : }
    1957             : 
    1958             : static inline float8
    1959    10318528 : point_dt(Point *pt1, Point *pt2)
    1960             : {
    1961    10318528 :     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     2559656 : point_sl(Point *pt1, Point *pt2)
    1981             : {
    1982     2559656 :     if (FPeq(pt1->x, pt2->x))
    1983      283396 :         return DBL_MAX;
    1984     2276260 :     if (FPeq(pt1->y, pt2->y))
    1985      282580 :         return 0.0;
    1986     1993680 :     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      600792 : point_invsl(Point *pt1, Point *pt2)
    1997             : {
    1998      600792 :     if (FPeq(pt1->x, pt2->x))
    1999      237548 :         return 0.0;
    2000      363244 :     if (FPeq(pt1->y, pt2->y))
    2001      233844 :         return DBL_MAX;
    2002      129400 :     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      672356 : statlseg_construct(LSEG *lseg, Point *pt1, Point *pt2)
    2096             : {
    2097      672356 :     lseg->p[0].x = pt1->x;
    2098      672356 :     lseg->p[0].y = pt1->y;
    2099      672356 :     lseg->p[1].x = pt2->x;
    2100      672356 :     lseg->p[1].y = pt2->y;
    2101      672356 : }
    2102             : 
    2103             : 
    2104             : /*
    2105             :  * Return slope of the line segment
    2106             :  */
    2107             : static inline float8
    2108     2559328 : lseg_sl(LSEG *lseg)
    2109             : {
    2110     2559328 :     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      977640 : lseg_interpt_lseg(Point *result, LSEG *l1, LSEG *l2)
    2292             : {
    2293             :     Point       interpt;
    2294             :     LINE        tmp;
    2295             : 
    2296      977640 :     line_construct(&tmp, &l2->p[0], lseg_sl(l2));
    2297      977640 :     if (!lseg_interpt_line(&interpt, l1, &tmp))
    2298      922928 :         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       54712 :     if (!lseg_contain_point(l2, &interpt))
    2305       40980 :         return false;
    2306             : 
    2307       13732 :     if (result != NULL)
    2308        4476 :         *result = interpt;
    2309             : 
    2310       13732 :     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             :  * Distance from a line to a point
    2353             :  */
    2354             : Datum
    2355         360 : dist_lp(PG_FUNCTION_ARGS)
    2356             : {
    2357         360 :     LINE       *line = PG_GETARG_LINE_P(0);
    2358         360 :     Point      *pt = PG_GETARG_POINT_P(1);
    2359             : 
    2360         360 :     PG_RETURN_FLOAT8(line_closept_point(NULL, line, pt));
    2361             : }
    2362             : 
    2363             : /*
    2364             :  * Distance from a point to a lseg
    2365             :  */
    2366             : Datum
    2367         288 : dist_ps(PG_FUNCTION_ARGS)
    2368             : {
    2369         288 :     Point      *pt = PG_GETARG_POINT_P(0);
    2370         288 :     LSEG       *lseg = PG_GETARG_LSEG_P(1);
    2371             : 
    2372         288 :     PG_RETURN_FLOAT8(lseg_closept_point(NULL, lseg, pt));
    2373             : }
    2374             : 
    2375             : /*
    2376             :  * Distance from a lseg to a point
    2377             :  */
    2378             : Datum
    2379         288 : dist_sp(PG_FUNCTION_ARGS)
    2380             : {
    2381         288 :     LSEG       *lseg = PG_GETARG_LSEG_P(0);
    2382         288 :     Point      *pt = PG_GETARG_POINT_P(1);
    2383             : 
    2384         288 :     PG_RETURN_FLOAT8(lseg_closept_point(NULL, lseg, pt));
    2385             : }
    2386             : 
    2387             : static float8
    2388         648 : dist_ppath_internal(Point *pt, PATH *path)
    2389             : {
    2390         648 :     float8      result = 0.0;   /* keep compiler quiet */
    2391         648 :     bool        have_min = false;
    2392             :     float8      tmp;
    2393             :     int         i;
    2394             :     LSEG        lseg;
    2395             : 
    2396             :     Assert(path->npts > 0);
    2397             : 
    2398             :     /*
    2399             :      * The distance from a point to a path is the smallest distance from the
    2400             :      * point to any of its constituent segments.
    2401             :      */
    2402        2016 :     for (i = 0; i < path->npts; i++)
    2403             :     {
    2404             :         int         iprev;
    2405             : 
    2406        1368 :         if (i > 0)
    2407         720 :             iprev = i - 1;
    2408             :         else
    2409             :         {
    2410         648 :             if (!path->closed)
    2411         288 :                 continue;
    2412         360 :             iprev = path->npts - 1; /* Include the closure segment */
    2413             :         }
    2414             : 
    2415        1080 :         statlseg_construct(&lseg, &path->p[iprev], &path->p[i]);
    2416        1080 :         tmp = lseg_closept_point(NULL, &lseg, pt);
    2417        1080 :         if (!have_min || float8_lt(tmp, result))
    2418             :         {
    2419         680 :             result = tmp;
    2420         680 :             have_min = true;
    2421             :         }
    2422             :     }
    2423             : 
    2424         648 :     return result;
    2425             : }
    2426             : 
    2427             : /*
    2428             :  * Distance from a point to a path
    2429             :  */
    2430             : Datum
    2431         324 : dist_ppath(PG_FUNCTION_ARGS)
    2432             : {
    2433         324 :     Point      *pt = PG_GETARG_POINT_P(0);
    2434         324 :     PATH       *path = PG_GETARG_PATH_P(1);
    2435             : 
    2436         324 :     PG_RETURN_FLOAT8(dist_ppath_internal(pt, path));
    2437             : }
    2438             : 
    2439             : /*
    2440             :  * Distance from a path to a point
    2441             :  */
    2442             : Datum
    2443         324 : dist_pathp(PG_FUNCTION_ARGS)
    2444             : {
    2445         324 :     PATH       *path = PG_GETARG_PATH_P(0);
    2446         324 :     Point      *pt = PG_GETARG_POINT_P(1);
    2447             : 
    2448         324 :     PG_RETURN_FLOAT8(dist_ppath_internal(pt, path));
    2449             : }
    2450             : 
    2451             : /*
    2452             :  * Distance from a point to a box
    2453             :  */
    2454             : Datum
    2455         264 : dist_pb(PG_FUNCTION_ARGS)
    2456             : {
    2457         264 :     Point      *pt = PG_GETARG_POINT_P(0);
    2458         264 :     BOX        *box = PG_GETARG_BOX_P(1);
    2459             : 
    2460         264 :     PG_RETURN_FLOAT8(box_closept_point(NULL, box, pt));
    2461             : }
    2462             : 
    2463             : /*
    2464             :  * Distance from a box to a point
    2465             :  */
    2466             : Datum
    2467      103488 : dist_bp(PG_FUNCTION_ARGS)
    2468             : {
    2469      103488 :     BOX        *box = PG_GETARG_BOX_P(0);
    2470      103488 :     Point      *pt = PG_GETARG_POINT_P(1);
    2471             : 
    2472      103488 :     PG_RETURN_FLOAT8(box_closept_point(NULL, box, pt));
    2473             : }
    2474             : 
    2475             : /*
    2476             :  * Distance from a lseg to a line
    2477             :  */
    2478             : Datum
    2479         320 : dist_sl(PG_FUNCTION_ARGS)
    2480             : {
    2481         320 :     LSEG       *lseg = PG_GETARG_LSEG_P(0);
    2482         320 :     LINE       *line = PG_GETARG_LINE_P(1);
    2483             : 
    2484         320 :     PG_RETURN_FLOAT8(lseg_closept_line(NULL, lseg, line));
    2485             : }
    2486             : 
    2487             : /*
    2488             :  * Distance from a line to a lseg
    2489             :  */
    2490             : Datum
    2491         320 : dist_ls(PG_FUNCTION_ARGS)
    2492             : {
    2493         320 :     LINE       *line = PG_GETARG_LINE_P(0);
    2494         320 :     LSEG       *lseg = PG_GETARG_LSEG_P(1);
    2495             : 
    2496         320 :     PG_RETURN_FLOAT8(lseg_closept_line(NULL, lseg, line));
    2497             : }
    2498             : 
    2499             : /*
    2500             :  * Distance from a lseg to a box
    2501             :  */
    2502             : Datum
    2503         160 : dist_sb(PG_FUNCTION_ARGS)
    2504             : {
    2505         160 :     LSEG       *lseg = PG_GETARG_LSEG_P(0);
    2506         160 :     BOX        *box = PG_GETARG_BOX_P(1);
    2507             : 
    2508         160 :     PG_RETURN_FLOAT8(box_closept_lseg(NULL, box, lseg));
    2509             : }
    2510             : 
    2511             : /*
    2512             :  * Distance from a box to a lseg
    2513             :  */
    2514             : Datum
    2515         160 : dist_bs(PG_FUNCTION_ARGS)
    2516             : {
    2517         160 :     BOX        *box = PG_GETARG_BOX_P(0);
    2518         160 :     LSEG       *lseg = PG_GETARG_LSEG_P(1);
    2519             : 
    2520         160 :     PG_RETURN_FLOAT8(box_closept_lseg(NULL, box, lseg));
    2521             : }
    2522             : 
    2523             : /*
    2524             :  * Distance from a line to a box
    2525             :  */
    2526             : Datum
    2527           4 : dist_lb(PG_FUNCTION_ARGS)
    2528             : {
    2529             : #ifdef NOT_USED
    2530             :     LINE       *line = PG_GETARG_LINE_P(0);
    2531             :     BOX        *box = PG_GETARG_BOX_P(1);
    2532             : #endif
    2533             : 
    2534             :     /* need to think about this one for a while */
    2535           4 :     ereport(ERROR,
    2536             :             (errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
    2537             :              errmsg("function \"dist_lb\" not implemented")));
    2538             : 
    2539             :     PG_RETURN_NULL();
    2540             : }
    2541             : 
    2542             : /*
    2543             :  * Distance from a box to a line
    2544             :  */
    2545             : Datum
    2546           4 : dist_bl(PG_FUNCTION_ARGS)
    2547             : {
    2548             : #ifdef NOT_USED
    2549             :     BOX        *box = PG_GETARG_BOX_P(0);
    2550             :     LINE       *line = PG_GETARG_LINE_P(1);
    2551             : #endif
    2552             : 
    2553             :     /* need to think about this one for a while */
    2554           4 :     ereport(ERROR,
    2555             :             (errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
    2556             :              errmsg("function \"dist_bl\" not implemented")));
    2557             : 
    2558             :     PG_RETURN_NULL();
    2559             : }
    2560             : 
    2561             : static float8
    2562         224 : dist_cpoly_internal(CIRCLE *circle, POLYGON *poly)
    2563             : {
    2564             :     float8      result;
    2565             : 
    2566             :     /* calculate distance to center, and subtract radius */
    2567         224 :     result = float8_mi(dist_ppoly_internal(&circle->center, poly),
    2568             :                        circle->radius);
    2569         224 :     if (result < 0.0)
    2570         120 :         result = 0.0;
    2571             : 
    2572         224 :     return result;
    2573             : }
    2574             : 
    2575             : /*
    2576             :  * Distance from a circle to a polygon
    2577             :  */
    2578             : Datum
    2579         224 : dist_cpoly(PG_FUNCTION_ARGS)
    2580             : {
    2581         224 :     CIRCLE     *circle = PG_GETARG_CIRCLE_P(0);
    2582         224 :     POLYGON    *poly = PG_GETARG_POLYGON_P(1);
    2583             : 
    2584         224 :     PG_RETURN_FLOAT8(dist_cpoly_internal(circle, poly));
    2585             : }
    2586             : 
    2587             : /*
    2588             :  * Distance from a polygon to a circle
    2589             :  */
    2590             : Datum
    2591           0 : dist_polyc(PG_FUNCTION_ARGS)
    2592             : {
    2593           0 :     POLYGON    *poly = PG_GETARG_POLYGON_P(0);
    2594           0 :     CIRCLE     *circle = PG_GETARG_CIRCLE_P(1);
    2595             : 
    2596           0 :     PG_RETURN_FLOAT8(dist_cpoly_internal(circle, poly));
    2597             : }
    2598             : 
    2599             : /*
    2600             :  * Distance from a point to a polygon
    2601             :  */
    2602             : Datum
    2603         252 : dist_ppoly(PG_FUNCTION_ARGS)
    2604             : {
    2605         252 :     Point      *point = PG_GETARG_POINT_P(0);
    2606         252 :     POLYGON    *poly = PG_GETARG_POLYGON_P(1);
    2607             : 
    2608         252 :     PG_RETURN_FLOAT8(dist_ppoly_internal(point, poly));
    2609             : }
    2610             : 
    2611             : Datum
    2612       22712 : dist_polyp(PG_FUNCTION_ARGS)
    2613             : {
    2614       22712 :     POLYGON    *poly = PG_GETARG_POLYGON_P(0);
    2615       22712 :     Point      *point = PG_GETARG_POINT_P(1);
    2616             : 
    2617       22712 :     PG_RETURN_FLOAT8(dist_ppoly_internal(point, poly));
    2618             : }
    2619             : 
    2620             : static float8
    2621       23188 : dist_ppoly_internal(Point *pt, POLYGON *poly)
    2622             : {
    2623             :     float8      result;
    2624             :     float8      d;
    2625             :     int         i;
    2626             :     LSEG        seg;
    2627             : 
    2628       23188 :     if (point_inside(pt, poly->npts, poly->p) != 0)
    2629         120 :         return 0.0;
    2630             : 
    2631             :     /* initialize distance with segment between first and last points */
    2632       23068 :     seg.p[0].x = poly->p[0].x;
    2633       23068 :     seg.p[0].y = poly->p[0].y;
    2634       23068 :     seg.p[1].x = poly->p[poly->npts - 1].x;
    2635       23068 :     seg.p[1].y = poly->p[poly->npts - 1].y;
    2636       23068 :     result = lseg_closept_point(NULL, &seg, pt);
    2637             : 
    2638             :     /* check distances for other segments */
    2639      171400 :     for (i = 0; i < poly->npts - 1; i++)
    2640             :     {
    2641      148332 :         seg.p[0].x = poly->p[i].x;
    2642      148332 :         seg.p[0].y = poly->p[i].y;
    2643      148332 :         seg.p[1].x = poly->p[i + 1].x;
    2644      148332 :         seg.p[1].y = poly->p[i + 1].y;
    2645      148332 :         d = lseg_closept_point(NULL, &seg, pt);
    2646      148332 :         if (float8_lt(d, result))
    2647        4392 :             result = d;
    2648             :     }
    2649             : 
    2650       23068 :     return result;
    2651             : }
    2652             : 
    2653             : 
    2654             : /*---------------------------------------------------------------------
    2655             :  *      interpt_
    2656             :  *              Intersection point of objects.
    2657             :  *              We choose to ignore the "point" of intersection between
    2658             :  *                lines and boxes, since there are typically two.
    2659             :  *-------------------------------------------------------------------*/
    2660             : 
    2661             : /*
    2662             :  * Return whether the line segment intersect with the line. If *result is not
    2663             :  * NULL, it is set to the intersection point.
    2664             :  */
    2665             : static bool
    2666     1580072 : lseg_interpt_line(Point *result, LSEG *lseg, LINE *line)
    2667             : {
    2668             :     Point       interpt;
    2669             :     LINE        tmp;
    2670             : 
    2671             :     /*
    2672             :      * First, we promote the line segment to a line, because we know how to
    2673             :      * find the intersection point of two lines.  If they don't have an
    2674             :      * intersection point, we are done.
    2675             :      */
    2676     1580072 :     line_construct(&tmp, &lseg->p[0], lseg_sl(lseg));
    2677     1580072 :     if (!line_interpt_line(&interpt, &tmp, line))
    2678        7748 :         return false;
    2679             : 
    2680             :     /*
    2681             :      * Then, we check whether the intersection point is actually on the line
    2682             :      * segment.
    2683             :      */
    2684     1572324 :     if (!lseg_contain_point(lseg, &interpt))
    2685     1510312 :         return false;
    2686       62012 :     if (result != NULL)
    2687             :     {
    2688             :         /*
    2689             :          * If there is an intersection, then check explicitly for matching
    2690             :          * endpoints since there may be rounding effects with annoying LSB
    2691             :          * residue.
    2692             :          */
    2693       61756 :         if (point_eq_point(&lseg->p[0], &interpt))
    2694        4140 :             *result = lseg->p[0];
    2695       57616 :         else if (point_eq_point(&lseg->p[1], &interpt))
    2696        8752 :             *result = lseg->p[1];
    2697             :         else
    2698       48864 :             *result = interpt;
    2699             :     }
    2700             : 
    2701       62012 :     return true;
    2702             : }
    2703             : 
    2704             : /*---------------------------------------------------------------------
    2705             :  *      close_
    2706             :  *              Point of closest proximity between objects.
    2707             :  *-------------------------------------------------------------------*/
    2708             : 
    2709             : /*
    2710             :  * If *result is not NULL, it is set to the intersection point of a
    2711             :  * perpendicular of the line through the point.  Returns the distance
    2712             :  * of those two points.
    2713             :  */
    2714             : static float8
    2715     1189640 : line_closept_point(Point *result, LINE *line, Point *point)
    2716             : {
    2717             :     Point       closept;
    2718             :     LINE        tmp;
    2719             : 
    2720             :     /*
    2721             :      * We drop a perpendicular to find the intersection point.  Ordinarily we
    2722             :      * should always find it, but that can fail in the presence of NaN
    2723             :      * coordinates, and perhaps even from simple roundoff issues.
    2724             :      */
    2725     1189640 :     line_construct(&tmp, point, line_invsl(line));
    2726     1189640 :     if (!line_interpt_line(&closept, &tmp, line))
    2727             :     {
    2728           0 :         if (result != NULL)
    2729           0 :             *result = *point;
    2730             : 
    2731           0 :         return get_float8_nan();
    2732             :     }
    2733             : 
    2734     1189640 :     if (result != NULL)
    2735         360 :         *result = closept;
    2736             : 
    2737     1189640 :     return point_dt(&closept, point);
    2738             : }
    2739             : 
    2740             : Datum
    2741         360 : close_pl(PG_FUNCTION_ARGS)
    2742             : {
    2743         360 :     Point      *pt = PG_GETARG_POINT_P(0);
    2744         360 :     LINE       *line = PG_GETARG_LINE_P(1);
    2745             :     Point      *result;
    2746             : 
    2747         360 :     result = (Point *) palloc(sizeof(Point));
    2748             : 
    2749         360 :     if (isnan(line_closept_point(result, line, pt)))
    2750         112 :         PG_RETURN_NULL();
    2751             : 
    2752         248 :     PG_RETURN_POINT_P(result);
    2753             : }
    2754             : 
    2755             : 
    2756             : /*
    2757             :  * Closest point on line segment to specified point.
    2758             :  *
    2759             :  * If *result is not NULL, set it to the closest point on the line segment
    2760             :  * to the point.  Returns the distance of the two points.
    2761             :  */
    2762             : static float8
    2763      600536 : lseg_closept_point(Point *result, LSEG *lseg, Point *pt)
    2764             : {
    2765             :     Point       closept;
    2766             :     LINE        tmp;
    2767             : 
    2768             :     /*
    2769             :      * To find the closest point, we draw a perpendicular line from the point
    2770             :      * to the line segment.
    2771             :      */
    2772      600536 :     line_construct(&tmp, pt, point_invsl(&lseg->p[0], &lseg->p[1]));
    2773      600536 :     lseg_closept_line(&closept, lseg, &tmp);
    2774             : 
    2775      600536 :     if (result != NULL)
    2776      316596 :         *result = closept;
    2777             : 
    2778      600536 :     return point_dt(&closept, pt);
    2779             : }
    2780             : 
    2781             : Datum
    2782         288 : close_ps(PG_FUNCTION_ARGS)
    2783             : {
    2784         288 :     Point      *pt = PG_GETARG_POINT_P(0);
    2785         288 :     LSEG       *lseg = PG_GETARG_LSEG_P(1);
    2786             :     Point      *result;
    2787             : 
    2788         288 :     result = (Point *) palloc(sizeof(Point));
    2789             : 
    2790         288 :     if (isnan(lseg_closept_point(result, lseg, pt)))
    2791          60 :         PG_RETURN_NULL();
    2792             : 
    2793         228 :     PG_RETURN_POINT_P(result);
    2794             : }
    2795             : 
    2796             : 
    2797             : /*
    2798             :  * Closest point on line segment to line segment
    2799             :  */
    2800             : static float8
    2801        2896 : lseg_closept_lseg(Point *result, LSEG *on_lseg, LSEG *to_lseg)
    2802             : {
    2803             :     Point       point;
    2804             :     float8      dist,
    2805             :                 d;
    2806             : 
    2807             :     /* First, we handle the case when the line segments are intersecting. */
    2808        2896 :     if (lseg_interpt_lseg(result, on_lseg, to_lseg))
    2809          16 :         return 0.0;
    2810             : 
    2811             :     /*
    2812             :      * Then, we find the closest points from the endpoints of the second line
    2813             :      * segment, and keep the closest one.
    2814             :      */
    2815        2880 :     dist = lseg_closept_point(result, on_lseg, &to_lseg->p[0]);
    2816        2880 :     d = lseg_closept_point(&point, on_lseg, &to_lseg->p[1]);
    2817        2880 :     if (float8_lt(d, dist))
    2818             :     {
    2819         964 :         dist = d;
    2820         964 :         if (result != NULL)
    2821         544 :             *result = point;
    2822             :     }
    2823             : 
    2824             :     /* The closest point can still be one of the endpoints, so we test them. */
    2825        2880 :     d = lseg_closept_point(NULL, to_lseg, &on_lseg->p[0]);
    2826        2880 :     if (float8_lt(d, dist))
    2827             :     {
    2828         696 :         dist = d;
    2829         696 :         if (result != NULL)
    2830         496 :             *result = on_lseg->p[0];
    2831             :     }
    2832        2880 :     d = lseg_closept_point(NULL, to_lseg, &on_lseg->p[1]);
    2833        2880 :     if (float8_lt(d, dist))
    2834             :     {
    2835         272 :         dist = d;
    2836         272 :         if (result != NULL)
    2837         168 :             *result = on_lseg->p[1];
    2838             :     }
    2839             : 
    2840        2880 :     return dist;
    2841             : }
    2842             : 
    2843             : Datum
    2844         256 : close_lseg(PG_FUNCTION_ARGS)
    2845             : {
    2846         256 :     LSEG       *l1 = PG_GETARG_LSEG_P(0);
    2847         256 :     LSEG       *l2 = PG_GETARG_LSEG_P(1);
    2848             :     Point      *result;
    2849             : 
    2850         256 :     if (lseg_sl(l1) == lseg_sl(l2))
    2851          52 :         PG_RETURN_NULL();
    2852             : 
    2853         204 :     result = (Point *) palloc(sizeof(Point));
    2854             : 
    2855         204 :     if (isnan(lseg_closept_lseg(result, l2, l1)))
    2856          60 :         PG_RETURN_NULL();
    2857             : 
    2858         144 :     PG_RETURN_POINT_P(result);
    2859             : }
    2860             : 
    2861             : 
    2862             : /*
    2863             :  * Closest point on or in box to specified point.
    2864             :  *
    2865             :  * If *result is not NULL, set it to the closest point on the box to the
    2866             :  * given point, and return the distance of the two points.
    2867             :  */
    2868             : static float8
    2869      103932 : box_closept_point(Point *result, BOX *box, Point *pt)
    2870             : {
    2871             :     float8      dist,
    2872             :                 d;
    2873             :     Point       point,
    2874             :                 closept;
    2875             :     LSEG        lseg;
    2876             : 
    2877      103932 :     if (box_contain_point(box, pt))
    2878             :     {
    2879          28 :         if (result != NULL)
    2880           4 :             *result = *pt;
    2881             : 
    2882          28 :         return 0.0;
    2883             :     }
    2884             : 
    2885             :     /* pairwise check lseg distances */
    2886      103904 :     point.x = box->low.x;
    2887      103904 :     point.y = box->high.y;
    2888      103904 :     statlseg_construct(&lseg, &box->low, &point);
    2889      103904 :     dist = lseg_closept_point(result, &lseg, pt);
    2890             : 
    2891      103904 :     statlseg_construct(&lseg, &box->high, &point);
    2892      103904 :     d = lseg_closept_point(&closept, &lseg, pt);
    2893      103904 :     if (float8_lt(d, dist))
    2894             :     {
    2895        4956 :         dist = d;
    2896        4956 :         if (result != NULL)
    2897          36 :             *result = closept;
    2898             :     }
    2899             : 
    2900      103904 :     point.x = box->high.x;
    2901      103904 :     point.y = box->low.y;
    2902      103904 :     statlseg_construct(&lseg, &box->low, &point);
    2903      103904 :     d = lseg_closept_point(&closept, &lseg, pt);
    2904      103904 :     if (float8_lt(d, dist))
    2905             :     {
    2906        5292 :         dist = d;
    2907        5292 :         if (result != NULL)
    2908           4 :             *result = closept;
    2909             :     }
    2910             : 
    2911      103904 :     statlseg_construct(&lseg, &box->high, &point);
    2912      103904 :     d = lseg_closept_point(&closept, &lseg, pt);
    2913      103904 :     if (float8_lt(d, dist))
    2914             :     {
    2915          24 :         dist = d;
    2916          24 :         if (result != NULL)
    2917           8 :             *result = closept;
    2918             :     }
    2919             : 
    2920      103904 :     return dist;
    2921             : }
    2922             : 
    2923             : Datum
    2924         180 : close_pb(PG_FUNCTION_ARGS)
    2925             : {
    2926         180 :     Point      *pt = PG_GETARG_POINT_P(0);
    2927         180 :     BOX        *box = PG_GETARG_BOX_P(1);
    2928             :     Point      *result;
    2929             : 
    2930         180 :     result = (Point *) palloc(sizeof(Point));
    2931             : 
    2932         180 :     if (isnan(box_closept_point(result, box, pt)))
    2933          20 :         PG_RETURN_NULL();
    2934             : 
    2935         160 :     PG_RETURN_POINT_P(result);
    2936             : }
    2937             : 
    2938             : 
    2939             : /* close_sl()
    2940             :  * Closest point on line to line segment.
    2941             :  *
    2942             :  * XXX THIS CODE IS WRONG
    2943             :  * The code is actually calculating the point on the line segment
    2944             :  *  which is backwards from the routine naming convention.
    2945             :  * Copied code to new routine close_ls() but haven't fixed this one yet.
    2946             :  * - thomas 1998-01-31
    2947             :  */
    2948             : Datum
    2949           4 : close_sl(PG_FUNCTION_ARGS)
    2950             : {
    2951             : #ifdef NOT_USED
    2952             :     LSEG       *lseg = PG_GETARG_LSEG_P(0);
    2953             :     LINE       *line = PG_GETARG_LINE_P(1);
    2954             :     Point      *result;
    2955             :     float8      d1,
    2956             :                 d2;
    2957             : 
    2958             :     result = (Point *) palloc(sizeof(Point));
    2959             : 
    2960             :     if (lseg_interpt_line(result, lseg, line))
    2961             :         PG_RETURN_POINT_P(result);
    2962             : 
    2963             :     d1 = line_closept_point(NULL, line, &lseg->p[0]);
    2964             :     d2 = line_closept_point(NULL, line, &lseg->p[1]);
    2965             :     if (float8_lt(d1, d2))
    2966             :         *result = lseg->p[0];
    2967             :     else
    2968             :         *result = lseg->p[1];
    2969             : 
    2970             :     PG_RETURN_POINT_P(result);
    2971             : #endif
    2972             : 
    2973           4 :     ereport(ERROR,
    2974             :             (errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
    2975             :              errmsg("function \"close_sl\" not implemented")));
    2976             : 
    2977             :     PG_RETURN_NULL();
    2978             : }
    2979             : 
    2980             : /*
    2981             :  * Closest point on line segment to line.
    2982             :  *
    2983             :  * Return the distance between the line and the closest point of the line
    2984             :  * segment to the line.  If *result is not NULL, set it to that point.
    2985             :  *
    2986             :  * NOTE: When the lines are parallel, endpoints of one of the line segment
    2987             :  * are FPeq(), in presence of NaN or Infinite coordinates, or perhaps =
    2988             :  * even because of simple roundoff issues, there may not be a single closest
    2989             :  * point.  We are likely to set the result to the second endpoint in these
    2990             :  * cases.
    2991             :  */
    2992             : static float8
    2993      601460 : lseg_closept_line(Point *result, LSEG *lseg, LINE *line)
    2994             : {
    2995             :     float8      dist1,
    2996             :                 dist2;
    2997             : 
    2998      601460 :     if (lseg_interpt_line(result, lseg, line))
    2999        7180 :         return 0.0;
    3000             : 
    3001      594280 :     dist1 = line_closept_point(NULL, line, &lseg->p[0]);
    3002      594280 :     dist2 = line_closept_point(NULL, line, &lseg->p[1]);
    3003             : 
    3004      594280 :     if (dist1 < dist2)
    3005             :     {
    3006      316272 :         if (result != NULL)
    3007      316136 :             *result = lseg->p[0];
    3008             : 
    3009      316272 :         return dist1;
    3010             :     }
    3011             :     else
    3012             :     {
    3013      278008 :         if (result != NULL)
    3014      277640 :             *result = lseg->p[1];
    3015             : 
    3016      278008 :         return dist2;
    3017             :     }
    3018             : }
    3019             : 
    3020             : Datum
    3021         320 : close_ls(PG_FUNCTION_ARGS)
    3022             : {
    3023         320 :     LINE       *line = PG_GETARG_LINE_P(0);
    3024         320 :     LSEG       *lseg = PG_GETARG_LSEG_P(1);
    3025             :     Point      *result;
    3026             : 
    3027         320 :     if (lseg_sl(lseg) == line_sl(line))
    3028          36 :         PG_RETURN_NULL();
    3029             : 
    3030         284 :     result = (Point *) palloc(sizeof(Point));
    3031             : 
    3032         284 :     if (isnan(lseg_closept_line(result, lseg, line)))
    3033          96 :         PG_RETURN_NULL();
    3034             : 
    3035         188 :     PG_RETURN_POINT_P(result);
    3036             : }
    3037             : 
    3038             : 
    3039             : /*
    3040             :  * Closest point on or in box to line segment.
    3041             :  *
    3042             :  * Returns the distance between the closest point on or in the box to
    3043             :  * the line segment.  If *result is not NULL, it is set to that point.
    3044             :  */
    3045             : static float8
    3046         480 : box_closept_lseg(Point *result, BOX *box, LSEG *lseg)
    3047             : {
    3048             :     float8      dist,
    3049             :                 d;
    3050             :     Point       point,
    3051             :                 closept;
    3052             :     LSEG        bseg;
    3053             : 
    3054         480 :     if (box_interpt_lseg(result, box, lseg))
    3055          96 :         return 0.0;
    3056             : 
    3057             :     /* pairwise check lseg distances */
    3058         384 :     point.x = box->low.x;
    3059         384 :     point.y = box->high.y;
    3060         384 :     statlseg_construct(&bseg, &box->low, &point);
    3061         384 :     dist = lseg_closept_lseg(result, &bseg, lseg);
    3062             : 
    3063         384 :     statlseg_construct(&bseg, &box->high, &point);
    3064         384 :     d = lseg_closept_lseg(&closept, &bseg, lseg);
    3065         384 :     if (float8_lt(d, dist))
    3066             :     {
    3067          96 :         dist = d;
    3068          96 :         if (result != NULL)
    3069          32 :             *result = closept;
    3070             :     }
    3071             : 
    3072         384 :     point.x = box->high.x;
    3073         384 :     point.y = box->low.y;
    3074         384 :     statlseg_construct(&bseg, &box->low, &point);
    3075         384 :     d = lseg_closept_lseg(&closept, &bseg, lseg);
    3076         384 :     if (float8_lt(d, dist))
    3077             :     {
    3078          12 :         dist = d;
    3079          12 :         if (result != NULL)
    3080           4 :             *result = closept;
    3081             :     }
    3082             : 
    3083         384 :     statlseg_construct(&bseg, &box->high, &point);
    3084         384 :     d = lseg_closept_lseg(&closept, &bseg, lseg);
    3085         384 :     if (float8_lt(d, dist))
    3086             :     {
    3087          12 :         dist = d;
    3088          12 :         if (result != NULL)
    3089           4 :             *result = closept;
    3090             :     }
    3091             : 
    3092         384 :     return dist;
    3093             : }
    3094             : 
    3095             : Datum
    3096         160 : close_sb(PG_FUNCTION_ARGS)
    3097             : {
    3098         160 :     LSEG       *lseg = PG_GETARG_LSEG_P(0);
    3099         160 :     BOX        *box = PG_GETARG_BOX_P(1);
    3100             :     Point      *result;
    3101             : 
    3102         160 :     result = (Point *) palloc(sizeof(Point));
    3103             : 
    3104         160 :     if (isnan(box_closept_lseg(result, box, lseg)))
    3105          20 :         PG_RETURN_NULL();
    3106             : 
    3107         140 :     PG_RETURN_POINT_P(result);
    3108             : }
    3109             : 
    3110             : 
    3111             : Datum
    3112           4 : close_lb(PG_FUNCTION_ARGS)
    3113             : {
    3114             : #ifdef NOT_USED
    3115             :     LINE       *line = PG_GETARG_LINE_P(0);
    3116             :     BOX        *box = PG_GETARG_BOX_P(1);
    3117             : #endif
    3118             : 
    3119             :     /* think about this one for a while */
    3120           4 :     ereport(ERROR,
    3121             :             (errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
    3122             :              errmsg("function \"close_lb\" not implemented")));
    3123             : 
    3124             :     PG_RETURN_NULL();
    3125             : }
    3126             : 
    3127             : /*---------------------------------------------------------------------
    3128             :  *      on_
    3129             :  *              Whether one object lies completely within another.
    3130             :  *-------------------------------------------------------------------*/
    3131             : 
    3132             : /*
    3133             :  *      Does the point satisfy the equation?
    3134             :  */
    3135             : static bool
    3136         696 : line_contain_point(LINE *line, Point *point)
    3137             : {
    3138         696 :     return FPzero(float8_pl(float8_pl(float8_mul(line->A, point->x),
    3139             :                                       float8_mul(line->B, point->y)),
    3140             :                             line->C));
    3141             : }
    3142             : 
    3143             : Datum
    3144         360 : on_pl(PG_FUNCTION_ARGS)
    3145             : {
    3146         360 :     Point      *pt = PG_GETARG_POINT_P(0);
    3147         360 :     LINE       *line = PG_GETARG_LINE_P(1);
    3148             : 
    3149         360 :     PG_RETURN_BOOL(line_contain_point(line, pt));
    3150             : }
    3151             : 
    3152             : 
    3153             : /*
    3154             :  *      Determine colinearity by detecting a triangle inequality.
    3155             :  * This algorithm seems to behave nicely even with lsb residues - tgl 1997-07-09
    3156             :  */
    3157             : static bool
    3158     2671828 : lseg_contain_point(LSEG *lseg, Point *pt)
    3159             : {
    3160     2671828 :     return FPeq(point_dt(pt, &lseg->p[0]) +
    3161             :                 point_dt(pt, &lseg->p[1]),
    3162             :                 point_dt(&lseg->p[0], &lseg->p[1]));
    3163             : }
    3164             : 
    3165             : Datum
    3166         288 : on_ps(PG_FUNCTION_ARGS)
    3167             : {
    3168         288 :     Point      *pt = PG_GETARG_POINT_P(0);
    3169         288 :     LSEG       *lseg = PG_GETARG_LSEG_P(1);
    3170             : 
    3171         288 :     PG_RETURN_BOOL(lseg_contain_point(lseg, pt));
    3172             : }
    3173             : 
    3174             : 
    3175             : /*
    3176             :  * Check whether the point is in the box or on its border
    3177             :  */
    3178             : static bool
    3179      291768 : box_contain_point(BOX *box, Point *point)
    3180             : {
    3181      528560 :     return box->high.x >= point->x && box->low.x <= point->x &&
    3182      415968 :         box->high.y >= point->y && box->low.y <= point->y;
    3183             : }
    3184             : 
    3185             : Datum
    3186       92152 : on_pb(PG_FUNCTION_ARGS)
    3187             : {
    3188       92152 :     Point      *pt = PG_GETARG_POINT_P(0);
    3189       92152 :     BOX        *box = PG_GETARG_BOX_P(1);
    3190             : 
    3191       92152 :     PG_RETURN_BOOL(box_contain_point(box, pt));
    3192             : }
    3193             : 
    3194             : Datum
    3195       95112 : box_contain_pt(PG_FUNCTION_ARGS)
    3196             : {
    3197       95112 :     BOX        *box = PG_GETARG_BOX_P(0);
    3198       95112 :     Point      *pt = PG_GETARG_POINT_P(1);
    3199             : 
    3200       95112 :     PG_RETURN_BOOL(box_contain_point(box, pt));
    3201             : }
    3202             : 
    3203             : /* on_ppath -
    3204             :  *      Whether a point lies within (on) a polyline.
    3205             :  *      If open, we have to (groan) check each segment.
    3206             :  * (uses same algorithm as for point intersecting segment - tgl 1997-07-09)
    3207             :  *      If closed, we use the old O(n) ray method for point-in-polygon.
    3208             :  *              The ray is horizontal, from pt out to the right.
    3209             :  *              Each segment that crosses the ray counts as an
    3210             :  *              intersection; note that an endpoint or edge may touch
    3211             :  *              but not cross.
    3212             :  *              (we can do p-in-p in lg(n), but it takes preprocessing)
    3213             :  */
    3214             : Datum
    3215         360 : on_ppath(PG_FUNCTION_ARGS)
    3216             : {
    3217         360 :     Point      *pt = PG_GETARG_POINT_P(0);
    3218         360 :     PATH       *path = PG_GETARG_PATH_P(1);
    3219             :     int         i,
    3220             :                 n;
    3221             :     float8      a,
    3222             :                 b;
    3223             : 
    3224             :     /*-- OPEN --*/
    3225         360 :     if (!path->closed)
    3226             :     {
    3227         180 :         n = path->npts - 1;
    3228         180 :         a = point_dt(pt, &path->p[0]);
    3229         420 :         for (i = 0; i < n; i++)
    3230             :         {
    3231         260 :             b = point_dt(pt, &path->p[i + 1]);
    3232         260 :             if (FPeq(float8_pl(a, b), point_dt(&path->p[i], &path->p[i + 1])))
    3233          20 :                 PG_RETURN_BOOL(true);
    3234         240 :             a = b;
    3235             :         }
    3236         160 :         PG_RETURN_BOOL(false);
    3237             :     }
    3238             : 
    3239             :     /*-- CLOSED --*/
    3240         180 :     PG_RETURN_BOOL(point_inside(pt, path->npts, path->p) != 0);
    3241             : }
    3242             : 
    3243             : 
    3244             : /*
    3245             :  * Check whether the line segment is on the line or close enough
    3246             :  *
    3247             :  * It is, if both of its points are on the line or close enough.
    3248             :  */
    3249             : Datum
    3250         320 : on_sl(PG_FUNCTION_ARGS)
    3251             : {
    3252         320 :     LSEG       *lseg = PG_GETARG_LSEG_P(0);
    3253         320 :     LINE       *line = PG_GETARG_LINE_P(1);
    3254             : 
    3255         320 :     PG_RETURN_BOOL(line_contain_point(line, &lseg->p[0]) &&
    3256             :                    line_contain_point(line, &lseg->p[1]));
    3257             : }
    3258             : 
    3259             : 
    3260             : /*
    3261             :  * Check whether the line segment is in the box or on its border
    3262             :  *
    3263             :  * It is, if both of its points are in the box or on its border.
    3264             :  */
    3265             : static bool
    3266         160 : box_contain_lseg(BOX *box, LSEG *lseg)
    3267             : {
    3268         172 :     return box_contain_point(box, &lseg->p[0]) &&
    3269          12 :         box_contain_point(box, &lseg->p[1]);
    3270             : }
    3271             : 
    3272             : Datum
    3273         160 : on_sb(PG_FUNCTION_ARGS)
    3274             : {
    3275         160 :     LSEG       *lseg = PG_GETARG_LSEG_P(0);
    3276         160 :     BOX        *box = PG_GETARG_BOX_P(1);
    3277             : 
    3278         160 :     PG_RETURN_BOOL(box_contain_lseg(box, lseg));
    3279             : }
    3280             : 
    3281             : /*---------------------------------------------------------------------
    3282             :  *      inter_
    3283             :  *              Whether one object intersects another.
    3284             :  *-------------------------------------------------------------------*/
    3285             : 
    3286             : Datum
    3287         320 : inter_sl(PG_FUNCTION_ARGS)
    3288             : {
    3289         320 :     LSEG       *lseg = PG_GETARG_LSEG_P(0);
    3290         320 :     LINE       *line = PG_GETARG_LINE_P(1);
    3291             : 
    3292         320 :     PG_RETURN_BOOL(lseg_interpt_line(NULL, lseg, line));
    3293             : }
    3294             : 
    3295             : 
    3296             : /*
    3297             :  * Do line segment and box intersect?
    3298             :  *
    3299             :  * Segment completely inside box counts as intersection.
    3300             :  * If you want only segments crossing box boundaries,
    3301             :  *  try converting box to path first.
    3302             :  *
    3303             :  * This function also sets the *result to the closest point on the line
    3304             :  * segment to the center of the box when they overlap and the result is
    3305             :  * not NULL.  It is somewhat arbitrary, but maybe the best we can do as
    3306             :  * there are typically two points they intersect.
    3307             :  *
    3308             :  * Optimize for non-intersection by checking for box intersection first.
    3309             :  * - thomas 1998-01-30
    3310             :  */
    3311             : static bool
    3312         640 : box_interpt_lseg(Point *result, BOX *box, LSEG *lseg)
    3313             : {
    3314             :     BOX         lbox;
    3315             :     LSEG        bseg;
    3316             :     Point       point;
    3317             : 
    3318         640 :     lbox.low.x = float8_min(lseg->p[0].x, lseg->p[1].x);
    3319         640 :     lbox.low.y = float8_min(lseg->p[0].y, lseg->p[1].y);
    3320         640 :     lbox.high.x = float8_max(lseg->p[0].x, lseg->p[1].x);
    3321         640 :     lbox.high.y = float8_max(lseg->p[0].y, lseg->p[1].y);
    3322             : 
    3323             :     /* nothing close to overlap? then not going to intersect */
    3324         640 :     if (!box_ov(&lbox, box))
    3325         416 :         return false;
    3326             : 
    3327         224 :     if (result != NULL)
    3328             :     {
    3329          56 :         box_cn(&point, box);
    3330          56 :         lseg_closept_point(result, lseg, &point);
    3331             :     }
    3332             : 
    3333             :     /* an endpoint of segment is inside box? then clearly intersects */
    3334         400 :     if (box_contain_point(box, &lseg->p[0]) ||
    3335         176 :         box_contain_point(box, &lseg->p[1]))
    3336          64 :         return true;
    3337             : 
    3338             :     /* pairwise check lseg intersections */
    3339         160 :     point.x = box->low.x;
    3340         160 :     point.y = box->high.y;
    3341         160 :     statlseg_construct(&bseg, &box->low, &point);
    3342         160 :     if (lseg_interpt_lseg(NULL, &bseg, lseg))
    3343          64 :         return true;
    3344             : 
    3345          96 :     statlseg_construct(&bseg, &box->high, &point);
    3346          96 :     if (lseg_interpt_lseg(NULL, &bseg, lseg))
    3347           0 :         return true;
    3348             : 
    3349          96 :     point.x = box->high.x;
    3350          96 :     point.y = box->low.y;
    3351          96 :     statlseg_construct(&bseg, &box->low, &point);
    3352          96 :     if (lseg_interpt_lseg(NULL, &bseg, lseg))
    3353           0 :         return true;
    3354             : 
    3355          96 :     statlseg_construct(&bseg, &box->high, &point);
    3356          96 :     if (lseg_interpt_lseg(NULL, &bseg, lseg))
    3357           0 :         return true;
    3358             : 
    3359             :     /* if we dropped through, no two segs intersected */
    3360          96 :     return false;
    3361             : }
    3362             : 
    3363             : Datum
    3364         160 : inter_sb(PG_FUNCTION_ARGS)
    3365             : {
    3366         160 :     LSEG       *lseg = PG_GETARG_LSEG_P(0);
    3367         160 :     BOX        *box = PG_GETARG_BOX_P(1);
    3368             : 
    3369         160 :     PG_RETURN_BOOL(box_interpt_lseg(NULL, box, lseg));
    3370             : }
    3371             : 
    3372             : 
    3373             : /* inter_lb()
    3374             :  * Do line and box intersect?
    3375             :  */
    3376             : Datum
    3377         200 : inter_lb(PG_FUNCTION_ARGS)
    3378             : {
    3379         200 :     LINE       *line = PG_GETARG_LINE_P(0);
    3380         200 :     BOX        *box = PG_GETARG_BOX_P(1);
    3381             :     LSEG        bseg;
    3382             :     Point       p1,
    3383             :                 p2;
    3384             : 
    3385             :     /* pairwise check lseg intersections */
    3386         200 :     p1.x = box->low.x;
    3387         200 :     p1.y = box->low.y;
    3388         200 :     p2.x = box->low.x;
    3389         200 :     p2.y = box->high.y;
    3390         200 :     statlseg_construct(&bseg, &p1, &p2);
    3391         200 :     if (lseg_interpt_line(NULL, &bseg, line))
    3392          44 :         PG_RETURN_BOOL(true);
    3393         156 :     p1.x = box->high.x;
    3394         156 :     p1.y = box->high.y;
    3395         156 :     statlseg_construct(&bseg, &p1, &p2);
    3396         156 :     if (lseg_interpt_line(NULL, &bseg, line))
    3397           8 :         PG_RETURN_BOOL(true);
    3398         148 :     p2.x = box->high.x;
    3399         148 :     p2.y = box->low.y;
    3400         148 :     statlseg_construct(&bseg, &p1, &p2);
    3401         148 :     if (lseg_interpt_line(NULL, &bseg, line))
    3402           0 :         PG_RETURN_BOOL(true);
    3403         148 :     p1.x = box->low.x;
    3404         148 :     p1.y = box->low.y;
    3405         148 :     statlseg_construct(&bseg, &p1, &p2);
    3406         148 :     if (lseg_interpt_line(NULL, &bseg, line))
    3407           0 :         PG_RETURN_BOOL(true);
    3408             : 
    3409             :     /* if we dropped through, no intersection */
    3410         148 :     PG_RETURN_BOOL(false);
    3411             : }
    3412             : 
    3413             : /*------------------------------------------------------------------
    3414             :  * The following routines define a data type and operator class for
    3415             :  * POLYGONS .... Part of which (the polygon's bounding box) is built on
    3416             :  * top of the BOX data type.
    3417             :  *
    3418             :  * make_bound_box - create the bounding box for the input polygon
    3419             :  *------------------------------------------------------------------*/
    3420             : 
    3421             : /*---------------------------------------------------------------------
    3422             :  * Make the smallest bounding box for the given polygon.
    3423             :  *---------------------------------------------------------------------*/
    3424             : static void
    3425       40346 : make_bound_box(POLYGON *poly)
    3426             : {
    3427             :     int         i;
    3428             :     float8      x1,
    3429             :                 y1,
    3430             :                 x2,
    3431             :                 y2;
    3432             : 
    3433             :     Assert(poly->npts > 0);
    3434             : 
    3435       40346 :     x1 = x2 = poly->p[0].x;
    3436       40346 :     y2 = y1 = poly->p[0].y;
    3437      481380 :     for (i = 1; i < poly->npts; i++)
    3438             :     {
    3439      441034 :         if (float8_lt(poly->p[i].x, x1))
    3440          48 :             x1 = poly->p[i].x;
    3441      441034 :         if (float8_gt(poly->p[i].x, x2))
    3442      240702 :             x2 = poly->p[i].x;
    3443      441034 :         if (float8_lt(poly->p[i].y, y1))
    3444      120268 :             y1 = poly->p[i].y;
    3445      441034 :         if (float8_gt(poly->p[i].y, y2))
    3446      120402 :             y2 = poly->p[i].y;
    3447             :     }
    3448             : 
    3449       40346 :     poly->boundbox.low.x = x1;
    3450       40346 :     poly->boundbox.high.x = x2;
    3451       40346 :     poly->boundbox.low.y = y1;
    3452       40346 :     poly->boundbox.high.y = y2;
    3453       40346 : }
    3454             : 
    3455             : /*------------------------------------------------------------------
    3456             :  * poly_in - read in the polygon from a string specification
    3457             :  *
    3458             :  *      External format:
    3459             :  *              "((x0,y0),...,(xn,yn))"
    3460             :  *              "x0,y0,...,xn,yn"
    3461             :  *              also supports the older style "(x1,...,xn,y1,...yn)"
    3462             :  *------------------------------------------------------------------*/
    3463             : Datum
    3464         262 : poly_in(PG_FUNCTION_ARGS)
    3465             : {
    3466         262 :     char       *str = PG_GETARG_CSTRING(0);
    3467             :     POLYGON    *poly;
    3468             :     int         npts;
    3469             :     int         size;
    3470             :     int         base_size;
    3471             :     bool        isopen;
    3472             : 
    3473         262 :     if ((npts = pair_count(str, ',')) <= 0)
    3474          16 :         ereport(ERROR,
    3475             :                 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
    3476             :                  errmsg("invalid input syntax for type %s: \"%s\"",
    3477             :                         "polygon", str)));
    3478             : 
    3479         246 :     base_size = sizeof(poly->p[0]) * npts;
    3480         246 :     size = offsetof(POLYGON, p) + base_size;
    3481             : 
    3482             :     /* Check for integer overflow */
    3483         246 :     if (base_size / npts != sizeof(poly->p[0]) || size <= base_size)
    3484           0 :         ereport(ERROR,
    3485             :                 (errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED),
    3486             :                  errmsg("too many points requested")));
    3487             : 
    3488         246 :     poly = (POLYGON *) palloc0(size);   /* zero any holes */
    3489             : 
    3490         246 :     SET_VARSIZE(poly, size);
    3491         246 :     poly->npts = npts;
    3492             : 
    3493         246 :     path_decode(str, false, npts, &(poly->p[0]), &isopen, NULL, "polygon", str);
    3494             : 
    3495         242 :     make_bound_box(poly);
    3496             : 
    3497         242 :     PG_RETURN_POLYGON_P(poly);
    3498             : }
    3499             : 
    3500             : /*---------------------------------------------------------------
    3501             :  * poly_out - convert internal POLYGON representation to the
    3502             :  *            character string format "((f8,f8),...,(f8,f8))"
    3503             :  *---------------------------------------------------------------*/
    3504             : Datum
    3505       47242 : poly_out(PG_FUNCTION_ARGS)
    3506             : {
    3507       47242 :     POLYGON    *poly = PG_GETARG_POLYGON_P(0);
    3508             : 
    3509       47242 :     PG_RETURN_CSTRING(path_encode(PATH_CLOSED, poly->npts, poly->p));
    3510             : }
    3511             : 
    3512             : /*
    3513             :  *      poly_recv           - converts external binary format to polygon
    3514             :  *
    3515             :  * External representation is int32 number of points, and the points.
    3516             :  * We recompute the bounding box on read, instead of trusting it to
    3517             :  * be valid.  (Checking it would take just as long, so may as well
    3518             :  * omit it from external representation.)
    3519             :  */
    3520             : Datum
    3521           0 : poly_recv(PG_FUNCTION_ARGS)
    3522             : {
    3523           0 :     StringInfo  buf = (StringInfo) PG_GETARG_POINTER(0);
    3524             :     POLYGON    *poly;
    3525             :     int32       npts;
    3526             :     int32       i;
    3527             :     int         size;
    3528             : 
    3529           0 :     npts = pq_getmsgint(buf, sizeof(int32));
    3530           0 :     if (npts <= 0 || npts >= (int32) ((INT_MAX - offsetof(POLYGON, p)) / sizeof(Point)))
    3531           0 :         ereport(ERROR,
    3532             :                 (errcode(ERRCODE_INVALID_BINARY_REPRESENTATION),
    3533             :                  errmsg("invalid number of points in external \"polygon\" value")));
    3534             : 
    3535           0 :     size = offsetof(POLYGON, p) + sizeof(poly->p[0]) * npts;
    3536           0 :     poly = (POLYGON *) palloc0(size);   /* zero any holes */
    3537             : 
    3538           0 :     SET_VARSIZE(poly, size);
    3539           0 :     poly->npts = npts;
    3540             : 
    3541           0 :     for (i = 0; i < npts; i++)
    3542             :     {
    3543           0 :         poly->p[i].x = pq_getmsgfloat8(buf);
    3544           0 :         poly->p[i].y = pq_getmsgfloat8(buf);
    3545             :     }
    3546             : 
    3547           0 :     make_bound_box(poly);
    3548             : 
    3549           0 :     PG_RETURN_POLYGON_P(poly);
    3550             : }
    3551             : 
    3552             : /*
    3553             :  *      poly_send           - converts polygon to binary format
    3554             :  */
    3555             : Datum
    3556           0 : poly_send(PG_FUNCTION_ARGS)
    3557             : {
    3558           0 :     POLYGON    *poly = PG_GETARG_POLYGON_P(0);
    3559             :     StringInfoData buf;
    3560             :     int32       i;
    3561             : 
    3562           0 :     pq_begintypsend(&buf);
    3563           0 :     pq_sendint32(&buf, poly->npts);
    3564           0 :     for (i = 0; i < poly->npts; i++)
    3565             :     {
    3566           0 :         pq_sendfloat8(&buf, poly->p[i].x);
    3567           0 :         pq_sendfloat8(&buf, poly->p[i].y);
    3568             :     }
    3569           0 :     PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
    3570             : }
    3571             : 
    3572             : 
    3573             : /*-------------------------------------------------------
    3574             :  * Is polygon A strictly left of polygon B? i.e. is
    3575             :  * the right most point of A left of the left most point
    3576             :  * of B?
    3577             :  *-------------------------------------------------------*/
    3578             : Datum
    3579         196 : poly_left(PG_FUNCTION_ARGS)
    3580             : {
    3581         196 :     POLYGON    *polya = PG_GETARG_POLYGON_P(0);
    3582         196 :     POLYGON    *polyb = PG_GETARG_POLYGON_P(1);
    3583             :     bool        result;
    3584             : 
    3585         196 :     result = polya->boundbox.high.x < polyb->boundbox.low.x;
    3586             : 
    3587             :     /*
    3588             :      * Avoid leaking memory for toasted inputs ... needed for rtree indexes
    3589             :      */
    3590         196 :     PG_FREE_IF_COPY(polya, 0);
    3591         196 :     PG_FREE_IF_COPY(polyb, 1);
    3592             : 
    3593         196 :     PG_RETURN_BOOL(result);
    3594             : }
    3595             : 
    3596             : /*-------------------------------------------------------
    3597             :  * Is polygon A overlapping or left of polygon B? i.e. is
    3598             :  * the right most point of A at or left of the right most point
    3599             :  * of B?
    3600             :  *-------------------------------------------------------*/
    3601             : Datum
    3602         196 : poly_overleft(PG_FUNCTION_ARGS)
    3603             : {
    3604         196 :     POLYGON    *polya = PG_GETARG_POLYGON_P(0);
    3605         196 :     POLYGON    *polyb = PG_GETARG_POLYGON_P(1);
    3606             :     bool        result;
    3607             : 
    3608         196 :     result = polya->boundbox.high.x <= polyb->boundbox.high.x;
    3609             : 
    3610             :     /*
    3611             :      * Avoid leaking memory for toasted inputs ... needed for rtree indexes
    3612             :      */
    3613         196 :     PG_FREE_IF_COPY(polya, 0);
    3614         196 :     PG_FREE_IF_COPY(polyb, 1);
    3615             : 
    3616         196 :     PG_RETURN_BOOL(result);
    3617             : }
    3618             : 
    3619             : /*-------------------------------------------------------
    3620             :  * Is polygon A strictly right of polygon B? i.e. is
    3621             :  * the left most point of A right of the right most point
    3622             :  * of B?
    3623             :  *-------------------------------------------------------*/
    3624             : Datum
    3625         196 : poly_right(PG_FUNCTION_ARGS)
    3626             : {
    3627         196 :     POLYGON    *polya = PG_GETARG_POLYGON_P(0);
    3628         196 :     POLYGON    *polyb = PG_GETARG_POLYGON_P(1);
    3629             :     bool        result;
    3630             : 
    3631         196 :     result = polya->boundbox.low.x > polyb->boundbox.high.x;
    3632             : 
    3633             :     /*
    3634             :      * Avoid leaking memory for toasted inputs ... needed for rtree indexes
    3635             :      */
    3636         196 :     PG_FREE_IF_COPY(polya, 0);
    3637         196 :     PG_FREE_IF_COPY(polyb, 1);
    3638             : 
    3639         196 :     PG_RETURN_BOOL(result);
    3640             : }
    3641             : 
    3642             : /*-------------------------------------------------------
    3643             :  * Is polygon A overlapping or right of polygon B? i.e. is
    3644             :  * the left most point of A at or right of the left most point
    3645             :  * of B?
    3646             :  *-------------------------------------------------------*/
    3647             : Datum
    3648         196 : poly_overright(PG_FUNCTION_ARGS)
    3649             : {
    3650         196 :     POLYGON    *polya = PG_GETARG_POLYGON_P(0);
    3651         196 :     POLYGON    *polyb = PG_GETARG_POLYGON_P(1);
    3652             :     bool        result;
    3653             : 
    3654         196 :     result = polya->boundbox.low.x >= polyb->boundbox.low.x;
    3655             : 
    3656             :     /*
    3657             :      * Avoid leaking memory for toasted inputs ... needed for rtree indexes
    3658             :      */
    3659         196 :     PG_FREE_IF_COPY(polya, 0);
    3660         196 :     PG_FREE_IF_COPY(polyb, 1);
    3661             : 
    3662         196 :     PG_RETURN_BOOL(result);
    3663             : }
    3664             : 
    3665             : /*-------------------------------------------------------
    3666             :  * Is polygon A strictly below polygon B? i.e. is
    3667             :  * the upper most point of A below the lower most point
    3668             :  * of B?
    3669             :  *-------------------------------------------------------*/
    3670             : Datum
    3671         196 : poly_below(PG_FUNCTION_ARGS)
    3672             : {
    3673         196 :     POLYGON    *polya = PG_GETARG_POLYGON_P(0);
    3674         196 :     POLYGON    *polyb = PG_GETARG_POLYGON_P(1);
    3675             :     bool        result;
    3676             : 
    3677         196 :     result = polya->boundbox.high.y < polyb->boundbox.low.y;
    3678             : 
    3679             :     /*
    3680             :      * Avoid leaking memory for toasted inputs ... needed for rtree indexes
    3681             :      */
    3682         196 :     PG_FREE_IF_COPY(polya, 0);
    3683         196 :     PG_FREE_IF_COPY(polyb, 1);
    3684             : 
    3685         196 :     PG_RETURN_BOOL(result);
    3686             : }
    3687             : 
    3688             : /*-------------------------------------------------------
    3689             :  * Is polygon A overlapping or below polygon B? i.e. is
    3690             :  * the upper most point of A at or below the upper most point
    3691             :  * of B?
    3692             :  *-------------------------------------------------------*/
    3693             : Datum
    3694         196 : poly_overbelow(PG_FUNCTION_ARGS)
    3695             : {
    3696         196 :     POLYGON    *polya = PG_GETARG_POLYGON_P(0);
    3697         196 :     POLYGON    *polyb = PG_GETARG_POLYGON_P(1);
    3698             :     bool        result;
    3699             : 
    3700         196 :     result = polya->boundbox.high.y <= polyb->boundbox.high.y;
    3701             : 
    3702             :     /*
    3703             :      * Avoid leaking memory for toasted inputs ... needed for rtree indexes
    3704             :      */
    3705         196 :     PG_FREE_IF_COPY(polya, 0);
    3706         196 :     PG_FREE_IF_COPY(polyb, 1);
    3707             : 
    3708         196 :     PG_RETURN_BOOL(result);
    3709             : }
    3710             : 
    3711             : /*-------------------------------------------------------
    3712             :  * Is polygon A strictly above polygon B? i.e. is
    3713             :  * the lower most point of A above the upper most point
    3714             :  * of B?
    3715             :  *-------------------------------------------------------*/
    3716             : Datum
    3717         196 : poly_above(PG_FUNCTION_ARGS)
    3718             : {
    3719         196 :     POLYGON    *polya = PG_GETARG_POLYGON_P(0);
    3720         196 :     POLYGON    *polyb = PG_GETARG_POLYGON_P(1);
    3721             :     bool        result;
    3722             : 
    3723         196 :     result = polya->boundbox.low.y > polyb->boundbox.high.y;
    3724             : 
    3725             :     /*
    3726             :      * Avoid leaking memory for toasted inputs ... needed for rtree indexes
    3727             :      */
    3728         196 :     PG_FREE_IF_COPY(polya, 0);
    3729         196 :     PG_FREE_IF_COPY(polyb, 1);
    3730             : 
    3731         196 :     PG_RETURN_BOOL(result);
    3732             : }
    3733             : 
    3734             : /*-------------------------------------------------------
    3735             :  * Is polygon A overlapping or above polygon B? i.e. is
    3736             :  * the lower most point of A at or above the lower most point
    3737             :  * of B?
    3738             :  *-------------------------------------------------------*/
    3739             : Datum
    3740         196 : poly_overabove(PG_FUNCTION_ARGS)
    3741             : {
    3742         196 :     POLYGON    *polya = PG_GETARG_POLYGON_P(0);
    3743         196 :     POLYGON    *polyb = PG_GETARG_POLYGON_P(1);
    3744             :     bool        result;
    3745             : 
    3746         196 :     result = polya->boundbox.low.y >= polyb->boundbox.low.y;
    3747             : 
    3748             :     /*
    3749             :      * Avoid leaking memory for toasted inputs ... needed for rtree indexes
    3750             :      */
    3751         196 :     PG_FREE_IF_COPY(polya, 0);
    3752         196 :     PG_FREE_IF_COPY(polyb, 1);
    3753             : 
    3754         196 :     PG_RETURN_BOOL(result);
    3755             : }
    3756             : 
    3757             : 
    3758             : /*-------------------------------------------------------
    3759             :  * Is polygon A the same as polygon B? i.e. are all the
    3760             :  * points the same?
    3761             :  * Check all points for matches in both forward and reverse
    3762             :  *  direction since polygons are non-directional and are
    3763             :  *  closed shapes.
    3764             :  *-------------------------------------------------------*/
    3765             : Datum
    3766        4198 : poly_same(PG_FUNCTION_ARGS)
    3767             : {
    3768        4198 :     POLYGON    *polya = PG_GETARG_POLYGON_P(0);
    3769        4198 :     POLYGON    *polyb = PG_GETARG_POLYGON_P(1);
    3770             :     bool        result;
    3771             : 
    3772        4198 :     if (polya->npts != polyb->npts)
    3773         136 :         result = false;
    3774             :     else
    3775        4062 :         result = plist_same(polya->npts, polya->p, polyb->p);
    3776             : 
    3777             :     /*
    3778             :      * Avoid leaking memory for toasted inputs ... needed for rtree indexes
    3779             :      */
    3780        4198 :     PG_FREE_IF_COPY(polya, 0);
    3781        4198 :     PG_FREE_IF_COPY(polyb, 1);
    3782             : 
    3783        4198 :     PG_RETURN_BOOL(result);
    3784             : }
    3785             : 
    3786             : /*-----------------------------------------------------------------
    3787             :  * Determine if polygon A overlaps polygon B
    3788             :  *-----------------------------------------------------------------*/
    3789             : Datum
    3790       19416 : poly_overlap(PG_FUNCTION_ARGS)
    3791             : {
    3792       19416 :     POLYGON    *polya = PG_GETARG_POLYGON_P(0);
    3793       19416 :     POLYGON    *polyb = PG_GETARG_POLYGON_P(1);
    3794             :     bool        result;
    3795             : 
    3796             :     Assert(polya->npts > 0 && polyb->npts > 0);
    3797             : 
    3798             :     /* Quick check by bounding box */
    3799       19416 :     result = box_ov(&polya->boundbox, &polyb->boundbox);
    3800             : 
    3801             :     /*
    3802             :      * Brute-force algorithm - try to find intersected edges, if so then
    3803             :      * polygons are overlapped else check is one polygon inside other or not
    3804             :      * by testing single point of them.
    3805             :      */
    3806       19416 :     if (result)
    3807             :     {
    3808             :         int         ia,
    3809             :                     ib;
    3810             :         LSEG        sa,
    3811             :                     sb;
    3812             : 
    3813             :         /* Init first of polya's edge with last point */
    3814        6952 :         sa.p[0] = polya->p[polya->npts - 1];
    3815        6952 :         result = false;
    3816             : 
    3817       84004 :         for (ia = 0; ia < polya->npts && !result; ia++)
    3818             :         {
    3819             :             /* Second point of polya's edge is a current one */
    3820       77052 :             sa.p[1] = polya->p[ia];
    3821             : 
    3822             :             /* Init first of polyb's edge with last point */
    3823       77052 :             sb.p[0] = polyb->p[polyb->npts - 1];
    3824             : 
    3825      383796 :             for (ib = 0; ib < polyb->npts && !result; ib++)
    3826             :             {
    3827      306744 :                 sb.p[1] = polyb->p[ib];
    3828      306744 :                 result = lseg_interpt_lseg(NULL, &sa, &sb);
    3829      306744 :                 sb.p[0] = sb.p[1];
    3830             :             }
    3831             : 
    3832             :             /*
    3833             :              * move current endpoint to the first point of next edge
    3834             :              */
    3835       77052 :             sa.p[0] = sa.p[1];
    3836             :         }
    3837             : 
    3838        6952 :         if (!result)
    3839             :         {
    3840        9196 :             result = (point_inside(polya->p, polyb->npts, polyb->p) ||
    3841        2928 :                       point_inside(polyb->p, polya->npts, polya->p));
    3842             :         }
    3843             :     }
    3844             : 
    3845             :     /*
    3846             :      * Avoid leaking memory for toasted inputs ... needed for rtree indexes
    3847             :      */
    3848       19416 :     PG_FREE_IF_COPY(polya, 0);
    3849       19416 :     PG_FREE_IF_COPY(polyb, 1);
    3850             : 
    3851       19416 :     PG_RETURN_BOOL(result);
    3852             : }
    3853             : 
    3854             : /*
    3855             :  * Tests special kind of segment for in/out of polygon.
    3856             :  * Special kind means:
    3857             :  *  - point a should be on segment s
    3858             :  *  - segment (a,b) should not be contained by s
    3859             :  * Returns true if:
    3860             :  *  - segment (a,b) is collinear to s and (a,b) is in polygon
    3861             :  *  - segment (a,b) s not collinear to s. Note: that doesn't
    3862             :  *    mean that segment is in polygon!
    3863             :  */
    3864             : 
    3865             : static bool
    3866         404 : touched_lseg_inside_poly(Point *a, Point *b, LSEG *s, POLYGON *poly, int start)
    3867             : {
    3868             :     /* point a is on s, b is not */
    3869             :     LSEG        t;
    3870             : 
    3871         404 :     t.p[0] = *a;
    3872         404 :     t.p[1] = *b;
    3873             : 
    3874         404 :     if (point_eq_point(a, s->p))
    3875             :     {
    3876          48 :         if (lseg_contain_point(&t, s->p + 1))
    3877           0 :             return lseg_inside_poly(b, s->p + 1, poly, start);
    3878             :     }
    3879         356 :     else if (point_eq_point(a, s->p + 1))
    3880             :     {
    3881         104 :         if (lseg_contain_point(&t, s->p))
    3882           0 :             return lseg_inside_poly(b, s->p, poly, start);
    3883             :     }
    3884         252 :     else if (lseg_contain_point(&t, s->p))
    3885             :     {
    3886           0 :         return lseg_inside_poly(b, s->p, poly, start);
    3887             :     }
    3888         252 :     else if (lseg_contain_point(&t, s->p + 1))
    3889             :     {
    3890           0 :         return lseg_inside_poly(b, s->p + 1, poly, start);
    3891             :     }
    3892             : 
    3893         404 :     return true;                /* may be not true, but that will check later */
    3894             : }
    3895             : 
    3896             : /*
    3897             :  * Returns true if segment (a,b) is in polygon, option
    3898             :  * start is used for optimization - function checks
    3899             :  * polygon's edges starting from start
    3900             :  */
    3901             : static bool
    3902      132184 : lseg_inside_poly(Point *a, Point *b, POLYGON *poly, int start)
    3903             : {
    3904             :     LSEG        s,
    3905             :                 t;
    3906             :     int         i;
    3907      132184 :     bool        res = true,
    3908      132184 :                 intersection = false;
    3909             : 
    3910      132184 :     t.p[0] = *a;
    3911      132184 :     t.p[1] = *b;
    3912      132184 :     s.p[0] = poly->p[(start == 0) ? (poly->npts - 1) : (start - 1)];
    3913             : 
    3914      653796 :     for (i = start; i < poly->npts && res; i++)
    3915             :     {
    3916             :         Point       interpt;
    3917             : 
    3918      521924 :         CHECK_FOR_INTERRUPTS();
    3919             : 
    3920      521924 :         s.p[1] = poly->p[i];
    3921             : 
    3922      521924 :         if (lseg_contain_point(&s, t.p))
    3923             :         {
    3924         524 :             if (lseg_contain_point(&s, t.p + 1))
    3925         312 :                 return true;    /* t is contained by s */
    3926             : 
    3927             :             /* Y-cross */
    3928         212 :             res = touched_lseg_inside_poly(t.p, t.p + 1, &s, poly, i + 1);
    3929             :         }
    3930      521400 :         else if (lseg_contain_point(&s, t.p + 1))
    3931             :         {
    3932             :             /* Y-cross */
    3933         192 :             res = touched_lseg_inside_poly(t.p + 1, t.p, &s, poly, i + 1);
    3934             :         }
    3935      521208 :         else if (lseg_interpt_lseg(&interpt, &t, &s))
    3936             :         {
    3937             :             /*
    3938             :              * segments are X-crossing, go to check each subsegment
    3939             :              */
    3940             : 
    3941         900 :             intersection = true;
    3942         900 :             res = lseg_inside_poly(t.p, &interpt, poly, i + 1);
    3943         900 :             if (res)
    3944         772 :                 res = lseg_inside_poly(t.p + 1, &interpt, poly, i + 1);
    3945             :         }
    3946             : 
    3947      521612 :         s.p[0] = s.p[1];
    3948             :     }
    3949             : 
    3950      131872 :     if (res && !intersection)
    3951             :     {
    3952             :         Point       p;
    3953             : 
    3954             :         /*
    3955             :          * if X-intersection wasn't found  then check central point of tested
    3956             :          * segment. In opposite case we already check all subsegments
    3957             :          */
    3958      130972 :         p.x = float8_div(float8_pl(t.p[0].x, t.p[1].x), 2.0);
    3959      130972 :         p.y = float8_div(float8_pl(t.p[0].y, t.p[1].y), 2.0);
    3960             : 
    3961      130972 :         res = point_inside(&p, poly->npts, poly->p);
    3962             :     }
    3963             : 
    3964      131872 :     return res;
    3965             : }
    3966             : 
    3967             : /*
    3968             :  * Check whether the first polygon contains the second
    3969             :  */
    3970             : static bool
    3971       56604 : poly_contain_poly(POLYGON *contains_poly, POLYGON *contained_poly)
    3972             : {
    3973             :     int         i;
    3974             :     LSEG        s;
    3975             : 
    3976             :     Assert(contains_poly->npts > 0 && contained_poly->npts > 0);
    3977             : 
    3978             :     /*
    3979             :      * Quick check to see if contained's bounding box is contained in
    3980             :      * contains' bb.
    3981             :      */
    3982       56604 :     if (!box_contain_box(&contains_poly->boundbox, &contained_poly->boundbox))
    3983       38212 :         return false;
    3984             : 
    3985       18392 :     s.p[0] = contained_poly->p[contained_poly->npts - 1];
    3986             : 
    3987      140592 :     for (i = 0; i < contained_poly->npts; i++)
    3988             :     {
    3989      130512 :         s.p[1] = contained_poly->p[i];
    3990      130512 :         if (!lseg_inside_poly(s.p, s.p + 1, contains_poly, 0))
    3991        8312 :             return false;
    3992      122200 :         s.p[0] = s.p[1];
    3993             :     }
    3994             : 
    3995       10080 :     return true;
    3996             : }
    3997             : 
    3998             : Datum
    3999         240 : poly_contain(PG_FUNCTION_ARGS)
    4000             : {
    4001         240 :     POLYGON    *polya = PG_GETARG_POLYGON_P(0);
    4002         240 :     POLYGON    *polyb = PG_GETARG_POLYGON_P(1);
    4003             :     bool        result;
    4004             : 
    4005         240 :     result = poly_contain_poly(polya, polyb);
    4006             : 
    4007             :     /*
    4008             :      * Avoid leaking memory for toasted inputs ... needed for rtree indexes
    4009             :      */
    4010         240 :     PG_FREE_IF_COPY(polya, 0);
    4011         240 :     PG_FREE_IF_COPY(polyb, 1);
    4012             : 
    4013         240 :     PG_RETURN_BOOL(result);
    4014             : }
    4015             : 
    4016             : 
    4017             : /*-----------------------------------------------------------------
    4018             :  * Determine if polygon A is contained by polygon B
    4019             :  *-----------------------------------------------------------------*/
    4020             : Datum
    4021       56364 : poly_contained(PG_FUNCTION_ARGS)
    4022             : {
    4023       56364 :     POLYGON    *polya = PG_GETARG_POLYGON_P(0);
    4024       56364 :     POLYGON    *polyb = PG_GETARG_POLYGON_P(1);
    4025             :     bool        result;
    4026             : 
    4027             :     /* Just switch the arguments and pass it off to poly_contain */
    4028       56364 :     result = poly_contain_poly(polyb, polya);
    4029             : 
    4030             :     /*
    4031             :      * Avoid leaking memory for toasted inputs ... needed for rtree indexes
    4032             :      */
    4033       56364 :     PG_FREE_IF_COPY(polya, 0);
    4034       56364 :     PG_FREE_IF_COPY(polyb, 1);
    4035             : 
    4036       56364 :     PG_RETURN_BOOL(result);
    4037             : }
    4038             : 
    4039             : 
    4040             : Datum
    4041         268 : poly_contain_pt(PG_FUNCTION_ARGS)
    4042             : {
    4043         268 :     POLYGON    *poly = PG_GETARG_POLYGON_P(0);
    4044         268 :     Point      *p = PG_GETARG_POINT_P(1);
    4045             : 
    4046         268 :     PG_RETURN_BOOL(point_inside(p, poly->npts, poly->p) != 0);
    4047             : }
    4048             : 
    4049             : Datum
    4050         288 : pt_contained_poly(PG_FUNCTION_ARGS)
    4051             : {
    4052         288 :     Point      *p = PG_GETARG_POINT_P(0);
    4053         288 :     POLYGON    *poly = PG_GETARG_POLYGON_P(1);
    4054             : 
    4055         288 :     PG_RETURN_BOOL(point_inside(p, poly->npts, poly->p) != 0);
    4056             : }
    4057             : 
    4058             : 
    4059             : Datum
    4060           4 : poly_distance(PG_FUNCTION_ARGS)
    4061             : {
    4062             : #ifdef NOT_USED
    4063             :     POLYGON    *polya = PG_GETARG_POLYGON_P(0);
    4064             :     POLYGON    *polyb = PG_GETARG_POLYGON_P(1);
    4065             : #endif
    4066             : 
    4067           4 :     ereport(ERROR,
    4068             :             (errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
    4069             :              errmsg("function \"poly_distance\" not implemented")));
    4070             : 
    4071             :     PG_RETURN_NULL();
    4072             : }
    4073             : 
    4074             : 
    4075             : /***********************************************************************
    4076             :  **
    4077             :  **     Routines for 2D points.
    4078             :  **
    4079             :  ***********************************************************************/
    4080             : 
    4081             : Datum
    4082      723650 : construct_point(PG_FUNCTION_ARGS)
    4083             : {
    4084      723650 :     float8      x = PG_GETARG_FLOAT8(0);
    4085      723650 :     float8      y = PG_GETARG_FLOAT8(1);
    4086             :     Point      *result;
    4087             : 
    4088      723650 :     result = (Point *) palloc(sizeof(Point));
    4089             : 
    4090      723650 :     point_construct(result, x, y);
    4091             : 
    4092      723650 :     PG_RETURN_POINT_P(result);
    4093             : }
    4094             : 
    4095             : 
    4096             : static inline void
    4097        1824 : point_add_point(Point *result, Point *pt1, Point *pt2)
    4098             : {
    4099        1824 :     point_construct(result,
    4100             :                     float8_pl(pt1->x, pt2->x),
    4101             :                     float8_pl(pt1->y, pt2->y));
    4102        1824 : }
    4103             : 
    4104             : Datum
    4105         324 : point_add(PG_FUNCTION_ARGS)
    4106             : {
    4107         324 :     Point      *p1 = PG_GETARG_POINT_P(0);
    4108         324 :     Point      *p2 = PG_GETARG_POINT_P(1);
    4109             :     Point      *result;
    4110             : 
    4111         324 :     result = (Point *) palloc(sizeof(Point));
    4112             : 
    4113         324 :     point_add_point(result, p1, p2);
    4114             : 
    4115         324 :     PG_RETURN_POINT_P(result);
    4116             : }
    4117             : 
    4118             : 
    4119             : static inline void
    4120        1656 : point_sub_point(Point *result, Point *pt1, Point *pt2)
    4121             : {
    4122        1656 :     point_construct(result,
    4123             :                     float8_mi(pt1->x, pt2->x),
    4124             :                     float8_mi(pt1->y, pt2->y));
    4125        1656 : }
    4126             : 
    4127             : Datum
    4128         324 : point_sub(PG_FUNCTION_ARGS)
    4129             : {
    4130         324 :     Point      *p1 = PG_GETARG_POINT_P(0);
    4131         324 :     Point      *p2 = PG_GETARG_POINT_P(1);
    4132             :     Point      *result;
    4133             : 
    4134         324 :     result = (Point *) palloc(sizeof(Point));
    4135             : 
    4136         324 :     point_sub_point(result, p1, p2);
    4137             : 
    4138         324 :     PG_RETURN_POINT_P(result);
    4139             : }
    4140             : 
    4141             : 
    4142             : static inline void
    4143        1324 : point_mul_point(Point *result, Point *pt1, Point *pt2)
    4144             : {
    4145        1324 :     point_construct(result,
    4146             :                     float8_mi(float8_mul(pt1->x, pt2->x),
    4147             :                               float8_mul(pt1->y, pt2->y)),
    4148             :                     float8_pl(float8_mul(pt1->x, pt2->y),
    4149             :                               float8_mul(pt1->y, pt2->x)));
    4150        1320 : }
    4151             : 
    4152             : Datum
    4153         192 : point_mul(PG_FUNCTION_ARGS)
    4154             : {
    4155         192 :     Point      *p1 = PG_GETARG_POINT_P(0);
    4156         192 :     Point      *p2 = PG_GETARG_POINT_P(1);
    4157             :     Point      *result;
    4158             : 
    4159         192 :     result = (Point *) palloc(sizeof(Point));
    4160             : 
    4161         192 :     point_mul_point(result, p1, p2);
    4162             : 
    4163         188 :     PG_RETURN_POINT_P(result);
    4164             : }
    4165             : 
    4166             : 
    4167             : static inline void
    4168         388 : point_div_point(Point *result, Point *pt1, Point *pt2)
    4169             : {
    4170             :     float8      div;
    4171             : 
    4172         388 :     div = float8_pl(float8_mul(pt2->x, pt2->x), float8_mul(pt2->y, pt2->y));
    4173             : 
    4174         380 :     point_construct(result,
    4175             :                     float8_div(float8_pl(float8_mul(pt1->x, pt2->x),
    4176             :                                          float8_mul(pt1->y, pt2->y)), div),
    4177             :                     float8_div(float8_mi(float8_mul(pt1->y, pt2->x),
    4178             :                                          float8_mul(pt1->x, pt2->y)), div));
    4179         368 : }
    4180             : 
    4181             : Datum
    4182          80 : point_div(PG_FUNCTION_ARGS)
    4183             : {
    4184          80 :     Point      *p1 = PG_GETARG_POINT_P(0);
    4185          80 :     Point      *p2 = PG_GETARG_POINT_P(1);
    4186             :     Point      *result;
    4187             : 
    4188          80 :     result = (Point *) palloc(sizeof(Point));
    4189             : 
    4190          80 :     point_div_point(result, p1, p2);
    4191             : 
    4192          72 :     PG_RETURN_POINT_P(result);
    4193             : }
    4194             : 
    4195             : 
    4196             : /***********************************************************************
    4197             :  **
    4198             :  **     Routines for 2D boxes.
    4199             :  **
    4200             :  ***********************************************************************/
    4201             : 
    4202             : Datum
    4203      160968 : points_box(PG_FUNCTION_ARGS)
    4204             : {
    4205      160968 :     Point      *p1 = PG_GETARG_POINT_P(0);
    4206      160968 :     Point      *p2 = PG_GETARG_POINT_P(1);
    4207             :     BOX        *result;
    4208             : 
    4209      160968 :     result = (BOX *) palloc(sizeof(BOX));
    4210             : 
    4211      160968 :     box_construct(result, p1, p2);
    4212             : 
    4213      160968 :     PG_RETURN_BOX_P(result);
    4214             : }
    4215             : 
    4216             : Datum
    4217         180 : box_add(PG_FUNCTION_ARGS)
    4218             : {
    4219         180 :     BOX        *box = PG_GETARG_BOX_P(0);
    4220         180 :     Point      *p = PG_GETARG_POINT_P(1);
    4221             :     BOX        *result;
    4222             : 
    4223         180 :     result = (BOX *) palloc(sizeof(BOX));
    4224             : 
    4225         180 :     point_add_point(&result->high, &box->high, p);
    4226         180 :     point_add_point(&result->low, &box->low, p);
    4227             : 
    4228         180 :     PG_RETURN_BOX_P(result);
    4229             : }
    4230             : 
    4231             : Datum
    4232         180 : box_sub(PG_FUNCTION_ARGS)
    4233             : {
    4234         180 :     BOX        *box = PG_GETARG_BOX_P(0);
    4235         180 :     Point      *p = PG_GETARG_POINT_P(1);
    4236             :     BOX        *result;
    4237             : 
    4238         180 :     result = (BOX *) palloc(sizeof(BOX));
    4239             : 
    4240         180 :     point_sub_point(&result->high, &box->high, p);
    4241         180 :     point_sub_point(&result->low, &box->low, p);
    4242             : 
    4243         180 :     PG_RETURN_BOX_P(result);
    4244             : }
    4245             : 
    4246             : Datum
    4247          80 : box_mul(PG_FUNCTION_ARGS)
    4248             : {
    4249          80 :     BOX        *box = PG_GETARG_BOX_P(0);
    4250          80 :     Point      *p = PG_GETARG_POINT_P(1);
    4251             :     BOX        *result;
    4252             :     Point       high,
    4253             :                 low;
    4254             : 
    4255          80 :     result = (BOX *) palloc(sizeof(BOX));
    4256             : 
    4257          80 :     point_mul_point(&high, &box->high, p);
    4258          80 :     point_mul_point(&low, &box->low, p);
    4259             : 
    4260          80 :     box_construct(result, &high, &low);
    4261             : 
    4262          80 :     PG_RETURN_BOX_P(result);
    4263             : }
    4264             : 
    4265             : Datum
    4266          40 : box_div(PG_FUNCTION_ARGS)
    4267             : {
    4268          40 :     BOX        *box = PG_GETARG_BOX_P(0);
    4269          40 :     Point      *p = PG_GETARG_POINT_P(1);
    4270             :     BOX        *result;
    4271             :     Point       high,
    4272             :                 low;
    4273             : 
    4274          40 :     result = (BOX *) palloc(sizeof(BOX));
    4275             : 
    4276          40 :     point_div_point(&high, &box->high, p);
    4277          40 :     point_div_point(&low, &box->low, p);
    4278             : 
    4279          40 :     box_construct(result, &high, &low);
    4280             : 
    4281          40 :     PG_RETURN_BOX_P(result);
    4282             : }
    4283             : 
    4284             : /*
    4285             :  * Convert point to empty box
    4286             :  */
    4287             : Datum
    4288          36 : point_box(PG_FUNCTION_ARGS)
    4289             : {
    4290          36 :     Point      *pt = PG_GETARG_POINT_P(0);
    4291             :     BOX        *box;
    4292             : 
    4293          36 :     box = (BOX *) palloc(sizeof(BOX));
    4294             : 
    4295          36 :     box->high.x = pt->x;
    4296          36 :     box->low.x = pt->x;
    4297          36 :     box->high.y = pt->y;
    4298          36 :     box->low.y = pt->y;
    4299             : 
    4300          36 :     PG_RETURN_BOX_P(box);
    4301             : }
    4302             : 
    4303             : /*
    4304             :  * Smallest bounding box that includes both of the given boxes
    4305             :  */
    4306             : Datum
    4307         100 : boxes_bound_box(PG_FUNCTION_ARGS)
    4308             : {
    4309         100 :     BOX        *box1 = PG_GETARG_BOX_P(0),
    4310         100 :                *box2 = PG_GETARG_BOX_P(1),
    4311             :                *container;
    4312             : 
    4313         100 :     container = (BOX *) palloc(sizeof(BOX));
    4314             : 
    4315         100 :     container->high.x = float8_max(box1->high.x, box2->high.x);
    4316         100 :     container->low.x = float8_min(box1->low.x, box2->low.x);
    4317         100 :     container->high.y = float8_max(box1->high.y, box2->high.y);
    4318         100 :     container->low.y = float8_min(box1->low.y, box2->low.y);
    4319             : 
    4320         100 :     PG_RETURN_BOX_P(container);
    4321             : }
    4322             : 
    4323             : 
    4324             : /***********************************************************************
    4325             :  **
    4326             :  **     Routines for 2D paths.
    4327             :  **
    4328             :  ***********************************************************************/
    4329             : 
    4330             : /* path_add()
    4331             :  * Concatenate two paths (only if they are both open).
    4332             :  */
    4333             : Datum
    4334         324 : path_add(PG_FUNCTION_ARGS)
    4335             : {
    4336         324 :     PATH       *p1 = PG_GETARG_PATH_P(0);
    4337         324 :     PATH       *p2 = PG_GETARG_PATH_P(1);
    4338             :     PATH       *result;
    4339             :     int         size,
    4340             :                 base_size;
    4341             :     int         i;
    4342             : 
    4343         324 :     if (p1->closed || p2->closed)
    4344         260 :         PG_RETURN_NULL();
    4345             : 
    4346          64 :     base_size = sizeof(p1->p[0]) * (p1->npts + p2->npts);
    4347          64 :     size = offsetof(PATH, p) + base_size;
    4348             : 
    4349             :     /* Check for integer overflow */
    4350          64 :     if (base_size / sizeof(p1->p[0]) != (p1->npts + p2->npts) ||
    4351             :         size <= base_size)
    4352           0 :         ereport(ERROR,
    4353             :                 (errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED),
    4354             :                  errmsg("too many points requested")));
    4355             : 
    4356          64 :     result = (PATH *) palloc(size);
    4357             : 
    4358          64 :     SET_VARSIZE(result, size);
    4359          64 :     result->npts = (p1->npts + p2->npts);
    4360          64 :     result->closed = p1->closed;
    4361             :     /* prevent instability in unused pad bytes */
    4362          64 :     result->dummy = 0;
    4363             : 
    4364         224 :     for (i = 0; i < p1->npts; i++)
    4365             :     {
    4366         160 :         result->p[i].x = p1->p[i].x;
    4367         160 :         result->p[i].y = p1->p[i].y;
    4368             :     }
    4369         224 :     for (i = 0; i < p2->npts; i++)
    4370             :     {
    4371         160 :         result->p[i + p1->npts].x = p2->p[i].x;
    4372         160 :         result->p[i + p1->npts].y = p2->p[i].y;
    4373             :     }
    4374             : 
    4375          64 :     PG_RETURN_PATH_P(result);
    4376             : }
    4377             : 
    4378             : /* path_add_pt()
    4379             :  * Translation operators.
    4380             :  */
    4381             : Datum
    4382         324 : path_add_pt(PG_FUNCTION_ARGS)
    4383             : {
    4384         324 :     PATH       *path = PG_GETARG_PATH_P_COPY(0);
    4385         324 :     Point      *point = PG_GETARG_POINT_P(1);
    4386             :     int         i;
    4387             : 
    4388        1008 :     for (i = 0; i < path->npts; i++)
    4389         684 :         point_add_point(&path->p[i], &path->p[i], point);
    4390             : 
    4391         324 :     PG_RETURN_PATH_P(path);
    4392             : }
    4393             : 
    4394             : Datum
    4395         324 : path_sub_pt(PG_FUNCTION_ARGS)
    4396             : {
    4397         324 :     PATH       *path = PG_GETARG_PATH_P_COPY(0);
    4398         324 :     Point      *point = PG_GETARG_POINT_P(1);
    4399             :     int         i;
    4400             : 
    4401        1008 :     for (i = 0; i < path->npts; i++)
    4402         684 :         point_sub_point(&path->p[i], &path->p[i], point);
    4403             : 
    4404         324 :     PG_RETURN_PATH_P(path);
    4405             : }
    4406             : 
    4407             : /* path_mul_pt()
    4408             :  * Rotation and scaling operators.
    4409             :  */
    4410             : Datum
    4411         324 : path_mul_pt(PG_FUNCTION_ARGS)
    4412             : {
    4413         324 :     PATH       *path = PG_GETARG_PATH_P_COPY(0);
    4414         324 :     Point      *point = PG_GETARG_POINT_P(1);
    4415             :     int         i;
    4416             : 
    4417        1008 :     for (i = 0; i < path->npts; i++)
    4418         684 :         point_mul_point(&path->p[i], &path->p[i], point);
    4419             : 
    4420         324 :     PG_RETURN_PATH_P(path);
    4421             : }
    4422             : 
    4423             : Datum
    4424          76 : path_div_pt(PG_FUNCTION_ARGS)
    4425             : {
    4426          76 :     PATH       *path = PG_GETARG_PATH_P_COPY(0);
    4427          76 :     Point      *point = PG_GETARG_POINT_P(1);
    4428             :     int         i;
    4429             : 
    4430         228 :     for (i = 0; i < path->npts; i++)
    4431         156 :         point_div_point(&path->p[i], &path->p[i], point);
    4432             : 
    4433          72 :     PG_RETURN_PATH_P(path);
    4434             : }
    4435             : 
    4436             : 
    4437             : Datum
    4438           4 : path_center(PG_FUNCTION_ARGS)
    4439             : {
    4440             : #ifdef NOT_USED
    4441             :     PATH       *path = PG_GETARG_PATH_P(0);
    4442             : #endif
    4443             : 
    4444           4 :     ereport(ERROR,
    4445             :             (errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
    4446             :              errmsg("function \"path_center\" not implemented")));
    4447             : 
    4448             :     PG_RETURN_NULL();
    4449             : }
    4450             : 
    4451             : Datum
    4452          60 : path_poly(PG_FUNCTION_ARGS)
    4453             : {
    4454          60 :     PATH       *path = PG_GETARG_PATH_P(0);
    4455             :     POLYGON    *poly;
    4456             :     int         size;
    4457             :     int         i;
    4458             : 
    4459             :     /* This is not very consistent --- other similar cases return NULL ... */
    4460          60 :     if (!path->closed)
    4461           4 :         ereport(ERROR,
    4462             :                 (errcode(ERRCODE_INVALID_PARAMETER_VALUE),
    4463             :                  errmsg("open path cannot be converted to polygon")));
    4464             : 
    4465             :     /*
    4466             :      * Never overflows: the old size fit in MaxAllocSize, and the new size is
    4467             :      * just a small constant larger.
    4468             :      */
    4469          56 :     size = offsetof(POLYGON, p) + sizeof(poly->p[0]) * path->npts;
    4470          56 :     poly = (POLYGON *) palloc(size);
    4471             : 
    4472          56 :     SET_VARSIZE(poly, size);
    4473          56 :     poly->npts = path->npts;
    4474             : 
    4475         168 :     for (i = 0; i < path->npts; i++)
    4476             :     {
    4477         112 :         poly->p[i].x = path->p[i].x;
    4478         112 :         poly->p[i].y = path->p[i].y;
    4479             :     }
    4480             : 
    4481          56 :     make_bound_box(poly);
    4482             : 
    4483          56 :     PG_RETURN_POLYGON_P(poly);
    4484             : }
    4485             : 
    4486             : 
    4487             : /***********************************************************************
    4488             :  **
    4489             :  **     Routines for 2D polygons.
    4490             :  **
    4491             :  ***********************************************************************/
    4492             : 
    4493             : Datum
    4494          84 : poly_npoints(PG_FUNCTION_ARGS)
    4495             : {
    4496          84 :     POLYGON    *poly = PG_GETARG_POLYGON_P(0);
    4497             : 
    4498          84 :     PG_RETURN_INT32(poly->npts);
    4499             : }
    4500             : 
    4501             : 
    4502             : Datum
    4503          28 : poly_center(PG_FUNCTION_ARGS)
    4504             : {
    4505          28 :     POLYGON    *poly = PG_GETARG_POLYGON_P(0);
    4506             :     Point      *result;
    4507             :     CIRCLE      circle;
    4508             : 
    4509          28 :     result = (Point *) palloc(sizeof(Point));
    4510             : 
    4511          28 :     poly_to_circle(&circle, poly);
    4512          28 :     *result = circle.center;
    4513             : 
    4514          28 :     PG_RETURN_POINT_P(result);
    4515             : }
    4516             : 
    4517             : 
    4518             : Datum
    4519          28 : poly_box(PG_FUNCTION_ARGS)
    4520             : {
    4521          28 :     POLYGON    *poly = PG_GETARG_POLYGON_P(0);
    4522             :     BOX        *box;
    4523             : 
    4524          28 :     box = (BOX *) palloc(sizeof(BOX));
    4525          28 :     *box = poly->boundbox;
    4526             : 
    4527          28 :     PG_RETURN_BOX_P(box);
    4528             : }
    4529             : 
    4530             : 
    4531             : /* box_poly()
    4532             :  * Convert a box to a polygon.
    4533             :  */
    4534             : Datum
    4535       12420 : box_poly(PG_FUNCTION_ARGS)
    4536             : {
    4537       12420 :     BOX        *box = PG_GETARG_BOX_P(0);
    4538             :     POLYGON    *poly;
    4539             :     int         size;
    4540             : 
    4541             :     /* map four corners of the box to a polygon */
    4542       12420 :     size = offsetof(POLYGON, p) + sizeof(poly->p[0]) * 4;
    4543       12420 :     poly = (POLYGON *) palloc(size);
    4544             : 
    4545       12420 :     SET_VARSIZE(poly, size);
    4546       12420 :     poly->npts = 4;
    4547             : 
    4548       12420 :     poly->p[0].x = box->low.x;
    4549       12420 :     poly->p[0].y = box->low.y;
    4550       12420 :     poly->p[1].x = box->low.x;
    4551       12420 :     poly->p[1].y = box->high.y;
    4552       12420 :     poly->p[2].x = box->high.x;
    4553       12420 :     poly->p[2].y = box->high.y;
    4554       12420 :     poly->p[3].x = box->high.x;
    4555       12420 :     poly->p[3].y = box->low.y;
    4556             : 
    4557       12420 :     box_construct(&poly->boundbox, &box->high, &box->low);
    4558             : 
    4559       12420 :     PG_RETURN_POLYGON_P(poly);
    4560             : }
    4561             : 
    4562             : 
    4563             : Datum
    4564          28 : poly_path(PG_FUNCTION_ARGS)
    4565             : {
    4566          28 :     POLYGON    *poly = PG_GETARG_POLYGON_P(0);
    4567             :     PATH       *path;
    4568             :     int         size;
    4569             :     int         i;
    4570             : 
    4571             :     /*
    4572             :      * Never overflows: the old size fit in MaxAllocSize, and the new size is
    4573             :      * smaller by a small constant.
    4574             :      */
    4575          28 :     size = offsetof(PATH, p) + sizeof(path->p[0]) * poly->npts;
    4576          28 :     path = (PATH *) palloc(size);
    4577             : 
    4578          28 :     SET_VARSIZE(path, size);
    4579          28 :     path->npts = poly->npts;
    4580          28 :     path->closed = true;
    4581             :     /* prevent instability in unused pad bytes */
    4582          28 :     path->dummy = 0;
    4583             : 
    4584         112 :     for (i = 0; i < poly->npts; i++)
    4585             :     {
    4586          84 :         path->p[i].x = poly->p[i].x;
    4587          84 :         path->p[i].y = poly->p[i].y;
    4588             :     }
    4589             : 
    4590          28 :     PG_RETURN_PATH_P(path);
    4591             : }
    4592             : 
    4593             : 
    4594             : /***********************************************************************
    4595             :  **
    4596             :  **     Routines for circles.
    4597             :  **
    4598             :  ***********************************************************************/
    4599             : 
    4600             : /*----------------------------------------------------------
    4601             :  * Formatting and conversion routines.
    4602             :  *---------------------------------------------------------*/
    4603             : 
    4604             : /*      circle_in       -       convert a string to internal form.
    4605             :  *
    4606             :  *      External format: (center and radius of circle)
    4607             :  *              "((f8,f8)<f8>)"
    4608             :  *              also supports quick entry style "(f8,f8,f8)"
    4609             :  */
    4610             : Datum
    4611         252 : circle_in(PG_FUNCTION_ARGS)
    4612             : {
    4613         252 :     char       *str = PG_GETARG_CSTRING(0);
    4614         252 :     CIRCLE     *circle = (CIRCLE *) palloc(sizeof(CIRCLE));
    4615             :     char       *s,
    4616             :                *cp;
    4617         252 :     int         depth = 0;
    4618             : 
    4619         252 :     s = str;
    4620         508 :     while (isspace((unsigned char) *s))
    4621           4 :         s++;
    4622         252 :     if ((*s == LDELIM_C) || (*s == LDELIM))
    4623             :     {
    4624         244 :         depth++;
    4625         244 :         cp = (s + 1);
    4626         492 :         while (isspace((unsigned char) *cp))
    4627           4 :             cp++;
    4628         244 :         if (*cp == LDELIM)
    4629         240 :             s = cp;
    4630             :     }
    4631             : 
    4632         252 :     pair_decode(s, &circle->center.x, &circle->center.y, &s, "circle", str);
    4633             : 
    4634         244 :     if (*s == DELIM)
    4635         244 :         s++;
    4636             : 
    4637         244 :     circle->radius = single_decode(s, &s, "circle", str);
    4638             :     /* We have to accept NaN. */
    4639         244 :     if (circle->radius < 0.0)
    4640           4 :         ereport(ERROR,
    4641             :                 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
    4642             :                  errmsg("invalid input syntax for type %s: \"%s\"",
    4643             :                         "circle", str)));
    4644             : 
    4645         712 :     while (depth > 0)
    4646             :     {
    4647         468 :         if ((*s == RDELIM) || ((*s == RDELIM_C) && (depth == 1)))
    4648             :         {
    4649         232 :             depth--;
    4650         232 :             s++;
    4651         472 :             while (isspace((unsigned char) *s))
    4652           8 :                 s++;
    4653             :         }
    4654             :         else
    4655           4 :             ereport(ERROR,
    4656             :                     (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
    4657             :                      errmsg("invalid input syntax for type %s: \"%s\"",
    4658             :                             "circle", str)));
    4659             :     }
    4660             : 
    4661         236 :     if (*s != '\0')
    4662           4 :         ereport(ERROR,
    4663             :                 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
    4664             :                  errmsg("invalid input syntax for type %s: \"%s\"",
    4665             :                         "circle", str)));
    4666             : 
    4667         232 :     PG_RETURN_CIRCLE_P(circle);
    4668             : }
    4669             : 
    4670             : /*      circle_out      -       convert a circle to external form.
    4671             :  */
    4672             : Datum
    4673        5852 : circle_out(PG_FUNCTION_ARGS)
    4674             : {
    4675        5852 :     CIRCLE     *circle = PG_GETARG_CIRCLE_P(0);
    4676             :     StringInfoData str;
    4677             : 
    4678        5852 :     initStringInfo(&str);
    4679             : 
    4680        5852 :     appendStringInfoChar(&str, LDELIM_C);
    4681        5852 :     appendStringInfoChar(&str, LDELIM);
    4682        5852 :     pair_encode(circle->center.x, circle->center.y, &str);
    4683        5852 :     appendStringInfoChar(&str, RDELIM);
    4684        5852 :     appendStringInfoChar(&str, DELIM);
    4685        5852 :     single_encode(circle->radius, &str);
    4686        5852 :     appendStringInfoChar(&str, RDELIM_C);
    4687             : 
    4688        5852 :     PG_RETURN_CSTRING(str.data);
    4689             : }
    4690             : 
    4691             : /*
    4692             :  *      circle_recv         - converts external binary format to circle
    4693             :  */
    4694             : Datum
    4695           0 : circle_recv(PG_FUNCTION_ARGS)
    4696             : {
    4697           0 :     StringInfo  buf = (StringInfo) PG_GETARG_POINTER(0);
    4698             :     CIRCLE     *circle;
    4699             : 
    4700           0 :     circle = (CIRCLE *) palloc(sizeof(CIRCLE));
    4701             : 
    4702           0 :     circle->center.x = pq_getmsgfloat8(buf);
    4703           0 :     circle->center.y = pq_getmsgfloat8(buf);
    4704           0 :     circle->radius = pq_getmsgfloat8(buf);
    4705             : 
    4706             :     /* We have to accept NaN. */
    4707           0 :     if (circle->radius < 0.0)
    4708           0 :         ereport(ERROR,
    4709             :                 (errcode(ERRCODE_INVALID_BINARY_REPRESENTATION),
    4710             :                  errmsg("invalid radius in external \"circle\" value")));
    4711             : 
    4712           0 :     PG_RETURN_CIRCLE_P(circle);
    4713             : }
    4714             : 
    4715             : /*
    4716             :  *      circle_send         - converts circle to binary format
    4717             :  */
    4718             : Datum
    4719           0 : circle_send(PG_FUNCTION_ARGS)
    4720             : {
    4721           0 :     CIRCLE     *circle = PG_GETARG_CIRCLE_P(0);
    4722             :     StringInfoData buf;
    4723             : 
    4724           0 :     pq_begintypsend(&buf);
    4725           0 :     pq_sendfloat8(&buf, circle->center.x);
    4726           0 :     pq_sendfloat8(&buf, circle->center.y);
    4727           0 :     pq_sendfloat8(&buf, circle->radius);
    4728           0 :     PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
    4729             : }
    4730             : 
    4731             : 
    4732             : /*----------------------------------------------------------
    4733             :  *  Relational operators for CIRCLEs.
    4734             :  *      <, >, <=, >=, and == are based on circle area.
    4735             :  *---------------------------------------------------------*/
    4736             : 
    4737             : /*      circles identical?
    4738             :  *
    4739             :  * We consider NaNs values to be equal to each other to let those circles
    4740             :  * to be found.
    4741             :  */
    4742             : Datum
    4743         258 : circle_same(PG_FUNCTION_ARGS)
    4744             : {
    4745         258 :     CIRCLE     *circle1 = PG_GETARG_CIRCLE_P(0);
    4746         258 :     CIRCLE     *circle2 = PG_GETARG_CIRCLE_P(1);
    4747             : 
    4748         258 :     PG_RETURN_BOOL(((isnan(circle1->radius) && isnan(circle1->radius)) ||
    4749             :                     FPeq(circle1->radius, circle2->radius)) &&
    4750             :                    point_eq_point(&circle1->center, &circle2->center));
    4751             : }
    4752             : 
    4753             : /*      circle_overlap  -       does circle1 overlap circle2?
    4754             :  */
    4755             : Datum
    4756       12800 : circle_overlap(PG_FUNCTION_ARGS)
    4757             : {
    4758       12800 :     CIRCLE     *circle1 = PG_GETARG_CIRCLE_P(0);
    4759       12800 :     CIRCLE     *circle2 = PG_GETARG_CIRCLE_P(1);
    4760             : 
    4761       12800 :     PG_RETURN_BOOL(FPle(point_dt(&circle1->center, &circle2->center),
    4762             :                         float8_pl(circle1->radius, circle2->radius)));
    4763             : }
    4764             : 
    4765             : /*      circle_overleft -       is the right edge of circle1 at or left of
    4766             :  *                              the right edge of circle2?
    4767             :  */
    4768             : Datum
    4769         256 : circle_overleft(PG_FUNCTION_ARGS)
    4770             : {
    4771         256 :     CIRCLE     *circle1 = PG_GETARG_CIRCLE_P(0);
    4772         256 :     CIRCLE     *circle2 = PG_GETARG_CIRCLE_P(1);
    4773             : 
    4774         256 :     PG_RETURN_BOOL(FPle(float8_pl(circle1->center.x, circle1->radius),
    4775             :                         float8_pl(circle2->center.x, circle2->radius)));
    4776             : }
    4777             : 
    4778             : /*      circle_left     -       is circle1 strictly left of circle2?
    4779             :  */
    4780             : Datum
    4781         256 : circle_left(PG_FUNCTION_ARGS)
    4782             : {
    4783         256 :     CIRCLE     *circle1 = PG_GETARG_CIRCLE_P(0);
    4784         256 :     CIRCLE     *circle2 = PG_GETARG_CIRCLE_P(1);
    4785             : 
    4786         256 :     PG_RETURN_BOOL(FPlt(float8_pl(circle1->center.x, circle1->radius),
    4787             :                         float8_mi(circle2->center.x, circle2->radius)));
    4788             : }
    4789             : 
    4790             : /*      circle_right    -       is circle1 strictly right of circle2?
    4791             :  */
    4792             : Datum
    4793         256 : circle_right(PG_FUNCTION_ARGS)
    4794             : {
    4795         256 :     CIRCLE     *circle1 = PG_GETARG_CIRCLE_P(0);
    4796         256 :     CIRCLE     *circle2 = PG_GETARG_CIRCLE_P(1);
    4797             : 
    4798         256 :     PG_RETURN_BOOL(FPgt(float8_mi(circle1->center.x, circle1->radius),
    4799             :                         float8_pl(circle2->center.x, circle2->radius)));
    4800             : }
    4801             : 
    4802             : /*      circle_overright    -   is the left edge of circle1 at or right of
    4803             :  *                              the left edge of circle2?
    4804             :  */
    4805             : Datum
    4806         256 : circle_overright(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(FPge(float8_mi(circle1->center.x, circle1->radius),
    4812             :                         float8_mi(circle2->center.x, circle2->radius)));
    4813             : }
    4814             : 
    4815             : /*      circle_contained        -       is circle1 contained by circle2?
    4816             :  */
    4817             : Datum
    4818         256 : circle_contained(PG_FUNCTION_ARGS)
    4819             : {
    4820         256 :     CIRCLE     *circle1 = PG_GETARG_CIRCLE_P(0);
    4821         256 :     CIRCLE     *circle2 = PG_GETARG_CIRCLE_P(1);
    4822             : 
    4823         256 :     PG_RETURN_BOOL(FPle(point_dt(&circle1->center, &circle2->center),
    4824             :                         float8_mi(circle2->radius, circle1->radius)));
    4825             : }
    4826             : 
    4827             : /*      circle_contain  -       does circle1 contain circle2?
    4828             :  */
    4829             : Datum
    4830         256 : circle_contain(PG_FUNCTION_ARGS)
    4831             : {
    4832         256 :     CIRCLE     *circle1 = PG_GETARG_CIRCLE_P(0);
    4833         256 :     CIRCLE     *circle2 = PG_GETARG_CIRCLE_P(1);
    4834             : 
    4835         256 :     PG_RETURN_BOOL(FPle(point_dt(&circle1->center, &circle2->center),
    4836             :                         float8_mi(circle1->radius, circle2->radius)));
    4837             : }
    4838             : 
    4839             : 
    4840             : /*      circle_below        -       is circle1 strictly below circle2?
    4841             :  */
    4842             : Datum
    4843         256 : circle_below(PG_FUNCTION_ARGS)
    4844             : {
    4845         256 :     CIRCLE     *circle1 = PG_GETARG_CIRCLE_P(0);
    4846         256 :     CIRCLE     *circle2 = PG_GETARG_CIRCLE_P(1);
    4847             : 
    4848         256 :     PG_RETURN_BOOL(FPlt(float8_pl(circle1->center.y, circle1->radius),
    4849             :                         float8_mi(circle2->center.y, circle2->radius)));
    4850             : }
    4851             : 
    4852             : /*      circle_above    -       is circle1 strictly above circle2?
    4853             :  */
    4854             : Datum
    4855         256 : circle_above(PG_FUNCTION_ARGS)
    4856             : {
    4857         256 :     CIRCLE     *circle1 = PG_GETARG_CIRCLE_P(0);
    4858         256 :     CIRCLE     *circle2 = PG_GETARG_CIRCLE_P(1);
    4859             : 
    4860         256 :     PG_RETURN_BOOL(FPgt(float8_mi(circle1->center.y, circle1->radius),
    4861             :                         float8_pl(circle2->center.y, circle2->radius)));
    4862             : }
    4863             : 
    4864             : /*      circle_overbelow -      is the upper edge of circle1 at or below
    4865             :  *                              the upper edge of circle2?
    4866             :  */
    4867             : Datum
    4868         256 : circle_overbelow(PG_FUNCTION_ARGS)
    4869             : {
    4870         256 :     CIRCLE     *circle1 = PG_GETARG_CIRCLE_P(0);
    4871         256 :     CIRCLE     *circle2 = PG_GETARG_CIRCLE_P(1);
    4872             : 
    4873         256 :     PG_RETURN_BOOL(FPle(float8_pl(circle1->center.y, circle1->radius),
    4874             :                         float8_pl(circle2->center.y, circle2->radius)));
    4875             : }
    4876             : 
    4877             : /*      circle_overabove    -   is the lower edge of circle1 at or above
    4878             :  *                              the lower edge of circle2?
    4879             :  */
    4880             : Datum
    4881         256 : circle_overabove(PG_FUNCTION_ARGS)
    4882             : {
    4883         256 :     CIRCLE     *circle1 = PG_GETARG_CIRCLE_P(0);
    4884         256 :     CIRCLE     *circle2 = PG_GETARG_CIRCLE_P(1);
    4885             : 
    4886         256 :     PG_RETURN_BOOL(FPge(float8_mi(circle1->center.y, circle1->radius),
    4887             :                         float8_mi(circle2->center.y, circle2->radius)));
    4888             : }
    4889             : 
    4890             : 
    4891             : /*      circle_relop    -       is area(circle1) relop area(circle2), within
    4892             :  *                              our accuracy constraint?
    4893             :  */
    4894             : Datum
    4895         256 : circle_eq(PG_FUNCTION_ARGS)
    4896             : {
    4897         256 :     CIRCLE     *circle1 = PG_GETARG_CIRCLE_P(0);
    4898         256 :     CIRCLE     *circle2 = PG_GETARG_CIRCLE_P(1);
    4899             : 
    4900         256 :     PG_RETURN_BOOL(FPeq(circle_ar(circle1), circle_ar(circle2)));
    4901             : }
    4902             : 
    4903             : Datum
    4904         256 : circle_ne(PG_FUNCTION_ARGS)
    4905             : {
    4906         256 :     CIRCLE     *circle1 = PG_GETARG_CIRCLE_P(0);
    4907         256 :     CIRCLE     *circle2 = PG_GETARG_CIRCLE_P(1);
    4908             : 
    4909         256 :     PG_RETURN_BOOL(FPne(circle_ar(circle1), circle_ar(circle2)));
    4910             : }
    4911             : 
    4912             : Datum
    4913        1052 : circle_lt(PG_FUNCTION_ARGS)
    4914             : {
    4915        1052 :     CIRCLE     *circle1 = PG_GETARG_CIRCLE_P(0);
    4916        1052 :     CIRCLE     *circle2 = PG_GETARG_CIRCLE_P(1);
    4917             : 
    4918        1052 :     PG_RETURN_BOOL(FPlt(circle_ar(circle1), circle_ar(circle2)));
    4919             : }
    4920             : 
    4921             : Datum
    4922         256 : circle_gt(PG_FUNCTION_ARGS)
    4923             : {
    4924         256 :     CIRCLE     *circle1 = PG_GETARG_CIRCLE_P(0);
    4925         256 :     CIRCLE     *circle2 = PG_GETARG_CIRCLE_P(1);
    4926             : 
    4927         256 :     PG_RETURN_BOOL(FPgt(circle_ar(circle1), circle_ar(circle2)));
    4928             : }
    4929             : 
    4930             : Datum
    4931         256 : circle_le(PG_FUNCTION_ARGS)
    4932             : {
    4933         256 :     CIRCLE     *circle1 = PG_GETARG_CIRCLE_P(0);
    4934         256 :     CIRCLE     *circle2 = PG_GETARG_CIRCLE_P(1);
    4935             : 
    4936         256 :     PG_RETURN_BOOL(FPle(circle_ar(circle1), circle_ar(circle2)));
    4937             : }
    4938             : 
    4939             : Datum
    4940         324 : circle_ge(PG_FUNCTION_ARGS)
    4941             : {
    4942         324 :     CIRCLE     *circle1 = PG_GETARG_CIRCLE_P(0);
    4943         324 :     CIRCLE     *circle2 = PG_GETARG_CIRCLE_P(1);
    4944             : 
    4945         324 :     PG_RETURN_BOOL(FPge(circle_ar(circle1), circle_ar(circle2)));
    4946             : }
    4947             : 
    4948             : 
    4949             : /*----------------------------------------------------------
    4950             :  *  "Arithmetic" operators on circles.
    4951             :  *---------------------------------------------------------*/
    4952             : 
    4953             : /* circle_add_pt()
    4954             :  * Translation operator.
    4955             :  */
    4956             : Datum
    4957         288 : circle_add_pt(PG_FUNCTION_ARGS)
    4958             : {
    4959         288 :     CIRCLE     *circle = PG_GETARG_CIRCLE_P(0);
    4960         288 :     Point      *point = PG_GETARG_POINT_P(1);
    4961             :     CIRCLE     *result;
    4962             : 
    4963         288 :     result = (CIRCLE *) palloc(sizeof(CIRCLE));
    4964             : 
    4965         288 :     point_add_point(&result->center, &circle->center, point);
    4966         288 :     result->radius = circle->radius;
    4967             : 
    4968         288 :     PG_RETURN_CIRCLE_P(result);
    4969             : }
    4970             : 
    4971             : Datum
    4972         288 : circle_sub_pt(PG_FUNCTION_ARGS)
    4973             : {
    4974         288 :     CIRCLE     *circle = PG_GETARG_CIRCLE_P(0);
    4975         288 :     Point      *point = PG_GETARG_POINT_P(1);
    4976             :     CIRCLE     *result;
    4977             : 
    4978         288 :     result = (CIRCLE *) palloc(sizeof(CIRCLE));
    4979             : 
    4980         288 :     point_sub_point(&result->center, &circle->center, point);
    4981         288 :     result->radius = circle->radius;
    4982             : 
    4983         288 :     PG_RETURN_CIRCLE_P(result);
    4984             : }
    4985             : 
    4986             : 
    4987             : /* circle_mul_pt()
    4988             :  * Rotation and scaling operators.
    4989             :  */
    4990             : Datum
    4991         288 : circle_mul_pt(PG_FUNCTION_ARGS)
    4992             : {
    4993         288 :     CIRCLE     *circle = PG_GETARG_CIRCLE_P(0);
    4994         288 :     Point      *point = PG_GETARG_POINT_P(1);
    4995             :     CIRCLE     *result;
    4996             : 
    4997         288 :     result = (CIRCLE *) palloc(sizeof(CIRCLE));
    4998             : 
    4999         288 :     point_mul_point(&result->center, &circle->center, point);
    5000         288 :     result->radius = float8_mul(circle->radius, HYPOT(point->x, point->y));
    5001             : 
    5002         288 :     PG_RETURN_CIRCLE_P(result);
    5003             : }
    5004             : 
    5005             : Datum
    5006          72 : circle_div_pt(PG_FUNCTION_ARGS)
    5007             : {
    5008          72 :     CIRCLE     *circle = PG_GETARG_CIRCLE_P(0);
    5009          72 :     Point      *point = PG_GETARG_POINT_P(1);
    5010             :     CIRCLE     *result;
    5011             : 
    5012          72 :     result = (CIRCLE *) palloc(sizeof(CIRCLE));
    5013             : 
    5014          72 :     point_div_point(&result->center, &circle->center, point);
    5015          64 :     result->radius = float8_div(circle->radius, HYPOT(point->x, point->y));
    5016             : 
    5017          64 :     PG_RETURN_CIRCLE_P(result);
    5018             : }
    5019             : 
    5020             : 
    5021             : /*      circle_area     -       returns the area of the circle.
    5022             :  */
    5023             : Datum
    5024         308 : circle_area(PG_FUNCTION_ARGS)
    5025             : {
    5026         308 :     CIRCLE     *circle = PG_GETARG_CIRCLE_P(0);
    5027             : 
    5028         308 :     PG_RETURN_FLOAT8(circle_ar(circle));
    5029             : }
    5030             : 
    5031             : 
    5032             : /*      circle_diameter -       returns the diameter of the circle.
    5033             :  */
    5034             : Datum
    5035          64 : circle_diameter(PG_FUNCTION_ARGS)
    5036             : {
    5037          64 :     CIRCLE     *circle = PG_GETARG_CIRCLE_P(0);
    5038             : 
    5039          64 :     PG_RETURN_FLOAT8(float8_mul(circle->radius, 2.0));
    5040             : }
    5041             : 
    5042             : 
    5043             : /*      circle_radius   -       returns the radius of the circle.
    5044             :  */
    5045             : Datum
    5046       12504 : circle_radius(PG_FUNCTION_ARGS)
    5047             : {
    5048       12504 :     CIRCLE     *circle = PG_GETARG_CIRCLE_P(0);
    5049             : 
    5050       12504 :     PG_RETURN_FLOAT8(circle->radius);
    5051             : }
    5052             : 
    5053             : 
    5054             : /*      circle_distance -       returns the distance between
    5055             :  *                                two circles.
    5056             :  */
    5057             : Datum
    5058         112 : circle_distance(PG_FUNCTION_ARGS)
    5059             : {
    5060         112 :     CIRCLE     *circle1 = PG_GETARG_CIRCLE_P(0);
    5061         112 :     CIRCLE     *circle2 = PG_GETARG_CIRCLE_P(1);
    5062             :     float8      result;
    5063             : 
    5064         112 :     result = float8_mi(point_dt(&circle1->center, &circle2->center),
    5065             :                        float8_pl(circle1->radius, circle2->radius));
    5066         112 :     if (result < 0.0)
    5067          48 :         result = 0.0;
    5068             : 
    5069         112 :     PG_RETURN_FLOAT8(result);
    5070             : }
    5071             : 
    5072             : 
    5073             : Datum
    5074          16 : circle_contain_pt(PG_FUNCTION_ARGS)
    5075             : {
    5076          16 :     CIRCLE     *circle = PG_GETARG_CIRCLE_P(0);
    5077          16 :     Point      *point = PG_GETARG_POINT_P(1);
    5078             :     float8      d;
    5079             : 
    5080          16 :     d = point_dt(&circle->center, point);
    5081          16 :     PG_RETURN_BOOL(d <= circle->radius);
    5082             : }
    5083             : 
    5084             : 
    5085             : Datum
    5086          36 : pt_contained_circle(PG_FUNCTION_ARGS)
    5087             : {
    5088          36 :     Point      *point = PG_GETARG_POINT_P(0);
    5089          36 :     CIRCLE     *circle = PG_GETARG_CIRCLE_P(1);
    5090             :     float8      d;
    5091             : 
    5092          36 :     d = point_dt(&circle->center, point);
    5093          36 :     PG_RETURN_BOOL(d <= circle->radius);
    5094             : }
    5095             : 
    5096             : 
    5097             : /*      dist_pc -       returns the distance between
    5098             :  *                        a point and a circle.
    5099             :  */
    5100             : Datum
    5101         500 : dist_pc(PG_FUNCTION_ARGS)
    5102             : {
    5103         500 :     Point      *point = PG_GETARG_POINT_P(0);
    5104         500 :     CIRCLE     *circle = PG_GETARG_CIRCLE_P(1);
    5105             :     float8      result;
    5106             : 
    5107         500 :     result = float8_mi(point_dt(point, &circle->center),
    5108             :                        circle->radius);
    5109         500 :     if (result < 0.0)
    5110          76 :         result = 0.0;
    5111             : 
    5112         500 :     PG_RETURN_FLOAT8(result);
    5113             : }
    5114             : 
    5115             : /*
    5116             :  * Distance from a circle to a point
    5117             :  */
    5118             : Datum
    5119       12484 : dist_cpoint(PG_FUNCTION_ARGS)
    5120             : {
    5121       12484 :     CIRCLE     *circle = PG_GETARG_CIRCLE_P(0);
    5122       12484 :     Point      *point = PG_GETARG_POINT_P(1);
    5123             :     float8      result;
    5124             : 
    5125       12484 :     result = float8_mi(point_dt(point, &circle->center), circle->radius);
    5126       12484 :     if (result < 0.0)
    5127           0 :         result = 0.0;
    5128             : 
    5129       12484 :     PG_RETURN_FLOAT8(result);
    5130             : }
    5131             : 
    5132             : /*      circle_center   -       returns the center point of the circle.
    5133             :  */
    5134             : Datum
    5135       12604 : circle_center(PG_FUNCTION_ARGS)
    5136             : {
    5137       12604 :     CIRCLE     *circle = PG_GETARG_CIRCLE_P(0);
    5138             :     Point      *result;
    5139             : 
    5140       12604 :     result = (Point *) palloc(sizeof(Point));
    5141       12604 :     result->x = circle->center.x;
    5142       12604 :     result->y = circle->center.y;
    5143             : 
    5144       12604 :     PG_RETURN_POINT_P(result);
    5145             : }
    5146             : 
    5147             : 
    5148             : /*      circle_ar       -       returns the area of the circle.
    5149             :  */
    5150             : static float8
    5151        5108 : circle_ar(CIRCLE *circle)
    5152             : {
    5153        5108 :     return float8_mul(float8_mul(circle->radius, circle->radius), M_PI);
    5154             : }
    5155             : 
    5156             : 
    5157             : /*----------------------------------------------------------
    5158             :  *  Conversion operators.
    5159             :  *---------------------------------------------------------*/
    5160             : 
    5161             : Datum
    5162       80052 : cr_circle(PG_FUNCTION_ARGS)
    5163             : {
    5164       80052 :     Point      *center = PG_GETARG_POINT_P(0);
    5165       80052 :     float8      radius = PG_GETARG_FLOAT8(1);
    5166             :     CIRCLE     *result;
    5167             : 
    5168       80052 :     result = (CIRCLE *) palloc(sizeof(CIRCLE));
    5169             : 
    5170       80052 :     result->center.x = center->x;
    5171       80052 :     result->center.y = center->y;
    5172       80052 :     result->radius = radius;
    5173             : 
    5174       80052 :     PG_RETURN_CIRCLE_P(result);
    5175             : }
    5176             : 
    5177             : Datum
    5178          32 : circle_box(PG_FUNCTION_ARGS)
    5179             : {
    5180          32 :     CIRCLE     *circle = PG_GETARG_CIRCLE_P(0);
    5181             :     BOX        *box;
    5182             :     float8      delta;
    5183             : 
    5184          32 :     box = (BOX *) palloc(sizeof(BOX));
    5185             : 
    5186          32 :     delta = float8_div(circle->radius, sqrt(2.0));
    5187             : 
    5188          32 :     box->high.x = float8_pl(circle->center.x, delta);
    5189          32 :     box->low.x = float8_mi(circle->center.x, delta);
    5190          32 :     box->high.y = float8_pl(circle->center.y, delta);
    5191          32 :     box->low.y = float8_mi(circle->center.y, delta);
    5192             : 
    5193          32 :     PG_RETURN_BOX_P(box);
    5194             : }
    5195             : 
    5196             : /* box_circle()
    5197             :  * Convert a box to a circle.
    5198             :  */
    5199             : Datum
    5200       12420 : box_circle(PG_FUNCTION_ARGS)
    5201             : {
    5202       12420 :     BOX        *box = PG_GETARG_BOX_P(0);
    5203             :     CIRCLE     *circle;
    5204             : 
    5205       12420 :     circle = (CIRCLE *) palloc(sizeof(CIRCLE));
    5206             : 
    5207       12420 :     circle->center.x = float8_div(float8_pl(box->high.x, box->low.x), 2.0);
    5208       12420 :     circle->center.y = float8_div(float8_pl(box->high.y, box->low.y), 2.0);
    5209             : 
    5210       12420 :     circle->radius = point_dt(&circle->center, &box->high);
    5211             : 
    5212       12420 :     PG_RETURN_CIRCLE_P(circle);
    5213             : }
    5214             : 
    5215             : 
    5216             : Datum
    5217       40056 : circle_poly(PG_FUNCTION_ARGS)
    5218             : {
    5219       40056 :     int32       npts = PG_GETARG_INT32(0);
    5220       40056 :     CIRCLE     *circle = PG_GETARG_CIRCLE_P(1);
    5221             :     POLYGON    *poly;
    5222             :     int         base_size,
    5223             :                 size;
    5224             :     int         i;
    5225             :     float8      angle;
    5226             :     float8      anglestep;
    5227             : 
    5228       40056 :     if (FPzero(circle->radius))
    5229           4 :         ereport(ERROR,
    5230             :                 (errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
    5231             :                  errmsg("cannot convert circle with radius zero to polygon")));
    5232             : 
    5233       40052 :     if (npts < 2)
    5234           4 :         ereport(ERROR,
    5235             :                 (errcode(ERRCODE_INVALID_PARAMETER_VALUE),
    5236             :                  errmsg("must request at least 2 points")));
    5237             : 
    5238       40048 :     base_size = sizeof(poly->p[0]) * npts;
    5239       40048 :     size = offsetof(POLYGON, p) + base_size;
    5240             : 
    5241             :     /* Check for integer overflow */
    5242       40048 :     if (base_size / npts != sizeof(poly->p[0]) || size <= base_size)
    5243           0 :         ereport(ERROR,
    5244             :                 (errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED),
    5245             :                  errmsg("too many points requested")));
    5246             : 
    5247       40048 :     poly = (POLYGON *) palloc0(size);   /* zero any holes */
    5248       40048 :     SET_VARSIZE(poly, size);
    5249       40048 :     poly->npts = npts;
    5250             : 
    5251       40048 :     anglestep = float8_div(2.0 * M_PI, npts);
    5252             : 
    5253      520528 :     for (i = 0; i < npts; i++)
    5254             :     {
    5255      480480 :         angle = float8_mul(anglestep, i);
    5256             : 
    5257      480480 :         poly->p[i].x = float8_mi(circle->center.x,
    5258             :                                  float8_mul(circle->radius, cos(angle)));
    5259      480480 :         poly->p[i].y = float8_pl(circle->center.y,
    5260             :                                  float8_mul(circle->radius, sin(angle)));
    5261             :     }
    5262             : 
    5263       40048 :     make_bound_box(poly);
    5264             : 
    5265       40048 :     PG_RETURN_POLYGON_P(poly);
    5266             : }
    5267             : 
    5268             : /*
    5269             :  * Convert polygon to circle
    5270             :  *
    5271             :  * The result must be preallocated.
    5272             :  *
    5273             :  * XXX This algorithm should use weighted means of line segments
    5274             :  *  rather than straight average values of points - tgl 97/01/21.
    5275             :  */
    5276             : static void
    5277          48 : poly_to_circle(CIRCLE *result, POLYGON *poly)
    5278             : {
    5279             :     int         i;
    5280             : 
    5281             :     Assert(poly->npts > 0);
    5282             : 
    5283          48 :     result->center.x = 0;
    5284          48 :     result->center.y = 0;
    5285          48 :     result->radius = 0;
    5286             : 
    5287         216 :     for (i = 0; i < poly->npts; i++)
    5288         168 :         point_add_point(&result->center, &result->center, &poly->p[i]);
    5289          48 :     result->center.x = float8_div(result->center.x, poly->npts);
    5290          48 :     result->center.y = float8_div(result->center.y, poly->npts);
    5291             : 
    5292         216 :     for (i = 0; i < poly->npts; i++)
    5293         168 :         result->radius = float8_pl(result->radius,
    5294             :                                    point_dt(&poly->p[i], &result->center));
    5295          48 :     result->radius = float8_div(result->radius, poly->npts);
    5296          48 : }
    5297             : 
    5298             : Datum
    5299          20 : poly_circle(PG_FUNCTION_ARGS)
    5300             : {
    5301          20 :     POLYGON    *poly = PG_GETARG_POLYGON_P(0);
    5302             :     CIRCLE     *result;
    5303             : 
    5304          20 :     result = (CIRCLE *) palloc(sizeof(CIRCLE));
    5305             : 
    5306          20 :     poly_to_circle(result, poly);
    5307             : 
    5308          20 :     PG_RETURN_CIRCLE_P(result);
    5309             : }
    5310             : 
    5311             : 
    5312             : /***********************************************************************
    5313             :  **
    5314             :  **     Private routines for multiple types.
    5315             :  **
    5316             :  ***********************************************************************/
    5317             : 
    5318             : /*
    5319             :  *  Test to see if the point is inside the polygon, returns 1/0, or 2 if
    5320             :  *  the point is on the polygon.
    5321             :  *  Code adapted but not copied from integer-based routines in WN: A
    5322             :  *  Server for the HTTP
    5323             :  *  version 1.15.1, file wn/image.c
    5324             :  *  http://hopf.math.northwestern.edu/index.html
    5325             :  *  Description of algorithm:  http://www.linuxjournal.com/article/2197
    5326             :  *                             http://www.linuxjournal.com/article/2029
    5327             :  */
    5328             : 
    5329             : #define POINT_ON_POLYGON INT_MAX
    5330             : 
    5331             : static int
    5332      164092 : point_inside(Point *p, int npts, Point *plist)
    5333             : {
    5334             :     float8      x0,
    5335             :                 y0;
    5336             :     float8      prev_x,
    5337             :                 prev_y;
    5338      164092 :     int         i = 0;
    5339             :     float8      x,
    5340             :                 y;
    5341             :     int         cross,
    5342      164092 :                 total_cross = 0;
    5343             : 
    5344             :     Assert(npts > 0);
    5345             : 
    5346             :     /* compute first polygon point relative to single point */
    5347      164092 :     x0 = float8_mi(plist[0].x, p->x);
    5348      164092 :     y0 = float8_mi(plist[0].y, p->y);
    5349             : 
    5350      164092 :     prev_x = x0;
    5351      164092 :     prev_y = y0;
    5352             :     /* loop over polygon points and aggregate total_cross */
    5353      757564 :     for (i = 1; i < npts; i++)
    5354             :     {
    5355             :         /* compute next polygon point relative to single point */
    5356      593548 :         x = float8_mi(plist[i].x, p->x);
    5357      593548 :         y = float8_mi(plist[i].y, p->y);
    5358             : 
    5359             :         /* compute previous to current point crossing */
    5360      593548 :         if ((cross = lseg_crossing(x, y, prev_x, prev_y)) == POINT_ON_POLYGON)
    5361          76 :             return 2;
    5362      593472 :         total_cross += cross;
    5363             : 
    5364      593472 :         prev_x = x;
    5365      593472 :         prev_y = y;
    5366             :     }
    5367             : 
    5368             :     /* now do the first point */
    5369      164016 :     if ((cross = lseg_crossing(x0, y0, prev_x, prev_y)) == POINT_ON_POLYGON)
    5370          60 :         return 2;
    5371      163956 :     total_cross += cross;
    5372             : 
    5373      163956 :     if (total_cross != 0)
    5374      126128 :         return 1;
    5375       37828 :     return 0;
    5376             : }
    5377             : 
    5378             : 
    5379             : /* lseg_crossing()
    5380             :  * Returns +/-2 if line segment crosses the positive X-axis in a +/- direction.
    5381             :  * Returns +/-1 if one point is on the positive X-axis.
    5382             :  * Returns 0 if both points are on the positive X-axis, or there is no crossing.
    5383             :  * Returns POINT_ON_POLYGON if the segment contains (0,0).
    5384             :  * Wow, that is one confusing API, but it is used above, and when summed,
    5385             :  * can tell is if a point is in a polygon.
    5386             :  */
    5387             : 
    5388             : static int
    5389      757564 : lseg_crossing(float8 x, float8 y, float8 prev_x, float8 prev_y)
    5390             : {
    5391             :     float8      z;
    5392             :     int         y_sign;
    5393             : 
    5394      757564 :     if (FPzero(y))
    5395             :     {                           /* y == 0, on X axis */
    5396        1132 :         if (FPzero(x))          /* (x,y) is (0,0)? */
    5397         128 :             return POINT_ON_POLYGON;
    5398        1004 :         else if (FPgt(x, 0))
    5399             :         {                       /* x > 0 */
    5400         668 :             if (FPzero(prev_y)) /* y and prev_y are zero */
    5401             :                 /* prev_x > 0? */
    5402          40 :                 return FPgt(prev_x, 0.0) ? 0 : POINT_ON_POLYGON;
    5403         628 :             return FPlt(prev_y, 0.0) ? 1 : -1;
    5404             :         }
    5405             :         else
    5406             :         {                       /* x < 0, x not on positive X axis */
    5407         336 :             if (FPzero(prev_y))
    5408             :                 /* prev_x < 0? */
    5409          16 :                 return FPlt(prev_x, 0.0) ? 0 : POINT_ON_POLYGON;
    5410         320 :             return 0;
    5411             :         }
    5412             :     }
    5413             :     else
    5414             :     {                           /* y != 0 */
    5415             :         /* compute y crossing direction from previous point */
    5416      756432 :         y_sign = FPgt(y, 0.0) ? 1 : -1;
    5417             : 
    5418      756432 :         if (FPzero(prev_y))
    5419             :             /* previous point was on X axis, so new point is either off or on */
    5420        1004 :             return FPlt(prev_x, 0.0) ? 0 : y_sign;
    5421      755428 :         else if ((y_sign < 0 && FPlt(prev_y, 0.0)) ||
    5422      386196 :                  (y_sign > 0 && FPgt(prev_y, 0.0)))
    5423             :             /* both above or below X axis */
    5424      480756 :             return 0;           /* same sign */
    5425             :         else
    5426             :         {                       /* y and prev_y cross X-axis */
    5427      274672 :             if (FPge(x, 0.0) && FPgt(prev_x, 0.0))
    5428             :                 /* both non-negative so cross positive X-axis */
    5429       99892 :                 return 2 * y_sign;
    5430      174780 :             if (FPlt(x, 0.0) && FPle(prev_x, 0.0))
    5431             :                 /* both non-positive so do not cross positive X-axis */
    5432       88296 :                 return 0;
    5433             : 
    5434             :             /* x and y cross axes, see URL above point_inside() */
    5435       86484 :             z = float8_mi(float8_mul(float8_mi(x, prev_x), y),
    5436             :                           float8_mul(float8_mi(y, prev_y), x));
    5437       86484 :             if (FPzero(z))
    5438           8 :                 return POINT_ON_POLYGON;
    5439       86476 :             if ((y_sign < 0 && FPlt(z, 0.0)) ||
    5440       49624 :                 (y_sign > 0 && FPgt(z, 0.0)))
    5441       48456 :                 return 0;
    5442       38020 :             return 2 * y_sign;
    5443             :         }
    5444             :     }
    5445             : }
    5446             : 
    5447             : 
    5448             : static bool
    5449        4062 : plist_same(int npts, Point *p1, Point *p2)
    5450             : {
    5451             :     int         i,
    5452             :                 ii,
    5453             :                 j;
    5454             : 
    5455             :     /* find match for first point */
    5456        4174 :     for (i = 0; i < npts; i++)
    5457             :     {
    5458        4150 :         if (point_eq_point(&p2[i], &p1[0]))
    5459             :         {
    5460             : 
    5461             :             /* match found? then look forward through remaining points */
    5462       12124 :             for (ii = 1, j = i + 1; ii < npts; ii++, j++)
    5463             :             {
    5464        8094 :                 if (j >= npts)
    5465          12 :                     j = 0;
    5466        8094 :                 if (!point_eq_point(&p2[j], &p1[ii]))
    5467          24 :                     break;
    5468             :             }
    5469        4054 :             if (ii == npts)
    5470        4030 :                 return true;
    5471             : 
    5472             :             /* match not found forwards? then look backwards */
    5473          56 :             for (ii = 1, j = i - 1; ii < npts; ii++, j--)
    5474             :             {
    5475          48 :                 if (j < 0)
    5476           8 :                     j = (npts - 1);
    5477          48 :                 if (!point_eq_point(&p2[j], &p1[ii]))
    5478          16 :                     break;
    5479             :             }
    5480          24 :             if (ii == npts)
    5481           8 :                 return true;
    5482             :         }
    5483             :     }
    5484             : 
    5485          24 :     return false;
    5486             : }
    5487             : 
    5488             : 
    5489             : /*-------------------------------------------------------------------------
    5490             :  * Determine the hypotenuse.
    5491             :  *
    5492             :  * If required, x and y are swapped to make x the larger number. The
    5493             :  * traditional formula of x^2+y^2 is rearranged to factor x outside the
    5494             :  * sqrt. This allows computation of the hypotenuse for significantly
    5495             :  * larger values, and with a higher precision than when using the naive
    5496             :  * formula.  In particular, this cannot overflow unless the final result
    5497             :  * would be out-of-range.
    5498             :  *
    5499             :  * sqrt( x^2 + y^2 ) = sqrt( x^2( 1 + y^2/x^2) )
    5500             :  *                   = x * sqrt( 1 + y^2/x^2 )
    5501             :  *                   = x * sqrt( 1 + y/x * y/x )
    5502             :  *
    5503             :  * It is expected that this routine will eventually be replaced with the
    5504             :  * C99 hypot() function.
    5505             :  *
    5506             :  * This implementation conforms to IEEE Std 1003.1 and GLIBC, in that the
    5507             :  * case of hypot(inf,nan) results in INF, and not NAN.
    5508             :  *-----------------------------------------------------------------------
    5509             :  */
    5510             : float8
    5511    10390496 : pg_hypot(float8 x, float8 y)
    5512             : {
    5513             :     float8      yx,
    5514             :                 result;
    5515             : 
    5516             :     /* Handle INF and NaN properly */
    5517    10390496 :     if (isinf(x) || isinf(y))
    5518        2280 :         return get_float8_infinity();
    5519             : 
    5520    10388216 :     if (isnan(x) || isnan(y))
    5521       15700 :         return get_float8_nan();
    5522             : 
    5523             :     /* Else, drop any minus signs */
    5524    10372516 :     x = fabs(x);
    5525    10372516 :     y = fabs(y);
    5526             : 
    5527             :     /* Swap x and y if needed to make x the larger one */
    5528    10372516 :     if (x < y)
    5529             :     {
    5530     5074088 :         float8      temp = x;
    5531             : 
    5532     5074088 :         x = y;
    5533     5074088 :         y = temp;
    5534             :     }
    5535             : 
    5536             :     /*
    5537             :      * If y is zero, the hypotenuse is x.  This test saves a few cycles in
    5538             :      * such cases, but more importantly it also protects against
    5539             :      * divide-by-zero errors, since now x >= y.
    5540             :      */
    5541    10372516 :     if (y == 0.0)
    5542     2656344 :         return x;
    5543             : 
    5544             :     /* Determine the hypotenuse */
    5545     7716172 :     yx = y / x;
    5546     7716172 :     result = x * sqrt(1.0 + (yx * yx));
    5547             : 
    5548     7716172 :     check_float8_val(result, false, false);
    5549             : 
    5550     7716172 :     return result;
    5551             : }

Generated by: LCOV version 1.13